Optimization method of magnetic flux distribution in composite core inductors

Abstract

This invention discloses a method for optimizing the magnetic flux distribution of a composite magnetic core inductor. The method comprises: constructing an initial finite element model and initializing a pseudo-density distribution field within the model's design domain; constructing and pre-training a graph neural network proxy model for predicting electromagnetic performance based on the initial finite element model; performing an iterative optimization process based on the pseudo-density distribution field until convergence; determining the converged pseudo-density distribution field as the final pseudo-density distribution field; and generating a final structural geometric model based on the pseudo-density distribution field. This method can efficiently and automatically generate an optimal magnetic core structure that takes into account multiple physical properties, improving inductor performance and R&D efficiency.

Description

Technical Field

[0001] The invention relates to a simulation optimization method, in particular to a magnetic flux distribution optimization method for a composite magnetic core inductor. Background Art

[0002] The rapid development of power electronics technology, particularly in cutting-edge fields such as new energy vehicles, data centers, 5G communications, and renewable energy, has placed unprecedentedly stringent demands on high efficiency, high power density, and miniaturization in power conversion systems. As a core component for power conversion and storage, magnetic core inductors often become a key bottleneck restricting overall system performance improvement and size reduction. Composite core inductors, by combining materials with different magnetic properties, offer an effective solution to the conflict between magnetic saturation under high DC bias and core losses under high-frequency AC inductors. However, the splicing of different materials and complex geometric structures inevitably lead to severe non-uniformity in the magnetic flux distribution within the core, which in turn causes a series of problems such as local saturation, concentrated hot spots, and dramatically increased losses. Therefore, precise and efficient optimization of the magnetic flux distribution in composite core inductors has become a core technology driving the advancement of high-frequency magnetic components and even the entire power electronics system. It has significant theoretical significance and engineering application value for improving energy efficiency and promoting device miniaturization.

[0003] Currently, the design and optimization of composite core inductors rely heavily on finite element analysis (FEA)-based simulation software, such as ANSYS Maxwell and COMSOL. The typical design process begins with engineers building an initial inductor geometry model based on their personal experience and theoretical knowledge. For example, this can be achieved by placing air gaps of specific shapes and sizes in the center or side legs of an E-shaped ferrite core, or by simply splicing it with materials such as iron powder cores. Subsequently, specific electrical (such as DC bias current and AC excitation current) and thermal boundary conditions are set in the simulation software to perform a high-fidelity electromagnetic field or thermal-magnetic coupling simulation. After the simulation is complete, the results are visualized and analyzed using post-processing software to visualize magnetic flux density contours, loss distribution, and temperature profiles. If localized magnetic flux crowding or hot spots are detected, engineers return to the geometric modeling phase to manually fine-tune the air gap shape, chamfer radius, or material size, and then restart the simulation. This modeling-simulation-analysis-modification cycle is repeated until a satisfactory design is achieved. In recent years, in order to improve efficiency, some studies have begun to try to combine finite element simulation with some standard optimization algorithms (such as genetic algorithms and particle swarm algorithms) to automatically optimize the parameters of a few preset geometric parameters (such as air gap length and core column width).

[0004] Despite the widespread adoption of these finite element simulation-based design methods, they still face deep, interrelated technical conflicts when addressing complex optimization problems such as high-performance composite core inductors. This makes the design process not only extremely time-consuming but also often fails to achieve true optimality. Existing core optimization methods fail to address optimization objectives at different frequencies and struggle to balance simulation efficiency with physical accuracy, resulting in long design cycles and poor performance. Summary of the Invention

[0005] The purpose of the invention is to provide a method for optimizing the magnetic flux distribution of a composite magnetic core inductor to solve the above-mentioned problems existing in the prior art.

[0006] The technical solution is a method for optimizing the magnetic flux distribution of a composite core inductor, which is used to generate the final structural geometry model, including:

[0007] Construct an initial finite element model and initialize the pseudo-density distribution field within the model design domain;

[0008] Based on the initial finite element model, a proxy model for predicting electromagnetic performance is generated and pre-trained; based on this, an iterative optimization process is performed until the pseudo-density distribution field converges;

[0009] The converged pseudo-density distribution field is determined as the final pseudo-density distribution field; and the final structural geometric model is generated accordingly.

[0010] Beneficial effect: the present invention can efficiently and automatically generate a magnetic core structure that takes into account multiple physical properties, thereby improving inductance performance and R&D efficiency. BRIEF DESCRIPTION OF THE DRAWINGS

[0011] Figure 1 A flowchart of the steps of a method for optimizing the magnetic flux distribution of a composite magnetic core inductor provided in an embodiment of the present application.

[0012] Figure 2 A flowchart of the steps for performing an iterative optimization process provided in an embodiment of the present application.

[0013] Figure 3 A flowchart of the steps for generating the final structural geometric model provided in an embodiment of the present application.

[0014] Figure 4 A flowchart of the steps for calculating the rotation effect provided in an embodiment of the present application. DETAILED DESCRIPTION

[0015] In order to enable those skilled in the art to better understand the solutions of the present invention, the technical solutions in the embodiments of the present invention will be clearly and completely described below in conjunction with the drawings in the embodiments of the present invention. Obviously, the embodiments described are only part of the embodiments of the present invention, not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without making creative efforts should fall within the scope of protection of the present invention.

[0016] It should be noted that the terms "include" and "have" and any variations thereof are intended to cover non-exclusive inclusions. For example, a process, method, system, product or apparatus that includes a series of steps or units is not necessarily limited to those steps or units explicitly listed, but may include other steps or units that are not explicitly listed or are inherent to these processes, methods, products or apparatus.

[0017] Research has revealed that existing topology optimization methods have inherent single-mode limitations when addressing multiple physical objectives. Specifically, a high-performance inductor must simultaneously meet two conflicting physical requirements: First, under strong DC bias currents, a low-reluctance, unobstructed magnetic path is required to avoid magnetic saturation, often requiring a relatively complex structure to adequately guide the magnetic flux. Second, under high-frequency AC ripple currents, the structure must be as smooth and simple as possible to suppress eddy currents caused by rapidly varying magnetic fields, thereby reducing core losses. Existing topology optimization algorithms typically use a fixed set of optimization rules (e.g., a fixed penalty factor) to guide material distribution. This fixed set of rules cannot intelligently distinguish between DC and high-frequency problems, and therefore can only seek a mediocre, suboptimal compromise between the two conflicting objectives. It is unable to boldly explore complex structures when optimizing the DC magnetic path, while cautiously preferring smooth boundaries when suppressing AC losses. This one-size-fits-all optimization approach is the fundamental reason why the final design fails to achieve optimal performance across all aspects.

[0018] Furthermore, there's an irreconcilable conflict between simulation efficiency and physical accuracy during the optimization process. Topology optimization requires thousands of iterative evaluations, and directly embedding high-precision finite element simulations into the optimization loop would be prohibitively time-consuming for engineering practice. To address this, researchers have explored the use of computationally inexpensive surrogate models. However, conventional surrogate models (such as polynomial regression and general neural networks) are black-box models that only learn superficial input-output relationships from data and are completely unaware of the underlying electromagnetic physics. When the optimization algorithm generates a completely new topology that differs significantly from the training data, these black-box models can easily make absurd predictions that violate physical laws, leading to misleading optimization processes and convergence to erroneous or even unphysical solutions. Therefore, existing technologies lack a surrogate model specifically designed for topology optimization that can achieve millisecond-level predictions while consistently ensuring that its predictions conform to the underlying electromagnetic physics. This constitutes a core technical bottleneck restricting optimization efficiency and reliability.

[0019] like Figure 1 As shown, a method for optimizing the magnetic flux distribution of a composite core inductor is proposed, comprising the following steps:

[0020] Construct an initial finite element model and initialize the pseudo-density distribution field within the model design domain.

[0021] In this embodiment, the initial finite element model is used to transform the actual physical problem (including geometry, material properties, and physical boundary conditions) into a digital model in a discretized mathematical form. The pseudo-density distribution field is a scalar field defined at each discrete element within the design domain, with a value range of [0, 1]; 0 typically represents air (or a material with extremely low magnetic permeability), 1 represents a solid core material, and values ​​between 0 and 1 represent an intermediate transition material. Specifically, pre-set design specifications are read, including the inductor's macroscopic geometric dimensions (e.g., the height and width of the E-shaped core), a material library (e.g., the BH curve and loss factor of N87 ferrite and Sendust powder core), and operating conditions (e.g., a DC bias of 20A and an AC ripple of 2A at 100kHz). A 3D model is created based on the geometric dimensions using commercial or in-house CAE software (e.g., ANSYS Maxwell, COMSOL Multiphysics). Within this model, a design domain where structural morphology variations are permitted is manually or automatically defined, such as specific corners or joint regions of the core. The entire model is meshed using tetrahedral or hexahedral grids to generate initial grid data. A uniform initial pseudo-density value, such as 0.5, is assigned to all grid cells within the design domain to initialize the pseudo-density distribution field. This embodiment transforms a specific engineering design problem into a structured mathematical optimization problem amenable to computer algorithms. By establishing a finite element model, complex electromagnetic field partial differential equations can be solved on a discrete grid.

[0022] In some embodiments, the initial pseudo-density distribution field does not need to be uniform. For example, based on the engineer's prior knowledge, higher initial density values ​​can be preset in some areas, while lower values ​​can be preset in other areas to help accelerate the convergence of the optimization.

[0023] Based on the initial finite element model, a surrogate model for predicting electromagnetic performance is generated and pre-trained.

[0024] In this embodiment, a surrogate model (SurrogateModel) is used to simulate computationally expensive high-fidelity simulations (such as finite element analysis). Its input is a pseudo-density distribution field, and its output is the corresponding magnetic field distribution and losses. This trades a small amount of expensive computation for a large amount of subsequent, inexpensive computation. Specifically, a data generation module is invoked. This module, based on the initial finite element model, performs a series of simulations to obtain input-output data pairs, forming an initial training set. A neural network is constructed and trained using the initial training set until its prediction error meets preset engineering requirements (for example, the prediction error of key performance indicators is less than 5%), thereby obtaining a surrogate model. Performing a complete high-frequency electromagnetic finite element simulation typically takes minutes to hours. Topology optimization algorithms, on the other hand, may require thousands or even tens of thousands of iterations. Directly embedding high-fidelity simulations into the optimization loop is not feasible. Once the surrogate model of this embodiment is trained, a single prediction can take only milliseconds, reducing the entire optimization process from months to hours or days, making it practical for engineering applications.

[0025] Based on the surrogate model, an iterative optimization process is performed until the pseudo-density distribution field converges.

[0026] In this embodiment, in each iterative loop, the current design solution (i.e., the pseudo-density distribution field) is evaluated and modified slightly toward a more optimal solution. Once the iterative process reaches a stable state, or converges, subsequent iterations no longer significantly change the design solution, or the preset performance target has been achieved. Specifically, at each step in the loop, a trained proxy model is invoked to perform performance predictions on the current pseudo-density distribution field. The topology optimization algorithm calculates the design sensitivity based on the predictions and updates the pseudo-density distribution field accordingly. Furthermore, this loop includes a calibration mechanism that periodically invokes high-fidelity simulations to verify and fine-tune the proxy model. Optionally, convergence criteria can be set in a complex manner. For example, convergence can be determined when the rate of change of the objective function (e.g., core loss) is less than 0.1% and the overall change in the pseudo-density distribution field is less than 1% over 10 consecutive iterations. Alternatively, to prevent infinite loops, an upper limit on the number of iterations, such as 2000, can be set.

[0027] The converged pseudo-density distribution field is determined as the final pseudo-density distribution field, and a final structural geometric model is generated based on the final pseudo-density distribution field, wherein the final structural geometric model can be recognized and processed by CAD / CAM software.

[0028] In this embodiment, this step is used to convert the optimized mathematical results into physical entities. Specifically, the final pseudo-density distribution field is thresholded. For example, all cells with a pseudo-density greater than 0.5 are defined as solid, and those less than 0.5 as air. Using an isosurface extraction algorithm such as MarchingCubes, a surface mesh with jagged edges is generated from this binarized field. This surface mesh is then processed using algorithms such as topological filtering and Laplace smoothing to obtain a smooth, free-standing structural geometry suitable for 3D printing or mold manufacturing.

[0029] Through the above steps, this embodiment can systematically and automatically generate a composite core inductor structure with optimized magnetic flux distribution. Compared with traditional design methods that rely on engineers' experience, both performance and R&D efficiency are improved.

[0030] like Figure 2 As shown, according to one aspect of the present application, an iterative optimization process is performed, comprising:

[0031] The proxy model is used to predict the electromagnetic performance of the current pseudo-density distribution field and obtain the predicted electromagnetic response data.

[0032] In this embodiment, the predicted electromagnetic response data includes the magnetic field vector and core loss scalar at each grid cell, forming a predicted field covering the entire design domain. Specifically, at the beginning of each iteration, the system feeds the current pseudo-density distribution field (in the first iteration, the initialized field) as an input tensor to the pretrained proxy model. Through a single forward computation of the neural network, predicted electromagnetic response data covering the entire design domain can be obtained in an extremely short time (e.g., milliseconds).

[0033] Based on the predicted electromagnetic response data, the pseudo-density distribution field is updated through the topology optimization algorithm for the next iteration.

[0034] Optionally, the topology optimization algorithm includes: extracting the frequency value f of the current analysis based on the predicted electromagnetic response data, and determining the penalty factor p according to the frequency value f: p(f)=p0+k·f α ; Among them, p0, k and α are preset parameters that characterize the high-frequency response characteristics of the material; in the topology optimization iteration, the interpolation relationship between the pseudo-density and the physical properties of the material is adjusted according to the penalty factor p.

[0035] In this embodiment, the topology optimization algorithm iteratively redistributes materials within a given design domain to maximize or minimize a certain performance indicator while satisfying specific constraints. Specifically, the design sensitivity of the objective function (for example, the predicted value of the total core loss) to the pseudo-density of each unit is calculated based on the predicted electromagnetic response data. The topology optimization algorithm includes: determining the value of the penalty factor p based on the frequency value f of the current analysis, wherein the penalty factor p is a monotonically increasing function of the frequency f; and applying the penalty factor p to adjust the interpolation relationship between the pseudo-density and the physical properties of the material during the topology optimization iteration. Preferably, the determination of the value of the penalty factor p is specifically implemented by a preset power law function relationship: p(f)=p0+k•f α Where f is the frequency of the current analysis, for example, 100 kHz; p0, k, and α are preset parameters that characterize the material's high-frequency response. For example, for a typical ferrite material, p0 can be set to 3.0, k to 1.0e-5, and α to 1.2. The optimizer (e.g., the moving asymptote method (MMA)) combines the design sensitivity information with the currently calculated penalty factor p to calculate the updated pseudo-density for each element and generate a new pseudo-density distribution for the next iteration. Traditional topology optimization cannot handle conflicting optimization objectives under different physical conditions (such as DC and high-frequency AC).

[0036] This embodiment addresses the core technical issue of conflicting optimization objectives by introducing a penalty factor p, whose value increases monotonically with frequency f, into the topology optimization algorithm. This provides a single optimization algorithm with dual-modal intelligent decision-making capabilities. Specifically, when the optimization algorithm addresses magnetic saturation issues associated with DC bias (where frequency f is assumed to be zero), a lower penalty factor p makes the algorithm more tolerant of intermediate density values ​​in material interpolation calculations, allowing and encouraging the algorithm to generate topological structures with fine, complex details to construct magnetic flux paths with minimal reluctance. Conversely, when the algorithm addresses eddy current losses associated with high-frequency ripple (where f is a high frequency value), a significantly higher penalty factor p strongly suppresses the presence of intermediate density values, forcing the algorithm to generate a simple structure with smoother boundaries and no small features. The resulting single geometric model has a macrostructure that meets DC bias requirements while simultaneously suppressing high-frequency losses at the microstructure level, achieving an optimal solution with comprehensive performance that is unattainable with conventional single-modal optimization algorithms.

[0037] The current pseudo-density distribution field is simulated periodically to obtain high-fidelity simulation results, and the proxy model is updated according to the high-fidelity simulation results for subsequent iterations.

[0038] In this embodiment, to ensure the accuracy of the proxy model throughout the iterative process, the system periodically performs calibration (e.g., every 50 iterations). Specifically, the proxy model is updated using an online learning strategy designed to mitigate catastrophic forgetting during the iterative optimization process. Catastrophic forgetting refers to the potential loss of the proxy model's understanding of the global design space acquired during early optimization when fine-tuned based on data from later optimization stages. This online learning strategy involves extracting the latest high-fidelity simulation results as high-timeliness calibration samples; retaining a subset of samples from the initial training data covering the global design space as a historical experience pool; and constructing a hybrid training set that balances local accuracy and global generalization by combining these high-timeliness calibration samples with samples from the historical experience pool. The proxy model is then updated based on this hybrid training set. For example, a hybrid training mini-batch consisting of one recent sample and nine randomly selected historical samples can be constructed, and the proxy model can be fine-tuned for several epochs using a small learning rate. In later iterations, the design solutions encountered by the surrogate model become increasingly similar and localized. If updated only with the latest data, it will quickly overfit to the current design region, leading to catastrophic forgetting. By introducing an experience replay mechanism, the model learns the latest knowledge (highly time-sensitive samples) while also forcing it to review older knowledge (the historical experience pool), enhancing the model's stability and global predictive capabilities throughout the optimization process.

[0039] This example describes an iterative optimization process that efficiently explores and converges to the optimal solution through a closed loop of proxy model prediction, topology optimization updates, and periodic calibration. In other words, after obtaining the initial finite element model and proxy model, the iterative optimization process continues until the pseudo-density distribution field converges. This allows for accurate and robust optimization of the core structure while ensuring computational efficiency, ultimately resulting in an excellent design that balances multiple physical properties.

[0040] In an optional embodiment, the optimizer in the topology optimization algorithm is implemented by using a moving asymptote method (MMA) as the optimizer to solve and update the pseudo-density distribution field.

[0041] The moving asymptote method does not directly solve the original, complex optimization problem, but instead constructs a simpler, convex, and separable approximate subproblem at each iteration point to solve it. It is particularly suitable for large-scale topology optimization problems. In each iteration of topology optimization, after obtaining the design sensitivity field (i.e., the gradient of the objective function with respect to the pseudo-density of each element), the MMA optimizer works as follows: The core of MMA is to calculate the pseudo-density of each element for each design variable (i.e., the pseudo-density ρ of each element) at the current iteration point k. i ) constructs a convex approximation function. The approximation function is based on the calculated design sensitivity field (first-order Taylor expansion) and introduces a pair of upward and downward moving asymptotes Li and U i To make corrections, a well-conditioned and easy-to-solve subproblem is formed. The intelligence of MMA lies in its moving asymptotes. After each iteration, the algorithm dynamically adjusts the asymptotes L used in the next iteration based on the oscillation of each design variable in the past few iterations. i and U i If a variable continues to oscillate, the algorithm will tighten its asymptotes and limit its range of variation, thereby suppressing oscillations and ensuring convergence; if a variable changes steadily, the algorithm will relax its asymptotes, allowing it to explore more boldly. The constructed convex subproblem can be quickly solved by an efficient dual solver to obtain the optimal update amount Δρ for each design variable in this iteration. i Apply the updated value to the current pseudo-density distribution field to complete this iteration.

[0042] The choice of MMA as the optimizer, rather than the simpler gradient descent method, is driven by the pursuit of stability and efficiency in the optimization process. Topology optimization is a highly nonlinear and complex problem with a large number of design variables. Simple gradient descent methods have difficult-to-control step sizes and are prone to falling into local optima or experiencing violent oscillations. MMA, through its unique, history-based moving asymptote strategy, enables more robust and efficient exploration, enabling it to consistently converge to a high-quality solution with fewer iterations.

[0043] like Figure 3 As shown, according to one aspect of the present application, generating a final structural geometric model includes:

[0044] The final pseudo-density distribution field is thresholded to distinguish the solid area from the air area and generate a binary structure.

[0045] In this embodiment, the final pseudo-density distribution field obtained after optimization convergence is thresholded, and a global threshold is set, for example, 0.5. All grid cells with pseudo-density values ​​greater than 0.5 are marked as solid material, and those less than or equal to 0.5 are marked as air, generating a binary structure.

[0046] Extract the initial geometric surface mesh from the binary structure (for example, using an isosurface extraction algorithm); perform topological filtering and Laplace smoothing on the initial geometric surface mesh to remove suspended islands and correct jagged boundaries to obtain the final structural geometric model.

[0047] In this embodiment, topological filtering refers to the screening and modification of structures within a geometric model based on topological rules (such as connectivity and volume) to remove isolated or excessively small components that are not practical for engineering purposes. Laplace smoothing eliminates high-frequency noise (jaggies) on the mesh surface by shifting each vertex to the average position of its neighboring vertices. Specifically, a marching cubes algorithm is used to traverse all cells of the binary structure, generating triangular facets at the interface between the solid and air (i.e., isosurfaces with a pseudo-density of 0.5) to extract the initial geometric surface mesh. This mesh accurately reflects the optimized macromorphology, but is typically rough and may contain defects. Before smoothing, the initial geometric surface mesh undergoes topological filtering: all independent mesh pieces not connected to the main structure (i.e., dangling islands) are identified and their volumes are calculated. All islands with a volume less than a preset threshold (e.g., 0.05% of the total design domain volume) are directly deleted. This effectively avoids the appearance of meaningless, unmanufacturable, and tiny fragments in the final model. Laplace smoothing is then performed on the filtered mesh. For example, 15 smoothing iterations can be performed, with the smoothing factor (which controls the step size of each movement) set to 0.2. This effectively eliminates the stair-stepping effect (jagged boundaries) introduced by the marching cubes algorithm, smoothing the model's surface and making it more compatible with fluid dynamics and mold manufacturing requirements. Ultimately, the final structural geometry is obtained, ready for use in downstream applications.

[0048] This embodiment solves common problems (such as aliasing and islanding) in the process of converting mathematical solutions into engineering entities, so that the output is not only a theoretically optimal solution, but also a high-quality industrial design with a robust structure, smooth surface, and high manufacturability in the real world.

[0049] According to one aspect of the present application, the proxy model is a graph neural network, the core of which is a physical propagation layer embedded with physical operators, where:

[0050] Mapping the mesh elements of the initial finite element model to nodes of the graph, where each node may contain initial features including its pseudo-density value and material identification;

[0051] The node feature information is processed and updated through the physical transmission layer. The specific implementation process of the physical transmission layer is as follows:

[0052] Aggregate the feature information of the neighboring nodes of each central node to form neighborhood information. Specifically, for any central node in the graph, the system collects the feature vectors of all its directly connected neighboring nodes;

[0053] Applying a learnable transformation to the neighborhood information, which simulates an electromagnetic curl operator, to calculate the rotational effect of the neighborhood magnetic field on the central node. The learnable transformation is preferably implemented as a small neural network (e.g., a two-layer multilayer perceptron) designed to approximate a discretized curl matrix operator. The transformation maps the input neighborhood information into a vector representing the rotational effect.

[0054] The feature information of the central node is updated by incorporating the rotation effect. For example, the original feature information of the central node can be fused with the calculated rotation effect vector through the update mechanism of a gated recurrent unit (GRU) to generate updated feature information for the node. Preferably, three to six physical propagation layers can be stacked to allow information to propagate further within the graph structure and capture physical interactions over a wider range.

[0055] This embodiment solves the technical bottleneck of the difficulty in balancing the simulation efficiency and physical accuracy of the proxy model by directly embedding the learnable transformation of the simulated electromagnetic curl operator into the propagation layer of the graph neural network. A proxy model that is both fast and physically aware is constructed. Conventional black-box neural networks require massive amounts of data to blindly learn the correlation between input and output, but the graph neural network of this embodiment, because its network structure itself contains the physical prior knowledge that magnetic field lines surround electric currents, makes its learning process no longer blind. When it learns data, its solution space is naturally constrained to a range that conforms to the laws of electromagnetic physics. This allows the proxy model to converge quickly with only a small amount of training data, and when faced with new topological structures generated during the optimization process that have never been seen in the training set, it can still make physically reasonable and accurate predictions, breaking the contradiction between efficiency and accuracy and achieving fast and reliable iterative optimization.

[0056] In an optional embodiment, the learnable transformation is specifically implemented as follows: the learnable transformation for simulating the electromagnetic curl operator is constructed as a three-layer multi-layer perceptron (MLP).

[0057] In this embodiment, the learnable transformation specifically refers to a function with trainable parameters, whose structure is designed to approximate a specific physical operator. Unlike fixed mathematical formulas, it can be fine-tuned through data training to more accurately match complex real-world physical effects. Specifically, construct the input layer: the input of the MLP is the aggregated neighborhood information tensor. For example, if a central node has N neighbors and the features of each node are D-dimensional, the dimension of the input tensor can be flattened to N×D. Design the hidden layer: The first hidden layer contains 128 neurons and uses SiLU as the activation function to extract preliminary nonlinear features from the original neighborhood information. The second hidden layer contains 64 neurons and also uses the SiLU activation function to further combine and abstract the features. Design the output layer: The output layer contains 3 neurons and uses a linear activation function. The output values ​​of these 3 neurons together constitute a three-dimensional vector, which represents the rotation effect of the neighborhood magnetic field on the central node. The whole process can be expressed as rotation effect_vector=MLP head (ReLU(W2ReLU(W1 neighborhood information + b1) + b2)), where W1, b1, W2, b2, etc. are the weights and bias parameters that the network needs to learn during training; MLP head A learnable transformer that simulates the physical curl operator.

[0058] This example designs the learnable transform as a refined MLP. The goal is not to force the network to learn arbitrarily complex functions, but rather to approximate local, well-behaved physical operators (curl). A shallow MLP with a moderate number of neurons has sufficient nonlinear fitting capabilities for this task, while avoiding the risk of overfitting and waste of computational resources associated with an overly large network. This is the key to achieving efficient and accurate physical information infusion.

[0059] In some other embodiments, the learnable transformation can also be implemented using a one-dimensional convolutional neural network (1D-CNN). In this case, the feature information of neighboring nodes can be viewed as a sequence, and the convolution kernel is used to extract their local spatial correlation, which can also effectively simulate the effect of discrete operators.

[0060] In order to make the learning direction of the neural network respect both the simulation data and the laws of physics, its training process is guided by a composite loss function. According to one aspect of the present application, the determination of the composite loss function includes:

[0061] The mean squared error between the predicted value of the graph neural network and the pre-stored high-fidelity simulation data is calculated to obtain the data fidelity loss term; here, the predicted value includes the magnetic field and loss, and the mean squared error of the two is weighted and summed during the calculation.

[0062] The physical residual loss term is obtained by calculating the mean squared value of the discrete divergence of the magnetic field predicted by the graph neural network at all nodes. Here, the discrete divergence calculation is represented on the graph as follows: for each node, the sum of the predicted magnetic flux flowing to all adjacent nodes is calculated. According to Gauss's law in electromagnetics, this sum should be zero.

[0063] The preset weight coefficient is used to adjust the weight of the physical residual loss term, and the physical residual loss term after weight adjustment is summed with the data fidelity loss term to obtain the final value of the composite loss function.

[0064] In a preferred embodiment, the composite loss function L total =L data +0.01L phys ; where L data is the data fidelity loss term; L phys The physics residual loss term is weighted by λ, which is set to 0.01. Setting this value is an engineering trade-off: a too high value may cause the network to adhere too rigidly to physical laws and ignore real data, while a too low value may make the physical constraints less effective. By constructing a dual optimization objective that combines data-driven and physical rule-driven approaches, the data fidelity loss term ensures that the model's predictions are consistent with reality (high-fidelity simulation), while the physical residual loss term acts as a regularizer. Even in the absence of real data labels, it still constrains the model's predictions, ensuring that its outputs are physically plausible and self-consistent.

[0065] This embodiment enhances the generalization and robustness of the proxy model by employing a composite loss function that incorporates a data fidelity loss term and a physical residual loss term when training the proxy model. This introduces a second layer of physical supervision to the neural network training that does not rely on labeled data. The data fidelity loss term aligns the model's predictions with real, high-fidelity simulation data, while the physical residual loss term acts as a physical rule checker. By calculating the discrete divergence of the predicted magnetic field, it independently verifies whether the predictions themselves violate the law of flux continuity (i.e., Gauss's law). This prevents the model from making absurd predictions that physically create or disappear magnetic field lines, even in regions of the design space where training data is sparse. Therefore, this composite loss function ensures that the proxy model remains physically consistent across the vast design space, a core guarantee for its use as a reliable evaluation tool throughout the optimization process.

[0066] According to one aspect of the present application, after processing through multiple physical propagation layers, the final node feature information obtained is a high-dimensional vector containing rich physical information. The graph neural network includes a parallel multi-head prediction module that decodes and outputs different physical quantities, wherein the multi-head prediction module includes:

[0067] The magnetic field prediction head is an independent, smaller neural network (e.g., a 2-layer MLP) that decodes the final node feature information output by the propagation layer into vector field data representing the magnetic field distribution.

[0068] In parallel with the magnetic field prediction head, the loss prediction head is another independent, smaller neural network used to decode the same final node feature information into scalar field data representing the core loss.

[0069] Predicting vector fields (magnetic fields) and scalar fields (losses) are two distinct tasks. If a single output network were used to accomplish both tasks simultaneously, they might interfere with each other, leading to a loss of focus. This embodiment utilizes a parallel multi-head architecture, allowing each head to focus on its own decoding task without interfering with each other. This improves the overall accuracy and training efficiency of predictions for different physical quantities. Using graph neural networks, we achieve a highly efficient and accurate proxy model specifically for electromagnetic topology optimization that far surpasses conventional machine learning methods.

[0070] According to one aspect of the present application, a method for generating data for pre-training of a proxy model adopts a two-stage sampling strategy to ensure that the initial training set has both static spatial diversity and dynamic evolutionary diversity. The method includes:

[0071] In the static spatial diversity sampling stage, Latin hypercube sampling is used to generate a low-bias initial pseudo-density distribution that covers the global design space.

[0072] In this embodiment, static spatial diversity refers to the uniformity and coverage of the initial design solutions (i.e., initial pseudo-density distribution) used for training across the entire possible design space. For datasets with high static spatial diversity, sample points are not clustered in similar areas, but rather are widely dispersed. Specifically, the pseudo-density distribution of the entire design domain is considered an ultra-high-dimensional space (the density of each grid cell is one dimension). The Latin Hypercube Sampling (LHS) algorithm stratifies each dimension with equal probability and ensures that only one sample point is extracted from each layer-dimensional combination. Compared to pure random sampling, LHS can achieve more uniform coverage of the entire design space with fewer sample points, effectively avoiding clustering and bias in the initial samples. For example, LHS can be used to generate 100 statistically unrelated, different initial pseudo-density distributions.

[0073] In the dynamic evolution diversity sampling stage, each initial pseudo-density distribution is taken as a starting point and an independent short-range topology optimization is performed to capture multiple parallel design evolution trajectories starting from different starting points.

[0074] In this embodiment, dynamic evolutionary diversity refers to the various structural evolution paths experienced during the optimization process, starting from different initial designs, within the dataset. A dataset with high dynamic evolutionary diversity not only contains a variety of snapshots but also a variety of change processes. Specifically, starting from the 100 initial pseudo-density distributions generated in the previous step, 100 independent optimization tasks are initiated in parallel. Each task performs a complete topology optimization for a fixed number of steps (e.g., 20 iterations). Due to their different starting points, these 100 optimization tasks will explore completely different regions of the design space, forming 100 unique design evolution trajectories.

[0075] All state points and corresponding simulation results on multiple design evolution trajectories are collected to construct the final initial training set.

[0076] In this example, the input (the pseudo-density distribution field at that time) and output (the corresponding high-fidelity simulation result) of each iteration (a total of 10020 = 2000 times) of the aforementioned 100 trajectories are collected as a single data point. All of these data points are combined to form a final, large-scale, and highly diverse initial training set for pre-training the surrogate model. A good surrogate model must not only recognize various different structures but also understand how one structure evolves into another. Traditional single-shot full optimization sampling only provides a single evolutionary path. The two-stage strategy, on the other hand, ensures static diversity through LHS and dynamic diversity through multi-path parallel short-range optimization. This allows the surrogate model to learn more essential and comprehensive knowledge about the design space, resulting in stronger predictive capabilities.

[0077] This embodiment solves the problem of poor model performance caused by sampling bias when constructing the initial training set for the surrogate model by adopting a two-stage strategy of Latin hypercube sampling and multi-path parallel short-range optimization. This generates a high-quality dataset that combines both static spatial diversity and dynamic evolutionary diversity. Specifically, Latin hypercube sampling in the first stage enables the starting points of the optimization to evenly and unbiasedly cover the entire design space, avoiding the narrow model vision caused by the clustering of initial samples. Multi-path short-range optimization in the second stage captures different evolutionary trajectories starting from different starting points, ensuring that the training data not only includes diverse structural snapshots but also diverse structural change processes. The surrogate model, having learned from rich and unbiased data, has a deeper and more comprehensive understanding of the macroscopic landscape and local details of the entire design space from the outset. It can enter the optimization loop with higher initial accuracy and stronger generalization ability, reducing the number of subsequent online calibrations and accelerating the convergence of the overall optimization process.

[0078] According to one aspect of the present application, a simulation is performed to obtain high-fidelity simulation results using an AC / DC decoupling analysis method. The AC / DC decoupling analysis is a simulation strategy for conditions where large DC bias and small-signal AC excitation exist simultaneously. It solves two physical processes of different natures separately to avoid the computational difficulties and precision loss caused by nonlinear coupling, and is particularly suitable for the analysis of components such as inductors in switching power supplies. When a high-fidelity simulation of a given pseudo-density distribution field is required, instead of performing a single transient or frequency domain simulation, the following decoupling process is used:

[0079] Perform a DC magnetic field analysis to determine the core saturation state and static magnetic flux distribution under DC bias current. For example, for an inductor designed to withstand a 20A DC current, perform a nonlinear static magnetic field simulation to calculate the permeability and static magnetic flux distribution at various locations within the core under this current.

[0080] Based on the material operating point determined by this static magnetic flux distribution, an AC small-signal analysis is performed to calculate the high-frequency core loss under ripple current excitation. This material operating point refers to the equivalent permeability and loss characteristics of the material exhibited for small-signal excitation due to saturation effects after the previous DC analysis. Based on this operating point, an AC excitation (for example, a ripple current of 2A at 100kHz) is applied and a frequency-domain simulation is performed to accurately calculate the resulting high-frequency core loss.

[0081] Combine the static flux distribution with the high-frequency core losses to create a complete high-fidelity simulation. This high-fidelity simulation is a multi-physics data package that serves as the most accurate ground truth for calibrating surrogate models or building training sets.

[0082] In many power electronics applications, inductors operate simultaneously under strong DC and weak AC. These two components have different mechanisms of influence on the magnetic core and influence each other, with strong nonlinearity. If a single transient simulation is used, the time step must be extremely small in order to accurately capture high-frequency ripples, and the total time required to simulate the DC establishment process is very long, resulting in extremely high computational costs. The AC / DC decoupling method adopted in this embodiment decomposes the problem and solves each part using the most appropriate method (using a static magnetic field solver for DC and a frequency domain solver for AC). This not only improves the computational efficiency of a single high-fidelity simulation, but also reflects the loss of the magnetic core under actual working conditions more accurately than a full transient simulation, providing higher-quality true value data for the training and calibration of the proxy model.

[0083] like Figure 4 As shown, according to one aspect of the present application, calculating the rotation effect of the neighborhood magnetic field on the central node includes:

[0084] The saturated and unsaturated regions are divided based on the local saturation of the central node, and a dual-mode operator structure is constructed: the curl weight matrix is ​​used to calculate the circulation characteristics in the unsaturated region, and the diffusion weight matrix is ​​used to calculate the potential field diffusion in the saturated region;

[0085] The local saturation is mapped to the modal fusion coefficient through a smooth switching function, and the circulation characteristics and potential field diffusion are weighted and fused according to the modal fusion coefficient to generate an adaptive weight matrix;

[0086] The adaptive weight matrix is ​​applied to the neighborhood magnetic field difference to calculate the rotation effect adapted to the saturation state.

[0087] According to one aspect of the present application, local saturation is achieved through a timing decoupling mechanism, including:

[0088] Maintaining a time series record queue containing historical local saturation values;

[0089] When the time series record queue length is greater than or equal to the preset threshold, the linear extrapolation formula is used to calculate the extrapolated saturation at the current moment; otherwise, the initial estimated value of the extrapolated saturation is set;

[0090] A preliminary dual-modal fusion calculation is performed based on the extrapolated saturation to obtain the predicted rotation effect;

[0091] The magnetic flux density field is reconstructed from the predicted rotation effect by an integration operation, and the actual local saturation is obtained by dividing the magnetic flux density amplitude of the reconstructed magnetic flux density field by the pre-stored material saturation magnetic flux density.

[0092] According to one aspect of the present application, when the deviation between the actual local saturation and the extrapolated saturation exceeds a preset deviation threshold, a local correction iteration is performed:

[0093] Re-calculating the bimodal fusion using the actual local saturation instead of the extrapolated saturation to obtain the re-predicted rotation effect (including recalculating the modal fusion coefficient using the actual saturation; reconstructing the adaptive weight matrix based on the updated fusion coefficient and recalculating the rotation effect);

[0094] Based on the re-predicted rotation effect, the corrected magnetic flux density field is obtained by integration again and the new local saturation is calculated;

[0095] The local correction process is repeated until the deviation is less than the deviation threshold or the maximum number of iterations is reached, and the final local saturation is output.

[0096] In one embodiment of the present application, the basic components of the bimodal operator are initialized. Before the optimization task begins, the system first performs a one-time initialization. Construct the standard curl modal weight matrix: read the grid topology data, and construct a discrete curl weight matrix W that represents the ideal curl operation based on the standard finite volume method.curl . Construct the saturated diffusion mode weight matrix: also read the grid topology data, construct the discrete Laplace operator, and obtain the diffusion weight matrix W diff This matrix is ​​used to simulate the diffusion behavior of the magnetic field in deeply saturated materials. Initialize the timing record structure: create a queue with a fixed capacity (for example, a length of 5) for each node (or each region) in the design domain to store its historical saturation values, and obtain an empty saturation history record S history At the same time, the switching control parameter set is set as 0.85 for the saturation threshold and 20 for the switching steepness.

[0097] Build a saturation timing extrapolation mechanism and fuse the bimodal weights. Each time the physical propagation layer is called: read the saturation history record S history If its length is greater than or equal to 2, the linear extrapolation formula S is used. ext = 2S[-1] - S[-2] to predict the saturation at the current moment, where S[-1] and S[-2] represent the saturation distribution field at the previous moment (most recent) and the one before that; if the length is insufficient, an initial guess value (such as 0.5) is assigned. The result is clipped to the interval [0, 1] to obtain the extrapolated saturation field S extrapolated Based on the extrapolated saturation field S extrapolated And switch control parameter set, through Sigmoid function α = 1 / (1 + exp(-20 (S extrapolated - 0.85))) Calculate the modal switching coefficient of each node and generate the modal switching coefficient field α field According to the mode switching coefficient field α field , the curl weight matrix W curl and the diffusion weight matrix W diff Perform weighted fusion: W adaptive = (1-α)W curl +αW diff Furthermore, the permeability modulation is introduced into the diffusion term, namely W adaptive =(1-α) W curl + α(1 / μ eff ) W diff_base , where μ eff It is the equivalent permeability determined by the BH curve according to the saturation; W diff_base is the basic diffusion weight matrix. Finally, the adaptive weight matrix W is obtained adaptive By extrapolating historical states, we bypass the endless loop of requiring the current saturation to calculate the current magnetic field, and vice versa. This allows us to transform a process that originally required complex nonlinear iterative solutions into an efficient, forward-propagation prediction process.

[0098] Perform prediction-correction iterative calculation, which is the micro-cycle inside this operator. Use the adaptive weight matrix W generated in the previous step adaptive , perform a matrix operation on the current input magnetic field intensity distribution to obtain the preliminary curl field Curl preliminary By solving Poisson's equation, we can get the initial curl field Curl preliminary Integrate and reconstruct the predicted magnetic flux density field B predicted . Compare it with the material's saturation flux density B sat Compared with the actual saturation S, the actual = |B predicted | / B sat , where B sat is the magnetic flux density; and its extrapolated saturation S is obtained extrapolated The prediction deviation between the fields ΔS= |S actual -S extrapolated Check the predicted deviation field ΔS. For nodes whose deviation exceeds a preset threshold (e.g. 0.1), the system will trigger a local correction microcycle (up to 3 iterations): At these nodes, the actual saturation S just calculated is used actual Replace the extrapolated value and re-execute the dual-modal fusion calculation to obtain a more accurate local solution. The final curl field is the corrected curl field Curl corrected The prediction-correction mechanism is the embodiment of the robustness of this embodiment. It acknowledges that there may be errors in the extrapolated prediction and establishes a low-cost error correction mechanism. In most cases, the extrapolation is accurate and the process can be passed quickly; only in a few areas where the saturation state changes drastically, a more accurate correction iteration is started. It is the optimal allocation of computing resources, taking into account both efficiency and accuracy. Update the timing information and output the results. The calculated actual saturation S actual The spatial average value is appended to the saturation history S history The end of the queue and remove the oldest record to maintain the queue length. corrected The relevant magnetic flux density, calculation status and other information are packaged into a structured adaptive curl calculation result package, which is used as the final output of the entire physical propagation layer for further processing by the next layer of the graph neural network.

[0099] Through the time-decoupled dual-modal curl operator of this embodiment, the proxy model has an unprecedented ability to handle strong nonlinear saturation effects, enabling it to maintain extremely high prediction accuracy and stability when facing the core design optimization problem under extreme working conditions.

[0100] In another embodiment of the present application, the process of constructing a physical propagation layer embedded with discretized physical operators, especially the core learnable transformation part, is to construct a parameterized discrete curl operator. Specifically, based on the adjacency relationship of the graph and the geometric coordinates of the nodes, the initial curl weight matrix W is constructed. curl_init The weights of this initial matrix are based on physical formulas (e.g., w j = (1 / |r ij | 2 )×r ij ×n ij ) calculation, which reflects the basic curl operation structure. The initial weight matrix is ​​used as input, processed by a small MLP (two-layer network), and its output range is constrained by a Sigmoid function to obtain the final learnable curl weight matrix W curl_learned This embodiment does not learn completely from scratch; it retains the basic structure of the curl operator derived from physical formulas (via initialization) while giving the network the ability to fine-tune and correct the basic operator through training data (via MLP). This ensures both physical rationality and the flexibility to adapt to errors in complex real-world data.

[0101] Construct a regularized loss that includes physical constraints. Specifically, in addition to the composite loss function, a physical regularized loss term L is further added to the model training. physics_reg The regularization term consists of three parts: antisymmetric regularization L antisym : Penalize the norm of the sum of the learned weight matrix and its transposed matrix to encourage the operator to maintain antisymmetry; norm constraint regularization L norm : Penalize the part of the weight matrix whose norm exceeds a preset upper bound (such as τ=3.0) to prevent numerical instability; local conservation regularization L conserve : Penalize the divergence of the curl field at each node to force it to satisfy the physical law that the divergence of the curl is zero. The total physical regularization loss is the weighted sum of these three terms, such as L physics_reg = 0.1 × L antisym + 0.05×L norm + 0.2×L conserve This embodiment constantly reminds the network during training that any transformations it learns must not violate the most basic laws of physics. This reduces the likelihood of the model making unphysical predictions, improving its generalization ability and the interpretability of its results.

[0102] In another embodiment of the present application, the iterative optimization process is optimized to ensure stable and efficient convergence. Specifically, an oscillation detection and suppression mechanism is implemented. During the optimization process, the system continuously monitors the objective function value sequence for the last 20 iterations. An oscillation indicator (OSC) is quantified by calculating the ratio of the standard deviation to the mean of this sequence. When this indicator exceeds a preset threshold (e.g., 0.1), the system determines that an oscillation state has occurred and automatically activates suppression measures. For example, the optimization algorithm step size is reduced by 50%, the current penalty factor p is temporarily increased by 10%, and a Gaussian smoothing filter is applied to the current pseudo-density distribution field. Preferably, a pool of historical optimal solutions containing the optimal solutions to the objective function from the last 50 iterations is maintained. When oscillation is detected, the most stable solution from this pool (e.g., the solution with the lowest variance over the next 10 iterations) is directly selected as the restart point. This is generally more effective than simple parameter reduction. An extension strategy and adaptive precision control are employed to accelerate convergence. In the early stages of optimization (e.g., the first 100 iterations), a low, fixed penalty factor (e.g., p = 3.0) is used to quickly form a rough outline of the structure. In the mid-term (e.g., 100-300 iterations), the penalty factor is gradually increased to a target value determined by frequency. In the early stages of optimization, the tolerance for the surrogate model's prediction error can be relaxed (e.g., 15%), allowing the algorithm to explore aggressively. As optimization progresses, the accuracy requirements are gradually tightened, and the frequency of high-fidelity calibration is increased, so that in the later stages of optimization, each iteration is based on high-precision predictions.

[0103] This embodiment forms an intelligent convergence strategy that combines coarse and fine tuning. The continuation strategy avoids optimization stalls caused by excessively high penalties in the early stages, while adaptive precision control focuses valuable computing resources (high-fidelity calibration) on the most valuable late stages of optimization. This reduces the typical number of convergence iterations.

[0104] In another embodiment of the present application, the process of determining the key parameters p0, k, and α in the frequency adaptive penalty factor p(f) is as follows: for the magnetic core material to be optimized (e.g., N87 ferrite), a material property test sample with a standard geometric structure (e.g., a toroidal core with an outer diameter of 30 mm) is prepared. A frequency sweep test is performed on this set of samples using an impedance analyzer. For example, 20 frequency points are selected at logarithmic intervals within the range of 1 kHz to 1 MHz, and the eddy current loss density P corresponding to each frequency point is measured and recorded under a constant magnetic flux density amplitude. eddy (f) The goal is to obtain reliable benchmark data reflecting the frequency-dependent variation of eddy current losses through rigorous experimental control. This data serves as the physical basis for all subsequent parameter fitting and model building.

[0105] Establish the relationship between frequency, critical microstructure size and optimal penalty factor. The critical microstructure size refers to the size threshold at which eddy current loss increases exponentially if features smaller than this size appear in the topological structure at a specific frequency. Specifically, establish the relationship between frequency f and critical microstructure size d through simulation or experiment. crit The relationship between the skin depth is obtained by analyzing the skin depth formula Δ=sqrt(2 / (ωμσ)) or by testing a group of samples with different microstructure sizes, where ω is the angular frequency, μ is the magnetic permeability, and σ is the electrical conductivity. For example, the test found that at 100kHz, d crit About 0.8mm. Based on this, the optimal penalty factor p is established opt The empirical relationship with frequency f, such as p opt (f) = 3 + log 10 (f / 1000)×2.5. At higher frequencies, the optimization algorithm will naturally avoid generating overly fine structures due to higher penalty values.

[0106] The obtained multiple (f, p opt ) data points, and the least squares method is used to fit the power law function form p(f)=p0+k•f α Specifically, we can first perform linear regression on the logarithm of the formula, and then optimize the p0 value through grid search and other methods to minimize the fitting error. For example, for N87 material, the final fitting results in a set of parameters: p0 = 3.0, k = 1.2×10 -5 , α = 1.15. This set of parameters was applied to a confirmatory topology optimization case to verify whether the optimization results (such as minimum feature size and loss reduction rate) met expectations, thus confirming the effectiveness of the parameters.

[0107] In a specific embodiment of this application, a technical solution for optimizing a specific composite magnetic core inductor is described. For example, this is used in a main inductor designed for a high-power 48V to 12V server power module, a step-down (Buck) converter. The physical structure utilizes an EI-type composite magnetic core. The E-type main core utilizes high-permeability ferrite material (e.g., TDK's N87) to provide the main magnetic circuit. The I-type baseplate utilizes high-saturation magnetic density Sendust powder core material (e.g., Magnetics' Kool Mµ) to provide an air gap and withstand strong DC bias. The right-angle region where the E-type core's center leg meets the I-type baseplate is a typical point of flux path abrupt change. During operation, magnetic flux lines here become highly congested, easily reaching local magnetic saturation, resulting in a sharp drop in inductance and a significant increase in core losses, thus impacting the efficiency and stability of the entire power supply. The goal is to automatically and optimally design the geometry of this right-angle region to address the flux congestion issue. This region is defined as the design domain for topology optimization. The inductor needs to withstand a 30A DC bias current, superimposed with an AC ripple current of 5A peak-to-peak and 200kHz frequency. The processing process is as follows:

[0108] An initial finite element model of the EI composite core was created in CAE software. The design domain was divided into approximately 50,000 tetrahedral elements, and the pseudo-density distribution field was initialized (for example, the initial value of all elements was set to 0.6). Based on the initial finite element model, an initial training set of approximately 2,500 data points was generated using a two-stage sampling strategy. A proxy model was pre-trained based on this initial training set. The core iterative optimization process was then executed. The following is a detailed description of the calculation process for a representative iteration step in this process: Assuming the current iteration is the kth iteration, the system will perform the following key calculations:

[0109] Step A: Proxy model prediction and calculation process of physical propagation layer.

[0110] The current proxy model is used to quickly predict the pseudo-density distribution field obtained at the end of the k-1th iteration. The core of this process lies in information processing at the physical propagation layer. Taking any grid cell in the design domain (i.e., the central node i in the graph) as an example, the computational process is as follows: Neighborhood Information Aggregation: The system first identifies all neighboring nodes directly adjacent to the central node i (e.g., j1, j2, j3, j4), extracts the current feature vectors of these neighboring nodes (each containing information such as pseudo-density and material type), and aggregates them into a unified neighborhood information tensor. A learnable curl operator transformation is then applied. The aggregated neighborhood information is input into a specially designed learnable transformation module. This module is a small neural network in structure, but its design goal is not to fit arbitrary functions but to simulate the curl operator in electromagnetics. In physics, the curl of the magnetic field intensity H is equal to the current density J, and the learnable transformation is a discretized simulation of this law. Through a series of matrix multiplications and nonlinear activations, a quantized vector of the vortex or rotational effect generated by the magnetic potential distribution of the neighboring nodes on the central node i is calculated. The computational process essentially involves learning an optimal, discretized curl matrix that most accurately reflects the physical behavior of the local magnetic field. The calculated rotation effect vector is then fused with the eigenvector of central node i (for example, via the update gate mechanism of a GRU unit). This yields updated eigenvalues ​​of central node i, enriched with physical information, after this layer of propagation. Once all nodes have completed this computation, the information has been physically propagated one layer throughout the entire graph (i.e., the entire design domain).

[0111] Step B: Topology optimization update and calculation of frequency adaptation penalty factor.

[0112] After obtaining the predicted electromagnetic response data from the surrogate model, the system updates the structure using a topology optimization algorithm. The core of this process lies in the calculation and application of a frequency-adaptive penalty factor. The goal is to reduce total losses, of which high-frequency AC losses are the primary component, corresponding to a frequency f of 200kHz. The penalty factor p is calculated as: p(200000)=3.0+1.0e-5(200000) 1.2 ≈7.97; where p0=3.0, k=1.0e-5, and α=1.2 are preset parameters for the material used in this example. The calculated higher p value (approximately 7.97) is substituted into the material interpolation equation. In subsequent sensitivity analysis and density update calculations, a high penalty value increases the cost of intermediate pseudo-density values ​​(non-0 or 1), forcing the optimization algorithm to favor solutions with clearer boundaries, smoother structures, and fewer small features when updating the pseudo-density distribution field. Such solutions have lower eddy current losses at high frequencies.

[0113] Step C: Periodic calibration and AC / DC decoupling simulation process.

[0114] Assuming the current iteration number k is a multiple of 50, the system will perform a high-fidelity calibration. The core of this process lies in the AC / DC decoupling analysis method. A nonlinear static magnetic field simulation is performed in the CAE software using a 30A DC bias current for the structure represented by the current pseudo-density distribution field. After the simulation, the static magnetic flux distribution and the equivalent permeability change due to saturation are obtained for each mesh element. Based on the material operating point (i.e., equivalent permeability) of each element under DC bias determined in the previous step, a frequency-domain simulation is performed with a 2A@200kHz AC small-signal excitation to accurately calculate the high-frequency core loss. The static magnetic flux distribution and high-frequency core loss data obtained in the two steps are combined to form a complete high-fidelity simulation result. These results are then used to fine-tune the surrogate model using an online learning strategy designed to suppress catastrophic forgetting. The optimization process converged after approximately 600 iterations. The resulting pseudo-density distribution field is post-processed to generate the final structural geometry. Compared to the original right-angle design, the optimized structure features a non-uniform, smoothly transitioned chamfer at the edge of the E-shaped core's center column. Simulations show that this structure reduces peak flux density and high-frequency core losses at the corners, improving the inductor's overall performance.

[0115] The preferred embodiments of the present invention are described in detail above. However, the present invention is not limited to the specific details in the above embodiments. Within the technical concept of the present invention, various equivalent transformations can be made to the technical solutions of the present invention, and these equivalent transformations all fall within the scope of protection of the present invention.

Claims

1. A method for optimizing the magnetic flux distribution of a composite core inductor, for generating a final structural geometric model, characterized in that: include: Construct an initial finite element model and initialize the pseudo-density distribution field within the model design domain; Based on the initial finite element model, a proxy model for predicting electromagnetic performance is generated and pre-trained; based on this, an iterative optimization process is performed until the pseudo-density distribution field converges; Determine the converged pseudo-density distribution field as the final pseudo-density distribution field; And generate the final structural geometric model accordingly; Perform an iterative optimization process, including: The surrogate model is used to predict the electromagnetic performance of the current pseudo-density distribution field and obtain the predicted electromagnetic response data; Based on the predicted electromagnetic response data, the pseudo-density distribution field is updated through the topology optimization algorithm for the next iteration; Periodically simulate the current pseudo-density distribution field to obtain high-fidelity simulation results, and update the proxy model accordingly for subsequent iterations; Topology optimization algorithms include: Based on the predicted electromagnetic response data, the frequency value f of the current analysis is extracted, and the penalty factor p is determined accordingly: p(f)=p0+k·f α ; Among them, p0, k and α are preset parameters that characterize the high-frequency response characteristics of the material; In the topology optimization iteration, the interpolation relationship between the pseudo-density and the material physical properties is adjusted according to the penalty factor p; Generate the final structural geometry model, including: The final pseudo-density distribution field is thresholded to distinguish the solid area from the air area and generate a binary structure; The initial geometric surface mesh is extracted from the binary structure, and topological filtering and Laplace smoothing are performed on it to obtain the final structural geometric model.

2. The method according to claim 1, characterized in that The proxy model is a graph neural network, where: Mapping the mesh elements of the initial finite element model to nodes of the graph; The node feature information is processed and updated through the physical propagation layer. The implementation of this physical propagation layer includes: aggregating the feature information of the neighboring nodes of each central node to form neighborhood information; applying a learnable transformation for simulating the electromagnetic curl operator to the neighborhood information to calculate the rotation effect of the neighborhood magnetic field on the central node; and updating the feature information of the central node based on the rotation effect.

3. The method according to claim 2, characterized in that The training process of the graph neural network uses a composite loss function, specifically: Calculate the mean square error between the predicted value of the graph neural network and the pre-stored high-fidelity simulation data to obtain the data fidelity loss term; Calculate the mean square value of the discrete divergence of the magnetic field predicted by the graph neural network at all nodes to obtain the physical residual loss term; The weight of the physical residual loss term is adjusted using the preset weight coefficient, and it is summed with the data fidelity loss term to obtain the final value of the composite loss function.

4. The method according to claim 2, characterized in that The graph neural network includes a parallel multi-head prediction module, which contains: The magnetic field prediction head is used to decode the final node feature information output by the propagation layer into vector field data representing the magnetic field distribution; The loss prediction head, in parallel with the magnetic field prediction head, is used to decode the same final node feature information into scalar field data representing the core loss.

5. The method according to claim 2, characterized in that Calculate the rotation effect of the neighboring magnetic field on the central node, including: The saturated and unsaturated regions are divided based on the local saturation of the central node, and a dual-mode operator structure is constructed: the curl weight matrix is ​​used to calculate the circulation characteristics in the unsaturated region, and the diffusion weight matrix is ​​used to calculate the potential field diffusion in the saturated region; The local saturation is mapped to the modal fusion coefficient through a smooth switching function, and the circulation characteristics and the potential field diffusion are weighted and fused to generate an adaptive weight matrix. The adaptive weight matrix is ​​applied to the neighborhood magnetic field difference to calculate the rotation effect adapted to the saturation state.

6. The method according to claim 5, characterized in that Local saturation is achieved through a timing decoupling mechanism, including: Maintaining a time series record queue containing historical local saturation values; When the time series record queue length is greater than or equal to the preset threshold, the linear extrapolation formula is used to calculate the extrapolated saturation at the current moment; otherwise, the initial estimated value of the extrapolated saturation is set; A preliminary dual-modal fusion calculation is performed based on the extrapolated saturation to obtain the predicted rotation effect; The magnetic flux density field is reconstructed from the predicted rotation effect by an integration operation, and the actual local saturation is obtained by dividing the magnetic flux density amplitude of the reconstructed magnetic flux density field by the pre-stored material saturation magnetic flux density.

7. The method according to claim 6, characterized in that When the deviation between the actual local saturation and the extrapolated saturation exceeds a preset deviation threshold, a local correction iteration is performed: The actual local saturation is used instead of the extrapolated saturation to recalculate the dual-modal fusion and obtain the re-predicted rotation effect; Based on the re-predicted rotation effect, the corrected magnetic flux density field is obtained by integration again and the new local saturation is calculated; The local correction process is repeated until the deviation is less than the deviation threshold or the maximum number of iterations is reached, and the final local saturation is output.

Images