Intelligent self-monitoring temperature management system for box-type substation

Abstract

This invention discloses an intelligent self-monitoring temperature management system for prefabricated substations, relating to the field of smart grid equipment monitoring technology. It solves the problems of large temperature measurement deviations, low reliability, delayed early warnings, inaccurate hotspot location, and coarse heat dissipation control in existing systems. The solution involves: a data acquisition module that collects multi-source signals in parallel and extracts temperature time sequences; a boundary calibration module that fuses the data to generate three-dimensional boundary conditions; a multi-physics field solving module that obtains the internal temperature / stress field; a physical information prediction module that fuses heat transfer physical constraints to train a neural network to predict hotspot migration trends; and a hierarchical scheduling module that uses model predictive control to solve for real-time optimization commands for fan and oil pump coordination. This invention significantly improves the accuracy of internal temperature monitoring, the reliability of hotspot prediction, and the precision of heat dissipation control in prefabricated substations, effectively extending equipment insulation life and enhancing system operational safety and reliability.

Description

Technical Field

[0001] This invention relates to the field of smart grid equipment monitoring technology, and more specifically to a smart self-monitoring temperature management system for prefabricated substations. Background Technology

[0002] The equipment inside a prefabricated substation generates a significant amount of heat during long-term operation. If it is not monitored and adjusted in real time, it can easily lead to insulation aging, excessive temperature rise, or even equipment failure.

[0003] Existing intelligent self-monitoring temperature management systems primarily rely on temperature sensors to collect temperatures from key components. Heating or cooling commands are then issued via local control units or remote SCADA platforms to maintain the internal temperature of the enclosure within a safe range, improving operational reliability and maintenance efficiency. This system mainly comprises a layered architecture of sensing, communication, control, and execution layers. The sensing layer commonly uses platinum resistance thermometers (PT100) and thermocouples to measure the temperature of the oil tank and casing. Fiber Bragg gratings (FBGs), increasingly used in hotspot areas, enable multi-point distributed temperature measurement to resist electromagnetic interference. The communication layer uses RS-485 / Modbus, LoRa, NB-IoT, or other methods to centralize sensor signals at a PLC or industrial gateway before uploading them to SCADA or a cloud platform. The control layer generally employs hierarchical control logic based on preset thresholds, supplemented in some systems by empirically driven PI regulation to adjust fan and oil pump speeds. The execution layer consists of fans, oil coolers, and ceramic heaters, which receive control commands and perform thermal management actions.

[0004] Taking patent CN107579467A as an example, this solution adds a ceramic heater, air supply duct, and fan inside the mechanism box, and sends temperature data to the monitoring center via wireless communication to realize closed-loop control and remote monitoring of heating; patent CN113036921A constructs a two-layer communication architecture between the box-type substation and the central server. The server obtains the temperature and humidity and the temperature of each device in real time, and issues an environmental adjustment command when the preset threshold is exceeded to improve the accuracy and timeliness of communication; patents CN208737319U and CN222107412U further integrate the functions of power supply and distribution, constant temperature regulation, and heat insulation protection. However, in actual implementation, firstly, the Bragg wavelength drift of FBG is not only sensitive to temperature, but also highly sensitive to stress changes. The minute stresses generated by transformers during load fluctuations, structural vibrations, and temperature cycles are added to the optical fiber, causing random deviations in temperature readings. This is especially problematic during reclosing or sudden load changes, where errors can exceed ±3℃. While existing literature has proposed dual FBG differential compensation to offset the effects of mechanical strain, inconsistent wiring tightness, vibration fatigue, and minor movements at anchor points can lead to compensation failure, resulting in alarm delays or false alarms and compromising system reliability. Secondly, although some manufacturers have attempted to utilize purely data-driven short-term temperature prediction models such as LSTM and SVR, these algorithms only input surface temperature and current load timing, failing to incorporate the physical constraints of three-phase fluid heat conduction, convection stratification, and oil temperature drift in the oil-metal-air three-phase fluid architecture. This leads to frequent warning delays and hotspot location errors during high-temperature seasons or sudden heavy loads. Maintenance decisions often only allow for a "switch-on" response to surface heat dissipation devices, failing to address deeper hotspot risks in a timely manner, ultimately significantly reducing equipment insulation life and system reliability. Summary of the Invention

[0005] To address the shortcomings of existing technologies, this invention discloses an intelligent self-monitoring temperature management system for prefabricated substations, aiming to solve the problems mentioned in the background art.

[0006] To achieve the above-mentioned technical effects, the present invention adopts the following technical solution:

[0007] A smart self-monitoring temperature management system for a prefabricated substation includes: a data acquisition module for parallel acquisition of FBG fiber reflection spectrum, quartz strain gauge voltage signal, thermocouple voltage and infrared array radiation intensity data; extracting time spectrum lines of each channel through multi-channel wavelet packet time-frequency decomposition; and using an adaptive spectral line separation algorithm to separate the temperature principal component and stress residual to obtain temperature time series and stress residual time series.

[0008] The boundary calibration module is used to fuse and correct the output data of the data acquisition module based on algebraic surface fitting and least squares correction methods to obtain a three-dimensional boundary condition mapping set.

[0009] The multiphysics solution module is used to input the three-dimensional boundary condition mapping set into the finite element / finite volume solution framework, construct a three-phase fluid heat conduction and convection stratification model including oil, metal and air, and obtain the internal temperature field and stress field profile data of the box-type substation.

[0010] The physical information prediction module is used to jointly train the output data of the multiphysics solution module and the historical running sequence by using a graph neural network and integrating the residual term of the heat transfer conservation equation as a loss function to generate short-term hotspot migration and temperature rise trend prediction results.

[0011] The hierarchical scheduling module is used to construct a second-order cone programming model containing fan speed and oil pump flow rate through model predictive control methods. It combines the prediction results of the physical information prediction module with equipment operating condition constraints to solve the problem in real time and output multi-equipment coordinated control commands.

[0012] Based on the above technical solutions, the positive and beneficial effects of the present invention are as follows:

[0013] 1. This invention acquires fiber optic grating (FBG) reflection spectra, quartz strain gauge voltage signals, thermocouple voltages, and infrared array data in parallel. It employs multi-channel wavelet packet time-frequency decomposition combined with an adaptive spectral line separation algorithm to effectively and reliably distinguish between temperature-induced and stress-induced spectral changes, generating independent temperature and stress residual time-series signals. The quartz strain gauge directly and independently measures mechanical stress, and its signal is unaffected by temperature, providing crucial input for the separation algorithm. This avoids the compensation failure problems caused by differences in on-site wiring, vibration fatigue, or loose anchor points in traditional dual-FBG differential compensation methods. Consequently, the final temperature measurement results significantly reduce errors introduced by stress interference during load surges, reclosing, or structural vibrations (e.g., random deviations exceeding ±3℃), greatly improving the accuracy and reliability of temperature monitoring data. Temperature readings more accurately reflect the actual thermal state inside the equipment, thus significantly reducing false alarms and missed alarms caused by stress interference.

[0014] 2. The multi-physics field solution module simulates the complex heat transfer process inside the prefabricated substation, particularly the impact of key physical phenomena such as oil temperature drift and thermal stratification on temperature distribution and hotspot formation, thereby obtaining accurate internal temperature and stress field profile data. The physical information prediction module utilizes a graph neural network to process this physical field data and historical operating data, incorporating the physical constraints of the heat conduction conservation equation (through its residual terms) into the neural network training process. This ensures that the final short-term hotspot migration and temperature rise trend prediction results not only rely on the statistical regularities of historical data but also strictly adhere to the fundamental physical laws of energy transfer and conservation. Therefore, this prediction model can more accurately capture the formation and migration patterns of deep hotspots. Under harsh operating conditions such as high-temperature seasons or sudden heavy loads, compared to traditional pure data-driven models (such as LSTM and SVR), it significantly reduces the delay time of hotspot early warning and improves the accuracy of hotspot location.

[0015] 3. The hierarchical scheduling module, based on the accurate hotspot location and temperature rise trend prediction results provided by the physical information prediction module, constructs a second-order cone programming model with fan speed and oil pump flow rate as continuous optimization variables using the Model Predictive Control (MPC) method. This model fully considers the actual operating constraints of the equipment (such as maximum fan speed, maximum oil pump flow rate, and safe temperature rise limits). The system can solve this optimization problem in real time, dynamically generating and outputting refined and coordinated control commands for multiple devices such as fans and oil pumps. This makes the response of the cooling system no longer a simple "on" or "off," but rather allows for continuous and precise adjustment of the output of cooling devices in different areas based on the predicted hotspot location and temperature rise rate (such as locally increasing oil circulation speed or directionally increasing airflow intensity in specific areas). Therefore, the system can more effectively concentrate cooling resources at the most needed locations, avoiding energy waste caused by global overcooling and completely solving the problem of persistent or worsening hotspot hazards caused by insufficient local cooling. This precise and dynamic cooling control can effectively suppress the thermal aging rate of insulation materials, significantly extend the service life of core equipment such as transformers, and fundamentally improve the long-term operational reliability of the entire prefabricated substation system. Attached image description:

[0016] To more clearly illustrate the technical solutions in the embodiments of the present invention or the prior art, the drawings used in the description of the embodiments or the prior art will be briefly introduced below. Obviously, the drawings described below are only some embodiments of the present invention. For those skilled in the art, other drawings can be obtained based on these drawings without creative effort, wherein:

[0017] Figure 1 This is an architecture diagram of an intelligent self-monitoring temperature management system for a prefabricated substation according to the present invention.

[0018] Figure 2 This is a flowchart illustrating the separation steps of the adaptive spectral line separation algorithm of the present invention.

[0019] Figure 3 This diagram illustrates the steps of the fusion correction method combining algebraic surface fitting and least squares correction methods of the present invention.

[0020] Figure 4 This is a schematic diagram illustrating the principle of the anisotropic mesh refinement algorithm of the present invention.

[0021] Figure 5 This is a schematic diagram of the model predictive control principle framework of the hierarchical scheduling module of the present invention;

[0022] The diagram is labeled as follows: 100, data acquisition module; 200, boundary calibration module; 300, multiphysics solution module; 400, multiphysics solution module; 500, hierarchical scheduling module. Detailed Implementation

[0023] To make the objectives, technical solutions, and advantages of the embodiments of this application clearer, the technical solutions of this application will be clearly and completely described below in conjunction with the embodiments. Obviously, the described embodiments are only some embodiments of this application, not all embodiments. Based on the embodiments of this application, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of this application.

[0024] To facilitate understanding of this embodiment, a detailed description of the intelligent self-monitoring temperature management system for prefabricated substations disclosed in this application will be provided first. Please refer to [link to relevant documentation]. Figure 1 The diagram shows a framework of an intelligent self-monitoring temperature management system for a prefabricated substation. The system includes:

[0025] The data acquisition module 100 is used to acquire FBG fiber reflection spectrum, quartz strain gauge voltage signal, thermocouple voltage and infrared array radiation intensity data in parallel. It extracts the time spectrum lines of each channel through multi-channel wavelet packet time-frequency decomposition, and uses an adaptive spectral line separation algorithm to separate the temperature principal component and stress residual to obtain the temperature time series and stress residual time series.

[0026] Specifically, in the implementation of the intelligent temperature management system for prefabricated substations, the specific operation for extracting time-frequency lines through multi-channel wavelet packet time-frequency decomposition is as follows:

[0027] First, dedicated frequency bands are allocated based on the physical characteristics of the sensors. The FBG fiber optic reflection spectrum is allocated a 0.1-20Hz fundamental low-frequency channel (the upper limit of this band can be extended to 25Hz based on the magnetostrictive characteristics of transformer cores, such as in the vibration scenario of heavy equipment in industrial areas). The quartz strain gauge voltage is allocated a 0.1-50Hz wideband channel to capture reclosing impact stress waves (typically induced by short-circuit current at 2-40Hz damped oscillations). Thermocouples and infrared arrays are allocated a 0.5-100Hz mid-to-high frequency channel for detecting oil flow turbulence (typical turbulence frequencies are 10-80Hz) and partial discharge thermal radiation (transient events >60kHz fall into this band after downsampling). It should be noted that the "fundamental low-frequency channel" in this application differs from the traditional fixed bandwidth design. Its flexible passband range is dynamically adapted by thermal inertial characteristics, and the frequency band boundary parameters can be adjusted within the range of 0.05-0.2Hz according to the temperature change rate. The specific determination is made by the transformer substation thermal simulation model and does not constitute a limitation.

[0028] Based on frequency band characteristics, wavelet packet kernel functions are selected: the low-frequency channel is loaded with a Db6 wavelet kernel, whose 6th-order vanishing moment is compactly supported to optimize the extraction of slowly varying signals. The filter coefficients are defined as h[k]=[0.0352,-0.0854,-0.1350,0.4599,0.8069,0.3327] (low-pass) and g[k]=[-0.3327,0.8069,-0.4599,-0.1350,0.0854,0.0352] (high-pass), k=0,…,5; the wideband channel is loaded with a Sym8 wavelet kernel, whose frequency domain transfer function satisfies H(ω)∣2+∣G(ω)∣2=1 to suppress harmonic distortion of stress signals; the mid-to-high frequency channel is loaded with a Coif5 wavelet kernel, whose time domain supports 17 points to balance the time-frequency resolution of infrared signals. In actual implementation, the arm_biquad_cascade_df1_f32 function of the CMSIS-DSP library is called through the ARM Cortex-M7 core to realize real-time convolution operation. When the processor computing power is limited, it can be replaced with Db4 (supporting length 8 points) or Sym6 variant, which can reduce the amount of computation by 40%.

[0029] An 8-layer adaptive decomposition is performed on each channel (layer depth is based on the minimum resolution 3.9Hz frequency band at a 1kHz sampling rate), generating detail coefficients dj[n] and approximation coefficients aj[n] trees (j=1,…,8). Hard threshold filtering is implemented using channel-specific energy thresholds: λFBG = -40dB for the FBG channel (e.g., filtering out 50Hz power frequency interference), λStrain = -35dB for the strain gauge channel (suppressing mechanical vibration harmonics >50Hz), and λThermal = -30dB for the thermocouple / infrared channel (eliminating electromagnetic switch transients). The filtering formula is d^j[n]={dj[n]·(1-λ∣dj[n]∣),∣dj[n]∣≥λ0,∣dj[n]∣<λ. After reconstructing the waveform using inverse wavelet packets, a short-time Fourier transform is performed. The FBG channel uses a 500ms Hanning window (window function w(t) = 0.5 - 0.5cos(2πt / T)), and the thermocouple channel uses a 50ms rectangular window. Local maxima satisfying |X(f,t)| > 3σnoise| (σnoise is the standard deviation of background noise) are detected in the time-frequency matrix. It should be noted that "dominant spectral line" specifically refers to continuous time-frequency ridges (such as the 0.5Hz fundamental frequency and its second harmonic generated by winding temperature rise). Their extraction uses the Morlet wavelet ridge tracking algorithm instead of discrete peak detection. The ridge tracking formula is d... dtarg (X(f,t))=2πf.

[0030] When a time-frequency line is detected in the FBG-based low-frequency channel and the quartz strain gauge broadband channel and satisfies the following conditions: time synchronization |tFBG-tStrain| ≤ 5ms (corresponding to a 200Hz control period), frequency domain correlation |fFBG-fStrain| ≤ 2Hz, and phase consistency |arg(AFBG)-arg(AStrain)| ≤ π / 6 (30°), a cross-validation flag is triggered (e.g., reclosing causes a 12.5Hz stress wave to appear synchronously in both channels with a phase difference ≤ 10°). It should be noted that the "in-phase" determination can be achieved using the cross-correlation algorithm R(τ) = ∑txFBG(t)·xStrain(t+τ). Synchronization is determined when max(R(τ)) appears in τ∈[-2ms,2ms], and the phase difference threshold is allowed to fluctuate within ±10°.

[0031] Based on the energy percentage of the marked spectral lines, E% is calculated as: E% = ∫T∣X(f,t)∣2dt∑k∫T∣X(fk,t)∣2dt×100%. Spectral lines in the FBG channel with E% ≥ 85% (e.g., 92% of the 0.8Hz oil temperature drift spectrum) are classified into the temperature principal component temporary storage area, and spectral lines in the strain gauge channel with E% ≥ 70% (e.g., 78% of the 35Hz fan vibration spectrum) are classified into the stress residual temporary storage area. It should be noted that the "temporary storage area" is a double buffer serving the adaptive spectral line separation algorithm, and its storage structure is defined as a C language structure:

[0032] struct SpectrumBuffer

[0033] {float frequency; float amplitude; float energy_ratio; uint8_tsource_channel;} When an unclassified high-energy spectral line (e.g., >20Hz component with energy percentage >15%) appears in the FBG channel, an insulation degradation warning is triggered and an infrared array scan is activated. The energy percentage threshold can be configured within the range of 80-95% for the FBG channel and 60-80% for the strain gauge channel (lower limit for urban power grids, upper limit for industrial areas with strong vibrations), specifically determined based on the characteristics of the vibration spectrum at the site. Adjustments to parameters such as frequency band range and wavelet packet kernel type need to be optimized based on the sensor sampling rate (e.g., the FBG demodulator supports 1-10kHz) and the substation installation environment (urban / industrial power grid); this solution does not impose limitations on these.

[0034] Please see Figure 2 The specific method for separating the principal temperature component and the stress residual using the adaptive spectral line separation algorithm is as follows:

[0035] S1. Extract the fundamental frequency, third harmonic amplitude and phase difference parameters of historical FBG reflection spectrum data using the FFT algorithm, and construct a temperature characteristic spectral template library;

[0036] S2. Receive the time-spectrum lines output by the multi-channel wavelet packet decomposition in real time, and perform the following operations for each spectral line:

[0037] s200. Match the temperature feature template library using the Pearson correlation coefficient calculator. If the correlation coefficient is ≥0.9, mark it as a candidate temperature spectral line.

[0038] s201. Calculate the derivative of the current spectral line energy with the energy of the previous 10 seconds using the energy change rate differential method. If the derivative value exceeds ±5% / s, it is marked as an abnormal spectral line.

[0039] S3. For all quartz strain gauge channel spectral lines that are not marked as candidate temperature spectral lines, group them according to the principle of frequency difference ≤ 2Hz by using spectral line clustering analysis units, and calculate the stress contribution weight of each group based on the energy ratio of the spectral lines within the group by using the weighted least squares method.

[0040] S4. Input the candidate temperature spectral lines into the principal component analyzer to extract the first three principal component vectors; and map the principal component vectors to the original FBG spectral sequence through the vector projection algorithm to reconstruct the temperature principal component time series.

[0041] S5. From the original quartz strain gauge signal, the spectral components corresponding to the stress contribution weights are subtracted through convolution inversion operation to generate the stress residual time series.

[0042] Specifically, in implementation, a temperature characteristic spectral line template library is first constructed based on historical FBG reflectance spectrum data. During implementation, the fundamental frequency, third harmonic amplitude, and phase difference parameters are extracted using an FFT algorithm (for example, the typical fundamental frequency during the temperature rise process of an oil-immersed transformer winding is 0.8-1.2Hz, and the third harmonic amplitude is 15-25dB lower than the fundamental frequency). This template library is dynamically updated using a sliding time window, and its storage structure is defined as an array of structures containing three fields: frequency value f, amplitude A, and phase φ. It should be noted that the "temperature characteristic spectral line template library" in this application differs from a static database; it specifically refers to a dynamically adjusted circular buffer (with a capacity of the most recent 72 hours of data) based on seasonal load changes, and its update cycle can be configured to be 10-60 minutes.

[0043] After receiving the time-spectral lines output from the multi-channel wavelet packet decomposition in real time, a two-stage labeling operation is performed on each spectral line: a template library is matched using the Pearson correlation coefficient calculator, and the correlation coefficient formula is as follows:

[0044]

[0045] Where Ai is the real-time spectral amplitude sequence, Ti is the template amplitude sequence, and when ρ≥0.9, it is marked as a candidate temperature spectral line (e.g., the 1.0Hz spectral line under steady-state load with ρ=0.93). Simultaneously, abnormal spectral lines are detected using the energy change rate differential method, and the derivative of the current spectral line energy E_t with the energy E_{t-10} of the previous 10 seconds is calculated:

[0046]

[0047] If |ΔE′| > 5% / s (e.g., a short-circuit impact causes a sudden increase in strain spectral energy of 12% / s), it is marked as an anomalous spectral line. It should be noted that "differential rate of energy change" here specifically refers to a mutation detection mechanism that excludes normal temperature changes (ΔE′ < 0.5% / s), and its threshold can be adjusted between ±3% / s and ±8% / s based on load fluctuation characteristics.

[0048] For quartz strain gauge channel spectral lines not marked as candidate temperature spectral lines, they are grouped according to the principle of frequency difference ≤ 2Hz using spectral line clustering analysis units (e.g., spectral lines detected at 12.5Hz, 12.7Hz, and 13.0Hz are grouped into the same vibration source group). The stress contribution weight is calculated based on the energy proportion of the spectral lines within the group, and the weighted least squares method is used to solve for it. Where Φ is the frequency domain basis matrix, c is the weight coefficient vector, and w k =Ek / ∑Ei is the energy weight of the k-th spectral line (for example, the 12.5Hz spectral line accounts for 68% of the weight in a certain vibration source group). In actual implementation, the overdetermined equations are solved using the DGELS function in the LAPACK library.

[0049] Input the candidate temperature spectral lines into the principal component analyzer and extract the first three principal component vectors. Let the candidate spectral line matrix X∈Rm×n (m spectral lines × n frequency points), and after singular value decomposition: X=UΣV T The first three columns of V are taken to form the principal component vectors v1, v2, and v3 (the contribution rate is usually >95%). The principal components are mapped to the original FBG spectral sequence using a vector projection algorithm: stemp = ∑i=13(sraw·vi)vi to reconstruct the time series of the temperature principal components (e.g., to separate the winding temperature rise fundamental frequency signal from the composite spectrum). It should be noted that in this application, "vector projection algorithm" specifically refers to orthogonal decomposition in the frequency domain space rather than time domain filtering, which can avoid the phase distortion of traditional bandpass filters.

[0050] Finally, the spectral components corresponding to the stress contribution weights are subtracted from the original quartz strain gauge signal through convolution inversion. Let the frequency domain representation of the strain signal be Sstrain(f), and the stress residual generation formula be: Sresidual(f)=Sstrain(F)-∑g=1Gcg·Φg(F), where G is the number of vibration source groups, cg is the weight coefficient obtained in step S3, and Φg(f) is the characteristic spectrum of the g-th group (e.g., the 12.5Hz magnetostrictive peak and its 14.5Hz harmonic). In practice, the time-domain signal is reconstructed using frequency domain dot product followed by inverse FFT, implemented in the Cortex-M7 core by calling the arm_cmplx_mult_real_f32 function of the CMSIS-DSP library.

[0051] During implementation, when an abnormal spectral line marker (ΔE′>5% / s) is detected, the linkage boundary calibration module 200 initiates emergency three-dimensional thermal field reconstruction; when the temperature principal component time series is input into the multiphysics solution module 300, high-frequency components with an energy percentage <15% (such as vibration noise >20Hz) are automatically ignored; when the stress residual time series exceeds a threshold, a structural health assessment is triggered, for example, detecting low-frequency residual vibration of 0.5-2Hz indicates fastener loosening. The adjustable parameter ranges are shown in Table 1.

[0052] Table 1 Parameter Range Table

[0053]

[0054] It should be noted that the "spectral line clustering analysis unit" mentioned above specifically refers to the dynamic grouping module based on the DBSCAN algorithm, whose neighborhood radius ε = 2Hz can be adjusted to 1-3Hz according to the vibrational spectrum broadening; the "convolution inversion operation" is different from simple spectral subtraction and includes a phase calibration step (compensating for the 0.2ms delay between sensor channels). In addition, the parameters such as the number of FFT points (1024 points recommended), the principal component analysis frequency band range (0.1-20Hz), and the weighted least squares regularization factor mentioned in this specification can be adjusted according to the embedded system resources; the reference time window of the "energy in the first 10 seconds" can be extended to 5-30 seconds to be compatible with devices with different thermal time constants, and the specific implementation does not constitute a limitation on the claims.

[0055] The boundary calibration module 200 is used to fuse and correct the output data of the data acquisition module 100 based on the algebraic surface fitting and least squares correction method to obtain a three-dimensional boundary condition mapping set.

[0056] It should be noted that the "three-dimensional boundary condition mapping set" in this application is different from the conventional temperature field cloud map. It specifically refers to a functional matrix that includes the normal temperature gradients and heat flux densities of the six box surfaces. Its spatial resolution reaches 10 mm / pixel, which can support the precise positioning of internal hotspots at the 5 cm level.

[0057] For details, please refer to Figure 3 The fusion correction method of algebraic surface fitting and least squares correction is as follows:

[0058] The internal space of the box is discretized into a three-dimensional grid. Each temperature / stress time series data point is mapped to the corresponding grid node through a coordinate mapping method to generate a spatial sampling vector set.

[0059] Based on the spatial sampling vector set, an algebraic surface of the initial temperature field is constructed using bicubic B-spline basis functions; and stress residual time series is introduced as a penalty term, and energy functional equations are constructed using the Lagrange multiplier method.

[0060] Solve the least squares solution of the energy functional equation and output the coordinate matrix of the surface control vertex; if the fitting residual of any grid node exceeds the preset threshold, trigger the local basis function order adaptive lifting mechanism to iteratively update the control vertex;

[0061] The final control vertex matrix is ​​analyzed, and the temperature gradient distribution function and heat flux density transfer coefficient of the six boundary surfaces of the box are extracted, outputting a three-dimensional boundary condition mapping set. The iterative update method of the local basis function order adaptive lifting mechanism is as follows:

[0062] T100. Calculate the residual deviation δ of this node using the historical residual distribution model. The calculation formula is as follows:

[0063] δ=(|ei | - μσ) / σ (1)

[0064] In formula (1), μ σ Let σ be the historical residual mean, and σ be the standard deviation.

[0065] T200, if δ≥2.0, then the Bernstein polynomial upgradeer is triggered to upgrade the local basis functions from bicubic to biquintic.

[0066] T300, calculate the rank-deficient direction of the original control vertex matrix through SVD decomposition, and based on the pseudo-inverse J of the Jacobian matrix. + The residual vector e generates the increment vector ΔV = J. + ·e, iteratively update the coordinates of the control vertex along the rank-deficient direction.

[0067] In implementation, the spatial discretization operation uses a 10 cm cubic mesh to divide the internal space of the enclosure. For the transformer winding area, an adaptive refinement mechanism is used to reduce the mesh size to 2 cm. Each mesh node is located using three-dimensional spatial coordinates. Temperature and stress data points are associated with neighboring nodes using an inverse squared distance weighted algorithm. The weight of fiber optic temperature measurement data is set to 1.0, the weight of thermocouple data is 0.8, and the weight of infrared data is fixed at 0.7. The association determination radius is set to 5 cm. Algebraic surface fitting uses bicubic B-spline basis functions, with the initial dimension of the control vertex matrix configured as 15×15×3. The residual time series data output by the quartz strain gauge participates in the surface curvature correction with a penalty weight coefficient of 0.7.

[0068] When the fitted residual of a grid node exceeds the median of historical residuals, the adaptive order upgrade mechanism is triggered as follows: First, the residual deviation is calculated. This value is obtained by extracting the average and standard deviation of the residual data of all nodes over the past 30 minutes, and then dividing the difference between the current residual and the average by the standard deviation. If the deviation reaches a threshold of 2.0, the basis function is replaced with a double quintic Bernstein polynomial within a spherical region centered on the node that exceeded the standard deviation and with a radius equal to three times the local grid size. Basis upgrade is achieved by inserting a midpoint node in the parameter domain. The number of newly added control vertices is generated at 40% of the original number of vertices, and the vertex coordinates are calculated according to the de Boor recursive algorithm. The control vertex correction process includes three operations:

[0069] Perform singular value decomposition on the original vertex matrix and extract the correction direction vector corresponding to the minimum singular value;

[0070] Incremental vectors are generated based on the pseudo-inverse of the Jacobian matrix and residual vectors.

[0071] The vertex coordinates are iteratively updated along the correction direction with a step size of 0.3 until the residual drops to 50% of the initial threshold or reaches the upper limit of 5 iterations.

[0072] The boundary condition output analyzes the final control vertex matrix. The calculation process includes: extracting the surface patches of the six outer surfaces of the box; calculating the temperature gradient vector of each surface (the vector product of the surface partial derivative and the surface normal); and determining the heat flux density based on the product of the material's thermal conductivity and the temperature gradient. The output data is in a structured array format, with each element containing a surface number, a three-dimensional gradient vector, and a heat flux density scalar.

[0073] The system performs calibration every 5 minutes. When a temperature gradient exceeding 8 degrees Celsius / cm is detected, the execution frequency is increased to once per minute. The coordinates and deviation values ​​of nodes with excessive residuals are written to the log file in real time. When more than 10% of the nodes exceed the limit within a single cycle, an E407 exception code is triggered.

[0074] The basic mesh size can be adjusted between 5 and 20 cm; the penalty weight coefficient can be configured between 0.2 and 1.5; the order radius coefficient supports adjustment of 2 to 5 times the mesh size; the deviation threshold setting range is recommended to be 1.5 to 3.0; the vertex update step size range is 0.1 to 0.5. The specific parameter values ​​are determined according to the transformer substation model and sensor layout, and are loaded through a JSON format configuration file during implementation.

[0075] Typical applications include monitoring the winding temperature zone of oil-immersed transformers (requiring a 2 cm dense grid) and monitoring the stress of epoxy resin insulation components (with an increased radius of 5 times the grid size).

[0076] During implementation, this module needs to be deployed on a quad-core ARM Cortex-A72 processor platform (clock speed ≥ 1.8GHz), with a dedicated memory area of ​​512 megabytes. Data transfer is completed through shared memory, with single-processing latency controlled within 50 milliseconds. In cases of sensor anomalies, an infrared data reconstruction mode can be enabled, with the specific switching threshold set to automatically activate when the deviation exceeds 5.0 for three consecutive cycles. Log recording uses a circular buffer mechanism, with historical residual data retained for a default of 30 minutes, which can be adjusted from 10 to 60 minutes depending on the system storage capacity. Taking the calibration of the copper busbar connection area of ​​a 110kV oil-immersed substation as an example:

[0077] When the residual at the node at coordinates (0.85m, 1.20m, 0.50m) reaches 4.8℃ (historical residual mean 1.6℃, standard deviation 1.2℃), the calculated deviation is (4.8-1.6) / 1.2 = 2.67 > 2.0, triggering a basis elevation operation within a spherical region with a radius of 60 mm (local mesh size 20 mm). After three iterations, the residual decreases to 1.1℃, and the output surface gradient vector is (6.2, 0.8, 3.4)℃ / cm. If the system is configured with a penalty weight coefficient of 0.5, the elevation radius needs to be adjusted to 4 times the mesh size to adapt to the new fitting characteristics.

[0078] Furthermore, it should be noted that in this application, "historical residuals" specifically refers to dynamic statistics within a sliding time window, distinct from static batch data; the "control vertex increment vector" is generated using a pseudo-inverse method, unlike the traditional gradient descent method, ensuring numerical stability under ill-conditioned matrices. The mesh refinement triggering condition can be adjusted according to the rate of change of the thermal gradient, and the shape of the upgraded neighborhood can be configured as a cube or an ellipsoid; the specific implementation form is not limited.

[0079] The multiphysics solver module 300 is used to input the three-dimensional boundary condition mapping set into the finite element / finite volume solver framework, construct a three-phase fluid heat conduction and convection stratification model including oil, metal, and air, and obtain the internal temperature and stress field profile data of the prefabricated substation. Specifically, taking a certain type of 110kV oil-immersed prefabricated substation as an example, its input data structure includes the temperature gradient distribution and heat flux density parameter set of six surfaces. For example, the data format of the front surface is ["front",(3.2,0.5,1.8),-680W / m 2 When the anisotropic mesh refinement algorithm is initiated, the temperature gradient field of the initial uniform mesh (basic size 5 cm) is first calculated. When the gradient value of the copper busbar connection area (coordinate example: X = 0.85 m, Y = 1.20 m, Z = 0.50 m) is detected to be >10℃ / cm, the octree recursive subdivision mechanism is triggered to refine the mesh of this area to 1.25 cm. At the same time, the geometric curvature of the box is scanned, and curvature adaptive subdivision is performed at the curved corners (such as the edge of the fan shroud)—the area with a curvature radius <15 cm is automatically refined by two levels of mesh. It should be noted that the "anisotropic mesh" in this application is different from the traditional uniform refinement. It performs differentiated subdivision along the temperature gradient direction and the curvature normal, which more accurately captures the boundary layer effect.

[0080] Based on the generated 1.6 million mesh elements, a dual discretization solution was performed: For solid regions such as the transformer winding copper busbar (material number #CU-1100) and the tank steel plate (#Steel-Q235), the finite element method was used to discretize the heat dissipation and conduction equations, with second-order Lagrange polynomials selected for the element shape functions and a 7-point Gaussian integral for the conduction matrix integration; For the transformer oil region (fluid domain #Oil-25#), the finite volume method was used to discretize the Navier-Stokes equations, with a second-order upwind difference scheme configured for the convection term discretization, and the time step size dynamically constrained by CFL conditions (CFL number = 0.8). An overlapping mesh layer was established at the oil-copper interface (e.g., the winding surface), and interface heat exchange conditions were set.

[0081] Where the heat transfer coefficient h tc The default value is 280W / (m²). 2·K), with the radiation coefficient σ∈ being taken as 0.08. In practice, the solution is obtained step by step through an operator splitting coupling mechanism: first, the fluid domain (oil phase) is solved by fixing the solid temperature field, and the convective heat transfer of the oil is calculated using the explicit Euler method; then, the fluid field is frozen to update the solid domain (metal phase), and the heat conduction equation is solved using the implicit Crank-Nickson method. This process is repeated three times in each time step until the residual interfacial heat flux is <1000W / m 2 .

[0082] In the air-side thermal field solution stage, the coupled temperature field between the solid and fluid is used as the boundary input for the SIMPLE algorithm accelerated by multi-grid processing. An extended computational domain containing 2 million cells (1.5 meters outside the box) is constructed, employing a V-cycle three-layer multi-grid structure: the pressure correction equation is solved in the coarse grid layer (40 cm cell size), the velocity field is calculated in the medium grid layer (20 cm), and the turbulence parameters are updated in the fine grid layer (5 cm). The k-epsilon turbulence model constants are adjusted to Cμ = 0.085, C1 = 1.42, and C2 = 1.68 to adapt to the flow characteristics around the box surface. Each iteration executes:

[0083] Solve the momentum equation to update the predicted velocities u* and v*;

[0084] Solving the pressure correction equation yields p′;

[0085] Corrected velocity field

[0086] The convergence criteria for updating turbine kinetic energy k and dissipation rate are based on a velocity field relative residual of <10⁻⁴, and the maximum number of iterations is set to 50. A case study of solving a 220kV transformer substation's wind turbine side shows that multi-grid acceleration reduces the convergence steps from 218 steps in a single-grid system to 76 steps, reducing computation time by 62%.

[0087] The default mesh size is 5 cm for solid domains (adjustable range 2-10 cm), and 3 cm for fluid domains (1-5 cm).

[0088] The CFL number for time step control is 0.8 (allowed 0.5-1.2), and the maximum step size is limited to 0.1 seconds;

[0089] The minimum iteration for interface coupling is 3 times per step (expandable to 10 times), with a residual threshold of 103 x 103 W / m. 2 (Can be reduced to 5×1025×102);

[0090] The constants Cμ for the turbulence model range from 0.08 to 0.09, and C1 ranges from 1.40 to 1.45. The specific values ​​are configured according to the power level of the transformer. For 35kV compact transformers, it is recommended to expand the fluid mesh to 5 cm, and for 220kV large transformers, shrink the solid mesh to 3 cm.

[0091] During implementation, this module needs to be deployed on a workstation equipped with dual NVIDIA Tesla V100 GPUs, with a single timestep solution time controlled within 120 seconds. Output data includes:

[0092] Node temperature distribution: binary floating-point array (single precision), spatial resolution corresponds to minimum grid size.

[0093] Thermal stress tensor field: A 6-component symmetric matrix storage format ([σxx,σyy,σzz,σxy,σyz,σxz][σxx,σyy,σzz,σxy,σyz,σxz]) is used as an example, with a typical node in the winding region (X=0.75m, Y=1.05m, Z=0.35m). The output temperature value is 86.3℃, and the stress tensor is (58.7,-12.4,33.9,4.2,0.8,-5.6) MPa. When the peak stress in the copper busbar region is detected to be >60MPa, warning code E508 is automatically triggered.

[0094] It should be noted that the "operator splitting coupling" in this application differs from traditional strongly coupled algorithms. Its decoupling iteration strategy reduces memory usage by 72%, making it particularly suitable for embedded deployment scenarios. The curvature response threshold in "anisotropic encryption" can be adjusted according to material properties: the curvature radius threshold for copper conductor regions is set to 15 cm, while for epoxy resin components, it is recommended to relax it to 25 cm. The "multi-grid acceleration structure" adopts a custom V-cycle configuration, and the coarse grid layer can be replaced with a geometric multigrid or an algebraic multigrid (AMG), with no specific implementation limitations.

[0095] During implementation, when the residual interfacial heat flux is greater than 5 × 10³ / 5 × 10³ W / m for 5 consecutive steps... 2 When convergence occurs, adaptive time step reduction (step size × 0.5) is initiated. If the multigrid iteration fails to converge after 50 iterations, the algorithm automatically switches to PISO and the sub-relaxation factor is increased to 0.85. In the event of a computational interruption, the most recently valid result is interpolated for output, and the faulty grid ID is recorded in the diagnostic log. It should be noted that in this application, "spatiotemporal resolved data" specifically refers to an output capability with a minimum time step of 0.01 seconds and a spatial resolution of 1 millimeter. The specific resolution can be adjusted to 0.1 seconds / 5 millimeters depending on hardware performance.

[0096] The thermal expansion coefficient of the solid domain and the Boussinesq approximation parameters of the fluid domain can be adjusted according to the actual oil type, and this application does not impose any restrictions. When deploying the module, it is necessary to ensure the spatial consistency between the boundary mapping set and the mesh topology. When the box structure changes, the mesh generation process must be re-executed.

[0097] For further details, please refer to Figure 4The refinement process of the anisotropic mesh refinement algorithm is as follows: Based on the three-dimensional temperature gradient tensor and curvature field output by the thermal field reconstruction module, the second-order temperature derivative matrix at each spatial node is decomposed to obtain the corresponding principal curvature direction vector and curvature intensity feature value. Then, based on the principal curvature direction vector and curvature intensity feature value, an anisotropic metric tensor is constructed. The principal curvature direction vector is used as the preferred direction for geometric stretching of the mesh unit. The unit size along each principal direction is determined according to the mapping relationship between the curvature intensity feature value and the preset local error threshold. The anisotropic refinement algorithm is used to densify the mesh along the principal curvature direction in regions where the temperature gradient increases sharply or the interface curvature changes abruptly, while keeping the mesh coarsened in low curvature regions.

[0098] Specifically, after the three-dimensional boundary condition mapping set output by the boundary calibration module 200 is input (taking a typical 110kV oil-immersed transformer as an example, its front surface gradient distribution is (3.2, 0.5, 1.8)℃ / cm), the anisotropic mesh refinement algorithm is first started: the second derivative matrix of the temperature field is calculated for each node in the initial uniform mesh (basic size 5 cm).

[0099]

[0100] The principal curvature direction vectors v1, v2, v3 and their corresponding eigenvalues ​​λ1, λ2, λ3 (arranged in descending order of absolute value) are obtained by performing eigenvalue decomposition using the ARPACK library. In the copper busbar connection area (coordinate example X = 0.85m, Y = 1.20m, Z = 0.50m), if |λ1| > 15 (corresponding to a temperature gradient > 12℃ / cm), then an anisotropic metric tensor is constructed. The mesh size h1 along the main direction v1 is calculated as h1 = min(1.25, 0.8 | λ1 |) cm (coefficient 0.8 adjustable from 0.5 to 1.2), while in the vertical directions v2, v3, h2 = h3 = 2h1. In practice, mesh refinement is achieved using the LibMesh library, generating a minimum of 1.25 mm cells in the winding region, while maintaining a 5 cm coarse mesh in the epoxy resin insulation region (radius of curvature > 25 cm). It should be noted that the "anisotropic metric tensor" in this application differs from traditional geometry-driven refinement. It integrates both temperature gradient and curvature physical quantities, enabling the refinement direction to accurately track the heat flow path—for example, under an overload condition, the angle between the refinement direction and the current direction at the copper busbar connection point is less than 5°.

[0101] Furthermore, the operator splitting coupling mechanism employs a second-order Strang splitting method: First, in the first half of the time step Δt of each joint solution of the temperature and stress fields, the solid temperature distribution is independently solved using the finite element discrete metal heat conduction equation; then, in the full step Δt, the solid boundary temperature of the first half is incorporated as a Dirichlet condition into the finite volume discrete oil convection equation to calculate the oil phase temperature field in an explicit upwind scheme, while simultaneously evaluating the heat transfer flux at the oil and metal interface through a heat transfer flux matching condition; finally, in the second half of the step, the finite element heat conduction operator is called again to update the metal temperature with an explicit boundary heat flux input. In implementation, the time step Δt is set to 0.1 seconds (CFL number = 0.8), and the solid heat conduction solver is called first within the first half step Δt / 2. The finite element discrete equation is applied to the copper winding region (material #CU-1100): The solid temperature field was obtained by using P2 continuous Galerkin discretization and integrating the stiffness matrix using the 7-point Gaussian formula, and then solving in parallel on the NVIDIA CUDA platform. Then, the oil convection equation is solved over the entire step length Δt: the winding surface temperature is taken into account. As the Dirichlet boundary condition input to the finite volume model, the oil region (#Oil-25#) was discretized using a second-order upwind scheme, and the momentum equation was supplemented with a Boussinesq buoyancy term ρgβ(T-T0). The heat transfer was evaluated at the oil-copper interface (e.g., at coordinates X = 0.82m, Y = 1.35m, Z = 0.48m) using heat flux matching conditions. The contact thermal resistance Rcontact is taken as 2.5 × 10⁻⁶. -4 m 2 K / W (adjustable range 1~5×10) -4 The radiation coefficient σ is set to 0.08. In the second step, the solid-state solver is called again to calculate the interfacial heat flux q. int Update the metal temperature to Neumann boundary conditions.

[0102] Specifically, the lower limit of |λ1| for triggering encryption requiring a feature value threshold is set to 12 (adjustable range 10-20), corresponding to a temperature gradient of approximately 8-15℃ / cm; the coefficient in the h1 calculation formula for the mesh size calculation factor is 0.8, with a recommended value of ≥0.7 for the copper conductor region and 1.0 for the insulation region; the default time step Δt is 0.1 seconds, automatically reduced to 0.05 seconds under overload conditions; and the recommended contact thermal resistance value for the copper-oil interface Rcontact is 2.5×10⁻⁶. -4 m 2 K / W, steel-oil interface value is 5×10 -4 Specific values ​​can be adjusted according to the transformer model: the time step of a 35kV compact box-type substation can be expanded to 0.2 seconds, and the contact thermal resistance of the UHV converter transformer needs to be reduced to 1×10. -4The splitter step size ratio is allocated by default as 1:2:1 (first half: full step: second half), and is adjusted to 1:1:1 when the oil phase temperature change rate is detected to be >5℃ / s.

[0103] In practice, the module was deployed on a dual Tesla V100 GPU platform, with a maximum single time step time of 120 seconds.

[0104] The output includes:

[0105] Node temperature field: stored in binary format, with spatial resolution up to the size of a encrypted mesh (e.g., 1.25 mm for winding areas).

[0106] Thermal stress tensor: Calculated based on the temperature gradient using the generalized Hooke's law, outputting the von Mises equivalent stress field. When the equivalent stress in the copper busbar region > 60 MPa (e.g., a detection value of 67.5 MPa at coordinates X = 0.75 m, Y = 1.05 m, Z = 0.35 m), an E508 warning code is triggered. It should be noted that in this application, "spatiotemporal resolution data" specifically refers to the output capability with a time step of 0.1 seconds and a spatial accuracy of 1 millimeter; the specific resolution can be adjusted according to hardware configuration. The principal curvature direction vector in this application specifically refers to the characteristic direction corresponding to the maximum eigenvalue of the Hessian matrix of the temperature field, distinct from the principal direction of geometric curvature—for example, in the winding bending region, the principal temperature direction coincides with the current path by 89%, while the principal geometric direction deviates by 34°. The heat flux matching condition adopts explicit flux transfer, distinct from implicit energy conservation constraints, reducing interface data exchange by 80%. The time allocation ratio for second-order Strang splitting can be dynamically adjusted according to the operating conditions: 1:2:1 for normal load and 1:1:1 for transient overload.

[0107] In case of anomalies, such as when the feature decomposition iteration exceeds 100 times, the system will automatically switch to power method iteration and output warning code E621; if the interface heat flux residual is >10 for 3 consecutive steps. 4 W / m 2 The time step is halved (Δt ← 0.5Δt) and interface sub-relaxation is activated (factor 0.7). It should be noted that the default encryption depth of "anisotropic refinement" in this application is level 3, which can be extended to level 5 for UHVDC transformers; the specific level is not limited.

[0108] The physical information prediction module 400 is used to jointly train the output data of the multiphysics solution module 300 and the historical running sequence by using a graph neural network and integrating the residual term of the heat transfer conservation equation as a loss function to generate short-term hotspot migration and temperature rise trend prediction results.

[0109] In this module implementation, G(V,E) constructs a discretized representation of the internal space of the transformer substation, where the node set V (e.g., |V| = 12,540 discrete locations) represents key monitoring points, and the edge set E maps the heat conduction paths. The adjacency matrix A = [A ij Define node connection weights, in specific implementation kij is the equivalent thermal conductivity (e.g., 48.7 W / m·K for copper-oil interface), and dij is the node spacing (e.g., 0.02 m between windings). Temperature vector. Store real-time monitoring values ​​(e.g., temperature of node #207) ), which serves as the initial input features for graph convolution. Ensure that physical quantities are directly embedded in the underlying layer of the neural network.

[0110] Specifically, the graph structure construction operation discretizes the internal space of the transformer into 12,540 nodes (taking a 110kV oil-immersed transformer as an example). The spatial coordinates of the nodes are determined through three-dimensional mesh mapping (e.g., the minimum node spacing in the winding region is 2cm, and the maximum spacing in the oil conservator region is 20cm). Adjacency matrix element A ij Calculated based on the actual heat conduction path:

[0111] Copper-to-copper connections (e.g., between winding turns):

[0112] Steel-oil interface (e.g., inner surface of the tank wall):

[0113] Air gap node: (Minimum threshold of 0.1 needs to be set); Initial feature vector Bind real-time sensor data (such as thermocouple at node #207 at 87.2℃, and fiber optic temperature measurement point at #581 at 112.3℃);

[0114] Spatial feature propagation is implemented using a 4-layer graph convolutional architecture, and the inter-layer parameter configurations are shown in Table 2:

[0115] Table 2 shows the convolutional architecture.

[0116]

[0117] In L-layer graph convolution (L=4 is recommended for engineering), the feature update formula for the l-th layer is... middle, Implement normalized adjacency weight (d) i =∑ k A ik (Assuming nodal degree), this design ensures the conservation of heat flux. Weight matrix W (l) ∈R fl×fl+1 and bias b (l)For trainable parameters, nonlinearity is introduced through the LeakyReLU (α = 0.1) activation function σ. This structure enables node i to fuse the thermal state information of its neighborhood N(i) (e.g., nodes within a radius of 0.05m), and can capture lateral thermal diffusion effects in densely wound regions.

[0118] Predicted output The temperature at the next time step is generated through the fully connected layer φ. Key physical constraints are achieved through the heat transfer residual term: in The term is the time derivative, Δ Lap (R t+Δt ) i For a discrete Laplace operator, Δt (typically 30 s) must be less than the system's thermal time constant. κ is the thermal diffusivity of the material (e.g., 7.5 × 10⁻⁸ m² / s for transformer oil), and Δt is the discrete Laplace operator. Lap Quantify the spatial heat flow intensity. This equation forces the prediction results to conform to Fourier's law of thermal conduction, suppressing numerical oscillations without physical basis. In one implementation, the physical residual calculation operates at a prediction step size Δt = 30s (when the load factor is <85%):

[0119] Dynamic assignment of material thermal diffusivity κ:

[0120] Insulating oil zone: κ oil = 7.5 × 10⁻⁸ m² / s (based on 25# mineral oil)

[0121] Copper winding region: κ cu = 1.1 × 10⁻⁴ m² / s (including enameled insulation layer)

[0122] Steel plate structural components: κ steel =3.4×10⁻⁵ m² / s

[0123] In equation (4), the regularization weight λ (recommended 1,200±300) is automatically adjusted through backpropagation, with a higher value (e.g., λ>1,500) used in the high-temperature region to strengthen physical constraints. During training, the gradient is calculated using automatic differentiation. L can reduce the false alarm rate by 83% in areas of abrupt changes in oil flow. In one implementation, the regularization weight λ in the loss function co-optimization mechanism is set differently according to regional characteristics:

[0124] Dense winding area (high hotspot risk): λ = 1,800 ± 200;

[0125] Oil flow channel region (convection-dominated): λ = 1,200 ± 100;

[0126] Enclosure area (heat dissipation boundary): λ = 900 ± 150;

[0127] The training process uses the Adam optimizer (learning rate 2×10⁻⁴) and converges in the following typical scenarios:

[0128] Steady-state temperature rise prediction: convergence in 1500 steps, Ldata < 0.3℃;

[0129] Short-circuit transient prediction: convergence in 5000 steps, Lphys < 15W;

[0130] This technical solution requires special handling in the following extreme scenarios:

[0131] Low-temperature start-up (<-25℃): Freeze the weights of layers 1-2 of the convolutional layer;

[0132] Lightning overvoltage: Switch to transient sub-model;

[0133] Oil chromatography anomaly: Injection of faulty heat source vector;

[0134] It should be noted that any implementation of box-type transformer temperature prediction based on physical constraint graph neural networks falls within the scope of protection of this patent. The "discrete Laplace operator" Δ in this application... Lap Unlike the traditional finite difference method, its innovation lies in reusing the adjacency weights A of the graph structure. ij Constructing a physically consistent gradient approximation improves computational accuracy by 38% at the oil-copper interface; "spatial variation" requires real-time matching of node material types (indexed by MaterialID) to avoid distortion of conservation laws caused by homogenization assumptions; the default number of graph convolution layers L is 3 layers, which can be increased to 5 layers for large box-type substations above 220kV (dl needs to be expanded to 128 simultaneously). The specific depth is dynamically adjusted according to the convergence of the validation set, and this application does not impose a mandatory limit.

[0135] The backpropagation process uses the Adam optimizer with a learning rate η = 10. -3 (The decay strategy is 0.8 × 0.8 every 20 cycles), gradient pruning threshold ||g|| 2 <10. Training data must cover typical operating conditions with a load rate of 60%–120%, and each epoch contains 200 time step sequences. During implementation, if the copper busbar area Fi > 5K / s is detected (indicating physical mismatch), the λ value will be automatically increased by 30%, and training will be restarted. Hardware deployment requirements: ≥16GB of GPU memory (NVIDIA Tesla V100), and real-time prediction latency <50ms.

[0136] Taking a certain type of 110kV oil-immersed box-type substation as an example, its input data structure is a temperature vector T containing 50,000 nodes. t ∈R 50000 and the adjacency matrix A∈R 50000×50000 Adjacency matrix element A ijThe assignment of values ​​must strictly follow the inverse ratio principle of thermal conductivity and distance: the formula for calculating the connection weight between the copper busbar node (material ID#CU-1100) and the adjacent oil node (#Oil-25#) is as follows: (Thermal conductivity is in W / m·K, distance is in meters), while the weight for the copper-air interface (#Air) is A. ij =0.05. Initial feature vectors of nodes during graph structure construction. Expanded to 8 dimensions: Where the spatial coordinates (x) i ,y i ,z i Normalized to the [0,1] interval, the physical property parameter κ i (thermal conductivity), ρ i (density), c p,i (Specific heat capacity) is injected based on the material library.

[0137] Graph convolution operations require three layers of processing (layer L=3 is the default configuration). In each layer l, the message passing operation follows this process: First, the node degree is calculated, and a normalization coefficient is applied to the copper busbar connection points (such as nodes with high connectivity (node ​​degree>20) at coordinates (0.85, 1.20, 0.50)m to prevent gradient explosion; then, a trainable weight matrix is ​​loaded (initialized using a He normal distribution with a standard deviation of 2 / 642 / 64); next, neighborhood information is aggregated: within a spherical neighborhood N(i) with a radius of 5 cm, the following steps are performed: Finally, the ReLU activation function is applied and a bias b is added. (l) (Default zero initialization) Output After three convolutional layers, the final features The predicted temperature is mapped through an affine layer (time step Δt defaults to 1.0 seconds, adjustable range 0.1-5.0 seconds).

[0138] After the predicted temperature field is generated, the heat transfer conservation residual is calculated by automatic differentiation. It should be noted that the discrete Laplace operator requires that the weights be positively correlated with the thermal conductivity and prohibits the use of simplified forms such as the reciprocal of the geometric distance.

[0139] In the composite loss function, the data fitting term is based on a historical database, covering typical operating conditions with load rates ranging from 60% to 120% (including 15 sets of overload transient data). The weighting coefficient λ for the physical regularization term is set as follows: Copper conductor region (MaterialID#CU) λ = 0.5; Transformer oil region (#Oil) λ = 0.3; Air region (#Air) λ = 0.1.

[0140] Backpropagation uses the Adam optimizer, and the physical residual gradient is calculated by automatic differentiation of the Jacobian matrix. When the copper bus node (e.g., coordinates (0.85, 1.20, 0.50) m) is detected for 5 consecutive iterations, the copper area λ is automatically increased to 0.7 and training is restarted.

[0141] In implementation, this module is deployed on a dual NVIDIA Tesla V100 GPU platform with a minimum of 16GB of VRAM. The training phase requires 200,000 time-step sequences (approximately 56 hours), with real-time prediction latency ≤50ms. The output prediction results are correlated with system behavior. It should be noted that "physical residual regularization" employs a partitioned weighting strategy: strong constraints in the copper region (λ = 0.5) ensure hotspot prediction accuracy, while weak constraints in the air region (λ = 0.1) avoid overfitting to background noise. If GPU resources are limited during real-time prediction, a simplified mode can be enabled: freezing the convolutional layer weights W(l) and only updating the affine layer parameters, reducing latency to 20ms (accuracy loss <8%). During implementation, the module implementation must ensure that the training data covers:

[0142] Ambient temperature range: -25℃ to +45℃;

[0143] Load fluctuation condition: Load step change ≥20% within 30 seconds;

[0144] Fault simulation scenarios: localized heat dissipation failure, oil pump jamming; minimum historical database size 50GB, sampling frequency 1Hz. The root mean square error (RMSE) between predicted and measured temperatures must be ≤1.5℃, and the copper busbar hotspot location deviation <3cm. Specific training cycles, data storage formats, and hardware acceleration architectures can be adjusted according to the substation scale; this application does not impose limitations.

[0145] The hierarchical scheduling module 500 is used to construct a second-order cone programming model containing fan speed and oil pump flow rate using model predictive control methods. It combines the prediction results from the physical information prediction module 400 with equipment operating condition constraints to solve the model in real time and output multi-equipment coordinated control commands. For details, please refer to [link to relevant documentation]. Figure 5 The model predictive control process of the hierarchical scheduling module 500 includes:

[0146] Obtain the short-term three-dimensional temperature field evolution sequence and node-level heat flux requirements, define the fan speed and oil pump flow rate as time-series decision variables, and construct a multi-step state transition equation;

[0147] During each control period, a second-order cone constraint method is used to limit the heat dissipation coupling capacity of the fan and pump to the predicted heat flux demand. The second-order cone constraint method normalizes the air volume and flow rate through an efficiency coefficient matrix.

[0148] The temperature over-limit penalty function and the weighted sum of equipment energy consumption are combined to form the MPC objective function, and together with the upper and lower limits of variables, the consistency cone constraint and the equipment dynamic response constraint, a second-order cone programming model is formed.

[0149] The cone programming problem is decomposed into a fan subproblem and a pump subproblem by the distributed alternating direction multiplier method, and parallel iterative calculations are performed in the embedded solver to update the control sequence of each device.

[0150] After parallel iteration convergence, the optimal fan speed and oil pump flow rate decision for the first moment are extracted and sent to the field execution unit via the real-time communication bus. The specific implementation is as follows:

[0151] Step S101: Obtain the short-term three-dimensional temperature field evolution sequence and node-level heat flux demand output by the physical information prediction module 400. The temperature field evolution sequence has a temporal resolution of 1 second and a spatial resolution of 1.25 millimeters. The heat flux demand includes data within a typical coordinate range of the winding region (e.g., the spatial range of the winding region is defined as a cubic region of X = 0.80-0.90 m, Y = 1.20-1.30 m, and Z = 0.50-0.55 m). It should be noted that the "node-level heat flux demand" in this application differs from traditional regional mean parameters. Its essence is a set of spatially distributed parameters obtained through high-resolution temperature field inversion, which can accurately characterize the local heat load intensity.

[0152] Step S102: Define the time-series decision variables for fan speed and oil pump flow rate. The fan speed range is limited to 800-2000 rpm (typical baseline 1500 rpm), and the oil pump flow rate range is limited to 20-100 liters / minute (typical baseline 60 liters / minute). These ranges can be adjusted according to the substation cooling system model. For example, for a compact 35kV prefabricated substation, the upper limit of fan speed can be reduced to 1800 rpm; for a large 220kV prefabricated substation, the upper limit of oil pump flow rate can be increased to 120 liters / minute. Specific limits are determined by the equipment technical specifications.

[0153] Step S103: Construct a multi-step state transition equation and implement second-order cone constraints. The state transition equation requires the injection of a spatially relevant equipment efficiency coefficient matrix (e.g., 0.07 for the winding region and 0.03 for the insulation region). The second-order cone constraints handle airflow and flow parameters using normalization coefficients (the normalization coefficient for the fan is 1 / 2000 of the reciprocal of the maximum speed, and the normalization coefficient for the oil pump is 1 / 100 of the reciprocal of the maximum flow rate). It should be noted that the "efficiency coefficient matrix" in this application differs from the traditional homogenized parameter model. Its essence is a spatially correlated parameter set established based on equipment thermal performance tests. For example, in the fan factory test, airflow-heat dissipation efficiency curve data from no fewer than 10 operating points must be collected; in the oil pump commissioning stage, heat exchange characteristic data from no fewer than 8 flow points must be obtained.

[0154] Step S104: The synthetic model predicts the control objective function and forms a second-order cone programming model. The objective function is composed of a temperature over-limit penalty function and a weighted sum of equipment energy consumption (the temperature penalty function threshold is set to 130℃, the fan energy consumption weight is 0.003, and the oil pump energy consumption weight is 0.001). When the system detects a historical temperature record exceeding 125℃, the temperature penalty function weight is automatically increased to 1.6 times the baseline value. The above temperature thresholds can be adjusted according to the insulation material grade; for example, a 120℃ threshold is used for Class A insulation materials, and a 150℃ threshold is used for Class H insulation materials. The specific values ​​are determined by the transformer design specifications.

[0155] Step S105: Solve the optimization problem using the distributed alternating direction multiplier method. The cone programming problem is decomposed into a fan control subtask (optimizing the speed sequence while maintaining a fixed oil pump flow rate) and an oil pump control subtask (optimizing the flow sequence while maintaining a fixed fan speed). Dual-thread parallel computation is initiated in an embedded solver (ECOS 2.0 solver is recommended). Boundary data exchange is performed in each iteration (exchange cycle 100 milliseconds). Convergence is determined by a duality gap of less than 0.01 or 50 iterations. It should be noted that the "distributed alternating direction multiplier method" in this application differs from traditional centralized solutions. Its essence is to decompose a 10,000-dimensional optimization problem into two 5,000-dimensional subproblems through a device decoupling strategy. For example, when implemented on the TI AM5728 embedded platform, the solution time is reduced from 120 milliseconds to 45 milliseconds.

[0156] Step S106: Extract the optimal control command and issue it for execution. After parallel iteration convergence, the optimal decision at the first moment (i.e., the control quantity at the 30th second) is output, and the execution command is issued via the CAN 2.0B bus protocol: the fan control command uses a three-byte frame format (first byte command code 0x3A, last two bytes speed value, e.g., 1850 rpm encoded as 0x073A); the oil pump control command uses a three-byte frame format (first byte command code 0x5B, last two bytes flow rate value, e.g., 85 liters / minute encoded as 0x0352). The above communication protocol can be adjusted according to the fieldbus type. For example, when using the Modbus RTU protocol, the fan control function code changes to 0x06, and the oil pump control function code changes to 0x10. The specific communication protocol is determined by the equipment interface standard.

[0157] Regarding dynamic response constraints, the fan speed change rate is limited to no more than 100 rpm, and the oil pump flow rate change rate is limited to no more than 10 liters / minute / second. When the predictive model detects that the winding temperature exceeds the 125°C threshold, an adaptive control parameter mechanism is triggered (the temperature penalty function weight is increased by 30%-100%). If the optimization iteration exceeds the limits, the system switches to a backup PID controller (proportional coefficient set to 0.8±0.1). If the bus communication packet loss rate exceeds 10%, the historical best control sequence is activated. The equipment efficiency coefficient matrix needs to be calibrated periodically (a calibration cycle of 3 months is recommended). The calibration process includes detecting the dust accumulation coefficient of the fan blades (analyzing dust coverage through image recognition) and measuring the oil viscosity parameters (obtained through an online viscosity sensor). The specific calibration cycle can be adjusted according to the operating environment.

[0158] While specific embodiments of the present invention have been described above, those skilled in the art should understand that these specific embodiments are merely illustrative. Those skilled in the art can omit, substitute, and modify the details of the above methods and systems in various ways without departing from the principles and essence of the present invention. For example, combining the above method steps to perform substantially the same function and achieve substantially the same result according to substantially the same method falls within the scope of the present invention. Therefore, the scope of the present invention is defined only by the appended claims.

Claims

1. An intelligent self-monitoring temperature management system for prefabricated substations; characterized in that: include: The data acquisition module is used to acquire FBG fiber reflection spectrum, quartz strain gauge voltage signal, thermocouple voltage and infrared array radiation intensity data in parallel. It extracts the time spectrum lines of each channel through multi-channel wavelet packet time-frequency decomposition, and uses an adaptive spectral line separation algorithm to separate the temperature principal component and stress residual to obtain the temperature time series and stress residual time series. The boundary calibration module is used to fuse and correct the output data of the data acquisition module based on algebraic surface fitting and least squares correction methods to obtain a three-dimensional boundary condition mapping set. The multiphysics solution module is used to input the three-dimensional boundary condition mapping set into the finite element / finite volume solution framework, construct a three-phase fluid heat conduction and convection stratification model including oil, metal and air, and obtain the internal temperature field and stress field profile data of the box-type substation. The physical information prediction module is used to jointly train the output data of the multiphysics solution module and the historical running sequence by using a graph neural network and integrating the residual term of the heat transfer conservation equation as a loss function to generate short-term hotspot migration and temperature rise trend prediction results. The hierarchical scheduling module is used to construct a second-order cone programming model containing fan speed and oil pump flow rate through model predictive control methods. It combines the prediction results of the physical information prediction module with equipment operating condition constraints to solve the problem in real time and output multi-equipment coordinated control commands.

2. The intelligent self-monitoring temperature management system for prefabricated substations according to claim 1, characterized in that: The extraction process of the time spectrum lines of each channel by the multi-channel wavelet packet time-frequency decomposition includes: Step 1: Allocate a 0.1-20Hz low-frequency channel for the FBG fiber reflection spectrum, a 0.1-50Hz wideband channel for the quartz strain gauge voltage signal, and a 0.5-100Hz mid-to-high frequency channel for the thermocouple voltage and infrared array radiation intensity, respectively. Step 2: Based on the frequency band range allocated to each channel, select the wavelet basis function. The selection method is as follows: the low-frequency channel is loaded with the Db6 wavelet packet kernel, the wideband channel is loaded with the Sym8 wavelet packet kernel, and the mid-to-high frequency channel is loaded with the Coif5 wavelet packet kernel. Step 3: Perform the following operations for each channel: 1) Perform 8-level adaptive decomposition using the wavelet packet kernel to generate a tree of detail coefficients and approximation coefficients; 2) Based on a preset channel-specific energy threshold, hard threshold filtering is applied to the detail coefficients and approximation coefficient trees to suppress power frequency harmonics and environmental noise; 3) The filtered time-domain waveform is recovered using the inverse wavelet packet reconstruction algorithm; 4) Perform a short-time Fourier transform on the recovered filtered time-domain waveform to extract the dominant spectral lines in the time-frequency matrix; Step 4: If in-phase energy peaks with a frequency difference ≤ 2Hz are detected in the time spectrum lines of the FBG base low-frequency channel and the quartz strain gauge broadband channel, then trigger the spectral line cross-verification marker. Step 5: Calculate the energy percentage of spectral lines based on the marked channels. Classify the spectral lines with an energy percentage ≥85% in the FBG channel as the temperature principal component temporary storage region, and classify the spectral lines with an energy percentage ≥70% in the quartz strain gauge channel as the stress residual temporary storage region.

3. The intelligent self-monitoring temperature management system for prefabricated substations according to claim 1, characterized in that: The adaptive spectral line separation algorithm separates the temperature principal component and the stress residual as follows: S1. Extract the fundamental frequency, third harmonic amplitude and phase difference parameters of historical FBG reflection spectrum data using the FFT algorithm, and construct a temperature characteristic spectral template library; S2. Receive the time-spectrum lines output by the multi-channel wavelet packet decomposition in real time, and perform the following operations for each spectral line: s200. Match the temperature feature template library using the Pearson correlation coefficient calculator. If the correlation coefficient is ≥0.9, mark it as a candidate temperature spectral line. s201. Calculate the derivative of the current spectral line energy with the energy of the previous 10 seconds using the energy change rate differential method. If the derivative value exceeds ±5% / s, it is marked as an abnormal spectral line. S3. For all quartz strain gauge channel spectral lines that are not marked as candidate temperature spectral lines, group them according to the principle of frequency difference ≤ 2Hz by using spectral line clustering analysis units, and calculate the stress contribution weight of each group based on the energy ratio of the spectral lines within the group by using the weighted least squares method. S4. Input the candidate temperature spectral lines into the principal component analyzer to extract the first three principal component vectors; and map the principal component vectors to the original FBG spectral sequence through the vector projection algorithm to reconstruct the temperature principal component time series. S5. From the original quartz strain gauge signal, the spectral components corresponding to the stress contribution weights are subtracted through convolution inversion operation to generate the stress residual time series.

4. The intelligent self-monitoring temperature management system for prefabricated substations according to claim 1, characterized in that: The fusion correction method of algebraic surface fitting and least squares correction is as follows: The internal space of the box is discretized into a three-dimensional grid. Each temperature / stress time series data point is mapped to the corresponding grid node through a coordinate mapping method to generate a spatial sampling vector set. Based on the spatial sampling vector set, an algebraic surface of the initial temperature field is constructed using bicubic B-spline basis functions. Furthermore, stress residual time series is introduced as a penalty term, and energy functional equations are constructed using the Lagrange multiplier method; Solve the least squares solution of the energy functional equation and output the coordinate matrix of the surface control vertex; if the fitting residual of any grid node exceeds the preset threshold, trigger the local basis function order adaptive lifting mechanism to iteratively update the control vertex; The final control vertex matrix is ​​analyzed, and the temperature gradient distribution function and heat flux density transfer coefficient of the six boundary surfaces of the box are extracted. The output is a three-dimensional boundary condition mapping set.

5. The intelligent self-monitoring temperature management system for a prefabricated substation according to claim 4, characterized in that: The iterative update method of the local basis function order adaptive lifting mechanism is as follows: T100. Calculate the residual deviation δ of this node using the historical residual distribution model. The calculation formula is as follows: δ=(|e i |-m σ ) / s (1) In formula (1), e i Let μ be the fitting residual of the i-th grid node. σ Let σ be the historical residual mean, and σ be the standard deviation. T200, if δ≥2.0, then the Bernstein polynomial upgradeer is triggered to upgrade the local basis functions from bicubic to biquintic. T300, calculate the rank-deficient direction of the original control vertex matrix through SVD decomposition, and based on the pseudo-inverse J of the Jacobian matrix. + The residual vector e generates the increment vector ΔV = J. + ·e, iteratively update the coordinates of the control vertex along the rank-deficient direction.

6. The intelligent self-monitoring temperature management system for prefabricated substations according to claim 1, characterized in that: The solution process of the multiphysics solution module includes: The three-dimensional boundary condition mapping set is input into the anisotropic mesh refinement algorithm to generate a non-uniform subdivision mesh based on the temperature gradient and geometric curvature. Based on the non-uniform subdivision grid, a dual discretization solution mechanism is adopted, using the finite element method to discretize the thermal conductivity equation of solid metal and the finite volume method to discretize the oil convection equation. Within each time step, an operator splitting coupling mechanism is used to alternately perform solid-fluid heat exchange calculations to maintain consistent heat flux at the oil-metal interface. The solid-fluid coupling results and air side boundary conditions are input into the SIMPLE algorithm accelerated by multigrid to iteratively solve the external thermal field of air convection. Output nodal temperature distribution data with spatiotemporal resolution and corresponding thermal stress tensor field data.

7. The intelligent self-monitoring temperature management system for prefabricated substations according to claim 6, characterized in that: The refinement process of the anisotropic mesh refinement algorithm is as follows: Based on the three-dimensional temperature gradient tensor and curvature field output by the thermal field reconstruction module, the second-order temperature derivative matrix at each spatial node is decomposed to obtain the corresponding principal curvature direction vector and curvature intensity feature value. Then, based on the principal curvature direction vector and curvature intensity feature value, an anisotropic metric tensor is constructed. The principal curvature direction vector is used as the preferred direction for geometric stretching of the mesh unit. The unit size along each principal direction is determined according to the mapping relationship between the curvature intensity feature value and the preset local error threshold. The anisotropic refinement algorithm is used to densify the mesh along the principal curvature direction in regions where the temperature gradient increases sharply or the interface curvature changes abruptly, and to keep the mesh coarsened in low curvature regions.

8. The intelligent self-monitoring temperature management system for a prefabricated substation according to claim 6, characterized in that: The operator splitting coupling mechanism adopts a second-order Strang splitting method: First, in the first half of the periodic time step Δt of each temperature field and stress field joint solution, the solid temperature distribution is solved independently through the finite element discrete metal heat conduction equation; then, in the full step Δt, the solid boundary temperature of the first half step is incorporated into the finite volume discrete oil convection equation as a Dirichlet condition, and the oil phase temperature field is calculated in an explicit upwind scheme. At the same time, the heat transfer flux at the oil and metal interface is evaluated through the heat transfer flux matching condition; finally, in the second half step, the finite element heat conduction operator is called again to update the metal temperature with explicit boundary heat flux input.

9. The intelligent self-monitoring temperature management system for a prefabricated substation according to claim 1, characterized in that: The short-term prediction process of the physical information prediction module includes: The nodal temperature vector output by the third multiphysics solver module and the corresponding adjacency matrix A = [A ij Construct a graph G(V,E) and set the initial features of each node. Set its corresponding temperature value In each graph convolutional layer l, message passing and node updates are performed, expressed as follows: In formula (2), N(i) is the set of neighbors of node i, and d i and d j W represents the degree values ​​of node i and node j, respectively. (l) and b (l) These are the weight matrix and bias vector of the graph convolution at layer l, respectively; A ij σ represents the connection weight between node i and node j in the graph; σ is the activation function. Let j be the feature vector of node j in the l-th layer; After L layers of graph convolution, the final features are... Input an affine mapping function to obtain the temperature prediction for the next time step. The heat transfer conservation residual is calculated using automatic differentiation based on finite difference, and the calculation formula is as follows: In formula (3), Δt is the predicted time step; κ represents the thermal diffusivity of transformer oil, metal, and air; Δ Lap (R t+Δt ) i This indicates that the temperature at node i is approximated by a discrete Laplace based on the adjacency matrix, which is used to reflect the heat flow generated by the temperature difference between adjacent nodes; Predicted deviation term With residual A composite loss function is constructed together, expressed as follows: In formula (4), λ is the residual regularization weight, which is applied to all W using the backpropagation algorithm. (l) and b (l) After iterative updates and joint fitting of spatial topology and physical conservation constraints, the results of short-term hotspot migration and temperature rise trend prediction are output.

10. The intelligent self-monitoring temperature management system for a prefabricated substation according to claim 1, characterized in that: The model prediction control process of the hierarchical scheduling module includes: Obtain the short-term three-dimensional temperature field evolution sequence and node-level heat flux requirements, define the fan speed and oil pump flow rate as time-series decision variables, and construct a multi-step state transition equation; During each control period, a second-order cone constraint method is used to limit the heat dissipation coupling capacity of the fan and pump to the predicted heat flux demand. The second-order cone constraint method normalizes the air volume and flow rate through an efficiency coefficient matrix. The temperature over-limit penalty function and the weighted sum of equipment energy consumption are combined to form the MPC objective function, and together with the upper and lower limits of variables, the consistency cone constraint and the equipment dynamic response constraint, a second-order cone programming model is formed. The cone programming problem is decomposed into a fan subproblem and a pump subproblem by the distributed alternating direction multiplier method, and parallel iterative calculations are performed in the embedded solver to update the control sequence of each device. After parallel iteration convergence, the optimal fan speed and oil pump flow rate decision at the first moment are extracted and sent to the field execution unit through the real-time communication bus.

Images