Multi-target electric field intensity intelligent optimization method and system based on individualized closed-loop nerve regulation

Abstract

The invention discloses a multi-target electric field intensity intelligent optimization method and system based on individualized closed-loop nerve regulation, and relates to the field of nerve regulation and brain electrical stimulation. The method comprises the following steps: constructing an electric field guide field matrix based on a finite element model, and sequentially injecting unit current under the condition of a fixed backflow electrode to establish a linear relation between the current and the electric field; according to the electric field intensity definition of the target area, constructing an optimized target function, constraining the target area to reach the set electric field intensity, and minimizing the electric field energy of a non-target area; processing non-convex constraint by adopting a convex-concave planning method, and introducing a slack variable and a penalty term to improve convergence; and setting a sparsity constraint to control the maximum number of activated electrodes, and finally outputting an optimized electrode flow injection scheme meeting intensity control and safety limitation. The method supports multi-target and multi-area combined control, has good stimulation focusing performance and practical feasibility, and is suitable for individualized transcranial electrical stimulation treatment and closed-loop nerve regulation and control systems.

Description

Technical Field

[0001] The present invention relates to a multi-objective electric field strength intelligent optimization method and system based on individualized closed-loop neural regulation. Background Art

[0002] Transcranial electrical stimulation (TES) is a noninvasive neuromodulation technique that injects weak currents into the brain via scalp electrodes, generating electric fields in specific brain regions to regulate neural function. Existing TES optimization methods are mostly based on the strategy of applying electric fields in specific directions to target brain regions. However, their application effectiveness is limited in subcortical regions or scenarios where the direction is unclear. Furthermore, current methods generally struggle to simultaneously control the electric field intensity in multiple target regions and lack unified modeling of practical physical constraints such as the number of electrodes and injected current. Summary of the Invention

[0003] To address the aforementioned issues in the existing technology, the present invention provides a multi-objective electric field intensity intelligent optimization method and system based on personalized closed-loop neuromodulation. These methods and systems can precisely control the electric field intensity in one or more target areas while meeting safe current constraints, while minimizing the electric field energy in non-target areas. This results in a more focused, safer, and more balanced neuromodulation effect. The present method and system are suitable for scenarios such as personalized brain stimulation, neurorehabilitation, and adjunctive treatment for psychiatric disorders.

[0004] A transcranial electrical stimulation control method for multi-objective electric field strength optimization according to an embodiment of the present invention comprises the following steps:

[0005] Step 1: Select a fixed electrode as the return path and inject a unit current into each of the remaining electrodes in sequence. Calculate the corresponding electric field response using the finite element method. Leveraging the linear relationship between the electric field and the injected current, a guiding field matrix is ​​constructed to quickly evaluate the electric field distribution generated by any combination of electrode currents.

[0006] Step 2: Based on the average value of the electric field intensity in the target area, construct a calculation expression for the electric field intensity in the target area. Discretize this expression to obtain a matrix-based linear relationship between the electrode current and the electric field intensity. Combined with the physical constraints, formulate the electrode current vector optimization problem.

[0007] Step 3: Using a convex-concave programming approach, the strength constraint is converted from an equation to an inequality and linearized. A sparsity constraint is introduced to control the number of activated electrodes, and slack variables and penalty terms are used to ensure good convergence and global search capabilities during the optimization process.

[0008] Step 4: Optimize iteratively until the objective function value and the electric field intensity in the target area no longer change significantly.

[0009] In step 1, the electric field simulation process includes: guiding field matrix construction, electrode arrangement, electric field simulation and boundary condition setting. First, under a fixed return electrode condition, unit current is injected into the remaining electrodes in sequence to generate a guiding field matrix. With the help of the linear relationship between electric field and current, this matrix can efficiently predict the electric field distribution characteristics under any current injection combination. In terms of electrode arrangement, a total of 64 electrodes are used, all arranged on the scalp surface according to the internationally accepted 10-10 EEG system.

[0010] In step 2, the electric field intensity modeling adopts the volume integral weighted average strategy, first defining the target area Ω t The average electric field intensity within the target area is calculated, and its integral average over the volume of the target area is calculated. Taking into account that the target area is usually small and the electric field distribution is approximately uniform, the expression can be simplified to a discrete form, and linear mapping is performed between the electrode current vector and the electric field matrix. This mapping involves the construction of a semi-positive symmetric matrix, which consists of a guide field matrix and a weight volume matrix. Among them, the weight matrix is ​​a diagonal matrix, the non-zero elements represent the volume of each voxel in the target area, and the remaining positions are set to zero, which is used to accurately weight the electric field influence of the target area. In addition, the modeling goal of the electric field optimization problem is to minimize the total electric field energy in the non-target area while ensuring that the electric field intensity in the target area Ω reaches the specified value. This optimization problem combines multiple practical constraints, and the modeling includes minimizing the total electric field energy outside the target area, constraining the average electric field intensity in the target area Ω to reach a preset value t; satisfying Kirchhoff current conservation, that is, the sum of the electrode currents is zero; limiting the total injection current to no more than the maximum current threshold 2I max ; Limit the current injected into a single electrode to no more than the threshold I ind,max The above constraints together constitute a non-convex optimization problem, and the model variables are the electrode current vector x∈R n , the objective function is the electric field energy term, and the constraints involve the electric field mean, current balance and limiting conditions.

[0011] In step 3, a convex-concave programming method is introduced to solve the non-convex optimization problem. This method converts the electric field intensity constraint in the target area into a convex-concave function difference and performs first-order Taylor linearization on the concave function. This method transforms the original non-convex problem into a series of convex optimization subproblems that are solved iteratively.

[0012] In step 3, to satisfy the CCP solution framework, the original equality constraint from step 2 is converted into an equivalent inequality constraint. Subsequently, in each iteration, a linear approximation is performed on the inequality constraint at the current solution point, constructing the corresponding convex subproblem. In each iteration, a sparsity control term is added to limit the maximum number of activated electrodes, and a slack variable and its penalty term are introduced to improve initial point robustness and convergence. The penalty coefficient increases exponentially with each iteration until it reaches an upper limit.

[0013] In step 4, the optimization solution uses a constrained eigenvalue initialization method to obtain the initial current distribution, which is then normalized using the total current and single current limit constraints. To enhance the global optimal search capability, each optimization is repeated multiple times in parallel with different initial points, and the convergence criterion is set as the change in the objective function and the electric field intensity below a threshold.

[0014] In step 4, if the preset target intensity t is too large, making the problem infeasible, the algorithm automatically switches to maximizing the electric field intensity in the target area. Furthermore, to support multi-objective electric field control, the weight matrices of multiple target areas can be merged to uniformly guide the optimization direction.

[0015] The advantages of the transcranial electrical stimulation control method and system with multi-objective electric field strength optimization proposed in the present invention include:

[0016] 1. The present invention solves the problems of unclear direction and / or limited effect of deep brain stimulation in traditional methods by directly controlling the electric field intensity of the target area rather than the specific directional component.

[0017] 2. The present invention can simultaneously apply electric field intensity control to multiple brain regions, effectively meeting the needs of multifunctional network intervention or cross-regional collaborative treatment in neuroregulation, and breaking through the bottleneck of existing methods that are difficult to optimize multiple targets. BRIEF DESCRIPTION OF THE DRAWINGS

[0018] Figure 1 This is a flowchart of a multi-objective electric field strength intelligent optimization method and system based on individualized closed-loop neural regulation according to an embodiment of the present invention.

[0019] Figure 2 A schematic diagram of a head model network according to one embodiment of the present invention.

[0020] Figure 3 Schematic diagram of an optimized injection scheme according to an embodiment of the present invention.

[0021] Figure 4 1 is a heat map of electric field distribution in multiple target areas according to one embodiment of the present invention.

[0022] Figure 5 Schematic diagram of the non-target area suppression effect according to an embodiment of the present invention. DETAILED DESCRIPTION

[0023] like Figure 1As shown, in a transcranial electrical stimulation control method based on multi-objective electric field strength optimization according to an embodiment of the present invention, by constructing a guidance field matrix, constructing a target electric field strength control expression, introducing a convex-concave programming solution strategy, and a sparsity constraint control module, precise control and personalized adjustment of the electric field distribution of multiple brain regions are achieved; at the same time, slack variables and penalty mechanisms are introduced to improve the convergence of non-convex optimization problems and the feasibility of solutions; and combined with the electrode activation number control strategy, the intelligent optimization of the electrode configuration scheme in the multi-channel TES system is completed, and finally an electrical stimulation scheme that meets safety, focus and individual adaptability is obtained.

[0024] The following describes in detail the transcranial electrical stimulation control method for multi-objective electric field strength optimization according to the present invention in conjunction with specific embodiments. The method includes:

[0025] Step S1: Select a fixed electrode as the return channel, and inject unit current into each of the other multiple electrodes in turn. Use the finite element method to calculate the corresponding electric field response, and then construct an electric field guidance field matrix to quickly predict the electric field distribution under any current combination. Figure 4 1 is a multi-target brain region distribution map according to one embodiment of the present invention.

[0026] According to a specific embodiment, the fixed return electrode is the Cz electrode in the scalp EEG 10-10 electrode system, and multiple electrodes are sequentially injected with unit current (such as 1 mA), and the electric field distribution response under current injection is calculated based on the finite element method (FEM) in a three-dimensional head model. Using the linear response characteristics of the electric field to the current input, the electric field response results of all electrodes under current injection are stacked in columns, and finally the electric field guidance field matrix A∈R is constructed. m×n , where m is the number of voxels of interest in the three-dimensional grid, n is the number of additional electrodes, R is a set of real numbers, R m×n represents a real matrix with m rows and n columns.

[0027] A specific embodiment of this step includes:

[0028] Step S1.1: Electrode placement includes selecting n (in one embodiment, n=64) electrodes and evenly distributing them on the scalp surface according to the internationally accepted EEG 10-10 placement system, covering the frontal lobe, parietal lobe, occipital lobe, temporal lobe and other regions.

[0029] Step S1.2: Boundary condition setting includes: To achieve accurate simulation, Neumann boundary conditions are set in the non-electrode area (current flux is zero), Dirichlet boundary conditions are set in the electrode-scalp contact area, and a linear electrode model (i.e., equipotential assumption) is adopted. On the head model mesh (a three-dimensional finite element mesh model generated by tissue segmentation, surface reconstruction, and volume meshing of the head structure magnetic resonance image, such as Figure 2 As shown in FIG, each tissue (including scalp, skull, cerebrospinal fluid, gray matter, white matter, and eyeball) is assigned its own bioconductivity attribute, and the conductivity of each tissue is set according to the built-in conductivity table of SimNIBS.

[0030] Step S1.3: Electric field response calculation includes: in the electric field guidance matrix construction stage, one electrode is selected as a fixed return current channel, and a unit current (1 mA) is injected into the remaining n-1 electrodes one by one, with only one injection electrode activated at a time. Figure 3 The schematic diagram of an example of the optimized injection scheme according to the present invention. Under each single injection condition, the MKLPARDISO solver is used to numerically solve the entire finite element model and output the electric field vector field E i (i is the electrode number). After n-1 repetitions, the unit injection responses of all electrodes form the electric field guidance matrix A = [E1, E2, ..., En_1], where Ei is the column vector, and each column represents the electric field response distribution of an injection electrode.

[0031] Step S1.4: Constructing a prediction model includes: using the linear superposition characteristics of the electric field, injecting a current vector x∈R into any electrode n , to quickly estimate the electric field distribution E in the target area total =Ax, which greatly improves the computational efficiency in the subsequent optimization steps. Where A represents the electric field guidance matrix, x is the column vector of the electrode injection current, n does not include the fixed return electrode, R n Represents a real column vector of length n, where each element is the current injection value of an electrode.

[0032] In one embodiment, the head model grid (Ernie model provided by SimNIBS 3.1, with a spatial resolution of 1 mm) 3 The types of tissues include white matter, gray matter, cerebrospinal fluid, skull, scalp, and eyeball, a total of 6 categories. The electrode positions are strictly aligned with the grid nodes. A total of n = 64 electrodes are arranged on the standard coordinates of the EEG 10-10 layout system of Fp1, Fp2, F3, F4, Cz, T7, T8, O1, and O2. The final output guidance field matrix size is m × 64, where m is the number of mask voxels of all target areas in the volume model.

[0033] Step S2: Based on the construction of the guiding field matrix, the electric field intensity expression is defined for the target area, and a constrained optimization objective function is established.

[0034] Step S2.1: Modeling the electric field intensity expression in the target area

[0035] Include: Set the target area voxel set to Ω t , the guidance field matrix is ​​A∈R m×n , the electrode current injection vector is x∈R n , the electric field at voxel i is E i =A i x, where A i is the row vector of the electric field mapping of the i-th voxel.

[0036] The average electric field strength can be approximated by the following expression:

[0037]

[0038] in

[0039] in represents the electric field response weight of the i-th non-target voxel, v is the voxel volume, that is, the physical space occupied by the unit in the grid, which is used to achieve weighted averaging. represents the transpose of the current injection vector x, i.e., a 1×n row vector.

[0040] Step S2.2: Modeling the electric field energy in the non-target area includes: In order to suppress the over-activation of the non-target area, the electric field weighting matrix Q of the non-target area is introduced, which is constructed by the following function:

[0041]

[0042] To reflect the total electric field energy of the entire brain area, where Ω non-target Represents the set of all voxels that do not belong to the target area in the 3D head model.

[0043] Step S2.3: Objective function construction includes: minimizing the electric field energy in the non-target area as the optimization goal, while constraining the electric field intensity in the target area to be equal to the set value t:

[0044]

[0045] Among them, formula (6) is used to minimize the total electric field energy outside the target area, and formula (7) is used to control the average electric field intensity in the target area Ωt to a given value t. represents the transpose of the current injection vector x, i.e., a 1×n row vector.

[0046] Step S2.4: Modeling additional constraints

[0047] include:

[0048]

[0049] ‖x‖1≤2I max (9)

[0050] -I ind,max ≤x i ≤I ind,max (10)

[0051] ‖x‖0≤N (11)

[0052] Among them, formula (8) is used for current conservation constraint to ensure that the total input current is zero; formula (9) is used for total current limitation to limit the total current intensity; formula (10) is used for single electrode limitation to limit the current intensity of a single electrode; formula (11) is used for sparse electrode activation number to control the total number of electrodes. Among them, x represents the electrode current injection vector, 1 is a column vector with all elements set to 1, and I max Indicates the upper limit of the maximum current intensity of a single pole, I ind,max represents the maximum allowed current of each electrode, ‖·‖0 represents the number of non-zero elements in the vector, N represents the upper limit of the number of electrodes allowed to be activated simultaneously, and ‖·‖1 represents the total current intensity.

[0053] In one embodiment, the electric field strength in the target area is set to 0.2 V / m(x T Qx=0.2 2 ), the total current limit is 4mA (‖x‖1≤2×2mA), and the upper limit of the single electrode current is 1mA (-1mA≤x i ≤1mA), the maximum number of active electrodes allowed is 10 (‖x‖0≤10).

[0054] Through the electric field intensity and energy modeling in this step, the optimization problem is standardized as a constrained quadratic optimization problem for subsequent solution using the convex-concave programming method.

[0055] Step S3: The convex-concave programming (CCP) method is used to solve the non-convex constrained optimization problem to ensure the accessibility of the electric field intensity constraint and improve the optimization convergence performance.

[0056] According to a specific embodiment, since the optimization model contains the target electric field strength equality constraint However, the constraint term is a non-convex quadratic form, which makes the original optimization problem non-convex as a whole. Conventional convex optimization algorithms cannot directly solve it. To this end, the inventors use the Convex-Concave Programming (CCP) method to transform the non-convex constraints into an iterative approximation convex problem solving framework through a linearization strategy. The application of linearized non-convex constraints to the specific electric field optimization problem is an innovative point of the present invention, and its beneficial effects include:

[0057] 1. Introducing the target intensity equality constraint into the TES electric field optimization problem;

[0058] 2. Based on the CCP framework, by introducing strategies such as slack variables, dynamic penalty terms, and sparse activation control, the original problem's engineering difficulties of infeasibility, difficulty in convergence, and inability to sparsify are resolved.

[0059] 3. Compared with traditional optimization methods, this method is more effective in balancing the stimulation intensity of multiple target areas and limiting the energy of non-target areas. Figure 5 Schematic diagram of the non-target area suppression effect according to an embodiment of the present invention.

[0060] This step specifically includes:

[0061] Step S3.1: Constraint relaxation process

[0062] Includes: constraining the target strength equation Convert to inequality relaxation form:

[0063]

[0064] This process can ensure that the solution space contains the solution corresponding to the original equality constraint and provide conditions for linearization processing.

[0065] Step S3.2: Linearization approximation processing includes: at the iteration point x (k) At the constraint function Perform a first-order Taylor expansion to obtain the linearized constraints:

[0066]

[0067] This inequality constitutes a convex constraint, which is used to replace the original non-convex constraint to solve the problem in each round of iteration.

[0068] Step S3.3: Introducing slack variables and penalty terms includes: In order to ensure the feasibility of the solution in the optimization initialization phase, the slack variable s≥0 and the penalty term ξ are introduced. k , update the optimization objective function to:

[0069]

[0070] At the same time, modify the linear constraints to:

[0071]

[0072] Among them, the penalty factor ξ k Dynamically increment as follows:

[0073] ξ k+1 =min(μξ k ,ξ max ) (16)

[0074] For example, the initial setting The maximum upper limit is ξ max =ξ0×10 4 , the update step factor is μ=2.

[0075] Step S3.4: The multi-point initialization and iterative optimization process includes: using a multi-point initialization strategy (e.g., 20 different initialization points), each of which is obtained by solving a constrained eigenvalue problem:

[0076]

[0077] The current distribution is then scaled based on the constraints. Scaling adjusts the initial current to ensure it meets the actual current limit while providing a reasonable starting point for subsequent optimization. CCP iterations are performed at each initial point, and the solution corresponding to the optimal target value is selected as the final current distribution solution.

[0078] Step S4: Through the electrode sparse activation constraint, the feasibility of the final stimulation scheme is controlled, and the solution termination criterion is set to ensure that the optimization process converges within a limited number of iterations.

[0079] In practical transcranial electrical stimulation systems, to meet the physical limitations of clinical or engineering implementation, such as the number of electrode channels, electrode switching accuracy, and stimulator hardware capacity, it is typically required that only a limited number of electrodes be used for each stimulation. Therefore, the present invention introduces a sparse activation control constraint into the optimization model to ensure that only at most N electrodes receive non-zero current in the final solution.

[0080] This step specifically includes:

[0081] Step S4.1: Modeling electrode sparse activation constraints

[0082] Include: Add e0 norm restriction in optimization constraints:

[0083] ‖x‖0≤N (18)

[0084] Where ‖x‖0 represents the number of nonzero elements in the current vector, i.e., the number of active electrodes; N is the maximum number of electrodes allowed to be activated (e.g., 10) (in this embodiment, the total number of electrodes is 64). This sparsity constraint is used to ensure that the optimization results are physically feasible, facilitating device implementation and electrode switching control.

[0085] Step S4.2: Solving the sparse optimization problem by branch and bound includes: Since the l0 norm is a non-convex function, the sparse optimization problem is essentially a combinatorial optimization problem. The present invention uses a branch and bound algorithm to solve the problem. The specific process includes:

[0086] 1. Divide the candidate electrode combination space into subsets during the main problem iteration;

[0087] 2. Solve the convex optimization problem for each subset separately (i.e., the operation of step S3);

[0088] 3. Eliminate suboptimal solutions through upper and lower bound strategies;

[0089] 4. Finally determine the current injection combination corresponding to the optimal sparse solution.

[0090] Step S4.3: Post-processing and sparse pruning of the solution include: To further enhance sparsity, the present invention adopts a soft threshold strategy to remove and merge components with amplitudes close to zero in the current vector after the objective function converges, simplify the stimulus configuration, and improve the robustness of actual system control.

[0091] Step S4.4: Setting the solution termination condition

[0092] In order to prevent the optimization from falling into invalid iterations, the present invention sets the following termination criteria:

[0093] The relative rate of change of the objective function is less than the set threshold:

[0094]

[0095] The relative rate of change of electric field intensity constraint is less than the set threshold:

[0096]

[0097] in

[0098] Step S4.5: Optimization output and visualization presentation include: outputting the final stimulation plan after optimization is completed, including the electrode number, current direction and magnitude; and calling the visualization module to display the electric field distribution in the target area and the inhibition effect in the non-target area, so as to facilitate clinicians or engineers to evaluate the feasibility of the plan.

[0099] Although the embodiments of the present application have been shown and described above, it can be understood that the above embodiments are exemplary and cannot be understood as limitations on the present application. Ordinary technicians in this field can change, modify, replace and modify the above embodiments within the scope of the present application.

Claims

1. A transcranial electrical stimulation control method with multi-objective electric field intensity optimization, characterized in that The following steps are involved: Step 1: Construct an electric field guidance matrix, which involves selecting a fixed electrode as a return current channel, injecting unit current into each of the remaining electrodes in sequence, and calculating the corresponding electric field response using the finite element method; Step 2: Based on the definition of the average electric field intensity in the target area, a relationship model between the electrode current vector and the electric field intensity is established, and an optimization objective function is constructed to minimize the electric field energy in the non-target area. At the same time, the target electric field intensity constraint and multiple physical restrictions are set. Step 3: Use the convex-concave programming method to linearize the non-convex strength constraint, introduce slack variables and penalty terms, and construct an iterative convex optimization subproblem; Step 4: Output step, which controls the number of electrode activations through sparsity constraints and sets convergence criteria, ultimately outputting an optimized current solution that meets the intensity target, current safety limit, and electrode sparsity.

2. The transcranial electrical stimulation control method according to claim 1, characterized in that Step S1 includes: Step S1.1: Electrode arrangement, including: selecting n electrodes and evenly distributing them on the scalp surface, covering the frontal lobe, parietal lobe, occipital lobe, and temporal lobe. Step S1.2: Setting boundary conditions, including: Neumann boundary conditions are set in the non-electrode region. Dirichlet boundary conditions are set in the electrode-scalp contact area, and a linear electrode model is used. In the head model mesh, each tissue, including the scalp, skull, cerebrospinal fluid, gray matter, white matter, and eyeball, is assigned its own bioconductivity attribute, and the conductivity of each tissue is set according to the built-in conductivity table of SimNIBS. The head model mesh is a three-dimensional finite element mesh model generated by performing tissue segmentation, surface reconstruction, and volume meshing on the magnetic resonance imaging of the head structure. Step S1.3: Determine the electric field response, including: in the electric field guidance matrix construction stage, select one electrode as a fixed return current channel, and inject unit current into the remaining n-1 electrodes one by one, activating only one injection electrode at a time, Under each individual injection condition, the entire finite element model is numerically solved using the MKL PARDISO solver to output the electric field vector field E i , where i is the electrode number. After n-1 repetitions, the unit injection responses of all electrodes form the electric field guidance field matrix A = [E1, E2, ..., En_1], where Ei is the column vector, and each column represents the electric field response distribution of an injection electrode. Step S1.4: Construct a prediction model, including: using the linear superposition characteristics of the electric field, injecting a current vector x∈R into any electrode n , to quickly estimate the electric field distribution E in the target area total =Ax, which greatly improves the computational efficiency in the subsequent optimization steps, where: A represents the electric field guidance matrix, x is the column vector of the electrode injection current, n does not include the fixed return electrode, R n represents a real column vector of length n, where each element is the current injection value of an electrode, Step S2 includes: Step S2.1: Modeling the electric field intensity expression of the target area, including: setting the voxel set of the target area to Ω t , the guidance field matrix is ​​A∈R m×n , the electrode current injection vector is x∈R n , the electric field at voxel i is E i =A i x, where A i is the electric field mapping row vector of the i-th voxel, The average electric field strength can be expressed approximately as: in in represents the electric field response weight of the i-th non-target voxel, v is the voxel volume, that is, the physical space occupied by the unit in the grid, which is used to achieve weighted averaging. represents the transpose of the current injection vector x, i.e. a 1×n row vector, Step S2.2: Modeling the electric field energy in the non-target area, including: In order to suppress the over-activation of the non-target area, introducing the electric field weighting matrix Q of the non-target area, with the function: Reflects the total electric field energy of the entire brain area, where Ω non-target represents the set of all voxels that do not belong to the target area in the 3D head model, Step S2.3: Construct an objective function, including: minimizing the electric field energy in the non-target area as the optimization goal, while constraining the electric field intensity in the target area to be equal to the set value t: Among them, formula (6) is used to minimize the total electric field energy outside the target area, and formula (7) is used to control the average electric field intensity in the target area Ωt to a given value t. represents the transpose of the current injection vector x, i.e. a 1×n row vector, Step S2.4: Model additional constraints, including: ‖x‖1≤2I max (9) -I ind,max ≤x i ≤I ind,max (10) ‖x‖0≤N (11) Among them, formula (8) is used for current conservation constraint to ensure that the total input current is zero; formula (9) is used for total current limitation to limit the total current intensity; formula (10) is used for single electrode limitation to limit the current intensity of a single electrode; formula (11) is used for sparse electrode activation number to control the total number of electrodes. x represents the electrode current injection vector, 1 is a column vector whose elements are all 1, I max Indicates the upper limit of the maximum current intensity of a single pole, I ind,max Indicates the maximum allowable current for each electrode. ‖·‖0 represents the number of non-zero elements in the vector, N represents the upper limit of the number of electrodes allowed to be activated simultaneously. ‖·‖1 represents the total current intensity, Step S3 includes: Step S3.1: Perform constraint relaxation, including: constraining the target strength equation Convert to inequality relaxation form: To ensure that the solution space contains the corresponding solution of the original equality constraint, and to provide conditions for linearization processing, Step S3.2: Perform linear approximation processing, including: at the iteration point x (k) At the constraint function Perform a first-order Taylor expansion to obtain the linearized constraints: This inequality constitutes a convex constraint, which is used to replace the original non-convex constraint to solve the problem in each round of iteration. Step S3.3: Introduce slack variables and penalty terms, including: In order to ensure the feasibility of the solution in the optimization initialization phase, introduce the slack variable s≥0 and the penalty term ξ k , update the optimization objective function to: At the same time, modify the linear constraints to: Among them, the penalty factor ξ k Dynamically increment as follows: x k+1 =min(μξ k ,x max ) (16), Step S3.4: Perform multi-point initialization and iterative optimization, where each point corresponds to a possible initial current distribution scheme, including: using a multi-point initialization strategy, solving a constrained eigenvalue problem for each initialization point to obtain: And it is scaled according to the constraints, including adjusting the initial current size ratio to make it conform to the physical limitations of the current and provide a reasonable starting point for subsequent optimization. CCP iteration is performed for each initialization point, and the corresponding optimization solution with the best objective function value is finally selected as the final current distribution solution. The initialization point refers to the initial current vector used to start the CCP iteration; the optimized solution refers to the optimal current distribution result obtained under all constraints. Step S4: includes: Step S4.1: Perform electrode sparse activation constraint modeling, including adding l0 norm constraints to the optimization constraints: ‖x‖0≤N (18) Where ‖x‖0 represents the number of non-zero elements in the current vector, that is, the number of active electrodes; N is the maximum number of electrodes allowed to be activated, Step S4.2: Solve the sparse optimization problem by branch and bound, including: S4.2.

1. In the main problem iteration, divide the candidate electrode combination space into subsets; S4.2.

2. Solve the convex optimization problem for each subset separately, i.e., perform the operation of step S3; S4.2.

3. Eliminate suboptimal solutions through upper and lower bounds; S4.2.

4. Determine the current injection combination corresponding to the optimal sparse solution, Step S4.3: Post-processing and sparse pruning of the solution, including: to further enhance sparsity, adopting a soft threshold strategy, removing and merging components of the current vector with amplitudes close to zero after the objective function converges, simplifying the stimulus configuration, and improving the robustness of the actual system control. Step S4.4: Set the termination conditions, including: To prevent the optimization from falling into invalid iterations, set the following termination criteria: The relative rate of change of the objective function is less than the set threshold: The relative rate of change of electric field intensity constraint is less than the set threshold: in Step S4.5: Optimize output and visualize the results, including: After the optimization is completed, the final stimulation plan is output, including the electrode number, current direction and size.

3. The transcranial electrical stimulation control method according to claim 2, characterized in that Step S4.5 further includes: Call the visualization module to display the electric field distribution in the target area and the suppression effect in the non-target area.

4. The transcranial electrical stimulation control method according to claim 2, wherein: The electrode arrangement adopts the international standard EEG 10-10 system, with n=64 electrodes.

5. A transcranial electrical stimulation control system with multi-objective electric field strength optimization, characterized in that include: Electric field guidance matrix construction module: used to construct the electric field guidance matrix, including selecting a fixed electrode as the return channel, injecting unit current into each of the remaining electrodes in turn, and using the finite element method to calculate the corresponding electric field response; Modeling module: used to establish the relationship model between the electrode current vector and the electric field intensity based on the definition of the average electric field intensity in the target area, and to construct the optimization objective function to minimize the electric field energy in the non-target area, while setting the target electric field intensity constraint and multiple physical restrictions; Convex-concave programming module: It uses convex-concave programming to linearize non-convex strength constraints, introduce slack variables and penalty terms, and construct iterative convex optimization subproblems; Output module: It is used to control the number of electrode activations through sparsity constraints and set convergence criteria, and ultimately output an optimized current solution that meets the intensity target, current safety limit and electrode sparsity.

6. The transcranial electrical stimulation control system according to claim 5, characterized in that the electric field The Guided Field Matrix building block includes components for performing the following operations: Step S1.1: Electrode arrangement, including: selecting n electrodes and evenly distributing them on the scalp surface, covering the frontal lobe, parietal lobe, occipital lobe, and temporal lobe. Step S1.2: Setting boundary conditions, including: Neumann boundary conditions are set in the non-electrode region. Dirichlet boundary conditions are set in the electrode-scalp contact area, and a linear electrode model is used. In the head model mesh, each tissue, including the scalp, skull, cerebrospinal fluid, gray matter, white matter, and eyeball, is assigned its own bioconductivity attribute, and the conductivity of each tissue is set according to the built-in conductivity table of SimNIBS. The head model mesh is a three-dimensional finite element mesh model generated by performing tissue segmentation, surface reconstruction, and volume meshing on the magnetic resonance imaging of the head structure. Step S1.3: Determine the electric field response, including: in the electric field guidance matrix construction stage, select one electrode as a fixed return current channel, and inject unit current into the remaining n-1 electrodes one by one, activating only one injection electrode at a time, Under each individual injection condition, the entire finite element model is numerically solved using the MKL PARDISO solver to output the electric field vector field E i , where i is the electrode number. After n-1 repetitions, the unit injection responses of all electrodes form the electric field guidance field matrix A = [E1, E2, ..., En_1], where Ei is the column vector, and each column represents the electric field response distribution of an injection electrode. Step S1.4: Construct a prediction model, including: using the linear superposition characteristics of the electric field, injecting a current vector x∈R into any electrode n , to quickly estimate the electric field distribution E in the target area total =Ax, which greatly improves the computational efficiency in the subsequent optimization steps, where: A represents the electric field guidance matrix, x is the column vector of the electrode injection current, n does not include the fixed return electrode, R n represents a real column vector of length n, where each element is the current injection value of an electrode, The Modeling module includes sections for performing the following operations: Step S2.1: Modeling the electric field intensity expression of the target area, including: setting the voxel set of the target area to Ω t , the guidance field matrix is ​​A∈R m×n , the electrode current injection vector is x∈R n , the electric field at voxel i is E i =A i x, where A i is the electric field mapping row vector of the i-th voxel, The average electric field strength can be expressed approximately as: in in represents the electric field response weight of the i-th non-target voxel, v is the voxel volume, that is, the physical space occupied by the unit in the grid, which is used to achieve weighted averaging. represents the transpose of the current injection vector x, i.e. a 1×n row vector, Step S2.2: Modeling the electric field energy in the non-target area, including: In order to suppress the over-activation of the non-target area, introducing the electric field weighting matrix Q of the non-target area, with the function: Reflects the total electric field energy of the entire brain area, where Ω non-target represents the set of all voxels that do not belong to the target area in the 3D head model, Step S2.3: Construct an objective function, including: minimizing the electric field energy in the non-target area as the optimization goal, while constraining the electric field intensity in the target area to be equal to the set value t: Among them, formula (6) is used to minimize the total electric field energy outside the target area, and formula (7) is used to control the average electric field intensity in the target area Ωt to a given value t. represents the transpose of the current injection vector x, i.e. a 1×n row vector, Step S2.4: Model additional constraints, including: ‖x‖1≤2i max (9) -I ind,max ≤x i ≤I ind,max (10) ‖x‖0≤N (11) Among them, formula (8) is used for current conservation constraint to ensure that the total input current is zero; formula (9) is used for total current limitation to limit the total current intensity; formula (10) is used for single electrode limitation to limit the current intensity of a single electrode; formula (11) is used for sparse electrode activation number to control the total number of electrodes. x represents the electrode current injection vector, 1 is a column vector whose elements are all 1, I max Indicates the upper limit of the maximum current intensity of a single pole, I ind,max Indicates the maximum allowable current for each electrode. ‖·‖0 represents the number of non-zero elements in the vector, N represents the upper limit of the number of electrodes allowed to be activated simultaneously. ‖·‖1 represents the total current intensity, The convex-concave programming module includes components for performing the following operations: Step S3.1: Perform constraint relaxation, including: constraining the target strength equation Convert to inequality relaxation form: To ensure that the solution space contains the corresponding solution of the original equality constraint, and to provide conditions for linearization processing, Step S3.2: Perform linear approximation processing, including: at the iteration point x (k) At the constraint function Perform a first-order Taylor expansion to obtain the linearized constraints: This inequality constitutes a convex constraint, which is used to replace the original non-convex constraint to solve the problem in each round of iteration. Step S3.3: Introduce slack variables and penalty terms, including: In order to ensure the feasibility of the solution in the optimization initialization phase, introduce the slack variable s≥0 and the penalty term ξ k , update the optimization objective function to: At the same time, modify the linear constraints to: Among them, the penalty factor ξ k Dynamically increment as follows: x k+1 =min(μξ k ,x max ) (16), Step S3.4: Perform multi-point initialization and iterative optimization, where each point corresponds to a possible initial current distribution scheme, including: using a multi-point initialization strategy, solving a constrained eigenvalue problem for each initialization point to obtain: And scale it according to the constraints, including adjusting the initial current size ratio to make it conform to the physical limitations of the current and provide a reasonable starting point for subsequent optimization. CCP iteration is performed for each initialization point, and the corresponding optimization solution with the best objective function value is finally selected as the final current distribution solution. The initialization point refers to the initial current vector used to start the CCP iteration; the optimized solution refers to the optimal current distribution result obtained under all constraints. The output module includes sections for performing the following operations: Step S4.1: Perform electrode sparse activation constraint modeling, including adding l0 norm constraints to the optimization constraints: ‖x‖0≤N (18) Where ‖x‖0 represents the number of non-zero elements in the current vector, that is, the number of active electrodes; N is the maximum number of electrodes allowed to be activated, Step S4.2: Solve the sparse optimization problem by branch and bound, including: S4.2.

1. In the main problem iteration, divide the candidate electrode combination space into subsets; S4.2.

2. Perform the convex optimization problem solution on each subset separately, i.e., perform the convex-concave programming module. S4.2.

3. Eliminate suboptimal solutions through upper and lower bounds; S4.2.

4. Determine the current injection combination corresponding to the optimal sparse solution, Step S4.3: Post-processing and sparse pruning of the solution, including: to further enhance sparsity, adopting a soft threshold strategy, removing and merging components of the current vector with amplitudes close to zero after the objective function converges, simplifying the stimulus configuration, and improving the robustness of the actual system control. Step S4.4: Set the termination conditions, including: To prevent the optimization from falling into invalid iterations, set the following termination criteria: The relative rate of change of the objective function is less than the set threshold: The relative rate of change of electric field intensity constraint is less than the set threshold: in Step S4.5: Optimize output and visualize the results, including: After the optimization is completed, the final stimulation plan is output, including the electrode number, current direction and size.

7. The transcranial electrical stimulation control system according to claim 6, characterized in that Step S4.5 further includes: calling a visualization module to display the electric field distribution in the target area and the suppression effect in the non-target area.

8. The transcranial electrical stimulation control system according to claim 6, characterized in that: The electrode arrangement adopts the international standard EEG 10-10 system, with n=64 electrodes.

Images