Millisecond-level cooperative inhibition method for flexible load of harbor district based on magnetic flux jump prediction

Abstract

The invention discloses a port area flexible load millisecond-level cooperative suppression method based on magnetic flux jump prediction, and the method comprises the steps: obtaining the external vibration and electrical data of a transformer, combining the Bayesian inversion of physical constraint, and reconstructing an internal three-dimensional magnetic flux field; based on time sequence evolution of the magnetic flux field, a dual-mode self-adaptive prediction algorithm is adopted, and a magnetic flux jump early warning vector containing jump time, amplitude and position is generated; and according to the early warning vector, through a multi-layer cooperative control strategy of time sequence decoupling, a control instruction for the flexible load of the harbor district is generated and executed. Advanced early warning and millisecond-level active suppression of the potential magnetic saturation risk are realized, and the operation reliability and equipment safety of a harbor district power distribution system are improved.

Description

TECHNICAL FIELD

[0001] The present application relates to the technical field of power system stability and control, and particularly relates to a millisecond-level collaborative suppression method for port flexible load based on magnetic flux jump prediction. BACKGROUND

[0002] As a key hub of international logistics and trade, the power distribution system of the port bears the important mission of ensuring the efficient and reliable operation of the port. Unlike conventional distribution networks, the port power distribution system has unique operating characteristics. On the one hand, the high-power ship load of the shore power system needs to be frequently connected and disconnected. On the other hand, the power demand of megawatt heavy lifting equipment such as shore cranes and yard cranes shows periodic and impulsive fluctuations during operation. High-frequency and large-amplitude power flow fluctuations directly impact the core equipment of the system, i.e., the power transformer, on the millisecond time scale. The dramatic current change rate can cause the magnetic flux inside the transformer core to jump sharply, which can easily lead to the local or overall saturation of the core. Transformer magnetic saturation can cause a series of negative effects, including excitation current distortion, dramatic increase in harmonic content, abnormal increase in core loss and vibration noise, and in severe cases, can cause protection device misoperation, equipment insulation accelerated aging, and even permanent damage, posing a threat to the reliability and safety of port power supply. Therefore, researching how to accurately perceive, predict and actively suppress transformer magnetic flux jump has important theoretical value and engineering significance for ensuring the safe and stable operation of key infrastructure in the port.

[0003] Currently, the existing technologies for monitoring and protection of transformer operating state mainly focus on several aspects. In terms of monitoring methods, online monitoring methods based on electrical quantities are commonly used, such as analyzing the voltage and current waveforms at the primary and secondary sides of the transformer to calculate indicators such as harmonic content and power factor to indirectly determine the operating state. At the same time, non-electrical quantity monitoring technologies are also widely used, such as monitoring the temperature distribution of the winding and core through infrared thermal imaging to detect overheating faults, or collecting box vibration signals through vibration acceleration sensors and using Fourier transform and other frequency domain analysis methods to diagnose mechanical faults such as winding looseness and core multi-point grounding. In terms of prediction, existing researches mainly focus on statistical methods or machine learning models based on historical load data to predict future minute-level or hour-level load curves to provide a basis for economic dispatch of the power grid. In terms of control and suppression, to cope with power flow impact, traditional methods mainly use reactive power compensation devices (such as SVC and STATCOM) for dynamic voltage support, or in extreme cases, use protection devices to implement load shedding or transformer tripping as passive protection strategies.

[0004] However, the existing technologies still face deep technical bottlenecks when dealing with the specific problem of millisecond-level magnetic flux jump of the transformer in the port. These problems are interrelated and collectively lead to the lag and inefficiency of suppression measures.

[0005] Current technologies generally face the challenge of observing the internal magnetic field state of transformers. Whether through external electrical quantities or single physical quantities such as vibration and temperature, traditional methods provide an indirect, macroscopic reflection of the complex electromagnetic processes within the transformer. They cannot penetrate the transformer's physical barriers to accurately reconstruct the three-dimensional spatial distribution of the internal magnetic flux density and its dynamic evolution in real time. Just as a doctor cannot obtain CT images based solely on external symptoms, diagnosis and decision-making lack the most direct and crucial physical basis. Without a precise picture of where, when, and to what extent the magnetic flux will saturate, subsequent predictions and control will lack data support and be poorly planned.

[0006] Due to the aforementioned observational challenges, existing technologies suffer from a time-scale mismatch in prediction. Magnetic flux jumps are rapid physical processes occurring on the order of milliseconds, while traditional methods based on electrical quantity statistics or load curves typically have response and prediction timescales on the order of seconds, minutes, or even longer. This time-scale difference means that existing prediction methods can only perform post-event analysis of magnetic flux jumps, failing to provide early warning. By the time severe current distortion or abnormal vibration is detected, the magnetic saturation process has often already occurred or even intensified, and the system has lost the optimal window for taking proactive suppression measures, leaving it only able to passively withstand the impact.

[0007] The limitations of observation and the lag in prediction lead to a dilemma in the coordinated suppression of existing control methods. Traditional control strategies, such as reactive power compensation or uniform load shedding, are single-layer, open-loop control modes. They cannot finely schedule various flexible load resources with different response characteristics within the port area (such as fast-responding but limited-capacity energy storage, slower-responding but large-scale charging piles, and quay cranes that can only be adjusted in one direction) based on the specific location, severity, and time of predicted saturation. The lack of coordinated control methods results in missed opportunities due to untimely responses and unnecessary economic losses or excessive interference with port operations due to single and crude adjustment methods. It is impossible to achieve a balance between speed, accuracy, and cost, and it is difficult to achieve accurate and efficient suppression of magnetic flux jumps. Summary of the Invention

[0008] The purpose of this invention is to propose a millisecond-level collaborative suppression method for flexible loads in port areas based on magnetic flux jump prediction, in order to solve the aforementioned problems existing in the prior art.

[0009] This invention proposes a millisecond-level collaborative suppression method for flexible loads in port areas based on magnetic flux jump prediction, comprising: The vibration signals and electrical operation data of the transformer's exterior are acquired, and the three-dimensional magnetic flux density distribution field characterizing the internal magnetic field state of the transformer is reconstructed based on these data. Receive the three-dimensional magnetic flux density distribution field, predict future magnetic flux jump events based on its temporal evolution law, and generate a magnetic flux jump early warning vector containing jump time, amplitude and location information; Based on the magnetic flux jump early warning vector and the real-time acquisition of the port area's flexible load operation status, hierarchical collaborative control commands are generated and issued for execution.

[0010] Compared with existing technologies, the beneficial effects are: This invention solves the technical problem of difficulty in real-time observation, rapid prediction and coordinated suppression of transformer magnetic flux jumps caused by frequent start-stop of high-power loads in port areas. It realizes advanced early warning and millisecond-level active suppression of potential magnetic saturation risks, and improves the operational reliability and equipment safety of the port area power distribution system. Attached Figure Description

[0011] Figure 1 This is a schematic diagram of a millisecond-level collaborative suppression method for flexible loads in port areas based on magnetic flux jump prediction. Figure 2 A schematic diagram of the three-dimensional magnetic flux density distribution field process for inverting and reconstructing the internal magnetic field state of a transformer; Figure 3 A schematic diagram of the steps for separating magnetostrictive vibration components; Figure 4 A flowchart illustrating the process of generating and issuing hierarchical collaborative control instructions for execution. Detailed Implementation

[0012] To enable those skilled in the art to better understand the present invention, the technical solutions of the present invention will be clearly and completely described below with reference to the accompanying drawings of the embodiments. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort should fall within the scope of protection of the present invention.

[0013] It should be noted that the terms "first," "second," etc., used in the specification and accompanying drawings of this invention are used to distinguish similar objects and are not necessarily used to describe a specific order or sequence. It should be understood that such data can be interchanged where appropriate so that the embodiments of the invention described herein can be implemented in orders other than those illustrated or described herein.

[0014] It should be noted that, for the purpose of clearly demonstrating the steps of this application, each step has been numbered in the specification. These numbers are for ease of explanation only and do not limit the execution order of the steps. In actual operation, depending on the technical requirements of the specific implementation scenario, the steps may be executed in a different order than that shown in the specification, and in some cases, parallel processing between steps can also be achieved.

[0015] Example 1 describes the basic system architecture and data acquisition environment required to implement the method of the present invention, providing a scenario and data foundation for the implementation of subsequent method steps.

[0016] In a specific application scenario, the system is deployed in the power distribution network of a large port area. The core of this network is one or more high-capacity power transformers, such as a dry-type transformer with a rated capacity of 20MVA and a voltage level of 10kV / 400V. This transformer supplies power to various flexible loads within the port area. Flexible loads are the objects of the collaborative control of this invention. They are diverse in type and have varying response characteristics. Exemplarily, they include: energy storage systems (ESS) with fast bidirectional power regulation capabilities, whose response time can reach the 5ms level; electric truck charging pile groups with a slightly slower response time, approximately 10ms; and large lifting equipment such as yard cranes and quay cranes with relatively slower response times, between 15ms and 20ms. The adjustable range of loads such as quay cranes is typically limited to power reduction, for example, [-0.3*P]. crane ,0].

[0017] To acquire the data required for subsequent calculations, the system is configured with corresponding sensors and data acquisition units. Specifically, to acquire multi-point vibration signals, high-frequency accelerometers are placed at four optimized locations on the transformer tank, with a sampling rate set to 50kHz to capture weak vibrations caused by the magnetostrictive effect. To ensure time synchronization of the multi-point signals, the acquisition system employs hardware clock synchronization and phase calibration technology. To acquire electrical operation data, the system synchronously acquires the instantaneous three-phase voltage values ​​u on the high-voltage or low-voltage side of the transformer using conventional voltage transformers and current transformers. abc (t) and instantaneous value of three-phase current i abc (t). In this embodiment, to facilitate subsequent analysis of the magnetization state, the collected three-phase current i abc (t) is further transformed into the dq-axis current component i in a rotating coordinate system through the Park transformation. d (t) and i q (t).

[0018] Example 2: A general flowchart for a millisecond-level collaborative suppression method for flexible loads in port areas based on magnetic flux jump prediction is proposed, as follows: Figure 1 As shown, the specific steps may include the following: Step 1: Obtain multi-point vibration signals and electrical operation data from the outside of the transformer. Based on the multi-point vibration signals and electrical operation data, reconstruct the three-dimensional magnetic flux density distribution field that characterizes the internal magnetic field state of the transformer.

[0019] Specifically, the multi-point vibration signals and electrical operation data were acquired by the system in Example 1. For example... Figure 2 As shown, the inversion and reconstruction process can preferably be a Bayesian inversion process based on physical constraints. This process correlates external vibrations with the internal magnetic field through a physical model of the magnetostrictive effect and integrates prior physical knowledge such as Maxwell's equations to analyze the internal three-dimensional magnetic flux density distribution field B(x,y,z,t) without intruding into the transformer body. In this embodiment, by establishing a mapping relationship between the easily measurable external vibration signal and the difficult-to-observe internal three-dimensional magnetic flux field, the technical problem of traditional methods being unable to perceive the internal magnetic field state of the transformer in real time and non-invasively is solved.

[0020] Step 2: Receive the three-dimensional magnetic flux density distribution field, predict possible magnetic flux jump events based on its temporal evolution, and generate a magnetic flux jump early warning vector containing jump time, amplitude, and location information.

[0021] Among them, the magnetic flux jump early warning vector is a key data structure containing the prediction results, which can be represented as: W=[ΔB max ,t jump ,x j ,y j ,z j ] T ;where ΔB max The predicted maximum amplitude of the magnetic flux jump is expressed in Tesla (T); t jump The predicted time of the jump is in milliseconds (ms). (x j ,y j ,z j The coordinates are the core spatial locations of the predicted flux jump. Specifically, this prediction process can be a dual-mode adaptive prediction process, intelligently switching between a high-precision, accurate prediction mode and a high-speed, rapid prediction mode based on real-time changes in the power flow status of the port area's power grid, thus balancing accuracy and timeliness. This step transforms the response to flux jumps from reactive remediation to proactive suppression, reserving a valuable millisecond-level response time window for subsequent coordinated control.

[0022] Step 3: Based on the magnetic flux jump early warning vector and the real-time acquisition of the port area's flexible load operation status, generate hierarchical collaborative control commands and issue them for execution.

[0023] For example, hierarchical collaborative control instructions may include instructions at the expected jump time t. jump The first layer of advanced phase compensation instructions and the second layer of key point reverse suppression instructions generated before arrival, and in t jump The third-layer system power optimization and reallocation command is generated after a certain time has elapsed. The operating status of the port area's flexible loads is the real-time power, adjustable capacity, response time, and other information of various loads in Example 1. Through a layered, time-decoupled control strategy, the differences in response characteristics of different flexible loads (such as fast response of energy storage and slow response of quay cranes) can be fully utilized to achieve the optimal balance between control costs and suppression effects, achieving precise collaborative suppression effects at the millisecond level.

[0024] Example 3 provides a preferred embodiment for inverting the spatial distribution of magnetic flux, describing how to invert the three-dimensional magnetic flux density distribution field inside a transformer from a complex multi-point vibration signal. In this embodiment, the inversion and reconstruction process can be specifically divided into two main stages: separation of magnetostrictive vibration components and Bayesian magnetic flux inversion based on physical constraints.

[0025] In the first stage, the magnetostrictive vibration components are separated.

[0026] Before performing the inversion and reconstruction to obtain the three-dimensional magnetic flux density distribution field characterizing the internal magnetic field state of the transformer, a step of separating the magnetostrictive vibration component is included. In this embodiment, since the vibration signal collected during the actual operation of the transformer is a mixture of multiple vibration sources, such as vibrations caused by the magnetostrictive effect of the iron core, vibrations caused by the electromagnetic force of the windings, mechanical vibrations caused by accessories such as cooling fans, and background noise, directly using the mixed signal for inversion would introduce a large error. Therefore, accurately separating the magnetostrictive vibration component in advance is a prerequisite for ensuring the accuracy of subsequent inversion.

[0027] like Figure 3 As shown, the separation step specifically includes: Step 1.1: Based on the instantaneous magnetization intensity extracted from the electrical operation data, and according to the physical principle that vibration is proportional to the square of magnetic flux density in the magnetostrictive effect, physical constraints are constructed to guide signal separation. Specifically, the dq-axis current component i obtained in Example 1 is used. d (t) and i q (t), from which the instantaneous magnetization M(t) = sqrt(i) can be calculated, which is approximately proportional to the magnetic flux density B. d (t) 2 +i q (t) 2 ) / N turns ;where N turns Let be the number of turns in the transformer winding, and be a known constant. Due to magnetostrictive vibration v... ms(t) and the square of the magnetic flux density B(t) 2 Proportional, i.e., v ms (t)∝M(t) 2 This physical relationship constitutes a strong constraint. The introduction of this constraint transforms the blind source separation problem into an optimization problem guided by prior physical knowledge, avoiding the misclassification of other vibration sources as target signals by traditional blind source separation algorithms and improving the accuracy of separation.

[0028] Step 1.2 involves applying physical constraints to the independent component analysis algorithm to perform constraint optimization and separation of the multi-point vibration signal, extracting the magnetostrictive vibration components. Specifically, this constraint optimization and separation process can be implemented using a modified constrained independent component analysis (cICA) algorithm. The objective function of this algorithm is no longer maximizing the non-Gaussianity of standard ICA, but rather optimizing an independence metric function embedded with physical constraints. J(W)=E{log(cosh(W T *V filtered ))}-E{log(cosh(W T *G*V filtered ))};where W is the separation matrix to be solved; V filtered is the preprocessed vibration signal matrix; E{} represents the expectation operation; G is the physical constraint matrix constructed according to step 1.1.

[0029] In the iterative solution of the separation matrix W, an adaptively decaying learning rate μ = 0.01 / (1 + 0.1k) (where k is the iteration number) can be used to balance convergence speed and stability. After separating multiple independent signal components, the kurtosis κ of each component is calculated, and components with kurtosis values ​​greater than a preset threshold (e.g., κ > 2, characterizing the super-Gaussian distribution) are identified as the desired magnetostrictive vibration components. Optionally, in addition to constrained independent component analysis, other constrained signal processing methods, such as Kalman filtering or particle filtering methods based on physical model constraints, can also be used to perform state estimation and separation of the mixed vibration signals.

[0030] Step 1.3 uses the magnetostrictive vibration components as the input signal for the subsequent inversion and reconstruction step. In other words, the output of this step is a vibration signal matrix V containing only magnetostrictive information. ms (t) serves as the observation data for the Bayesian inversion process in stage two.

[0031] Furthermore, before performing constraint optimization separation, to improve the separation effect, the original multi-point vibration signal can be preprocessed to filter out strong interference from the power grid frequency and its harmonics. A preferred approach is to use an adaptive notch filter array. Unlike a notch filter with fixed parameters, this adaptive notch filter can adjust the amplitude A of each harmonic based on real-time spectrum analysis.n Dynamically adjust the quality factor Q of the notch filter. n For example, through the formula: Q n =35*(1+0.1*log(A n / A noise Adjustments will be made to A; noise This is an estimate of the background noise level. Adaptive adjustment allows the notch filter to accurately filter out strong harmonic interference while preserving as much useful signal detail as possible in the frequency band near the harmonics, providing a higher signal-to-noise ratio input for subsequent cICA algorithms.

[0032] Phase Two involves performing Bayesian flux inversion based on physical constraints.

[0033] In obtaining the magnetostrictive vibration component V ms After (t), the step of inverting and reconstructing the three-dimensional magnetic flux density distribution field of the transformer's internal magnetic field state is performed. Specifically, this is a Bayesian inversion process based on physical constraints, such as... Figure 2 As shown, this process places the inverse problem of vibration-magnetic flux inversion under a probabilistic reasoning framework for solution, in order to integrate physical laws to constrain the solution space and obtain a stable solution with true physical meaning.

[0034] The inversion process specifically includes: Step 2.1: Based on Maxwell's equations and the saturation magnetic flux density constraint of the core material, a physical prior probability distribution is constructed to constrain the solution space. In this embodiment, this composite physical prior probability distribution is constructed to address the non-uniqueness of the solution caused by the number of sensors being far less than the number of grid points in the magnetic field to be solved. These physical constraints limit the solution space to a range that conforms to the fundamental physical laws of electromagnetic fields, improving the stability and physical accuracy of the inversion.

[0035] Constructing the physical prior probability distribution for constraining the solution space involves generating a probability distribution composed of the following components: a magnetic flux divergence constraint based on finite difference operators, characterizing the continuity of magnetic flux, i.e., mathematically requiring ▽•B to approach 0; a magnetic flux curl constraint based on Ampere's law and the current density derived from electrical operating data, characterizing the curl of the magnetic field generated by the current distribution, i.e., mathematically requiring ▽×B to match the current density J; and a saturation soft constraint defined by a probability function, used to characterize the nonlinear magnetization saturation characteristics of the core material, for example, through a soft penalty function p(B)∝exp(-α*max(0,|B|-B). sat ) 2 This is achieved by ) where |B| is the magnitude of the magnetic flux density; B satα is the saturation magnetic flux density of the core material, which can be taken as 1.7T for commonly used silicon steel sheets; α is a coefficient that controls the constraint strength, for example, it can be taken as 100; and the spatial smoothness constraint established by the Markov random field model characterizes the characteristic that the magnetic flux density should have a continuous and gradual change in space, so as to avoid abrupt changes or oscillations that do not conform to physical laws.

[0036] Step 2.2 initializes particles representing multiple candidate three-dimensional magnetic flux density distribution fields, forming an initial particle set. Specifically, this process generates N p Particles (e.g., N) p =1000), each particle B i 0 (i is the particle index, and the superscript 0 indicates the initial time) are all complete three-dimensional vector fields, representing a possible magnetic flux density distribution. These initial particles can be sampled from a truncated Gaussian distribution with the magnetic flux distribution under rated conditions as the mean and satisfying saturation constraints.

[0037] Step 2.3: Based on the multi-point vibration signals, an observation likelihood function describing the relationship between vibration and magnetic flux is established. Combined with the physical prior probability distribution, the weights of each particle in the initial particle set are iteratively updated to obtain the updated particle set. The observation likelihood function p(V ms |B i k (k is the number of iterations) is a physical model V based on the magnetostrictive effect. ms ∝B 2 Constructed, quantified for a given candidate magnetic flux distribution B i k Under these circumstances, the actual vibration signal V was observed. ms The probability of the new weight w. The iterative formula for weight updates follows Bayes' theorem, that is, the probability of the new weight w. i k Proportional to the old weight w i k-1 Observational likelihood p(V) ms k |B i k ) and the physical prior probability p(B) constructed in step 2.1 i k The product of ) . During the iteration process, to avoid particle degeneration (i.e., a few particles having excessively large weights), the effective number of particles N can be monitored in real time. eff =1 / Σ(w i 2 When N eff Below a preset threshold (e.g., N) p When / 3), the resampling step is triggered.

[0038] Step 2.4: Extract the maximum a posteriori (MAP) estimate from the weighted probability distribution of the updated particle set and define this estimate as a three-dimensional magnetic flux density distribution field for subsequent magnetic flux jump event prediction. Specifically, this MAP estimate can be obtained by calculating the weighted mean of all particles after the last iteration: B MAP =Σ(w i K *B i K ); where K is the total number of iterations. B MAP This is the final output of the three-dimensional magnetic flux density distribution field B(x,y,z,t) in this embodiment. Furthermore, to evaluate the reliability of the inversion results, this embodiment can also calculate the posterior covariance matrix from the updated particle set to obtain the uncertainty or variance of the magnetic flux density estimate at each spatial grid point. Based on this uncertainty, a three-dimensional confidence map can be generated. For regions with low confidence, spatial statistical methods such as Kriging interpolation can be used for post-processing correction to improve the smoothness and accuracy of the final output magnetic flux density distribution field.

[0039] Example 4 provides a preferred embodiment of dual-mode adaptive flux jump prediction, describing how to dynamically switch prediction modes according to power grid operating conditions to balance the speed and accuracy of prediction, and provide high-quality early warning information for subsequent millisecond-level suppression control.

[0040] In this embodiment, the step of predicting potential future magnetic flux jump events and generating a magnetic flux jump early warning vector is specifically a dual-mode adaptive prediction process, which includes: Step 1: Monitor and judge power flow switching events.

[0041] Real-time monitoring of the power flow status of the port area power grid, characterized by electrical operation data, is used to determine whether a power flow switching event has occurred. In this embodiment, a two-factor dynamic adaptive switching threshold is designed to make the switching decision of the prediction mode more intelligent, rather than using a fixed, empirical threshold. This design comprehensively considers the severity of macroscopic disturbances in the power grid (current change rate) and the microscopic vulnerability of the transformer (saturation risk). It can switch to a fast prediction mode in a timely manner when the system truly needs a rapid response, and maintain an accurate prediction mode during stable system operation to save computing resources and provide a high-precision benchmark. The judgment steps specifically include: The real-time comprehensive current change rate is calculated from the electrical operation data. Specifically, the dq-axis current component i obtained in Example 1 can be used. d (t) and i q To improve the accuracy of numerical differentiation, the instantaneous rate of change di can be calculated using the five-point difference method. d / dt and di q / dt, and further calculate the comprehensive current change rate |di / dt|=sqrt((di d / dt) 2 +(di q / dt) 2 ); where |di / dt| is the comprehensive current change rate, in amperes per second (A / s), which reflects the intensity of power flow fluctuations in the power grid.

[0042] The maximum value of the local saturation index, characterizing the current saturation risk of the transformer, is calculated from the three-dimensional magnetic flux density distribution field. This calculation specifically includes: obtaining the three-dimensional magnetic flux density distribution field characterizing the current magnetic field state (obtained by the method in Example 3), and calculating its first-order time derivative at each spatial grid point; for each spatial grid point, multiplying the ratio of its magnetic flux density amplitude to the core saturation magnetic flux density with a dynamic exponential function related to the time derivative of the magnetic flux density at that point to obtain the local saturation index for that point; iterating through the local saturation indices of all spatial grid points, and determining the maximum value as the maximum local saturation index. A preferred calculation formula is: SI(x,y,z)=|B(x,y,z)| / B sat *exp(k*|dB / dt|); where SI(x,y,z) is the local saturation index at the spatial point (x,y,z), which is a dimensionless value; |B(x,y,z)| is the magnetic flux density amplitude at that point; B sat Let be the saturation magnetic flux density of the core material, such as 1.7T; k is the dynamic influence coefficient, for example, 0.15; |dB / dt| is the absolute value of the time-varying rate of change of the magnetic flux density amplitude at that point. The calculated maximum local saturation index SI is... max It characterizes the risk level of the most vulnerable and closest point to saturation point in the entire transformer core.

[0043] An adaptive switching threshold is dynamically generated based on the statistical frequency of historical trend switching events and the maximum value of the local saturation index. For example, this adaptive switching threshold is thresh. adapt It can be determined by the base threshold thresh base (e.g., 500A / s) can be dynamically adjusted. On one hand, the system can statistically analyze the switching frequency f over a past period (e.g., 60 seconds). switch When f switch A higher threshold (e.g., greater than 0.5Hz) indicates that the power grid is in a state of continuous fluctuation, and the threshold can be appropriately lowered (e.g., thresh). adapt =0.8*thresh base To improve sensitivity; when f switch When the threshold is low, it can be increased to prevent false triggering. On the other hand, this threshold is also based on the maximum local saturation index SI. maxMake real-time adjustments, for example, using the formula thresh adapt =thresh base *(2-SI max This formula makes SI max As the transformer increases and approaches saturation, the switching threshold will decrease accordingly, making it more sensitive to power flow fluctuations when the transformer is more vulnerable.

[0044] The comprehensive current change rate is compared with the adaptive switching threshold in real time, and a final judgment is made on whether a power flow switching event has occurred based on the comparison result. Specifically, if |di / dt|>thresh adapt If the event occurs, a power flow switching event is determined to have occurred; otherwise, no power flow switching event is determined to have occurred.

[0045] Step two: Execute the corresponding prediction mode based on the judgment result.

[0046] Step 2.1: When it is determined that no power flow switching event has occurred, the system operates in precise prediction mode. In this mode, based on the three-dimensional magnetic flux density distribution field, the first prediction algorithm is used to perform refined evolution prediction over a longer time scale. The first prediction algorithm used in precise prediction mode is specifically implemented as follows: based on the current value B(x,y,z,t) of the three-dimensional magnetic flux density distribution field as the initial condition, a magnetic diffusion physical equation describing the diffusion process of the magnetic field in the conductive medium is constructed; the physical space of the transformer is discretized using the finite element method, for example, dividing the iron core into thousands of tetrahedral elements; and the Crank-Nicholson method (CN method) time-stepping scheme, which combines stability and accuracy in numerical calculations, is applied to numerically solve the magnetic diffusion physical equation to obtain a refined evolution prediction result of the three-dimensional magnetic flux density distribution field over a longer time scale (e.g., the next 160ms). The magnetic diffusion physical equation can be expressed as: dB / dt=D m *▽ 2 B-σ*v×(▽×B); where dB / dt is the time derivative of the magnetic flux density; D m The magnetic diffusion coefficient is 1 / (μσ), where μ is the magnetic permeability and σ is the electrical conductivity; ▽ 2 is the Laplace operator; v is the velocity field, which can be used to characterize effects such as rotating magnetic fields. This model has a large computational load but high prediction accuracy. Its results can be used not only for early warning under normal operating conditions but also as a benchmark for the following fast models.

[0047] Step 2.2: When a power flow switching event is detected, the system switches to a fast prediction mode. In this mode, based on the three-dimensional magnetic flux density distribution field, a second prediction algorithm is used to perform incremental fast prediction on a shorter time scale. The steps for incremental fast prediction in fast prediction mode include: acquiring the three-dimensional magnetic flux density distribution field at the current moment and the most recent few sampling moments (e.g., the most recent 5 sampling points) to form a time series dataset; based on the time series dataset, using a high-order finite difference method (such as the five-point difference method) to ensure accuracy, calculating the first time derivative dB / dt and the second time derivative dt of the three-dimensional magnetic flux density distribution field at the current moment. 2 B / dt 2 A second-order Taylor expansion prediction model is applied. Based on the current magnetic flux density distribution field B(t) and the calculated first and second-order time derivatives, the magnetic flux evolution within a short future time window (e.g., 5ms to 20ms) is calculated and used as the output of the fast prediction mode. The second-order Taylor expansion prediction model is: B(t+Δt)=B(t)+(dB / dt)Δt+0.5(d 2 B / dt 2 )Δt 2 Where B(t+Δt) is the predicted value of the magnetic flux density at the future time Δt. This model has low computational cost and fast response speed, and can meet the emergency prediction needs when there are severe fluctuations in the power flow.

[0048] Furthermore, to improve the efficiency and accuracy of the fast prediction mode, the mode can be calculated incrementally. When the system is running in accurate prediction mode, it can periodically (e.g., every 160ms) maintain a historical reference flux library. This library can use methods such as Principal Component Analysis (PCA) to compress the magnetic flux field obtained from accurate predictions, forming a reference descriptor which is then stored. When switching to fast mode, the system only needs to calculate the increment ΔB of the current magnetic flux field relative to the latest reference, and only perform Taylor expansion prediction on this increment, then superimpose the predicted increment back onto the reference. This incremental calculation method reduces the computational load and is one of the techniques for achieving millisecond-level fast prediction.

[0049] As an alternative implementation, the second prediction algorithm, in addition to the Taylor expansion model, can also employ a data-driven lightweight time series prediction model, such as an autoregressive integral moving average model or a simplified recurrent neural network variant (such as GRU). These models can learn dynamic patterns from historical flux evolution data and achieve rapid prediction during the inference phase.

[0050] Step 3: Output the warning vector. The prediction results under the current operating mode, whether in precise or rapid prediction mode, are formatted into a unified magnetic flux jump warning vector W for use by the subsequent cooperative suppression control module.

[0051] Example 5 provides a preferred embodiment of hierarchical collaborative suppression control, describing how to decouple a single suppression control problem into three-layer collaborative control commands executed at different time scales for different physical processes, thereby achieving accurate, effective, and economical suppression of magnetic flux jumps.

[0052] In this embodiment, the need to suppress flux jump events has different urgency in time, and the responsiveness of the port area's flexible loads (such as response speed, adjustment cost, and adjustment direction) also varies. By decoupling the control tasks in the time dimension, the most suitable load resources can be matched for different control objectives, achieving overall optimization.

[0053] like Figure 4 As shown, the steps for generating and issuing hierarchical collaborative control instructions include: The expected time of the flux jump is extracted from the early warning vector. Based on the time of the jump, multiple control nodes are set on the time axis, and control instructions with specific physical objectives are generated for different nodes: at the first and second preset time points before the time of the jump, a first-level control instruction for advanced phase compensation and a second-level control instruction for key point reverse suppression are generated in sequence; at the third preset time point after the time of the jump, a third-level control instruction for system power optimization and redistribution is generated; the first, second, and third-level control instructions are integrated to form a hierarchical collaborative control instruction and issued to the corresponding flexible load.

[0054] It is understandable that parsing the expected jump time from the magnetic flux jump warning vector and constructing hierarchical collaborative control commands based on the jump time is a time-decoupled generation process. Specifically, the key timestamp information t is extracted from the warning vector W output in Example 4. jump The t jump It will serve as a unified time reference for the generation and execution of subsequent three-layer control instructions, ensuring precise timing coordination of control actions at each layer.

[0055] The generation process of this time-series decoupling specifically includes: Before the expected transition occurs, a first-level control command for advanced phase compensation is generated. This first-level control, by actively injecting reactive power into the advanced phase, can pre-establish a magnetomotive force on the grid side that is opposite to the trend of the impending magnetic flux change, thus acting as a pre-cancellation and gently suppressing the initial development of magnetic flux disturbances. This step preferably dispatches loads with the fastest dispatch response speed and four-quadrant reactive power regulation capability, such as the energy storage system (ESS) described in Example 1. The execution time is typically set at t. jump An earlier window, such as t jump -5ms.

[0056] The steps for generating the first-level control command for lead phase compensation include: obtaining the expected jump amplitude ΔB from the flux jump warning vector. max Based on the jump amplitude, determine the reactive power compensation amount to offset part of the magnetic flux change; set an injection timing with a specific leading phase for the reactive power compensation amount so that the reactive power injection leads the magnetic flux disturbance in phase; generate a reactive power command containing the reactive power compensation amount and the injection timing as the first-level control command.

[0057] A preferred method for calculating reactive power compensation is: Q1(t) = K p *ΔB max *sin(ωt+π / 4+φ0); where Q1(t) is the instantaneous reactive power command that varies with time, in kilovars (kVar); K p This is the proportional gain coefficient, which can be set to, for example, 150 kVar / T. Its value can be calibrated according to the transformer parameters; ΔB max ω is the predicted maximum flux jump amplitude obtained from the warning vector; ω is the electrical angular frequency of the power grid; π / 4 is the key lead phase angle, ensuring that the reactive power injection effect leads the peak value of the flux disturbance in time; φ0 is the excitation impedance angle of the transformer itself, which can be obtained by arctan(X). m *ω / R m ) is calculated to obtain, where X m and R m These are the magnetizing reactance and magnetizing resistance of the transformer, respectively.

[0058] For a specific example, based on ΔB in the magnetic flux jump early warning vector max and t jump , in t jump Lead compensation is initiated at -5ms. The reactive power demand corresponding to the change in magnetic flux is calculated: Q. required =V 2 ·ΔB max / (X m ·B rated ), where V is voltage, X m For the magnetizing reactance. Considering the phase relationship between reactive power and magnetic flux, the reactive power compensation amount with a phase lead of π / 4 is: Q1(t) = Q1(t) = K p *ΔB max *sin(ωt+π / 4+φ0), where K p =150kVar / T is the proportional gain coefficient, φ0=arctan(X m *ω / R m The excitation impedance angle is denoted by ). This is based on the reactive power capacity Q of the energy storage system. ess_max Limit the amplitude: Q 1_final =min(Q1,0.8*Qess_max Generate PWM modulation parameters: m a =Q 1_final / (sqrt(3)*V dc *I rated ), outputting reactive power command u1={Q 1_final ,ma,φ inject}

[0059] Furthermore, before the expected flux jump occurs, a second layer of control commands is generated for critical point reverse suppression. When the flux jump is severe and the first layer's gentle compensation is insufficient to suppress it, a portion of the active load is rapidly disconnected, instantly reducing the total current flowing through the transformer and directly weakening the saturating magnetomotive force at its source. This suppression method is more direct and effective in dealing with severe flux jump events. The execution time is usually set immediately adjacent to t. jump At a time, for example t jump -2ms.

[0060] The steps for generating second-layer control instructions for keypoint inverse suppression include: Step 3.1: Identify the expected saturation danger zone inside the transformer from the magnetic flux jump warning vector. In this embodiment, this identification step is achieved by constructing a comprehensive risk index distribution map, specifically including: obtaining the local saturation index distribution reflecting the current saturation state (calculated by the method in Embodiment 4), the jump probability distribution parsed from the magnetic flux jump warning vector, and the real-time hot spot temperature distribution characterizing the thermal stress of the equipment; calculating the comprehensive risk index for each spatial point by weighted fusion of the values ​​of the local saturation index, jump probability, and hot spot temperature; and determining the set of spatial points whose index values ​​exceed a preset danger threshold (e.g., 0.85) as the expected saturation danger zone based on the distribution of the comprehensive risk index. The formula for calculating the comprehensive risk index can be: R(x,y,z)=α*SI(x,y,z)+β*P jump (x,y,z)+γ*T thermal (x,y,z); where R(x,y,z) is the comprehensive risk index; SI is the local saturation index; P jump T is the probability of a jump; thermal The normalized hotspot temperature; α, β, and γ are weighting coefficients, which can be 0.5, 0.3, and 0.2 respectively, and α+β+γ=1.

[0061] Step 3.2: For saturation-hazardous areas, a reverse magnetic flux distribution is designed in the magnetic field domain to counteract the saturation effect. Based on the principle of energy equivalence conversion, the magnetic field energy required to maintain this reverse magnetic flux distribution is converted into the active power regulation amount to be executed in the power grid domain. Based on the active power regulation amount, loads with rapid disconnection capabilities (such as some charging piles) are selected from the flexible load operation status of the port area, and a second-level control command is generated. Specifically, the reverse magnetic flux distribution B... counter It can be designed as a spatial Gaussian distribution centered on the transition center, with amplitudes opposite to the predicted transition amplitudes. The magnetic field energy E required to maintain this reverse magnetic flux... risk The relationship between P2 and the active power P2 that needs to be adjusted can be achieved through energy conservation and power relationships. For example, P2 is proportional to E. risk Switching frequency f with power grid flow switch The product of.

[0062] Specifically, read the set of saturated danger zones Ω risk Calculate the magnetic flux energy E in this region. risk =i Ω_risk B 2 / (2μ)dV. Design the reverse magnetic flux distribution B. counter =-0.7*ΔB predicted *exp(-r 2 / r0 2 (), where r is the distance to the transition center, and r0 = 0.1m is the radius of influence. Based on the principle of energy equivalence, the required active power adjustment is calculated: P2 = 2 * E risk *f switch , where f switch The frequency of power flow switching is decomposed into a fast response component P. fast (accounting for 70%) and slow response component P slow (30%). Identify loads that can be quickly disconnected: Select the set L of loads with a response time <10ms from the load response capability matrix. fast Generate a rapid load cutoff command u2={L fast ,P fast ,t cut =t jump -2ms}.

[0063] Furthermore, after the expected transition occurs, a third-level control command is generated for optimized power redistribution in the system. This third-level control addresses the aftereffects of the second-level fast control, namely, the system power imbalance. This occurs after the flux transition shock has passed (e.g., at t...). jumpWith a +5ms startup, the output of other adjustable loads in the port area (including slower-responding quay cranes and yard cranes) is reallocated in a way that minimizes total cost and ensures the fairest adjustment of the load, in order to make up for the power cut off in the second layer and restore the system to a stable operating state.

[0064] The steps for generating third-level control instructions for system power optimization and reallocation include: setting the active power adjustment P2 determined by the second-level control instructions as the system's total power balance target (i.e., ΣΔP) for the current power reallocation process. i =-P2); Based on the cost weights and operating limits contained in the flexible load operation status of the port area, the objective function of the multi-objective optimization scheduling model is constructed. The objective function of the optimization scheduling model is to minimize the weighted adjustment cost of each load and balance the relative adjustment amplitude of each load. Under the triple constraints of satisfying the total power balance target of the system, the adjustable capacity of a single load and the ramp rate, the optimization scheduling model is solved to obtain the optimal power adjustment sequence, which constitutes the third layer of control commands.

[0065] The objective function of a preferred multi-objective optimization scheduling model is: minJ=Σ(w i *|ΔP i |)+λ*Σ(ΔP i / P i_rated ) 2 Where J is the comprehensive cost function to be minimized; ΔP i w is the power adjustment variable for the i-th flexible load; i The unit adjustment cost weight for this load reflects its economic efficiency; P i_rated λ is the rated power of the load; λ is the equalization coefficient (e.g., 0.1), used to penalize behavior where the regulation amplitude is too large relative to its rated capacity, so as to achieve fair allocation of regulation tasks. Constraints: (1) Power balance ΣΔP i =-P2;(2) Capacity constraint|ΔP i |ΔP i (3) Climbing constraint |dP i / dt|.

[0066] The optimization problem can be solved using efficient algorithms such as interior point method and sequential quadratic programming. Preferably, the hierarchical dimensionality reduction optimization algorithm described in Example 7 can also be used.

[0067] Optionally, before integrating and issuing the three-layer hierarchical collaborative control instructions, a timing coordination and conflict resolution step is also included. This step ensures that the logically separate three-layer control instructions can be executed in a coordinated and conflict-free manner in terms of physical execution, which is a key engineering step to ensure the stable and reliable execution of the entire collaborative control strategy.

[0068] This step specifically includes: detecting whether there are timing conflicts among the first, second, and third layer control commands that target the same flexible load (especially the energy storage system) and whose execution time windows overlap; when a timing conflict is detected, adjusting the execution timing of the conflicting commands according to a preset priority rule, i.e., the leading phase compensation command has the highest priority, followed by the critical point reverse suppression command, and the system power optimization and redistribution command has the lowest priority; the adjustment includes at least delaying the execution time of low-priority commands to avoid peak loads, or superimposing the amplitudes of multiple commands targeting the same load (for example, if the energy storage system receives reactive power commands and active power commands at the same time, its total apparent power limit needs to be considered comprehensively), in order to generate a final coordinated and conflict-free hierarchical collaborative control command.

[0069] Example 6 provides a preferred embodiment of closed-loop correction of control effect, which introduces a feedback correction mechanism to address the control effect deviation caused by factors such as model inaccuracy, load response uncertainty or unexpected power grid disturbances, so that the final magnetic flux suppression effect is always maintained within the expected target range.

[0070] After the port area's flexible loads are dispatched and the hierarchical coordinated control commands are executed, a closed-loop correction step is also included, which specifically includes: Step one: Real-time acquisition of measured magnetic flux changes after command execution. Specifically, within a very short time window, such as 2ms, after the execution of the three-layer control commands in Example 5 (especially the most powerful second-layer command), the system needs to quickly evaluate the actual effect of the control action. This measured magnetic flux change ΔB measured The magnetic flux can be obtained in several ways. One way is to directly measure the change in magnetic flux by using a detection coil pre-deployed at a specific location in the transformer core. Another non-invasive method is to indirectly estimate the actual change in magnetic flux inside the core based on the signal changes of the leakage flux sensor outside the transformer, combined with a pre-calibrated leakage flux-main flux mapping model.

[0071] Step two involves comparing the measured magnetic flux change with the expected jump amplitude contained in the magnetic flux jump warning vector to generate a control deviation signal that quantifies the control effect. Here, the benchmark for comparison is not the original, uncontrolled expected jump amplitude ΔB. predicted Instead, it is the desired, suppressed target value of magnetic flux change ΔB. expected The target value can be the original predicted value multiplied by the expected inhibition rate (e.g., if the target inhibition rate is 80%, then ΔB). expected For ΔB predicted (20% of the total). The control deviation signal e can be obtained through e=|ΔB measured |-|ΔB expected The calculation is performed using the formula. If e > 0, it indicates insufficient inhibition; if e < 0, over-inhibition may exist.

[0072] Step 3: When the amplitude of the control deviation signal exceeds the preset correction threshold, a supplementary power adjustment command is generated based on the deviation signal. The correction threshold can be set according to the transformer's safe operating margin, for example, 0.1T. When the absolute value of |e| exceeds 0.1T, the system determines that closed-loop correction needs to be initiated. At this time, the system generates a supplementary active power adjustment command ΔP based on the deviation signal e. ess A preferred method for generating instructions is to use an integral control law to eliminate potential steady-state errors. The calculation formula is: ΔP ess =K i *Ks is preferred; where K is the preferred choice. i This is the integral gain coefficient, and its unit can be kW / (T·s), for example, it can be 500. This coefficient value determines the speed and strength of the correction response and can be tuned according to the dynamic characteristics of the system.

[0073] Step four: The unit with the highest dynamic response capability among the flexible loads is scheduled to execute a supplementary power adjustment command to quickly compensate for control deviations. Specifically, this supplementary power adjustment command ΔP ess The request will be prioritized for execution by the fastest-responding energy storage system (ESS) in Example 1. The energy storage system can quickly adjust its active power output within milliseconds. If e>0 (insufficient suppression), it will increase the discharge power or decrease the charging power; if e<0 (oversuppression), it will decrease the discharge power or increase the charging power, thus quickly and accurately compensating for control deviations.

[0074] This closed-loop correction step enhances the robustness and adaptability of the entire suppression system by constructing a complete prediction-control-feedback-recorrection closed-loop system.

[0075] Example 7 provides a preferred embodiment for solving the optimized scheduling model, which adopts a hierarchical dimensionality reduction optimization strategy to avoid the problems of dimensionality curse and excessive computation time that may be encountered when directly solving high-dimensional, non-convex multi-objective optimization problems.

[0076] The steps to solve the optimal scheduling model are specifically a hierarchical dimensionality reduction optimization process, which includes: Step 1: Analyze the Hessian matrix of the objective function of the optimization scheduling model relative to each flexible load power adjustment variable, and identify its sparsity pattern to quantify the coupling strength between variables. Specifically, the objective function J (as in Example 5) is related to each flexible load power adjustment variable ΔP1, ΔP2, ..., ΔP N The function is J. The Hessian matrix is ​​an N×N square matrix, where the element in the i-th row and j-th column is determined by the objective function J with respect to ΔP. i and ΔP jThe second-order mixed partial derivative d 2 J / (dΔP i dΔP j ) consists. In this embodiment, the absolute value magnitude of this second-order mixed partial derivative reflects the mutual influence or coupling strength between the power adjustments of load i and load j during the optimization process. For example, if two loads are on the same feeder, their power adjustments may mutually affect the line voltage or losses, which is reflected as a strong coupling relationship in the objective function, corresponding to a relatively large non-diagonal element value in the Hessian matrix.

[0077] Step 2: Based on the coupling strength, apply the spectral clustering algorithm to automatically group the power adjustment variables, and decompose the original high-dimensional optimization problem into multiple associated low-dimensional subproblems. Specifically, the spectral clustering algorithm takes the Hessian matrix (Hessian matrix) composed of the coupling strength or its derived association matrix as input. By analyzing its eigenvectors, it automatically divides the variables with strong coupling (i.e., loads with high physical or electrical correlation) into the same group and divides the variables with sparse coupling into different groups at the level of graph theory. Through this step, the original high-dimensional (N-dimensional) optimization problem involving all N loads and difficult to directly solve is decomposed into M (M < N) subproblems with lower dimensions, strong internal coupling, and weak external coupling. Each subproblem corresponds to a load group.

[0078] Step 3: Perform recursive optimization in multiple optimization levels composed of low-dimensional subproblems and fuse the solutions of each level to obtain the optimal power adjustment sequence. This recursive optimization process can be constructed as a multi-level optimization architecture. At the top level, a rough and decoupled optimization can be performed on each subproblem to quickly obtain an initial solution. Take the solution at the top level as the constraint or initial value at the next level, and consider some key cross-group coupling relationships at the next level for more refined optimization. This process can be carried out recursively, approaching the global optimal solution layer by layer. Through methods such as weighted average or Pareto front selection, the optimal solutions of each level or each subproblem are fused to form a globally consistent optimal power adjustment sequence.

[0079] This hierarchical dimension reduction optimization process, by first grouping and then stratifying, decomposes large-scale complex problems into multiple small problems that are easy to handle. It can control the calculation time within milliseconds while ensuring the solution quality, meeting the rapidity requirements of the third-layer power reallocation.

[0080] As an optional implementation method or when the load scale is small, the optimal scheduling model can also be solved using traditional centralized optimization algorithms, such as the interior point method or sequential quadratic programming. These methods can theoretically find the exact optimal solution, but when the number of flexible loads in the port area is huge and the model dimension is very high, their computational efficiency may be lower than the hierarchical dimension reduction optimization method of this embodiment.

[0081] Example 8 provides an alternative multiphysics fusion feature extraction embodiment, describing another implementation method for the magnetic flux spatial distribution inversion step. Specifically, this embodiment abandons the attempt to accurately reconstruct the computationally intensive and extremely difficult three-dimensional physical field. Instead, it directly extracts feature indicators strongly correlated with magnetic saturation by fusing key information from multiple physical fields such as vibration, heat, magnetism, and electricity generated by the transformer during operation. This achieves perception of the transformer's internal state with lower computational cost and higher robustness.

[0082] In this embodiment, the aforementioned magnetic flux spatial distribution inversion step—acquiring multi-point vibration signals and electrical operation data outside the transformer, and reconstructing a three-dimensional magnetic flux density distribution field characterizing the internal magnetic field state of the transformer based on the multi-point vibration signals and electrical operation data—is replaced by a step of extracting multi-physics field fusion features, specifically including: Step one: In addition to collecting multi-point vibration signals and electrical operation data, simultaneously collect distributed temperature signals and external leakage magnetic field signals of the transformer. Specifically, in addition to the vibration and electrical sensors in Example 1, this example adds: 1) distributed fiber optic temperature measurement points, for example, 16 temperature measurement points are arranged at key parts of the transformer windings and core to obtain a refined temperature field distribution T(x,y,z,t); 2) external magnetic flux probes (leakage magnetic field sensors), for example, 4 probes are non-invasively installed at specific locations outside the transformer tank to measure the external leakage magnetic field B intensified by the internal magnetic field saturation of the transformer. leak (t).

[0083] Step two involves performing bispectral analysis on the multi-point vibration signal to extract the vibration nonlinear characteristics representing second harmonic coupling caused by magnetic saturation. Bispectral analysis is a high-order spectral analysis method that can detect nonlinearity and second-phase coupling phenomena in signals. In this embodiment, when the transformer core enters the saturation region, the nonlinearity of the magnetostrictive effect intensifies, leading to second harmonic coupling in the vibration signal, i.e., a nonlinear energy transfer from the fundamental frequency to the second harmonic frequency. The bispectral S of the vibration signal is calculated... vv By taking (f1,f2) and integrating it over a specific coupling region (e.g., the region where f1-2f2 are fixedly coupled), the characteristic Ψ that quantifies the nonlinear coupling strength can be obtained. v This feature Ψ v It is a strong indicator of saturation because its theoretical value should be zero in a linear system.

[0084] Step 3: Extract the corresponding thermal, magnetic, and electrical saturation characterization quantities from the distributed temperature signal, external leakage flux signal, and electrical operation data, respectively. Specifically: Thermal saturation characterization quantity Ψ tIt can be the ratio of the Laplace operator of the temperature field to the rate of change of temperature over time, i.e., Ψ t =▽ 2 T / (dT / dt) reflects the abnormal hot spot diffusion mode caused by the sharp increase in eddy current loss due to saturation inside the iron core.

[0085] Magnetic saturation characterization quantity Ψ m It can be expressed as the product of the square of the leakage magnetic signal amplitude and the direction of its rate of change over time, i.e., Ψ. m =|B leak | 2 *sign(dB leak / dt), this index captures the nonlinear and asymmetric enhancement of the leakage magnetic field due to saturation.

[0086] Electrical saturation characterization quantity Ψ e The total harmonic distortion (THD) of the input current can be used as a metric. i The product of functions related to the power factor, i.e., Ψ e =THD i *(1+|cos(φ)|), because core saturation will cause severe distortion of the excitation current waveform and generate a large number of harmonics.

[0087] Step four involves fusing the vibrational nonlinear characteristics with thermal, magnetic, and electrical saturation parameters to generate a saturation feature vector. This saturation feature vector then replaces the three-dimensional magnetic flux density distribution field to perform subsequent magnetic flux transition event prediction. Specifically, this data fusion process can integrate the extracted feature quantities Ψ... v ,Ψ t ,Ψ m ,Ψ e Combined into a multidimensional saturated eigenvector Ψ(t) = [Ψ v ,Ψ t, Ψ m ,Ψ e The vector Ψ(t) no longer describes the complete magnetic flux physical field, but rather provides a condensed feature representation of the saturation state from multiple physical dimensions. Subsequent prediction modules (such as variations of Embodiment 4 or Embodiment 9) will directly use the time series of the feature vector Ψ(t) as input to predict magnetic flux jump events.

[0088] Example 9 provides a prediction method based on deep temporal graph networks, describing an alternative implementation of the dual-mode adaptive prediction step. Specifically, it utilizes cutting-edge deep learning models in artificial intelligence to automatically learn complex, nonlinear magnetic flux evolution patterns from massive historical data, achieving accurate prediction of magnetic flux jumps. This example addresses the issue of purely data-driven models potentially violating physical laws by introducing a physical constraint loss function.

[0089] In this embodiment, the steps of receiving a three-dimensional magnetic flux density distribution field and predicting possible future magnetic flux jump events based on its temporal evolution, generating a magnetic flux jump early warning vector containing jump time, amplitude, and location information, can also be a prediction method based on a deep temporal graph network, including: A weighted graph representing the physical magnetic circuit topology of a transformer is constructed. The nodes of the graph are defined as physical regions such as the transformer's core pillars and yokes, and the connection weights between nodes are assigned based on the magnetic reluctance between these physical regions. Specifically, the node set V of the graph G=(V,E,A) can contain several core physical regions (e.g., 12 nodes), such as the transformer's three core pillars and two yokes. The edge set E between nodes indicates that these regions are physically adjacent and exchange magnetic flux. The adjacency matrix A contains elements A... ij The value assigned is the ease with which magnetic flux can be exchanged between two regions, and can be proportional to the reciprocal of the magnetic reluctance.

[0090] Furthermore, the temporal evolution data of the three-dimensional magnetic flux density distribution field is input into a spatiotemporal graph convolutional neural network (ST-GCN) for spatiotemporal feature extraction on the weighted graph to obtain the future evolution trend of the magnetic flux. The ST-GCN is a deep learning model that can simultaneously handle dependencies in both time and space dimensions. In the spatial dimension, graph convolution operations (such as Chebyshev graph convolution) allow the state updates of each node (physical region) to aggregate information from its neighboring nodes, thus simulating the physical diffusion process of magnetic flux in the magnetic circuit topology. In the temporal dimension, structures such as gated recurrent units (GRUs) or causal dilation convolutions can be used to capture the temporal dependencies of magnetic flux evolution. The input data is obtained by averaging or integrating the three-dimensional magnetic flux field obtained in Example 3 within each graph node region to obtain the time-series signal of each node.

[0091] During the training phase, the spatiotemporal graph convolutional neural network is constrained by a loss function that incorporates physical constraints from Maxwell's equations to ensure the physical authenticity of its output. Specifically, the network's total loss function L... total Data-driven prediction error loss L data and physical constraint loss L physics Weighted composition: L total =L data +λ1*L physics Among them, L data This is the standard mean squared error loss. Physical constraint loss L physics Then the predicted magnetic flux field B output by the network will be... pred Substitute the values ​​into the discrete form of Maxwell's equations for calculation, for example, L... physics =||▽·Bpred || 2 +||▽×B pred -μJ|| 2 By minimizing L during backpropagation total This allows the purely data-driven model to follow the basic laws of electromagnetic field physics during the learning process, avoiding the generation of physically impossible predictions.

[0092] The network also employs a dual-branch prediction architecture. One branch predicts the smooth evolution of magnetic flux, while the other predicts the probability and magnitude of abrupt changes in magnetic flux. The outputs of the two branches are fused to generate a magnetic flux abrupt change warning vector. Specifically, after the spatiotemporal features are extracted by ST-GCN, these features are fed into two parallel prediction heads (output layers). The prediction head of the smooth evolution branch outputs the predicted magnetic flux value B for multiple future time steps. smooth (t+τ). The prediction head of the mutation prediction branch outputs the mutation probability P. jump (t+τ) and the magnitude of the mutation ΔB jump The final prediction result is obtained through B. pred =B smooth +P jump *ΔB jump The mixture is then fused. The fused result is then processed into the standard magnetic flux jump warning vector W defined in Example 2.

[0093] Example 10: A specific end-to-end numerical calculation case. This example provides an example that runs through the entire process of sensing-prediction-control-correction of this invention, including specific parameters and calculation steps, so that those skilled in the art can clearly understand how the various links of this invention work together.

[0094] This case study is set against the backdrop of the port area power distribution system described in Example 1. The core research object is a 20MVA dry-type transformer, whose core material has a saturation magnetic flux density B. sat The current is 1.7T. At a certain moment t0, a large gantry crane suddenly starts its main hoisting mechanism, causing the total current of the system to rise sharply.

[0095] The first step is to conduct state awareness and risk assessment. At time t0, the monitoring system deployed on the transformer continues to operate.

[0096] The magnetic flux inversion method is running in real time. Within milliseconds of the quay crane starting, the system inverts the current magnetic flux density amplitude |B(t0)| at a T-joint in the transformer core (assumed to be the weakest point (x0, y0, z0)) to be 1.6T, with a time change rate |dB / dt| as high as 5.0T / s. Based on the above data, the system immediately calculates the local saturation index SI at this point using the method of Example 4. Calculation process: SI(x0,y0,z0)=|B(t0)| / B sat *exp(k*|dB / dt|)=1.6 / 1.7*exp(0.15*5.0) Local saturation index of a certain iron core. The prediction head outputs the row average or product of the magnetic flux prediction values ​​for multiple future time steps.

[0097] At this point, the maximum local saturation index SI of the entire transformer is... max The value is 1.992, which is far beyond the danger threshold of 0.95, indicating that a severe dynamic magnetic saturation is about to occur at this point.

[0098] The second step involves making prediction mode switching and jump prediction switching decisions. Specifically, the system synchronously calculates the comprehensive current change rate |di / dt| caused by the quay crane startup, obtaining a value of 650 A / s. Simultaneously, the system adjusts the current based on the current SI... max The threshold for switching prediction modes is dynamically calculated.

[0099] Calculate the adaptive threshold: thresh adapt =thresh base *(2-SI max )=500A / s*(2-1.992)=500*0.008=4A / s.

[0100] Since the measured value of |di / dt| = 650 A / s is much greater than the adaptively adjusted switching threshold of 4 A / s, the system determines that a drastic power flow switching event has occurred. The system immediately switches from precise prediction mode to fast prediction mode.

[0101] The system retrieves the magnetic flux density values ​​at four sampling points up to and including time t0 (assuming the values ​​at point (x0, y0, z0) are [1.45, 1.49, 1.54, 1.58, 1.60]T respectively). It then calculates the current first-order time derivative dB / dt / dt using higher-order differences and the second-order time derivative dt / dt. 2 B / dt 2 ≈ / interval derivative 2 A second-order Taylor expansion model is applied for prediction, aiming to predict the magnetic flux density 18 ms later: B(t0+18ms)=B(t0)+(dB / dt)*0.018s+0.5*(d 2 B / dt 2 )*(0.018s) 2 ≈*(0.018s)0.5*(Predicting the future 1.60]T8=4A / s. The prediction head outputs the magnetic flux prediction value for multiple future time steps, either row average or...

[0102] Since the predicted value of 1.703T exceeds the saturation magnetic flux density of 1.7T, the system determines that a magnetic flux jump will occur. The system generates and broadcasts a magnetic flux jump warning vector to the downstream control module: W=[ΔB max =0.28T,t jump =18ms,x j =x0,y j =y0,z j =z0] T ; Among them, the predicted jump amplitude ΔB max This is the spatial maximum value obtained based on a more complex full-field prediction; in this example, it is 0.28T; the transition time t jump After being set to 18ms.

[0103] The third step is to generate and execute hierarchical collaborative control commands. Upon receiving the warning vector W, the control module immediately... jump =18ms is used as the time base to generate three-layer control commands.

[0104] First-level control (t=t) jump -5ms=13ms): Perform advanced reactive power compensation, specifically, according to the formula: Q1(t) = K p *ΔB max Substitute *sin(ωt+π / 4+φ0) into K p =150kVar / T, ΔB max =0.28T, the peak reactive power compensation is calculated to be approximately 42kVar. At 13ms, an instruction is sent to the 5MW / 10MWh energy storage system (ESS) in the port area, requiring it to inject 42kVar of advanced reactive power according to the calculated curve.

[0105] Second-level control (t=t) jump -2ms=16ms): Key point reverse suppression is performed. Specifically, based on the identified saturation danger zone (the spatial range centered on (x0,y0,z0)), the required active power regulation of P2=250kW is calculated through the principle of energy equivalence to effectively suppress saturation. At 16ms, a command is sent to a group of electric truck charging piles (total power 500kW) with fast response speed in the port area, requiring them to urgently cut off the 250kW charging load.

[0106] Third-level control (t=t) jump +5ms=23ms): Perform system power optimization and redistribution, specifically, rebalance the 250kW power deficit caused by the second-level control in the system. The system solves a multi-objective optimization problem. Assume the adjustable loads are: the slow-responding field bridge A (adjustment cost w)_ A = 0.5 yuan / kWh), refrigerated container load B (adjustment cost w) _ B = 0.8 yuan / kWh). The solution is: the power of the field bridge A is reduced by 150kW, and the power of the refrigerated container load B is reduced by 100kW. At 23ms, power adjustment commands are sent to both the field bridge A and the refrigerated container load B.

[0107] The fourth step is to perform closed-loop calibration. Approximately 2ms after all control commands have been executed, i.e., at 25ms, the system initiates a closed-loop calibration check.

[0108] The actual maximum magnetic flux jump ΔB at point (x0, y0, z0) after suppression, as measured by the detection coil, was obtained. measured It is 0.06T.

[0109] Assuming the system's target inhibition rate is 80%, the expected residual jump is: ΔB expected =ΔB max *(1-80%)=0.28T*0.2=0.056T.

[0110] Calculate the deviation e = |ΔB measured |-|ΔB expected |=0.06T-0.056T=0.004T.

[0111] Since the deviation e = 0.004T, which is less than the correction threshold of 0.1T, the system determines that the suppression control effect meets the standard and there is no need to start supplementary power adjustment.

[0112] Through the aforementioned series of precisely coordinated actions on a millisecond-level timescale, the power flow surge that was initially predicted to cause severe magnetic saturation was successfully suppressed, the transformer's magnetic flux was ultimately controlled within a safe range, and the stable operation of the port area's power grid was ensured. This case demonstrates the entire closed-loop workflow of this invention, from high-precision sensing and intelligent prediction to hierarchical collaborative control, verifying its technical feasibility and advancement.

[0113] Example 11: A preferred method for extracting multi-scale modal features of vibration signals. This example provides a deep signal processing method between two stages: separation of magnetostrictive vibration components and Bayesian flux inversion based on physical constraints. Preferably, the separated time-domain vibration signal V is not used directly. ms Instead of using (t) as the input to the inversion algorithm, it is first transformed to a modal feature space with denser information and more robust features to improve the stability and accuracy of the subsequent inversion process.

[0114] In one specific implementation, when the magnetostrictive vibration component V msBefore (t) (obtained by the method in stage one of Example 3) is used as the input signal for subsequent inversion steps, the following steps are also included: Step one involves performing multi-resolution analysis on each channel of the magnetostrictive vibration component to obtain its joint time-frequency distribution. Specifically, continuous wavelet transform (CWT) can be used to analyze the signal V. ms (t) is processed. Preferably, the Morlet wavelet basis ψ(t) = exp(-t) is used. 2 The formula ∂² / ∂t (i*ω₀*t) represents time, i is the sampling point number, and w₀ is the angular frequency. It exhibits good localization characteristics in both time and frequency, making it suitable for analyzing non-stationary vibration signals. By performing wavelet transforms on multiple scales 'a' (scale a is inversely proportional to frequency f), a two-dimensional time-frequency coefficient matrix CWT(a,b) can be obtained; where b is the time shift. This matrix accurately describes how the energy of each frequency component in the vibration signal evolves over time.

[0115] Step two: Based on the joint time-frequency distribution, identify the dominant vibration modes. Since the vibration of a transformer is mainly dominated by its inherent mechanical modes, the energy in the time-frequency distribution is also mainly concentrated in the frequency bands corresponding to these modes. To automatically and objectively identify these dominant modes, singular value decomposition (SVD) can be performed on the time-frequency coefficient matrix CWT obtained in step one. SVD can decompose this matrix into U*S*V... T The form is given by the matrix S, where the diagonal elements (singular values) of the diagonal matrix S directly correspond to the energy or importance of each modal component. The larger the singular value, the greater the contribution of the corresponding mode to the total vibration.

[0116] Step 3: Construct modal feature vectors based on the dominant vibration modes. Specifically, sort the singular values ​​obtained after SVD decomposition and select the modes corresponding to the top N largest singular values ​​as the dominant modes (for example, selecting N=6 dominant modes, whose total energy may account for more than 95% of the total vibration energy). For the N dominant modes, extract their instantaneous amplitude and instantaneous phase information at each moment. Combine the 2N feature quantities (N amplitudes and N phases) of the N modes to form a 2N-dimensional modal feature vector Φ(f,t).

[0117] Step four: Use the modal feature vector Φ(f,t) as the input signal for the inversion and reconstruction step. That is, in the subsequent Bayesian inversion algorithm of stage two in embodiment three, the observation equation will no longer be V. ms It is not the relationship between Φ(f,t) and B, but the relationship between Φ(f,t) and B.

[0118] In this embodiment, the introduction of this multi-scale modal feature extraction process has multiple effects. Specifically, by selecting the dominant mode through SVD, it is equivalent to performing energy-based intelligent filtering, effectively eliminating low-energy, non-dominant vibrational components caused by measurement noise, background interference, etc.; and transforming the high-dimensional time-domain signal V... ms Transforming (t) into a low-dimensional modal feature vector Φ(f,t) achieves dimensionality reduction and information compression of the data, which can reduce the computational complexity of subsequent inversion algorithms. On the other hand, using modes with more explicit physical meanings instead of the original vibration signals as inputs allows the inversion model to better capture the core physical correlation between vibration and magnetic field, improving the accuracy and robustness of the inversion.

[0119] Example 12: A dynamic reference library method for enhancing dual-mode prediction. This example provides a detailed systematic implementation scheme for enhancing the efficiency and accuracy of dual-mode adaptive prediction methods. It designs and maintains a dynamic historical reference flux library, which solves the problem of fast prediction modes' dependence on high-quality references and enables incremental fast prediction.

[0120] In one specific implementation, the dual-mode adaptive prediction process further includes the following steps: Step 1: When the system is running in accurate prediction mode, a historical reference flux library is periodically generated and maintained. This process specifically includes: generating and compressing the reference. Specifically, when the system is running in accurate prediction mode, it periodically (e.g., every 160ms) obtains a high-precision, complete three-dimensional flux density distribution field B(x,y,z,t). Considering that directly storing massive amounts of three-dimensional field data is impractical, this embodiment uses Principal Component Analysis (PCA) or Singular Value Decomposition (SVD) to compress the three-dimensional field data. Specifically, after vectorizing the three-dimensional flux field, SVD is performed, retaining only the top k principal components (e.g., k=20) that can capture the vast majority of energy (e.g., 95%), compressing the high-dimensional physical field into a low-dimensional reference descriptor D containing core information. ref The descriptor D ref Along with statistical information such as timestamps, mean, and variance, the data is stored in the benchmark database. Furthermore, the benchmark database is dynamically updated. Specifically, to avoid redundancy and maintain its timeliness, the system maintains a sliding window containing only the most recent M benchmarks (e.g., M=10). When a new benchmark is generated, the system calculates its similarity (sim) to the latest benchmark in the database (e.g., using vector cosine similarity). Only when sim is below a preset threshold (e.g., 0.85) is the operating condition considered changed, requiring the new benchmark to be added to the database (and the oldest benchmark removed); otherwise, the operating condition is considered stable, and no update is needed.

[0121] Step two: When the system switches to fast prediction mode, incremental calculations are performed based on the historical reference flux library. Specifically, after the switching decision in Example four occurs, the system retrieves the latest timestamp of the reference flux distribution B from the reference library. ref (This benchmark is achieved through the stored descriptor D) ref (As obtained from reconstruction). The system no longer performs a Taylor expansion on the complete magnetic flux field B(t), but instead calculates the increment of the current magnetic flux field relative to the reference, ΔB = B(t) - B. ref Since the magnetic flux change caused by the power flow switching is usually localized and sparsity-dependent, the increment ΔB is typically smaller and simpler than B(t) itself.

[0122] Step three involves performing prediction on the increment and applying accuracy compensation. A fast prediction algorithm (such as a second-order Taylor expansion) is applied to the increment ΔB, rather than the complete B(t), significantly reducing the computational cost. This yields the predicted increment ΔB. pred Then, it is superimposed back onto the reference B. ref The final magnetic flux prediction result B is obtained. pred =B ref +ΔB pred .

[0123] Furthermore, to compensate for the difference between the current time t and the reference time t ref The error introduced by the time difference can be compensated by introducing a phase or trend compensation term. A preferred compensation method is: ΔB compensated =ΔB+(tt ref )*dB ref / dt;ΔB compensated For trend compensation, dB ref / dt is the time rate of change of the reference magnetic flux field, which can be estimated by linear fitting to the two or three most recent references in the reference library.

[0124] The dynamic benchmark library mechanism described in this embodiment, through the collaborative working method of precise model library construction and fast model library usage, as well as a series of means such as compressed storage, incremental calculation, and trend compensation, provides a high-quality, near real-time reference benchmark for fast prediction models while ensuring rapid response, thus ensuring that the dual-model prediction system can have both speed and accuracy.

[0125] Example 13: A Preferred Implementation of Distributed Synchronous Execution of Control Commands. This example provides a preferred system-level implementation scheme for the execution stage after the control commands are generated. It describes the preferred technical architecture adopted by the present invention to achieve microsecond-level time synchronization and millisecond-level low-latency communication between multiple physically distributed flexible load execution units. In a specific system implementation, after the generation of hierarchical collaborative control commands or closed-loop correction commands, their efficient and accurate execution faces two major technical challenges: how to transmit the commands from the main controller to the load execution units distributed throughout the port area with extremely low and deterministic latency; and how to ensure that execution units in different locations execute their respective control actions at exactly the same time. To address these challenges, this example adopts the following steps: Step 1: Perform fast communication based on IEC61850 GOOSE messages.

[0126] To ensure the rapid and reliable issuance of control commands, the system in this embodiment preferably uses GOOSE (General Object-Oriented Substation Event) messages that conform to the IEC61850 standard for intelligent substation automation for communication.

[0127] Specifically, when the main controller generates a set of control commands (e.g., reactive power command u1 for the energy storage system and load disconnection command u2 for the charging pile group), it encapsulates these commands in GOOSE messages. Unlike traditional application layer communication methods based on TCP / IP, GOOSE messages are event-driven messages based on the publish / subscribe model and transmitted directly at the data link layer (OSI layer 2). They are broadcast directly from the publisher (main controller) to all subscribers (each load execution unit) in the network via Ethernet multicast.

[0128] In this embodiment, GOOSE message communication bypasses the complex TCP / IP protocol stack, eliminating the encapsulation and decapsulation processes at the network and transport layers, thus reducing communication latency. Experimental data shows that in a port area LAN Ethernet environment, the end-to-end transmission latency of GOOSE messages can be stably controlled within 3 milliseconds. GOOSE messages have built-in retransmission and heartbeat mechanisms, ensuring high communication reliability and providing a foundation for achieving millisecond-level control responses.

[0129] Step two: Perform precise timing coordination based on GPS-synchronized clocks.

[0130] To achieve precise timing coordination among multiple flexible load execution units distributed in different physical locations within the port area, each load execution unit in this embodiment is equipped with a synchronization clock module based on the Global Positioning System (GPS).

[0131] Specifically, the GPS synchronization clock module, by receiving satellite signals, can provide the local load controller with a UTC (Coordinated Universal Time) standard time reference with sub-microsecond accuracy. Upon receiving the GOOSE message from the previous step, each load execution unit does not execute immediately, but instead parses the encapsulated control commands and the crucial target execution timestamp t from the message. exec This timestamp is generated by the main controller based on the alert's transition time t. jump Precisely set (e.g., t) exec =t jump -2ms).

[0132] The controller continuously compares the target execution timestamp with the local high-precision GPS clock, and only triggers the corresponding control action (such as changing the output by controlling the inverter or disconnecting the load by the circuit breaker) when the local clock reaches the target timestamp.

[0133] Step 3: Implement an integrated distributed synchronous execution workflow.

[0134] The complete, collaborative instruction execution workflow is as follows: The main controller generates a time stamp t containing specific control parameters and the target execution timestamp at time t1. exec The control command is broadcast via the IEC61850 GOOSE network. After a delay of less than 3ms, all relevant load execution units receive the command at approximately t1+3ms. Each execution unit waits independently, based on its local GPS synchronization clock. At a precise t... exec Even if these execution units are physically located hundreds of meters apart, they can still achieve synchronous control actions with an error within 1 microsecond.

[0135] By organically combining the above-mentioned communication and time synchronization architecture, this invention solves the key engineering problem of realizing millisecond-level collaborative control of multiple agents (flexible loads) in a distributed system, and realizes that the timing decoupling strategy can be accurately executed in the physical world, ensuring the final magnetic flux suppression effect.

[0136] As an alternative implementation method, in scenarios where GPS signals are unavailable or where cost requirements are more stringent, time synchronization can also be achieved using a local area network-based precision time protocol (PTP, IEEE1588), which can also achieve sub-microsecond synchronization accuracy.

[0137] Example 14: In a preferred embodiment, during the separation stage of the magnetostrictive vibration components in the magnetic flux spatial distribution inversion process driven by vibration modes, specifically, it can also be as follows: Read the time-aligned vibration signal matrix V(t), obtain the spectrum F(ω) through Fast Fourier Transform, and identify the actual frequency f of the 50Hz fundamental frequency and its harmonic components. n (Considering grid frequency fluctuations of ±0.2Hz). Based on the amplitude A of each harmonic... n Design the quality factor Q of the notch filter. n =35*(1+0.1*log(A n / A noise )), where A noise To determine the background noise level, construct a cascaded notch filter array: H(z) = Π(1-2r) n *cos(2πf n / f s )z -1 +r n 2 z -2 ) / (1-2r n *cos(2πf n / f s )z -1 +r n 2 z -2 ), where r n =1-π*f n / (Q n *f s After processing, the power frequency vibration signal V is obtained. filtered (t).

[0138] Using the dq axis current component i d (t) and i q (t) Calculate the instantaneous magnetization M(t) = sqrt(i d 2 +i q 2 ) / N turns , where N turns This represents the number of turns in the winding. According to the magnetostrictive physical relationship ε=λ... s *B 2 Construct a constraint matrix G, where G ij =λ s *M i (t)*M j (t) represents the magnetostrictive coupling relationship between time i and time j, λ s M is the magnetostrictive coefficient. i (t) represents the magnetization intensity at time i. Define the independence metric function: J(W)=E{log(cosh(W T *V filtered))}-E{log(cosh(W T *G·*V filtered ))}, generate the ICA constraint matrix G and the independence objective function J(W).

[0139] Initialize the separation matrix W0 as the identity matrix, and then process the power frequency vibration signal V. filtered V was obtained by performing whitening pretreatment. white Perform constrained ICA iteration: In the k-th iteration, calculate the gradient: ▽J=E{tanh(W T *V white )*V white T}-G*E{tanh(W T *G*V white )*V white T Update the separation matrix W(k+1) = Wk - μ*▽J, where the learning rate μ = 0.01 / (1 + 0.1k) adaptively decays. Use the kurtosis criterion: κ=E{s 4} / E{s 2} 2 -3 Identify the magnetostrictive component (κ>2 indicates a supergaussian distribution), iterate until ||W k+1 -W k ||<10 -4 Output magnetostrictive vibration component matrix V ms (t)=W final *V white W final These are the coefficients of the final separation matrix.

[0140] For the magnetostrictive vibration component matrix V ms Multi-resolution analysis is performed on each channel. The Morlet wavelet basis ψ(t) = exp(-t) is used. 2 / 2)·exp(i*ω0*t) performs continuous wavelet transform CWT(a,b)=-t ms (t)*ψ*((tb) / a)dt / sqrt(a), where scale a corresponds to frequency f=f s / (2πa), f s The fundamental frequency is used. Thirty-two logarithmically spaced frequency points are selected within the 100-5000Hz range, and the instantaneous amplitude |CWT(f,t)| and phase arg(CWT(f,t)) of each frequency are extracted. The dominant mode is identified through singular value decomposition. [U,S,V]=SVD(CWT), retaining the first 6 principal modes, construct the modality feature vector: Φ(f,t)=[|CWT1|,arg(CWT1),...,|CWT6|,arg(CWT6)] T .

[0141] Example 15, another preferred embodiment, in the inversion process of magnetic flux spatial distribution driven by vibration modes, the Bayesian magnetic flux inversion step based on physical constraints can specifically be as follows: Based on the modal eigenvector Φ(f,t) and the magnetostrictive vibration component matrix V ms Construct the nonlinear observation equation V ms =λ s *(B 2 (x sensor )-B0 2 )+λ d *B*dB / dt+η, where x sensor For 4 sensor positions, λ d Let be the d-axis magnetostriction, and η be the bias coefficient. Performing a Taylor expansion of the magnetic flux density B near the operating point B0: V ms ≈s s *B0*(B-B0)+λ d *B0·dB / dt+O((B-B0) 2 Construct a linearized transfer matrix H, whose elements H ij =2λ s *B0*G ij (x i ,ξ j ) represents a spatial point ξ j The magnetic flux relative to the sensor position x i The contribution of vibration, where the Green's function G ij Obtained through finite element pre-calculation. Output linearized observation equation V. ms =H·B+η and the transfer matrix H.

[0142] Physical constraints are constructed based on Maxwell's equations. Divergence constraint: ∂B / ∂D is calculated using the finite difference operator D, requiring ||D·B|| < ε. div , where ε div =10 -6 T / m. Curl constraint: According to Ampere's law ▽×B=μ0·J, where the current density J is calculated from the dq-axis current components, J=(i d ·cos(θ)-i q ·sin(θ)) / A core A core Let B be the cross-sectional area of ​​the iron core. Saturation constraint: Soft constraint p(B) ∝ exp(-α·max(0,|B|-B)). sat ) 2), where α=100 controls the constraint strength. Spatial smoothness: modeled by Markov random fields p(B i |B neighbors )∝exp(-β·Σ||B i -B j || 2 ), β=0.5, B i Let be the magnetic flux at time i. Construct the prior probability density p(B) = pi, combining all constraints. div (B)·p curl (B)·p sat (B)·p smooth (B), p div (B) represents the divergence constraint factor, p curl (B) is the curl constraint factor, p sat (B) is the saturation constraint factor, p smooth (B) is the smoothing constraint factor.

[0143] Generate N p =1000 initial particles {B i 0 Each particle is sampled through a truncated Gaussian distribution: B i 0 ~N(B nominal ,Σ0)·I(|B| sat ), where B nominal =1.2T is the rated magnetic flux density, and the covariance Σ0=0.1 2 • I, I() are indicator functions. Calculate the importance weight w for each particle. i 0 =p(V ms |B i 0 )·p(B i 0 ) / q(B i The proposed distribution q is a Laplace distribution to enhance robustness. Normalized weights w i =w i / Σw j Calculate the effective number of particles N eff =1 / Σw i 2 When N eff <N p Resampling is triggered at / 3, outputting the initial particle set {B}. i 0 ,w i 0}

[0144] ​Perform K=20 iterations. In the k-th iteration, the prediction step uses a random walk B. i k =B i k-1 +σ w ·ξ i The process noise ξ i ~N(0,I), σ w =0.01·(1-k / K) adaptively decreases. Constraint projection: Perform gradient projection B on particles that violate physical constraints. i k =B i k -γ·▽C(B i k ), where C(B) = ||D·B|| 2 +max(0,|B|-B sat ) 2 To constrain the violation rate, the step size γ = 0.1. Update step: Based on the new observation V ms k Calculate the likelihood p(V) ms k )|B i k )=exp(-||V ms k -H·(B i k ) 2 || / (2σ n 2 )), where the noise variance σ n 2 Estimating online using innovative sequences. Updating weights w. i k =w i k-1 ·p(V ms k |B i k ), output update particle set {B i k ,w i k}

[0145] Extract the maximum a posteriori estimate from the updated particle set. Calculate the weighted mean B. MAP =Σw i K ·B i K As a point estimate, the posterior covariance Cov(B) = Σw is calculated using the particle distribution. i K ·(B iK -B MAP )·(B i K -B MAP ) T Extracting the diagonal elements yields the uncertainty σ at each point. B (x,y,z)=sqrt(diag(Cov)). This identifies the high-confidence region Ω. confident ={(x,y,z)|σ B <0.05T} and low confidence region Ω uncertain ={(x,y,z)|σ B >0.15T}. Kriging interpolation B is used for low-confidence regions. kriging =Σλ i ·B i , where λ i The weights are used to calculate the three-dimensional magnetic flux density distribution field B(x,y,z,t) and the confidence plot C. con f(x,y,z)=1-σ B / B sat .

[0146] This invention solves the problem of observing the internal magnetic field state of a transformer through a series of mutually supporting technical features, achieving real-time, non-invasive imaging of the internal three-dimensional magnetic flux density distribution field. This invention uses easily measurable multi-point vibration signals from the transformer's exterior as input corresponding to specific data in the technical field, and employs a constrained independent component analysis algorithm to apply constraints embedded with the physical laws of magnetostriction to the signal separation process, extracting the magnetostrictive vibration component V directly related to the internal magnetic field state. ms Building upon this foundation, a Bayesian inversion framework based on particle filtering is employed, placing the inverse problem solution process within a probabilistic reasoning system. Furthermore, a strongly physical prior constraint framework, comprised of Maxwell's equations and saturation characteristics, is introduced to constrain the solution space. This deep integration of algorithmic features and physical rule characteristics ensures that even with limited observational data, the inverted output three-dimensional magnetic flux density field B(x,y,z,t) (corresponding to specific data in the technical field) possesses high physical fidelity. Therefore, this invention enables the system to, for the first time, accurately determine where saturation will occur, providing high-resolution physical state input for all subsequent advanced predictions and accurate control, thus changing the previous reliance on indirect, macroscopic electrical quantities for rough judgments.

[0147] This invention solves the time scale mismatch problem by constructing a unique dual-mode adaptive prediction architecture, realizing a paradigm shift from post-event analysis to pre-event millisecond-level early warning. In this invention, the system is designed with intelligent switching decision logic. The input parameters of this logic are closely related to technical data, simultaneously monitoring the comprehensive current change rate (representing the severity of external disturbances) and the maximum value of the local saturation index (representing the vulnerability of the transformer's internal structure), and dynamically adjusting the switching threshold based on these, making the selection of the prediction mode both sensitive and reliable. When the system is stable, a precise prediction algorithm based on the magnetic diffusion equation is used to solve the problem in a refined manner over a longer time scale, maintaining a high-precision dynamic benchmark for the system. Once a violent power flow switch is detected, it immediately switches to a fast prediction algorithm based on higher-order Taylor expansion, outputting the prediction result within milliseconds. This combination of fast and slow architecture, through the mutual support of algorithmic features, achieves a good balance between computational efficiency and prediction accuracy. For port scenarios, when the quay crane equipment is about to start high-power operation, this technology can accurately output an early warning vector containing the moment of change, amplitude and position within tens of milliseconds before the actual saturation occurs, thus providing the control system with an effective time window to take active suppression measures.

[0148] This invention designs a three-layer cooperative control strategy with time-decoupling. The deep coupling and interaction between its control rules and algorithmic features effectively solves the dilemma of cooperative suppression failure, achieving rapid and precise suppression of magnetic flux jumps with minimal control cost and system impact. The input to this strategy includes a precise timestamp t. jump The warning vector. Based on this, the control algorithm decomposes the originally single control task into three layers that are staggered in time and complementary in function. The first layer, at t jump In the first 5ms, the energy storage system with the fastest dispatch response performs low-cost proactive reactive power compensation for flexible pre-suppression; the second layer, in t jump In the first 2ms, based on the predicted saturation danger zone, the next fastest load (such as charging piles) is scheduled to perform small-scale, high-precision active power cut-off for targeted strong suppression; the third layer, in t jump In the last 5ms, a multi-objective optimization algorithm balancing adjustment costs and load balance is initiated, scheduling all resources, including those with slower response times, for system-level economic recovery. By combining abstract control rules with the specific technical characteristics (response time, adjustment costs) of flexible loads (energy storage, charging piles, quay cranes) in the port area, a scheduling method is implemented that achieves functional mutual support. This ensures that each suppression action is executed with the optimal resource combination at the best time, improving the system's internal control performance and efficiency. Compared to traditional single load reduction, this invention reduces interference with port operations, improving port operational efficiency while ensuring grid safety.

[0149] Finally, it should be noted that the above embodiments are merely illustrative of the technical solutions of the present invention and not intended to limit it. Those skilled in the art should understand that modifications or equivalent substitutions can be made to the specific embodiments of the present invention, but such modifications or alterations are all within the scope of protection of the pending claims.

Claims

1. A millisecond-level collaborative suppression method for flexible loads in port areas based on magnetic flux jump prediction, characterized in that, include: The vibration signals and electrical operation data of the transformer's exterior are acquired, and the three-dimensional magnetic flux density distribution field characterizing the internal magnetic field state of the transformer is reconstructed based on these data. Receive the three-dimensional magnetic flux density distribution field, predict future magnetic flux jump events based on its temporal evolution law, and generate a magnetic flux jump early warning vector containing jump time, amplitude and location information; Based on the magnetic flux jump early warning vector and the real-time acquisition of the port area's flexible load operation status, hierarchical collaborative control commands are generated and issued for execution.

2. The millisecond-level collaborative suppression method for flexible loads in port areas based on magnetic flux jump prediction as described in claim 1, characterized in that, The inversion and reconstruction of the three-dimensional magnetic flux density distribution field includes: Based on Maxwell's equations and the saturation magnetic flux density constraint of the iron core material, a physical prior probability distribution is constructed to constrain the solution space; Initialize particles representing multiple candidate three-dimensional magnetic flux density distribution fields to form an initial particle set; An observational likelihood function describing the relationship between vibration and magnetic flux is established based on multi-point vibration signals. Combined with the physical prior probability distribution, the weights of each particle in the initial particle set are iteratively updated to obtain the updated particle set. The maximum a posteriori estimate is extracted from the weighted probability distribution of the updated particle set and identified as the three-dimensional magnetic flux density distribution field.

3. The millisecond-level collaborative suppression method for flexible loads in port areas based on magnetic flux jump prediction as described in claim 2, characterized in that, Constructing the physical prior probability distribution for constraining the solution space includes: The magnetic flux divergence constraint established based on the finite difference operator characterizes the continuity of magnetic flux; Based on Ampere's law and the current density derived from electrical operating data, the flux curl constraint is used to characterize the curl of the magnetic field generated by the current distribution. The nonlinear magnetization saturation characteristics of iron core materials are characterized by saturation soft constraints defined by probability functions. Spatial smoothness constraints established by Markov random field models characterize the continuous and gradual variation of magnetic flux density in space.

4. The millisecond-level collaborative suppression method for flexible loads in port areas based on magnetic flux jump prediction as described in claim 1, characterized in that, Before reconstructing the three-dimensional magnetic flux density distribution field characterizing the internal magnetic field state of the transformer, the process also includes separating the magnetostrictive vibration components, including: Based on the instantaneous magnetization intensity extracted from electrical operation data, physical constraints are constructed to guide signal separation; Physical constraints are applied to the independent component analysis algorithm to perform constraint optimization separation of multi-point vibration signals and extract the magnetostrictive vibration components.

5. The millisecond-level collaborative suppression method for flexible loads in port areas based on magnetic flux jump prediction as described in claim 1, characterized in that, Generate hierarchical collaborative control instructions and issue them for execution, including: The expected time of flux jump occurrence is extracted from the flux jump warning vector; based on the time of jump occurrence, a hierarchical collaborative control command is constructed through a time-decoupled generation process, specifically including: Before the expected transition occurs, the first layer of control commands for advanced phase compensation and the second layer of control commands for key point reverse suppression are generated sequentially. After the expected transition occurs, a third-level control command is generated for optimized power reallocation of the system. The first, second, and third layer control commands are integrated to form a hierarchical collaborative control command, which is then sent to the corresponding flexible load.

6. The millisecond-level collaborative suppression method for flexible loads in port areas based on magnetic flux jump prediction as described in claim 5, characterized in that, Generate the first-level control commands for lead phase compensation, including: Obtain the expected jump amplitude from the magnetic flux jump warning vector; Based on the jump amplitude, determine the reactive power compensation amount used to offset part of the magnetic flux change; Set an injection timing sequence with a specific leading phase for reactive power compensation, so that reactive power injection leads magnetic flux disturbance in phase; Generate reactive power commands that include reactive power compensation amount and injection timing, as the first-level control commands.

7. The millisecond-level collaborative suppression method for flexible loads in port areas based on magnetic flux jump prediction as described in claim 5, characterized in that, Generate a second layer of control instructions for keypoint inverse suppression, including: Identify the expected saturation danger zone inside the transformer from the magnetic flux jump warning vector; For the saturation danger zone, a reverse magnetic flux distribution is designed to counteract the saturation effect; Based on the principle of energy equivalent conversion, the magnetic field energy required to maintain this reverse magnetic flux distribution is converted into the active power regulation amount to be performed in the power grid domain. Based on the active power adjustment, loads with rapid disconnection capabilities are selected from the flexible load operation status of the port area, and a second-level control command is generated.

8. The millisecond-level collaborative suppression method for flexible loads in port areas based on magnetic flux jump prediction as described in claim 5, characterized in that, Generate third-level control instructions for system power optimization and reallocation, including: The active power adjustment amount determined by the second-level control command is set as the system's total power balance target for the current power redistribution process; Based on the cost weights and operational limits included in the flexible load operation status of the port area, an objective function for a multi-objective optimization scheduling model is constructed. Under the triple constraints of satisfying the system's total power balance target, the adjustable capacity of a single load, and the ramp rate, the optimal scheduling model is solved to obtain a set of optimal power adjustment sequences, which constitute the third layer of control commands.

9. The millisecond-level collaborative suppression method for flexible loads in port areas based on magnetic flux jump prediction as described in claim 1, characterized in that, Predict future magnetic flux jump events and generate magnetic flux jump early warning vectors, including: Real-time monitoring of the power flow status of the port area power grid, as represented by electrical operation data, to determine whether a power flow switching event has occurred; When it is determined that the event has not occurred, the system operates in precise prediction mode, which uses the first prediction algorithm to perform refined evolution prediction on a longer time scale based on the three-dimensional magnetic flux density distribution field. When the judgment occurs, the system switches to the fast prediction mode, and uses the second prediction algorithm to perform incremental fast prediction on a shorter time scale based on the three-dimensional magnetic flux density distribution field. The prediction results under the current operating mode are output as a magnetic flux jump warning vector.

10. A millisecond-level collaborative suppression method for flexible loads in port areas based on magnetic flux jump prediction, as described in claim 9, is characterized in that... Determining whether a power flow switching event has occurred includes: The real-time comprehensive current change rate is calculated from electrical operation data; The maximum value of the local saturation index, which characterizes the current saturation risk of the transformer, is calculated from the three-dimensional magnetic flux density distribution field. Based on the statistical frequency of historical trend switching events and the maximum value of local saturation index, an adaptive switching threshold is dynamically generated. The comprehensive current change rate is compared with the adaptive switching threshold in real time, and a judgment is made on whether a power flow switching event has occurred based on the comparison result.

Images