Multi-objective optimization of beam transport in medical heavy ion accelerators using an improved non-dominated sorting differential evolution algorithm (NSDE)

To address the issues of high-dimensional coupling parameters easily falling into local optima and multi-objective conflicts in the beam transport of medical heavy ion accelerators, this paper proposes an improved non-dominated sorting differential evolution (NSDE) algorithm. The algorithm employs inverse learning for initialization and introduces an adaptive mechanism to adjust the mutation factor and crossover probability online, balancing exploration and exploitation. Additionally, it incorporates local enhancement based on crowding distance in particle swarm optimization (PSO) to refine non-dominated elite solutions. Large-scale experiments based on FLUKA Monte Carlo coupled simulation (nine-dimensional decision variables) have shown that the improved NSDE has increased the beam transport efficiency from the baseline of 92.42% to 99.21% (an improvement of 6.79%), while also achieving continuous improvements in key physical indicators such as the beam spot size at the end point, system power consumption, and energy retention rate. The research indicates that the proposed method exhibits significant advantages in enhancing optimization quality and maintaining robustness, making it suitable for accelerator engineering optimization that demands stringent real-time performance and multi-objective accuracy.

1 Introduction

With the widespread application of heavy ion therapy in tumor treatment, the performance optimization of medical heavy ion accelerators has become a key link in ensuring clinical treatment effects (Sokol and Durante, 2023; Yamada et al., 2022). During the process of beam transmission from the accelerator outlet to the treatment terminal, its transmission efficiency and beam spot quality directly affect the accuracy of dose distribution and treatment safety. Since there are a large number of parameters involved in the transmission system, such as magnet current, pipeline geometry, vacuum state, etc., there is a significant coupling relationship between these high-dimensional parameters, making optimization extremely difficult. In addition, there is often a performance conflict between high beam transmission efficiency and small beam spot size, and traditional methods find it difficult to take both into account. Seeking efficient and robust multi-objective optimization strategies to improve beam performance has important engineering value and clinical significance.

In order to solve the multi-objective conflict problem of medical heavy ion accelerator beam in high-dimensional coupled parameter space, the existing system has the following defects: traditional optimization methods such as NSGA-II (Ma et al., 2023) and MOEA/D (Gao et al., 2024) are prone to fall into local optimality when dealing with high-dimensional parameters, and it is difficult to effectively explore complex parameter space; and these methods have insufficient adaptive adjustment capabilities for parameters, resulting in slow convergence speed and insufficient diversity of solution sets when balancing multi-objective performance, making it difficult to balance the conflicts of key indicators such as beam transmission efficiency and beam spot size. To address the above problems, this paper proposes an improved NSDE algorithm that integrates reverse learning initialization, adaptive parameter control and PSO local refinement. The initial population with diversity is generated by reverse mapping and non-dominated crowding screening, which enhances the global search capability; in the main iteration process, the mutation factor and crossover probability are dynamically adjusted according to the fitness, achieving a balance between global exploration and local convergence; the PSO algorithm is applied to the Pareto elite solution for local refinement every 20 generations, which improves the convergence accuracy and stability of the algorithm in high-dimensional multi-objective optimization and solves the local traps and multi-objective conflicts of the existing system. After 200 iterative experiments, the improved algorithm achieved a beam transmission efficiency of 99.21% ± 0.49%, a minimum beam spot size of 5.13 mm ± 0.10 mm, a minimum power consumption of 153.68 kW ± 1.42 kW, and an energy retention rate of 94.12% ± 0.39%. Its performance is better than NSDE, NSGA-II (non-dominated sorting genetic algorithm II), MOEA/D (multi-objective evolutionary algorithm based on decomposition), MOPSO (multi-objective particle swarm optimization algorithm) and other methods, and the beam transmission efficiency has the most significant improvement compared with the original medical heavy ion accelerator, specifically an increase of 6.79%, which proves the effectiveness of this method.

Paper Contribution:

1. The algorithm is integrated with reverse learning initialization, adaptive differential parameter adjustment and PSO local refinement mechanism to improve the global exploration and local convergence capabilities of the algorithm in high-dimensional multi-objective optimization.

2. A multi-objective model including beam efficiency, beam spot size, power consumption and energy retention rate is constructed, and high-fidelity Monte Carlo simulation verification is performed in combination with FLUKA.

3. Under the 9-dimensional decision variable, the beam transmission efficiency is significantly improved, and the four objectives are comprehensively superior to NSGA-II, MOPSO, etc., which has practical engineering guidance value for the design optimization of heavy ion accelerators.

2 Related works

In the research of medical heavy ion accelerators, in order to improve the beam transmission efficiency and treatment accuracy of medical heavy ion accelerators, many particle therapy centers around the world have conducted in-depth research on accelerator beam transmission systems. Pivi M T F and other scholars introduced the system layout and beam transmission improvement work of the MedAustron particle therapy facility in Austria in carbon ion and proton beam therapy, especially the introduction of a synchronous control mechanism in the rotator system to improve the beam transmission quality and treatment consistency in the patient’s body (Pivi, 2024). Han M C and other scholars have ensured the stability and safety of the beam transmission process through precise clinical workflow and quality control processes (Han et al., 2024). Shayanmoghadam A A and other scholars have studied the energy deposition characteristics of heavy ion beams in inertial fusion and the transmission effects under different ion types and energy conditions with the help of the GEANT4 simulation platform, providing an important theoretical reference for beam transmission efficiency modeling (Shayanmoghadam et al., 2025). Liang X et al. reviewed the cutting-edge progress of heavy ion therapy in modern radiotherapy from the perspective of clinical radiation oncology, emphasizing the core role of beam quality control in treatment safety and efficacy (Liang et al., 2025). The above research has promoted the development of medical heavy ion accelerators from multiple aspects, including system architecture, clinical process, simulation modeling and treatment mechanism. Most existing methods focus on physical structure improvement and treatment strategy optimization, but lack the ability to coordinate multi-parameter and multi-objective control of beams in complex transmission channels. In particular, when there is a significant conflict between performance indicators such as transmission efficiency and beam spot size, there is still a lack of effective high-dimensional global optimization strategies. Introducing stronger intelligent optimization algorithms to solve the problems of high-dimensional coupling, multi-objective trade-offs, and local optimal traps in the transmission process and improve the overall performance of the beam system has become an important research direction.

With the emergence of intelligent algorithms, researchers have gradually introduced intelligent optimization algorithms into this field to improve transmission efficiency and system performance. Ge Y et al. proposed a multi-objective optimization strategy based on NSGA-II, which effectively reduced the variable dimension and significantly improved the efficiency of the accelerator neutron source system by combining multivariate statistical methods (Ge et al., 2024). In order to solve the problem of strong divergence of laser-accelerated proton beams, Yan Y and other scholars designed a high-efficiency beam collection system that combines permanent magnetic quadrupoles with GA (Genetic Algorithm). The system significantly improved the transmission efficiency through pre-focusing and achieved an efficiency gain of up to 6 times in the experiment, verifying the feasibility of GA in structural parameter configuration (Yan et al., 2024). Yan W et al. applied DE (Differential Evolution) to the design of beam optical systems for X-ray sources. They demonstrated through examples that DE has good global search capabilities and solution quality in high-dimensional complex spaces, improving the efficiency of the traditional design process (Yan et al., 2022). Liu Y et al. integrated DE with an expert system to create an intelligent control system, which significantly shortened the convergence time and enhanced the solution quality in synchrotron radiation nanobeamline tests, demonstrating the potential of this method in rapid beam adjustment and multi-target collaboration (Liu Y. et al., 2024). The above research has achieved certain results in parameter optimization, optical system design and transmission system structure configuration, but some strategies have not been able to fully avoid falling into local optimality, and in the case of multi-objective performance conflicts, it is difficult to balance the diversity and convergence of the solution set. At the same time, there has been no targeted and in-depth research on the optimization of the beam transmission efficiency of medical heavy ion accelerators.

In recent years, multi-objective evolutionary algorithms have shown significant advantages in the field of high-dimensional nonlinear optimization. However, in order to address the complex characteristics of beam transmission in medical heavy ion accelerators, it is necessary to further enhance the diversity, global search and local refinement capabilities of the algorithm. Reverse learning initialization generates symmetric and information-rich candidate solutions at the beginning of the population, which can better balance the spatial coverage and quickly construct a diverse Pareto approximation set (Cao et al., 2023; Ahmad et al., 2022; Qtaish et al., 2025). The adaptive mechanism monitors the evolutionary state of the population in real time, dynamically adjusts the mutation factor and crossover probability, and realizes smooth switching between global exploration and local development to avoid premature maturity (Wang et al., 2022; Huang et al., 2025). The PSO algorithm is combined with the local stage of non-dominated elite solutions, which enables the algorithm to perform efficient local search and accelerate convergence near the Pareto frontier (Ramirez-Ochoa et al., 2022; Lin et al., 2023). The above method has been proven to reduce computational overhead, improve solution quality, and effectively get rid of local optimality in applications in other fields. Existing studies have failed to apply it to the optimization of beam transmission efficiency of medical heavy ion accelerators, and are prone to falling into local optimality and high-dimensional multi-objective conflicts.

Research gaps:

• Multi-objective Bayesian active learning for MeV ultrafast electron diffraction: In the MeV-UED experiment, a multi-objective Bayesian optimization algorithm was used to achieve efficient tuning of beam parameters, reduce experimental time, and provide a new method for optimizing beam transmission efficiency (Ji et al., 2024).

• Application of multi-objective Bayesian optimization in accelerators: A multi-objective Bayesian optimization scheme was proposed, which can efficiently find the complete Pareto frontier of accelerator optimization problems, reduce the number of required observations, and provide a key step for online multi-objective optimization (Roussel et al., 2021).

• Application of Bayesian optimization in beam injection process: Studies have shown that the Bayesian optimization method performs well in beam injection tuning, outperforming traditional manual tuning and Nelder-Mead optimization algorithms, especially when dealing with up to nine tuning parameters, with higher efficiency and stability (Xu et al., 2023).

The literature meta-analysis is shown in Table 1.

Table 1

Table 1. Literature meta-analysis.

The literature reviewed in this paper covers multiple aspects such as physical architecture improvement, clinical process optimization, simulation modeling and intelligent algorithm application of medical heavy ion accelerator beam transmission system, highlighting the important role of intelligent multi-objective optimization methods in improving transmission efficiency and system performance. For example, methods based on NSGA-II, genetic algorithm (GA), differential evolution (DE) and multi-objective Bayesian optimization have shown good global search capabilities and solution quality in high-dimensional complex space, significantly improving efficiency and convergence speed; however, some methods still have problems such as falling into local optimality, insufficient solution diversity and inflexible parameter adjustment. This paper innovatively introduces an improved NSDE algorithm that combines reverse learning initialization, adaptive parameter control and PSO local reinforcement, effectively solving high-dimensional coupling, multi-objective conflict and local optimal traps, significantly improving beam transmission efficiency and system robustness, and filling the gap in the current research on intelligent optimization algorithms in the field of medical heavy ion accelerator beam transmission efficiency optimization.

3 Optimization design of beam transmission efficiency of medical heavy ion accelerators

3.1 Multi-objective optimization modeling

3.1.1 Decision variable definition

In order to achieve global coordinated optimization of multiple physical processes in the beam transmission path of medical heavy ion accelerators, this paper incorporates nine physical parameters that are strongly related to beam dynamics into the decision variable set. The expression formula of the decision vector $\mathrm{x}$ is shown in Equation 1.

$\mathrm{x}={\left[{\mathrm{x}}_1,{\mathrm{x}}_2,{\mathrm{x}}_3,{\mathrm{x}}_4,{\mathrm{x}}_5,{\mathrm{x}}_6,{\mathrm{x}}_7,{\mathrm{x}}_8,{\mathrm{x}}_9\right]}^{\mathrm{T}}(1)$

Among them, ${\mathrm{x}}_1$ represents the quadrupole magnet gradient, which adjusts the beam spot convergence; ${\mathrm{x}}_2$ represents the deflection magnet current, which directly controls the particle trajectory deflection angle. ${\mathrm{x}}_3$ represents Incident beam lateral emittance; ${\mathrm{x}}_4$ represents the divergence angle, and ${\mathrm{x}}_5$ represents Relative energy dispersion. ${\mathrm{x}}_6$ represents the magnet spacing, which regulates the coupling behavior between magnetic elements, ${\mathrm{x}}_7$ represents Adjustable beam limiter opening, ${\mathrm{x}}_8$ represents the steering angle, and ${\mathrm{x}}_9$ represents the pipeline pressure.

3.1.2 Multi-objective function design

For the global collaborative optimization of the beam performance of medical heavy ion accelerator, this paper constructs a four-dimensional multi-objective function vector $\mathrm{F}\left(\mathrm{x}\right)$ as shown in Formula 2.

$\mathrm{F}\left(\mathrm{x}\right)={\left[{\mathrm{f}}_1\left(\mathrm{x}\right),{\mathrm{f}}_2\left(\mathrm{x}\right),{\mathrm{f}}_3\left(\mathrm{x}\right),{\mathrm{f}}_4\left(\mathrm{x}\right)\right]}^{\mathrm{T}}(2)$

${\mathrm{f}}_1\left(\mathrm{x}\right)$ represents the beam transmission efficiency objective function, ${\mathrm{f}}_2\left(\mathrm{x}\right)$ represents the beam spot size objective function, ${\mathrm{f}}_3\left(\mathrm{x}\right)$ represents the system power consumption objective function, and ${\mathrm{f}}_4\left(\mathrm{x}\right)$ represents the energy retention rate objective function.

The formula of the beam transmission efficiency objective function is shown in Equation 3.

${\mathrm{f}}_1\left(\mathrm{x}\right)=-\frac{{\mathrm{N}}_{\text{out}}\left(\mathrm{x}\right)}{{\mathrm{N}}_{\text{in}}}(3)$

${\mathrm{N}}_{\text{in}}$ represents the number of incident particles, and ${\mathrm{N}}_{\text{out}}\left(\mathrm{x}\right)$ represents the number of particles received at the terminal cross section.

The formula of the beam spot size objective function is shown in Equation 4.

${\mathrm{f}}_2\left(\mathrm{x}\right)=\sqrt{\frac{1}{{\mathrm{N}}_{\text{out}}}{\displaystyle \sum _{\mathrm{i}=1}^{{\mathrm{N}}_{\text{out}}}}\left[{\left({\mathrm{x}}_{\mathrm{i}}-\bar{\mathrm{x}}\right)}^2+{\left({\mathrm{y}}_{\mathrm{i}}-\bar{\mathrm{y}}\right)}^2\right]}(4)$

${\mathrm{x}}_{\mathrm{i}}$ and ${\mathrm{y}}_{\mathrm{i}}$ represent the coordinates of the particle on the terminal plane, and $\bar{\mathrm{x}}$ and $\bar{\mathrm{y}}$ represent the center of mass position.

The formula of the system power consumption objective function is shown in Equation 5.

${\mathrm{f}}_3\left(\mathrm{x}\right)={\displaystyle \sum _{\mathrm{j}=1}^{\text{nA}}}{\mathrm{R}}_{\text{Aj}}{\mathrm{I}}_{\text{Aj}}^2+{\displaystyle \sum _{\mathrm{k}=1}^{\text{nB}}}{\mathrm{R}}_{\text{Bk}}{\mathrm{I}}_{\text{Bk}}^2(5)$

${\mathrm{R}}_{\text{Aj}}$ and ${\mathrm{R}}_{\text{Bk}}$ represent the coil resistance of the quadrupole magnet and the deflection magnet, respectively, and ${\mathrm{I}}_{\text{Aj}}^2$ and ${\mathrm{I}}_{\text{Bk}}^2$ represent the corresponding excitation currents.

The formula of the energy retention rate objective function is shown in Equation 6.

${\mathrm{f}}_4\left(\mathrm{x}\right)=-\frac{1}{{\mathrm{N}}_{\text{out}}}{\displaystyle \sum _{\mathrm{i}=1}^{{\mathrm{N}}_{\text{out}}}}\left(\frac{{\mathrm{E}}_{\mathrm{i}}}{{\mathrm{E}}_{\text{in}}}\right)(6)$

${\mathrm{E}}_{\mathrm{i}}$ represents the residual kinetic energy of the particle when it reaches the terminal, and ${\mathrm{E}}_{\text{in}}$ represents the incident energy.

The heat map of the objective function correlation coefficient is shown in Figure 1. For the calculation of correlation coefficients in multi-objective function design, this paper uses a data set obtained by uniformly sampling the values of each objective function in the decision variable space. First, multiple sample points x are generated according to a uniform distribution within the value range of the decision variable, and then the corresponding multi-objective function values F(x) are calculated. Based on these sampled objective function values, the linear correlation between the objective functions is quantified by calculating the Pearson correlation coefficient, ensuring that the correlation analysis reflects the overall trend of the objective function in the entire decision space, rather than being limited to the local optimal point.

Figure 1

Figure 1. Objective function correlation coefficient heat map.

In Figure 1, there is a significant positive correlation between beam transmission efficiency and energy retention rate, with a coefficient of 0.95. It can be seen that in the optimization process, improving beam transmission efficiency can be accompanied by an increase in energy retention rate, which is crucial to accelerator performance. There is a strong negative correlation between beam transmission efficiency and beam spot size, with a coefficient of −0.80, indicating that when beam transmission efficiency increases, the beam spot size decreases. There is a strong positive correlation between power consumption and beam spot size, with a coefficient of 0.75, indicating that in the process of increasing the beam spot size, the system needs more power to maintain beam stability. There is a negative correlation between power consumption and energy retention rate, with a coefficient of −0.25, which reveals that increased power consumption can have a certain negative impact on energy retention, which is caused by factors such as energy loss or thermal effect.

3.2 NSDE algorithm

Based on the NSDE algorithm (Kuo et al., 2023; Farda and Thammano, 2023), this paper optimizes the beam transmission efficiency of medical heavy ion accelerators. NSDE combines the non-dominated sorting mechanism on the framework of the classic differential evolution algorithm to ensure the balance between global search and local search.

3.2.1 Initializing the population

The NSDE algorithm first randomly generates an initial population consisting of multiple individuals, where each individual represents a decision variable vector. The generation formula $x_i^{\left(0\right)}$ for the initial individuals is shown in Equation 7.

$x_i^{\left(0\right)}=x_{{\mathit{min}}_i}+\left(x_{{\mathit{max}}_i}-x_{{\mathit{min}}_i}\right)·{\alpha }_i(7)$

$x_{{\mathit{min}}_i}$ and $x_{{\mathit{max}}_i}$ represent the lower and upper bounds of the decision variables, respectively, and ${\alpha }_i$ represents a random number. The index i represents the ith decision variable, and the ith dimension in the decision variable vector corresponding to each individual in the population ranges from 1 to the total dimension of the decision variables.

3.2.2 Differential mutation and crossover operation

Differential mutation is the core step of NSDE, which generates new candidate solutions by introducing differential strategies between individuals in the population. For each individual, the formula of differential mutation $w_i^{\left(k\right)}$ is shown in Equation 8.

$w_i^{\left(k\right)}=x_{\gamma 1}^{\left(k\right)}+\beta ·\left(x_{\gamma 2}^{\left(k\right)}-x_{\gamma 3}^{\left(k\right)}\right)(8)$

$x_{\gamma 1}^{\left(k\right)}$, $x_{\gamma 2}^{\left(k\right)}$, and $x_{\gamma 3}^{\left(k\right)}$ represent three different individuals randomly selected from the current population, and $\beta$ represents the differential mutation factor.

For the crossover operation, it is used to generate new candidate solutions, which enhances the diversity of the population by exchanging the genetic information of the parent individual and the mutant individual. For each graph, the expression formula of the crossover operation is shown in Equation 9.

$u_i^{\left(k\right)}=\left\{\begin{array}{c}{w}_i^{\left(k\right)}if{\delta }_j\le C_{\delta }\\ x_i^{\left(k\right)}if{\delta }_j>C_{\delta }\end{array}\right.(9)$

${\delta }_j$ represents a random number between 0 and 1, and $C_{\delta }$ represents the crossover probability.

3.2.3 Non-dominated sorting and crowding sorting

Non-dominated sorting is a key step in solving multi-objective optimization problems in NSDE. Individuals in each generation are sorted according to their performance on multiple objectives. In a non-dominated sorting, an individual is considered to dominate another individual if it is not inferior to the other individual in all objectives and is superior to the other individual in at least one objective. The formula for the individual dominance relationship is shown in Equation 10.

$x_i\prec x_j\iff \forall k\in \left\{1,2,\dots ,m\right\},f_k\left(x_i\right)\le f_k\left(x_j\right),\exists k{f}_k\left(x_i\right)<f_k\left(x_j\right)(10)$

$\prec$ represents the dominance relationship and $f_k\left(x_i\right)$ represents the objective function value.

Each individual is assigned a non-dominated level according to its dominance relationship and is sorted by crowding within each level. NSDE divides individuals in the population into multiple levels through non-dominated sorting, and prioritizes the solution set on the Pareto frontier.

In order to avoid multiple individuals being in the same non-dominated level and resulting in insufficient convergence of the solution, NSDE introduces a crowding sorting mechanism. In the non-dominated level, the calculation formula for the crowding of each individual is shown in Equation 11.

$d_i={\displaystyle \sum _{k=1}^m}\left(\frac{f_k\left(x_{i+1}\right)-f_k\left(x_{i-1}\right)}{f_k^{\mathit{max}}-f_k^{\mathit{min}}}\right)(11)$

$f_k\left(x_{i+1}\right)$ and $f_k\left(x_{i-1}\right)$ represent the function values of the neighboring individuals on the kth target after sorting, $f_k^{\mathit{max}}$ and $f_k^{\mathit{min}}$ represent the maximum and minimum values of the target function, respectively.

3.2.4 Selection operation

In each iteration of the generation, NSDE uses non-dominated sorting and crowding sorting to select the next-generation of population. For each individual, if $x_i$ dominates $x_j$, then $x_i$ is selected. If there is no dominance relationship between $x_i$ and $x_j$, the individual with greater crowding degree is selected. For the NSDE algorithm, the algorithm stops when the maximum number of iterations is reached or the change in the solution is less than the preset convergence threshold.

In order to further improve the multi-objective optimization capability of NSDE, this paper introduces reverse learning to initialize the population, and adopts an adaptive mechanism to adaptively adjust the mutation factor and crossover probability in the algorithm, and combines the PSO algorithm for fine tuning to locally strengthen the non-dominated elite solution. The improved NSDE algorithm framework is shown in Figure 2.

Figure 2

Figure 2. Improved NSDE algorithm framework.

In Figure 2, the overall optimization framework of the improved NSDE algorithm covers the initial population generation, main iteration process, adaptive control mechanism and PSO local reinforcement module. The improved NSDE algorithm generates a double-sized initial population through random initialization and reverse learning mapping, and combines non-dominated sorting and crowding evaluation to select an initial solution set with better diversity. Then, differential mutation and crossover operations are performed in each generation to generate candidate solutions, and the parent population is combined for non-dominated sorting and a new generation of population is selected. During the iteration process, the algorithm dynamically adjusts the mutation factor and crossover probability to achieve a balance between perturbation and convergence, and activates PSO local optimization every several generations to fine-tune the elite solutions in the non-dominated frontier to enhance the distribution quality of the Pareto frontier and the local accuracy of the solution. The entire process continues to iterate until the maximum number of iterations or convergence conditions are met, and finally a high-quality multi-objective optimal solution set is output. In Figure 2, T represents the current number of iterations, and Tmax represents the maximum number of iterations.

3.3 Improved NSDE algorithm

3.3.1 Reverse learning initialization

For the NSDE algorithm, there is a problem of insufficient population diversity in the initial population generation process. This paper introduces a reverse learning initialization strategy to improve it (Yang and Qiu, 2023; Liu R. et al., 2024). By adjusting the initial population, reverse learning enhances the population’s coverage of the Pareto optimal solution during the optimization process, improving the effect of multi-objective optimization, especially the balance optimization between objectives such as beam transmission efficiency and beam spot size.

In this paper, in order to further improve the diversity and search effect of the initial population, the reverse learning initialization strategy is applied to all decision variables. All decision variables of each individual are reverse mapped, which helps to evenly expand the population distribution in the entire decision space, improve the coverage and diversity of the Pareto frontier, and avoid the problem of insufficient diversity caused by local variable adjustment. The key idea of reverse learning initialization is to increase the diversity of the population by reverse mapping the randomly generated initial population individuals, and to enable the generated initial solution set to better cover the Pareto frontier (Cao and Huang, 2024; Sun et al., 2023). During the initialization process of each individual, it is reversely mapped to another possible solution space through certain transformations to promote the improvement of population diversity. The process of reverse learning is as follows:

3.3.1.1 Population initialization

Before reverse mapping, an initial population is generated, and the decision variables of each individual are randomly distributed between the upper and lower bounds of the decision space.

3.3.1.2 Reverse mapping

For each individual, the formula for the reverse learning transformation $x_i^{RL}$ is shown in Equation 12.

$x_i^{RL}=x_{\mathit{max}}-\left(x_i-x_{\mathit{min}}\right)(12)$

$x_i^{RL}$ represents the new individual generated by the transformation.

The symmetric inverse mapping method used in Formula 12 faces certain challenges when dealing with restricted spaces or non-uniformly distributed feasible regions. Since the values of decision variables are strictly limited to the upper and lower bounds, directly applying the symmetric inverse mapping will cause the individuals generated after mapping to fall into the infeasible region, or fail to fully cover the feasible solution space of complex shapes. In response to this situation, a boundary correction strategy is introduced to map the out-of-bounds points back to the boundary to ensure that the individuals after inverse mapping are still in the legal search space. For non-uniformly distributed feasible regions, simple symmetric mapping cannot effectively enhance population diversity because it does not consider the impact of constraints on the solution space morphology. In this way, a sampling adjustment strategy based on probability density is combined to improve the applicability and optimization performance of reverse learning under complex constraints, and better promote the uniform coverage of the Pareto frontier by the initial population.

In this paper, a sampling adjustment strategy based on probability density is proposed to improve the mapping effect in order to solve the problem of non-uniformly distributed feasible domain in the reverse learning process. First, the sample density of the initial population in each region in the decision space is statistically analyzed to establish a regional probability density function to characterize the distribution characteristics of the solution in the feasible domain. Then, in the reverse mapping process, the sampling positions of the generated new individuals are adjusted according to the probability density function, and sampling is performed preferentially in low-density areas to enhance the coverage and diversity of the population in sparse areas. In terms of algorithm implementation, kernel density estimation is used to calculate the probability density function, and then the importance sampling technique is combined to adjust and resample the reverse mapping results to ensure that the individuals are evenly distributed after mapping and meet the boundary constraints.

In this paper, for the out-of-bounds individuals generated in the reverse learning process, a boundary correction strategy is adopted to ensure that all individuals are in the legal search space. When a decision variable exceeds its upper and lower bounds after reverse mapping, its value is immediately truncated to the corresponding boundary value, that is, if the variable is less than the lower bound, it is assigned to the lower bound; if the variable is greater than the upper bound, it is assigned to the upper bound. This truncation mapping method is simple and effective, avoiding the infeasibility problem caused by out-of-bounds solutions, ensuring the effectiveness and diversity of the initial population, while taking into account the complete coverage of the search space and preventing holes or non-uniform distribution in the solution space.

3.3.1.3 Merge the initial solution set

Merge the individual set after the reverse learning transformation with the original population to form a new population with a size of 2N. The new population $X_{total}$ is shown in Formula 13.

$X_{total}=X\cup X^{RL}(13)$

$X^{RL}$ represents the individual set after reverse learning transformation, and X represents the original population.

3.3.1.4 Selecting high-quality individuals

After merging the populations, the non-dominated sorting and crowding evaluation methods are used to select N high-quality individuals from the merged initial population to form the final initial population. The steps for selecting high-quality individuals are to first perform non-dominated sorting to ensure that non-dominated solutions are selected from the population, and then sort the individuals in the same non-dominated level according to the crowding degree, select individuals with higher crowding degree, and maintain the diversity of the solution set.

3.3.2 Adaptive mechanism

In order to improve the global search capability of the NSDE algorithm and avoid falling into the local optimum when solving high-dimensional coupling problems, this paper introduces an adaptive mechanism (Wongsa et al., 2024; Cheng et al., 2024) to dynamically adjust the mutation factor and crossover probability. The adaptive mechanism adjusts according to the search status of the current population during the algorithm operation, optimizes the parameter configuration, and improves the convergence of the algorithm and the quality of the solution (Li and Tam, 2024; Zhong and Yu, 2024).

In this paper, the reverse learning initialization strategy provides the algorithm with a wider and evenly distributed initial solution space by expanding the diversity of the initial population, effectively avoiding the risk of the initial search falling into the local area; and the adaptive mechanism dynamically adjusts the mutation factor and crossover probability according to the fitness and search status of the population during the iteration process, further promoting the maintenance of population diversity and the refined optimization of local search. The two complement each other, reverse learning ensures a good starting point diversity, and the adaptive mechanism ensures the balance between diversity and convergence in the entire evolutionary process, jointly improving the search efficiency and solution quality of the improved NSDE algorithm in complex multi-objective optimization.

In this paper, the design of the adaptive formula is inspired by the control idea based on the S-type function. Its core principle is that the S-type function can smoothly map the fitness state of the population to the parameter adjustment range, achieving a regulation effect of slow growth in the global search stage and rapid convergence in the local convergence stage (Li et al., 2022). This nonlinear mapping helps to achieve a dynamic balance between exploration and utilization at different evolutionary stages, avoiding oscillation or premature convergence of the solution due to abrupt parameter changes. At the same time, the controllability and flexibility of the S-type function facilitate the combination of performance feedback mechanisms to flexibly adjust the mutation factor and crossover probability, thereby improving the convergence and robustness of the improved NSDE algorithm in complex multi-objective optimization.

The curves of the variation factor and speed factor of the S-type function changing with the number of iterations are shown in Figure 3.

Figure 3

Figure 3. Curves of variation factor and speed factor of S-type function with the number of iterations.

In Figure 3, as the number of iterations increases, both parameters gradually and smoothly decrease from a higher initial value to a lower level, reflecting the transition of the algorithm from the global search stage to the local fine search stage. The high initial value stage ensures population diversity, enhances global exploration capabilities, and prevents falling into local optimality; while the decrease in parameters in the later stage enhances the accuracy of local search, improves convergence speed and solution stability. This adaptive control mechanism effectively balances exploration and utilization, and helps to improve the convergence and robustness of the improved NSDE algorithm in complex multi-objective optimization.

For each individual, there is a fitness, which corresponds to the multi-objective function value corresponding to the individual. Individuals with higher fitness represent better current solutions, while individuals with lower fitness require more searches to improve the quality of the solution. The adaptive adjustment of the mutation factor is dynamically adjusted according to the change of individual fitness. When the fitness difference of individuals in the population is large, a larger mutation factor is allowed to promote global search, while when the fitness difference is small, the mutation factor is reduced to enhance local search. The adjustment formula of the mutation factor is shown in Equation 14.

${\beta }_i\left(t\right)={\beta }_{\mathit{min}}+\left({\beta }_{\mathit{max}}-{\beta }_{\mathit{min}}\right)·\left(\frac{1}{1+\mathit{exp}\left(\frac{\Delta f_i}{\eta }\right)}\right)(14)$

Among them, ${\beta }_{\mathit{min}}$ and ${\beta }_{\mathit{max}}$ represent the minimum and maximum values of the variation factor, $\Delta f_i$ represents the change in individual fitness,It is calculated by the difference between the current fitness of the individual and the fitness of the previous generation. and $\eta$ represents the adjustment factor.

In the adaptive adjustment of crossover probability, when the proportion of individuals with higher fitness in the population is large, the crossover probability should be lower to maintain the current high-quality solution. When the proportion of individuals with lower fitness is larger, the crossover probability should be higher to promote diversity and a wider exploration of the solution space. The adaptive adjustment formula of the crossover probability is shown in Formula 15. In Formula 15, all fitness values are scaled to the [0,1] interval by Min–Max normalization at the beginning of each iteration.

$C_i\left(t\right)=C_{\mathit{min}}+\left(C_{\mathit{max}}-C_{\mathit{min}}\right)·\left(\frac{1}{1+\mathit{exp}\left(-\theta ·\frac{f\left(x_i\right)}{\mathrm{\Sigma }f\left(x\right)}\right)}\right)(15)$

Among them, $C_{\mathit{min}}$ and $C_{\mathit{max}}$ represent the minimum and maximum values of the crossover probability, respectively, and $\theta$ represents the adjustment speed control factor. $\mathrm{\Sigma }f\left(x\right)$ represents the sum of the fitness of the population, and $f\left(x_i\right)$ represents the fitness of the individual.

Although the adaptive mutation factor driven by fitness difference performs well in balancing global exploration and local utilization, the selection of key parameters such as adjustment factor η and speed control factor θ has a great impact on the performance of the algorithm. The sensitivity of these parameters means that their optimal values vary greatly under different problem scenarios and optimization objectives, which directly affects the dynamic adjustment effect of the mutation factor and crossover probability. This paper uses the adaptive parameter update mechanism of the system to automatically adjust the values of η and θ in combination with the characteristics of the problem, avoiding premature convergence of the algorithm or reduced search efficiency due to fixed parameter settings. With the help of cross-validation, grid search and online adjustment strategies based on performance feedback, dynamic adaptation of parameters is achieved, improving the robustness and applicability of the improved NSDE algorithm in diverse and complex problems.

When adaptively updating η and θ, the overall diversity and convergence speed indicators (IGD, HV change rate) of the population are first evaluated every fixed generation (10 generations), and the current indicators are incrementally compared with the previous evaluation results; if the diversity decline rate exceeds the preset threshold τ₁ (ΔIGD/Δgeneration <−0.001), η and θ are automatically increased (multiplied by a factor of 1+μ, μ ≈ 0.1) to encourage stronger mutation and higher crossover exploration; conversely, if the convergence speed is too fast, resulting in premature solution (ΔHV/Δgeneration <τ₂, 0.0005), η and θ are gradually reduced (multiplied by a factor of 1−μ) to converge to the fine solution region. All updates are truncated within the range of [η_min,η_max] and [θ_min,θ_max] ([0.4,0.9] and [0.2,0.8]) to prevent excessive parameter fluctuations. At the same time, an exponentially weighted moving average is used to smooth historical updates to avoid severe jitters caused by single noise or abnormal evaluations, thus achieving a balance between suppressing overfitting and ensuring algorithm stability.

In the early stages of iteration, η and θ are maintained at a high level (close to their respective upper bounds) to promote diverse exploration. As generations grow, when the diversity index stabilizes or the convergence rate slows down, η and θ will gradually decay linearly or exponentially toward the middle and low range (decreasing by 1%–2% per generation) to strengthen local refinement. If a sudden drop in diversity or premature maturation occurs in the later stages, η and θ will be temporarily rebounded (increased by 5%–10%) to reinject exploration capabilities, and then continue to decay, forming a dynamic evolution curve of “exploration-convergence-re-exploration-refinement” to ensure an organic balance between global and local search.

3.3.3 PSO algorithm assisted elite local reinforcement

During the main search process of NSDE, the population converges to a sparse non-dominated frontier in the high-dimensional complex solution space, resulting in the underdevelopment of some elite solutions, which affects the accuracy and diversity of the final solution. In order to further improve the optimization effect of beam transmission efficiency and beam spot size, this paper introduces the PSO algorithm (Gad, 2022; Xu et al., 2025; Shami et al., 2022) on the non-dominated elite subset of each generation to perform local enhanced search, so as to achieve the refinement and expansion of the multi-objective Pareto frontier.

This paper chooses Particle Swarm Optimization (PSO) instead of other fast search algorithms (such as Bayesian optimization) as the local reinforcement strategy, mainly based on the following considerations. Both PSO and NSDE are population-based swarm search algorithms, capable of seamless data structure integration and directly acting on non-dominated elite subsets without additional solution encoding or complex mapping. Moreover, by introducing an improved PSO guided by crowding distance, particles can not only move towards the global optimum but also actively expand the solution distribution in the sparse region of the Pareto front, thereby improving the diversity and uniformity of the multi-objective solution set, which is highly consistent with the multi-objective optimization framework presented in this paper. Furthermore, PSO’s velocity-position collaborative update mechanism enables rapid fine-tuning of elite solutions during the local fine-tuning search phase, improving the local accuracy of the non-dominated front. In contrast, methods such as Bayesian optimization have high computational costs and are not easily parallelized in the local reinforcement of high-dimensional, multi-objective, and sparse Pareto fronts. Therefore, PSO in the improved NSDE algorithm of this paper balances local accuracy and overall diversity, achieving efficient and balanced local optimization.

In this article, the “particles” in PSO refer to candidate solutions used for local search in the algorithm, which move in the solution space mainly through speed and position updates, emphasizing group coordination and dynamic adjustment. The “individuals” in NSDE are solutions in the population, focusing on differential mutation and crossover operations in the global search and population evolution process. Although both represent points in the solution space, there are significant differences in algorithm mechanisms and role positioning.

After each generation of NSDE iteration, fast non-dominated sorting is performed to stratify the solutions in the population and select the first frontier as the elite subset. The expression formula of the elite subset $E_t$ is shown in Equation 16.

$E_t=\left\{x_i\in P_t|x_j\in P_t,x_j\prec x_i\right\}(16)$

$E_t$ represents the elite subset and $P_t$ represents the population.

As a type of optimization algorithm, the PSO algorithm guides the search in the solution space through the synergy between the individual position and the historical optimal solution of the group, and is suitable for fine-tuning the elite solution. The update rules of the particle position and velocity at the iteration k+1 time are shown in Formulas 17, 18.

$v_i^p\left(k+1\right)=\iota ·v_i^p\left(k\right)+{\vartheta }_1·{\lambda }_1·\left(p_i*-x_i^p\left(k\right)\right)+{\vartheta }_2·{\lambda }_2·\left(g^{*}-x_i^p\left(k\right)\right)(17)$

$x_i^p\left(k+1\right)=x_i^p\left(k\right)+v_i^p\left(k+1\right)(18)$

Among them, $\iota$ represents the inertia weight, ${\vartheta }_1$ and ${\vartheta }_2$ represent individual and group learning factors, and ${\lambda }_1$ and ${\lambda }_2$ represent random variables between 0 and 1 that follow a uniform distribution. $p_i*$ represents the individual historical optimal solution of the particle, and $g^{*}$ represents the global optimal solution in the current elite subset. i represents the number of the ith particle in the elite subset, which is used to identify the individual particle currently being updated; p represents the particle group number, which is used to identify the group to which the particle belongs in different PSO local optimization stages. Elite individual guidance is to adjust the particle speed and position through individual optimal solutions and global optimal solutions, guiding the search to gather in the space of better solutions, which is an important mechanism for local strengthening of PSO.

PSO is originally a single-objective optimization method. This paper introduces the crowding distance criterion to guide particles to shift to the target space with lower density and expand the diversity of the solution set. After each update of the particle position, the non-dominated frontier is reconstructed according to the objective function value, and the crowding degree of the updated elite particle is calculated. If the following Formula 19 is satisfied, the new position is retained, otherwise it is retreated.

In this paper, “particles” specifically refer to individuals used to search the solution space in the PSO local reinforcement algorithm. They are different from the “population individuals” in the main algorithm NSDE, and have nothing to do with the particles used to represent physical entities or physical processes in physics or Monte Carlo simulation. The particle here is a parameter vector representing a candidate solution. It is regarded as a point in the multidimensional optimization space and moves in the solution space through position and velocity updates to find the optimal solution or Pareto frontier. In this paper, PSO particles are mathematical individuals in the algorithm that are used to solve optimization problems, rather than particles used to simulate beams, energy deposition or Monte Carlo processes in a physical sense. PSO particles are derived from the elite subset generated by the NSDE algorithm and are further fine-tuned by PSO local search to refine and optimize the distribution and diversity of the solution set.

$\text{rank}\left(x_i^p\left(k+1\right)\right)<\text{rank}\left(x_i^p\left(k\right)\right)\text{\hspace{0.17em}or\hspace{0.17em}}{d}_i^{\left(k+1\right)}>d_i^{\left(k\right)}(19)$

In this paper, local reinforcement is decoupled and integrated with the main algorithm in a periodic manner. After each T-generation NSDE iteration, a PSO local optimization is triggered. The process is as follows:

1. Perform non-dominated sorting on the current population and select the first frontier.

2. The first solution with the largest crowding distance in the first $N_c$ frontiers is extracted as the initial particle.

3. Perform PSO local search for K generations to obtain the optimized solution.

4. Merge the optimized solution with the current population and perform elite retention.

The selection of T value should be determined based on the convergence rate of the algorithm and the balance of computing resources. First, the HV improvement curve and single PSO cost at different T (such as 5, 10, 20, 50) are measured through preliminary experiments to find the “inflection point”. That is, the position where the gain margin decreases rapidly when T is reduced, and then T is set near the inflection point (if the performance improvement is obvious and the increase in calculation is acceptable when T = 20, then T ≈ 20 is taken) to ensure the effect of periodic local reinforcement and avoid unnecessary computing burden caused by excessively frequent PSO.

In view of the application of PSO algorithm in multi-objective optimization, this paper introduces crowding distance as an important indicator to guide particle search to enhance the diversity and distribution uniformity of solution set. Crowding distance reflects the sparsity of the solution in the target space. Solutions with larger crowding distance are located in sparser areas and have higher diversity value. Based on this, in the process of PSO speed update, the algorithm preferentially guides particles to move to areas with larger crowding distance to avoid excessive aggregation of particles in local dense areas and promote extensive exploration of solution space. This mechanism adjusts the weight of speed update to make particles more inclined to explore sparse and potential areas in solution space, thereby improving the coverage and uniformity of Pareto frontier solution set.

In the particle selection and position update link, crowding distance is also used to determine whether the newly generated position is better than the current solution. After each particle position update, the algorithm recalculates the non-dominated sorting of the solution and its corresponding crowding distance. If the crowding distance of the solution corresponding to the new position is better than or equal to the current solution, the new position is retained, otherwise the particle falls back to the old position. This selection strategy ensures that the diversity is not destroyed during the particle search process, while strengthening the development of sparse areas. By introducing the crowding distance, PSO effectively takes into account both local accuracy and overall diversity in a multi-objective environment, and improves the quality and distribution balance of elite solutions in the improved NSDE algorithm.

To avoid repeated oscillations in local dense areas and ensure the stability of the convergence process, this algorithm performs a backoff operation after each PSO velocity-position update if the crowding distance corresponding to the new position fails to surpass (or equal to) the old position: the particle directly returns to the last accepted position and marks the position as “reject update” to prevent the same position from being repeatedly tried in a short period of time. In order to prevent excessive backoff from causing global search stagnation, a maximum backoff threshold R_max (3 times) is set. When a particle is rejected for R_max consecutive updates, the particle is forced to re-randomly initialize its velocity vector or fine-tune its learning factor to help it jump out of the local dense area. This backoff + threshold mechanism can not only ensure that the crowding distance screening strictly maintains the diversity of the solution set, but also avoid search pauses caused by excessive rejection of updates, and strike a balance between convergence accuracy and stability.

In order to prevent particles from oscillating back and forth in a narrow area or gathering in the local optimum too early due to too frequent PSO updates during the high-dimensional local fine-tuning stage, the algorithm also introduces two mechanisms to the PSO module: first, a linearly decreasing inertia weight and adaptive learning factor are added to the speed update. When local area oscillation is detected (such as the position change of particles for several consecutive generations is less than ε = 10⁻⁴), the individual learning factor c₁ will be instantly reduced and the group learning factor c₂ will be increased, weakening the dependence on its own historical optimal point and strengthening the attraction to the global diversity area; second, the speed vector is limited to the upper and lower limits, and a small amount of random perturbation is applied to the over-aggregated particle group after each K generations of PSO is completed, thereby injecting a new search direction, which not only ensures the fine adjustment of the elite solution, but also avoids falling into microscopic oscillations or premature convergence.

3.4 Algorithm parameter setting and tuning strategy

The improved NSDE algorithm and PSO local enhancement module in this paper are implemented based on the MATLAB R2021b environment. The core part of the algorithm is based on the optimization toolbox that comes with MATLAB, and the multi-objective differential evolution (NSDE) and PSO local search modules are built by self-compiled functions to realize non-dominated sorting, crowding calculation, adaptive parameter update and other functions. The Pareto boundary processing and performance index (such as crowding distance, deviation) calculation modules of multi-objective optimization refer to the open source library PlatEMO (Platform for Evolutionary Multi-Objective Optimization), on which secondary development and function expansion are carried out to meet the special needs of multi-objective optimization problems of medical heavy ion accelerator beams. The interface processing part of FLUKA and OPERA simulation data is completed by mixed calls of MATLAB and Python scripts to ensure the efficiency and accuracy of objective function evaluation during the optimization process. All experiments were completed on a workstation equipped with an Intel i7-12700 CPU and 32 GB of memory.

In the PSO algorithm, the inertia weight is set as a nonlinear descent function to control the balance between convergence speed and disturbance intensity. The formula for the inertia weight is shown in Equation 20:

${\mathrm{\iota }}_{\mathrm{t}}={\mathrm{\iota }}_{\max }-\left({\mathrm{\iota }}_{\max }-{\mathrm{\iota }}_{\min }\right)·{\left(\frac{\mathrm{t}}{\mathrm{T}}\right)}^{\mathrm{o}}(20)$

$\mathrm{o}$ is the control coefficient, and its value is 1.5.

The parameter settings in this paper are formulated in coordination with the problem characteristics and the algorithm structure, taking into account both global exploration and local convergence capabilities. The differential mutation factor and crossover probability adopt a linear adaptive mechanism to achieve high-disturbance exploration in the early stage and convergence in the later stage. The reverse learning initialization uses symmetric mapping to enhance the initial population diversity and alleviate the local optimal trap. The local strengthening part of PSO improves the local accuracy of non-dominated solutions through nonlinear reduction of inertia weight and elite individual guidance. The number of particles, local step size and trigger frequency are tuned based on the balance between computing resources and optimization gain. Each parameter is screened through multiple rounds of experimental comparison and multi-objective deviation measurement function to ensure that the final configuration has consistent optimization performance in terms of performance indicator convergence, diversity and stability. The parameter settings are shown in Table 2.

Table 2

Table 2. Parameter settings.

In this paper, the selection of the number of particles Nc selected for local reinforcement should take into account both the computational resource limitations and the optimization effect requirements. Nc should maintain a certain ratio with the population size N to ensure that local reinforcement can cover enough elite solutions while avoiding excessive computational overhead. In this paper, Nc is set to the interval of 10%–20% of the population size, and is dynamically adjusted according to the complexity of the problem and the width of the Pareto front: when the problem scale is large or the Pareto front is wide, Nc is appropriately increased to ensure the coverage and diversity of local search; on the contrary, for smaller problems or narrower fronts, a smaller Nc can meet the needs. In addition, combined with the actual number and distribution density of non-dominated front solutions during the optimization process, an adaptive adjustment strategy is used to dynamically determine Nc, so that it can more flexibly adapt to different problems and algorithm iteration stages, and improve the efficiency and effect of the local reinforcement module.

4 Experimental design for optimizing beam transmission efficiency of medical heavy ion accelerator

4.1 Experimental data

The experimental data in this paper are derived from the Monte Carlo particle transmission simulation results based on the FLUKA platform. The experiment uses a typical medical heavy ion accelerator transmission pipeline model, which includes 5 quadrupole magnet focusing units and two deflection magnet units, with a total pipeline length of 12 m. The static magnetic field distribution is calculated using OPERA 4.2 and imported into FLUKA for coupled simulation. In the experiment, the incident heavy ions are $12{\mathrm{C}}^{6+}$, the energy range is 200–430MeV/u, the number of particles in a single simulation is 1*10⁵, the statistical error is less than 1%, and the initial distribution of divergence angle, beam intensity, etc., is uniformly sampled within the range of decision variables. In data collection, a simulation is run for each decision vector configuration to record the terminal cross-section particle coordinates and kinetic energy. A total of 5,218 sets of data are collected, including the number of terminal particles, particle plane coordinates, particle kinetic energy, energy loss in the pipeline, particle loss position, etc.

For the 5,218 groups of initial data sets based on uniform sampling, this paper mainly uses this data set for subsequent preprocessing steps such as normalization, correlation analysis, and principal component analysis (PCA) as the initial sample basis for the optimization algorithm. In the actual multi-objective optimization process, as the population is iteratively updated, the new decision vector needs to recalculate the corresponding objective function value through Monte Carlo simulation to ensure the dynamic accuracy of the objective evaluation. The initial data set is only used for statistical analysis such as auxiliary scaling and feature extraction, and is not directly used for the target calculation of optimization iterations, to ensure the real-time and accuracy of the objective function during the optimization process.

The schematic diagram of the transmission pipeline layout of the medical heavy ion accelerator is shown in Figure 4.

Figure 4

Figure 4. Schematic diagram of the transmission pipeline layout of the medical heavy ion accelerator.

Figure 4 shows that five quadrupole magnets (blue rectangles, marked as Q1 to Q5) and two deflection magnets (red rectangles, marked as B1 and B2) are arranged in sequence along the 14-m-long pipeline. The quadrupole magnets are mainly used for beam spot focusing adjustment, and the deflection magnets are used to control the particle trajectory deflection. The magnet units are unevenly distributed as a whole, reflecting the layout relationship of the magnetic components in the actual pipeline, which helps to understand the spatial distribution of key physical devices in the beam transmission path and their synergy.

4.2 Data preprocessing

4.2.1 Missing value and outlier processing

For missing values in the original data caused by simulation termination or complete particle loss, this paper uses the nearest neighbor interpolation method to estimate and fill in the missing values accordingly. For outliers, the experiment uses the interquartile range to detect outliers. Data points that exceed the upper and lower bounds are marked and removed to ensure that the distribution of input variables and objective function values meets the optimization requirements. The range formula for outlier detection is shown in Equation 21.

$\left[\mathrm{Q}1-1.5*\left(\mathrm{Q}3-\mathrm{Q}1\right),\mathrm{Q}3+1.5*\left(\mathrm{Q}3-\mathrm{Q}1\right)\right](21)$

$\mathrm{Q}1$ represents the first quartile and $\mathrm{Q}3$ represents the third quartile.

In terms of missing value processing, this article further distinguishes the sources of missing values. For missing values caused by simulation termination (such as numerical divergence or program interruption), considering that the corresponding configuration may not be physically feasible, it is directly eliminated and not included in the subsequent statistical analysis. For missing values caused by complete particle loss (particles fail to reach the terminal cross section), this situation reflects poor transmission performance, and the corresponding beam transmission efficiency objective function value is recorded as zero. Other objective functions (such as terminal kinetic energy and spot size) are no longer counted to prompt the optimization algorithm to avoid such invalid solution configurations. This processing method ensures the consistency of objective function evaluation and avoids interference in the optimization search process due to abnormal data.

For missing value imputation, the K-nearest neighbor interpolation method is used, where K is set to 5. This means the missing value is filled with the weighted average of the 5 nearest valid samples to maintain local data distribution characteristics. For outlier detection, the interquartile range (IQR) method is used, with upper and lower limits of Q1−1.5 × IQR and Q3 + 1.5 × IQR, respectively, where IQR = Q3 − Q1, and Q1 and Q3 are the first and third quartiles, respectively. Samples outside this range are marked as outliers and removed to ensure a reasonable distribution of input variables and objective function values, and to avoid outlier data interfering with the optimization process.

4.2.2 Normalization and feature decorrelation processing

After processing missing values and outliers, this paper uses Min-Max normalization to linearly normalize each decision variable and target indicator and map them to [0,1], as shown in Formula 22.

${\widehat{x}}_{ij}=\frac{x_{ij}-{\mathit{min}}_j\left(x.j\right)}{{\mathit{max}}_j\left(x.j\right)-{\mathit{min}}_j\left(x.j\right)}(22)$

Among them, ${\widehat{\mathrm{x}}}_{\text{ij}}$ represents the normalized value, ${\mathrm{x}}_{\text{ij}}$ represents the original decision variable, ${\min }_{\mathrm{j}}\left(\mathrm{x}.\mathrm{j}\right)$ and ${\max }_{\mathrm{j}}\left(\mathrm{x}.\mathrm{j}\right)$ represent the minimum and maximum values of the decision variable respectively. The normalization processing of the target indicator is the same as that of the decision variable.

After the above processing, the principal component analysis is now performed on the decision variable, the principal components with a cumulative contribution rate of less than 1% are eliminated, and the variables are reconstructed.

4.3 Evaluation indicators

In this paper, in order to comprehensively evaluate the performance of the improved multi-objective optimization algorithm NSDE, the following six common indicators are selected: GD (generational distance), C-metric, IGD (inverted generational distance), SPV (Spacing Variance) of the spacing between adjacent solutions in the solution set, SI (Spread Indicator) and HV (hypervolume).

The formula of GD is shown in Equation 23.

$\text{GD}\left(\mathrm{P},{\mathrm{P}}^{*}\right)=\frac{1}{\left|\mathrm{P}\right|}\sqrt{{\displaystyle \sum _{\mathrm{i}=1}^{\left|\mathrm{P}\right|}}\left({}_{\mathrm{y}\in {\mathrm{P}}^{*}}{}^{\min }‖{\mathrm{F}}_{\mathrm{i}}-\mathrm{y}‖^2\right)}(23)$

$\mathrm{P}$ represents the non-dominated solution set obtained by the algorithm, ${\mathrm{P}}^{*}$ represents the reference Pareto optimal solution set. ${\mathrm{F}}_{\mathrm{i}}$ represents the target vector, and $\mathrm{y}$ represents any point in the reference optimal solution. “Reference optimal solution” refers to an optimal solution on the reference Pareto frontier, which is used to evaluate the closeness of the algorithm solution set to the ideal optimal solution, such as distance calculation in GD or IGD; “reference Pareto optimal solution set” refers to a complete or representative Pareto optimal solution set that is known or approximated by a high-precision algorithm, which is used as a benchmark for performance evaluation. The relationship between the two is: the reference optimal solution is a single solution point in the reference Pareto optimal solution set, and both are “ideal solutions” used for comparison and measurement in the evaluation process.

The formula of C-metric is shown in Equation 24.

$\mathrm{C}\left({\mathrm{\rho }}^{*},{{\mathrm{\omega }}^{-}}^{*}\right)=\frac{\left|\left\{{\mathrm{\omega }}^{-}\in {{\mathrm{\omega }}^{-}}^{*}\left|\exists \mathrm{\rho }\in \right.{\mathrm{\rho }}^{*}:\mathrm{\rho }\prec {\mathrm{\omega }}^{-}\right\}\right|}{\left|{{\mathrm{\omega }}^{-}}^{*}\right|}(24)$

${\mathrm{\rho }}^{*}$ represents a non-dominated solution set, and ${{\mathrm{\omega }}^{-}}^{*}$ represents another non-dominated solution set.

The formula of IGD is shown in Equation 25.

$\text{IGD}\left(\mathrm{P},{\mathrm{P}}^{*}\right)=\frac{1}{\left|{\mathrm{P}}^{*}\right|}{\displaystyle \sum _{\mathrm{j}=1}^{\left|{\mathrm{P}}^{*}\right|}}{}_{\mathrm{F}\in \mathrm{P}}{}^{\min }‖\mathrm{F}-{\mathrm{y}}_{\mathrm{j}}‖(25)$

The formula for SPV is shown in Equation 26.

$\text{SPV}=\sqrt{\frac{1}{\mathrm{N}-1}{\displaystyle \sum _{\mathrm{i}=1}^{\mathrm{N}}}{\left({\text{dis}}_{\mathrm{i}}-\overline{\text{dis}}\right)}^2}(26)$

${\text{dis}}_{\mathrm{i}}$ represents the distance between the solution and the nearest neighbor, and $\overline{\text{dis}}$ represents the mean of all minimum distances.

The formula of SI is shown in Equation 27.

$\text{SI}=\frac{{\displaystyle {\sum }_{\mathrm{i}=1}^{\mathrm{N}-1}}\left({\text{dis}}_{\mathrm{i}}-\overline{\text{dis}}\right)+{\text{dis}}_{\min }+{\text{dis}}_{\max }}{{\text{dis}}_{\min }+{\text{dis}}_{\max }+\left(\mathrm{N}-1\right)\overline{\text{dis}}}(27)$

${\text{dis}}_{\min }$ represents the distance between the first boundary solution and the minimum endpoint of the reference solution set, and ${\text{dis}}_{\max }$ represents the distance between the last boundary solution and the maximum endpoint of the reference solution set.

The formula of HV is shown in Equation 28.

$\text{HV}\left(\mathrm{P}\right)=\text{vol}\left({}_{\mathrm{F}\in \mathrm{P}}{}^{\mathrm{U}}\left[{\mathrm{f}}_1,{\mathrm{\varrho }}_1\right]*\left[{\mathrm{f}}_2,{\mathrm{\varrho }}_2\right]*\dots *\left[{\mathrm{f}}_{\mathrm{m}},{\mathrm{\varrho }}_{\mathrm{m}}\right]\right)(28)$

Among them, ${\mathrm{\varrho }}_{\mathrm{m}}$ represents the reference point, corresponding to the inferior solution boundary in the target space, and vol (.) represents the volume of the set in the target space. The “reference point” refers to a fixed point in the target space used to determine the volume boundary when calculating the hypervolume (HV). It is usually selected at a worse value of each objective function to measure the target space volume range covered by the non-dominated solution set.

4.4 Experimental design

This study employs an improved NSDE algorithm, incorporating inverse learning initialization, an adaptive mechanism, and a PSO local enhancement module, as the primary optimization framework. To comprehensively validate its performance, this paper compares the algorithm with seven control algorithms, including five mainstream multi-objective optimization algorithms: NSDE, NSGA-II, MOEA/D, MOPSO, and MOCOA (multi-objective Coati optimization algorithm), and two algorithm variants used for ablation experiments: NSDE using only inverse learning initialization, and NSDE using both inverse learning initialization and an adaptive mechanism. All algorithms are run on the same nine-dimensional decision variable space, the same four optimization objectives (beam transmission efficiency, beam spot size, system power consumption, and energy retention rate), and a consistent preprocessed dataset. Each experiment is independently repeated 30 times to reduce the impact of randomness; each optimization process terminates after 200 iterations or when the convergence threshold is reached, and the resulting optimal Pareto front solution set is recorded for subsequent performance analysis.

To ensure fairness in comparison, all algorithms are executed in the same hardware environment, a 16-core CPU (Central Processing Unit) parallel computing node, and the same FLUKA simulation interpolation model is called to obtain the objective function value. The experimental indicators include HV, IGD, GD, SPV, SI and C-metric. The experiment compares the mean and standard deviation of each algorithm on the above indicators, and analyzes the contribution of reverse learning, adaptive mechanism, and PSO enhancement in the improved module to the algorithm convergence speed, solution set diversity, and single-objective optimization performance. This can fully verify the superiority and robustness of the improved NSDE in the optimization of beam transmission efficiency of medical heavy ion accelerators.

To ensure the fairness and reproducibility of the experimental results, this paper uniformly sets and clearly records the main hyperparameters of the comparison algorithms. The NSGA-II algorithm uses simulated binary crossover (SBX) and polynomial mutation operators, with the crossover probability set to 0.9, the mutation probability set to 1/decision variable dimension 1/9 ≈ 0.111, the crossover distribution index set to 20, and the mutation distribution index set to 20. The number of neighbors of the MOEA/D algorithm is 20, the crossover probability is 0.9, the mutation probability is 1/9, and the weight vector is generated with uniform distribution. In the MOPSO algorithm, the number of particles is consistent with the population size (120), the initial value of the inertia weight is 0.9, the final value is 0.4, and the learning factors c1 and c2 are both 1.8. The MOCOA algorithm uses standard parameter configuration with a population size of 120. The maximum number of iterations for all baseline algorithms is set to 200, and the population size is consistent (120).

To enhance the reproducibility of the comparative experiments and the transparency of the baseline definition, the source and function of the “original medical heavy ion accelerator design” (“original values”) referred to in this paper are as follows. This original parameter set is not derived from a publicly operating particle therapy facility, but rather is a representative initial engineering configuration determined by our research group during the accelerator transmission line engineering design and experimental phases through traditional “empirical-simulation” manual parameter tuning (i.e., the engineering starting point before systematic multi-objective global optimization). These baseline parameters serve as inputs to the paper-based Monte Carlo coupled simulations (FLUKA and OPERA) and initial uniform sampling (1 × 10^5 particles/sample) verification process, reflecting the engineering performance level before full optimization by our intelligent algorithm (evaluated under the same simulation platform, the same random seed sequence, and the same computing environment: Intel i7-12700, 32 GB memory). In this paper, this serves as a unified and reproducible benchmark for comparison, ensuring that the conclusions are a fair measure relative to the original engineering starting point under the same simulation model and experimental conditions.

5 Optimization results display

5.1 Pareto frontier quality display diagram of different algorithms

In order to verify the optimization performance of different algorithms, the Pareto frontier distribution in the target space obtained by comparing the improved NSDE algorithm with the NSDE algorithm, NSGA-II, MOEA/D, MOPSO, and MOCOA is used to see whether it is closer to the ideal point, and the Pareto frontier quality of different algorithms is visualized. The results are shown in Figure 5. In Figure 5, beam transmission efficiency, beam spot size, and system power consumption are respectively represented as three-dimensional coordinates, and energy retention is represented by color mapping. The yellow pentagrams in each subgraph represent ideal points and serve as a benchmark for comparing the approximation and diversity of solution sets among different algorithms.

Figure 5

Figure 5. Pareto front quality display of different algorithms. (a) Results of the improved NSDE algorithm; (b) Results of the NSDE algorithm. (c) Results of the NSGA-II algorithm; (d) Results of the MOEA/D algorithm. (e) Results of the MOPSO algorithm; (f) Results of the MOCOA algorithm.

In Figure 5, the scattered point cloud of the improved NSDE is more closely clustered towards the ideal point. The color distribution is biased towards the high energy retention rate area, indicating that it has a better overall performance when taking into account the four objectives, while the solution sets of other algorithms deviate from the efficiency dimension and perform unevenly in the power consumption-beam spot size balance.

The improved NSDE can obtain a frontier closer to the ideal point, mainly due to the synergy of the three improved strategies. In the initial population stage, the reverse learning initialization collects more widely distributed candidate solutions, allowing the algorithm to start from a better diversity basis and avoid falling into the local optimum in the early stage. The adaptive mechanism dynamically adjusts the mutation factor and crossover probability according to the fitness of the population, so that the entire target space can be fully explored in the early stage, while the later stage focuses on fine convergence and achieves a more accurate compromise in the multi-target conflict area. PSO local reinforcement fine-tunes the speed-position of the non-dominated elite solution, further improving the approximation of the Pareto frontier in terms of efficiency and energy retention. The three complement each other and jointly construct an optimization process that can both globally search and locally refine, and finally make the solution set of the improved NSDE show the minimization of f2 and f3 and the maximization of f1 and f4 in the scatter plot.

From the perspective of multi-objective optimization, there is a typical triangular conflict between beam transmission efficiency, beam spot size, and system power consumption. Increasing the focusing intensity can improve efficiency but increase power consumption and beam spot divergence, and vice versa. During the search process, the improved NSDE continuously screens the Pareto hierarchy through non-dominated sorting, and combines crowding distance to maintain the diversity of the solution set, so that the algorithm neither loses the exploration of extreme compromise solutions nor falls into single path convergence. This mechanism complements the global jump ability of differential evolution based on vector difference, allowing the solution set to cross the local valley and quickly converge to the multi-objective optimal area. At the same time, the gradient adjustment in the subtle area during the PSO enhancement phase helps the algorithm to eliminate the rough boundary solutions, achieving a good balance between stability and diversity. NSGA-II and MOEA/D are slightly insufficient in terms of adaptability and local accuracy, while MOPSO lacks systematic screening of non-dominated sorting, resulting in its frontier being overall backward or unevenly distributed in the graph. Through the integration of the above technologies, the improved NSDE in this paper maintains efficient global exploration while carefully cultivating the details of the Pareto frontier, presenting the best frontier quality.

5.2 Convergence and solution set diversity evaluation of multi-objective algorithms

In the experiment, for multi-objective algorithms, the measurement indicators are divided into three aspects, including convergence, solution set diversity, and comprehensive indicators. This paper uses convergence index GD, C-metric, comprehensive index IGD, solution set diversity index SPV, distribution width SI, HV to evaluate the multi-objective algorithm, and the results are shown in Figure 6. In Figure 6, C-metric shows the comparison results between the improved NSDE algorithm and other algorithms.

Figure 6

Figure 6. Convergence and solution diversity evaluation of multi-objective algorithms. (a) GD evaluation results. (b) C-metric evaluation results. (c) IGD evaluation results. (d) SPV evaluation results. (e) SI evaluation results. (f) HV evaluation results.

In Figure 6, the improved NSDE achieves the best results in both GD and IGD, with GD being 0.028 and NSDE being only 0.035. The IGD of the improved NSDE is 0.035, while that of NSDE is 0.045, indicating that the solution set of the improved NSDE is closer to the reference Pareto frontier. NSGA-II achieves 0.042 and 0.052 on GD and IGD, respectively, while MOCOA performs the worst, reaching 0.052 and 0.062, respectively. In the C-metric comparison, the dominance rate of the improved NSDE over NSDE reaches 0.82, and the dominance rates over NSGA-II, MOEA/D, MOPSO, and MOCOA are 0.75, 0.78, 0.70, and 0.65, respectively. It can be seen that the frontier of the improved NSDE is superior to other methods in most objectives. Overall, the improved NSDE is significantly better than the control algorithms in terms of convergence speed and frontier quality.

In terms of solution set diversity, the SPV of the improved NSDE is 0.009, the smallest among all algorithms, and NSDE reaches 0.012. The SI of the improved NSDE is 0.10, and the MOEA/D is 0.16. On HV, the improved NSDE reaches 0.48, while NSDE is only 0.42. The data shows that the improved NSDE can cover the target space in a more uniform and wider way while maintaining a high-quality solution set, and obtains comprehensive advantages in the diversity and coverage of multi-objective optimization.

The performance improvement of the improved NSDE in this paper is due to the fact that the reverse learning initialization constructs a wide coverage of candidate solutions in the initial population stage, preventing the algorithm from concentrating on the local area too early, and avoiding the problems of high GD and IGD in the first few generations. The adaptive mechanism dynamically adjusts the mutation factor and crossover probability, balancing exploration and development according to the population convergence, so that the algorithm can continue to advance to the Pareto frontier in the middle and late stages, further reducing GD/IGD. PSO local enhancement fine-tunes the speed-position of the elite solution, improves local convergence accuracy, and promotes the continuous growth of HV.

Multi-objective optimization is essentially a compromise between conflicting objectives. Although NSDE, NSGA-II, and MOEA/D have certain mechanisms in non-dominated sorting and crowding maintenance, they lack dynamic adaptation of search strategies. When facing highly coupled four-dimensional objectives, the solution set is prone to “clustering” or “sparseness”, resulting in high SPV and SI and low HV. MOPSO and MOCOA also have difficulty maintaining a good solution set distribution due to the lack of strong non-dominated screening or local fine-tuning. The improved NSDE relies on adaptation and PSO reinforcement to continuously supplement the missing areas on the frontier, while balancing the global and local through the crowding distance. This ensures exploration in all directions and maintains uniform coverage of the border and middle areas, leading in the three diversity indicators of SPV, SI and HV.

5.3 Convergence curve analysis

The experiment was iterated 200 times, and the convergence performance of different algorithms was statistically analyzed. The results of the convergence curve analysis are shown in Figure 7. In Figure 7, for the convenience of representation, Z1, Z2, and Z3 represent reverse learning initialization, adaptive mechanism, and PSO local reinforcement, respectively. In Figure 6, due to the multi-objective nature, the experiment uses the HV indicator for statistics.

Figure 7

Figure 7. Convergence curve analysis results.

The convergence curve in Figure 7 shows that after 140 iterations, the HVs of all algorithms have basically converged. The improved NSDE (Z1+Z2+Z3) has stabilized at 0.48 after 90 iterations, the improved NSDE algorithm (Z1+Z2) has reached convergence after 100 iterations, and the improved NSDE algorithm (Z1) has reached convergence after 110 iterations. NSDE algorithm, NSGA-II, MOEA/D, MOPSO, and MOCOA all reach their respective stable platforms after 120–140 generations, and their values are generally no higher than 0.42.

Rapid improvement of HV requires that the solution set quickly approaches the Pareto frontier and maintains good diversity. Although NSDE and NSGA-II have non-dominated sorting and crowding maintenance mechanisms, they lack dynamic adjustment of search parameters and local refinement, resulting in a slowdown in HV improvement between 60 and 80 generations. MOEA/D and MOPSO converge to an average HV of 0.36–0.4 because decomposition or particle swarm strategies are difficult to take into account both global and local considerations in a high-dimensional four-objective space. The improved NSDE continuously controls mutation and crossover through an adaptive mechanism, maintaining a high upward slope in the mid-term 50–80 generations. The PSO enhancement performs local refinement on the non-dominated elite solution in the late 80–100 generations, eliminating the holes on the boundary and quickly pushing HV to 0.48. The optimization dynamics of “wide-area exploration first, local refinement later” in this paper enables the improved NSDE to surpass similar algorithms in multi-objective conflict management.

5.4 Beam transmission efficiency, end beam spot size, system power consumption, and energy retention rate under different algorithm optimizations

The simulation results of the beam transmission process of the original medical heavy ion accelerator without optimized design are shown in Figure 8.

Figure 8

Figure 8. Simulation results of the beam transmission process of the original medical heavy ion accelerator without optimized design.

In Figure 8, the beam spot radius shown in the figure adopts the root mean square radius (RMS, denoted as σ_x) in the transverse x direction, and its calculation definition is shown in Equation 29.

${\sigma }_x=\sqrt{<{\left(x-x_c\right)}^2>}(29)$

In Formula 29, x is the x-coordinate of a single particle on the terminal cross section (unit: mm), $x_c$ =<x> is the center of mass position, and the angle sign represents the arithmetic mean of all particles reaching the terminal. This statistic reflects the second-order broadening of the beam in the transverse x-direction and is commonly used to evaluate beam optics and transmission efficiency.

Figure 8 shows the simulation results of the beam transmission process of the original medical heavy ion accelerator without optimized design. The left axis in the figure shows the change of the beam spot radius along the transmission line direction (0–14 m), and the right axis shows the beam transmission efficiency. It can be observed that the beam spot radius shows a certain fluctuation during the transmission process, reflecting the phenomenon of alternating focusing and divergence in the optical system, indicating that the magneto-optical parameters have not been stably matched. At the same time, the beam transmission efficiency gradually decreases with the increase of transmission distance, and finally reaches about 92.42% at the end. It shows that there are particle losses caused by factors such as beam divergence, stray field interference and beam limiting structure in the original design. This simulation verifies the physical basis for the insufficient beam transmission efficiency of the original system, and also supports the necessity of optimized design from the perspective of particle dynamics.

The original medical heavy ion accelerator has a beam transmission efficiency of only 92.42%, not due to the imbalance of a single physical factor, but because the traditional “experience-simulation” iteration method of magneto-optical parameters is difficult to meet the global requirements of multiple objectives (high efficiency, tight beam spot, low power consumption, high energy retention) at the same time. Although the quadrupole magnet gradient, deflection magnet current, pipeline geometry, vacuum parameters, etc., are optimized for multiple objectives at the beginning of the design, this method of manual parameter adjustment + local search has a long iteration cycle, is sensitive to high-dimensional coupling, and is easy to fall into local optimality. If you are not careful, side effects such as beam divergence, stray field enhancement, or mechanical beam limitation will occur, and ultimately it can only compromise at a level of just over 92%. It is precisely because the traditional optimization process is time-consuming and difficult to ensure global optimality that the value of high-dimensional multi-objective global optimization methods based on intelligent algorithms is highlighted. Within limited computing resources and time budget, it can efficiently jump out of local traps, achieve coordinated balance of all physical objectives, and improve beam transmission efficiency.

The results of beam transmission efficiency, end beam spot size, system power consumption, and energy retention rate under different algorithm optimizations are shown in Figure 9. Each objective value shown in Figure 9 represents the statistical performance (mean ± standard deviation) of the optimal or near-optimal value of the single objective indicator in the Pareto frontier solution set obtained by the corresponding algorithm in the multi-objective optimization process, rather than a single extreme optimal solution. These indicators reflect the representative performance of a set of Pareto optimal solutions obtained by different algorithms under multi-objective trade-offs, reflecting the overall optimization effect and stability of the algorithm in the objective space, rather than the extreme value of a single solution point.

Figure 9

Figure 9. Beam transmission efficiency, end beam spot size, system power consumption, energy retention rate under different algorithm optimizations. (a) Beam transmission efficiency under different algorithm optimizations. (b) End beam spot size under different algorithm optimizations. (c) System power consumption under different algorithm optimizations. (d) Energy retention rate under different algorithm optimizations.

The error of ±0.49% is based on the standard deviation of the beam transmission efficiency obtained from multiple independent experiments. Specifically, the statistically improved NSDE algorithm converges to the optimal solution set in different random initial populations and multiple optimization processes, and calculates the mean and standard deviation of the transmission efficiency. This result reflects the stability and robustness of the algorithm under multiple optimizations, and reflects the consistency of the improved NSDE in optimizing the transmission efficiency under different random conditions. In Figure 9, it can be seen that the beam transmission efficiency of all algorithms is higher than that of the original medical heavy ion accelerator. The beam transmission efficiency of the improved NSDE algorithm is 99.21% ± 0.49%, which is significantly higher than all the comparison algorithms and 6.79% higher than the original medical heavy ion accelerator. NSDE is only 97.47% ± 0.63%, and MOCOA performs the worst, reaching 92.89% ± 0.85%. In terms of the end beam spot size, the improved NSDE achieved a minimum value of 5.13 mm ± 0.10 mm, NSDE reached 5.68 mm ± 0.12 mm, and NSGA-II reached 5.74 mm ± 0.15 mm. This shows that the improved NSDE can maintain a more compact beam spot focus while achieving higher transmission efficiency.

In terms of system power consumption, the improved NSDE reached 153.68 kW ± 1.42 kW, while NSDE reached 168.23 kW ± 1.85 kW. MOCOA reached 174.76 kW ± 2.04 kW, indicating that the improved NSDE algorithm significantly reduced power consumption while optimizing beam quality. The improved NSDE achieved an energy retention rate of 94.12% ± 0.39%, while NSDE reached 90.46% ± 0.47%. MOCOA is only 87.97% ± 0.65%, indicating that the improved NSDE is effective in reducing particle dissipation loss.

Among the comparison algorithms, although NSGA-II and MOEA/D have multi-objective decomposition and non-dominated sorting mechanisms, they lack strategies for initial diversity and parameter adaptation, and are prone to fall into local areas in high-dimensional coupling problems, resulting in their efficiency and energy retention rate not being further improved. Although MOCOA integrates multiple operators, its fixed parameters make it incapable of compromising multiple target conflicts, and its efficiency and energy retention are both at the lowest values. The improved NSDE uses a triple strategy to achieve targeted compensation for the weaknesses of each algorithm, and has achieved comprehensive breakthroughs in the four indicators of efficiency, focusing accuracy, energy consumption, and energy retention.

To further explore the differences between the different algorithms, a significance test was conducted (two-group comparison: improved NSDE vs. each control algorithm), reporting the t-value, degrees of freedom (df), p-value, effect size Cohen’s d, and 95% confidence intervals for the differences. The test data were based on the sample mean and sample standard deviation of 30 independent runs for each algorithm; to account for potential unequal variances, a Welch t-test (two-tailed) was used, and the Welch degrees of freedom and 95% confidence intervals based on this df were reported; Cohen’s d was also calculated (using the pooled SD estimate of the two groups’ sample variance means) to measure the effect size; to control for multiple comparisons, Bonferroni correction was used (5 comparisons, with a corrected significance threshold α = 0.05/5 = 0.01). The null hypothesis H0: The population mean of improved NSDE and a control algorithm is equal in beam transmission efficiency; the alternative hypothesis H1: The population means are unequal (two-tailed test). If p < 0.01 (Bonferroni correction), the difference is considered statistically significant. The results of the significance test are shown in Table 3.

Table 3

Table 3. Significance test results.

As shown in Table 3, the improved NSDE significantly outperforms all control algorithms in beam transmission efficiency. Compared with the original NSDE, the t-value is 11.94, the degrees of freedom are 54.69, p < 0.001, and Cohen’s d = 3.08, indicating a very strong effect size. Compared with NSGA-II, t = 18.05, df = 51.12, and Cohen’s d = 4.66; compared with MOEA/D, t = 23.01, df = 48.80, and Cohen’s d = 5.94; compared with MOPSO, t = 26.80, df = 48.07, and Cohen’s d = 6.92; compared with MOCOA, t = 35.28, df = 46.36, and Cohen’s d = 9.11. All comparisons were significant at the Bonferroni-corrected significance threshold of α = 0.01 (p < 0.001), and Cohen’s d values were all much greater than 0.8, indicating that the differences were not only statistically significant but also significant in terms of practical optimization effects. This demonstrates that the improved NSDE algorithm is more stable and effective than the traditional multi-objective algorithm and the original NSDE in improving beam transmission efficiency.

5.5 Beam transmission efficiency under different decision variable dimensions and noise

To investigate the impact of different decision variable dimensions and noise on beam transmission efficiency, experiments with 9–15 different dimensions and 10–30 dB Gaussian white noise were designed for verification. The results are shown in Figure 10. In Figure 10, the intensity of the Gaussian white noise is quantified using the signal-to-noise ratio (SNR). The dimensions include: 1) quadrupole magnet gradient; 2) deflection magnet current; 3) incident beam transverse emittance; 4) divergence angle; 5) relative energy dispersion; 6) magnet spacing; 7) adjustable beam limiter aperture; 8) deflection angle; 9) pipeline pressure; 10) magnet correction current; 11) vacuum parameters; 12) beam energy stability; 13) magnetic field nonlinearity correction coefficient; 14) beam eccentricity; and 15) beam spot shape parameters.

Figure 10

Figure 10. Beam transmission efficiency under different decision variable dimensions and noise. (a) Beam transmission efficiency under different decision variable dimensions. (b) Beam transmission efficiency under different noise intensities.

In Figure 10a, as the dimension of decision variables increases from 9 to 15, the transmission efficiency of all algorithms shows a downward trend. Taking NSDE as an example, the efficiency drops from 97.47% to 91.26%, and NSGA-II drops from 96.16% to 90.33%. MOEA/D, MOPSO, and MOCOA drop from 95.34% to 89.28%, 94.62%–88.61%, and 92.89%–87.21%, respectively. After the improvement, NSDE drops from 99.21% to 93.88%, but the drop is small. The increase in dimension increases the search complexity, resulting in a decrease in the performance of all algorithms, but the improved NSDE is more robust to this.

In Figure 10b, when the Gaussian white noise signal-to-noise ratio increases from 10 dB to 30 dB, the transmission efficiency increases with the increase of the signal-to-noise ratio. NSDE increased from 90.14% to 98.44%, and NSGA-II increased from 88.46% to 97.13%. MOEA/D, MOPSO, and MOCOA increased from 87.39% to 96.31%, 85.61%–95.59%, and 84.18%–93.86%, respectively. The improved NSDE increased from 92.33% to 100.00%. All algorithms benefit from a higher signal-to-noise ratio, but the improved NSDE always maintains the highest efficiency in the low noise range of 20–30 dB, and the efficiency improvement curve is smoother, indicating that it is less sensitive to noise.

The increase in variable dimensions brings a larger search space and stronger parameter coupling. Traditional algorithms such as NSGA-II and MOEA/D rely on fixed non-dominated sorting and crossover mutation, lack dynamic search strategies, and are prone to falling into local areas in high-dimensional space, resulting in a more obvious decline in efficiency. The improved NSDE combines reverse learning initialization and adaptive mechanisms to ensure that it can still effectively jump out of the local optimum under high-dimensional conditions, with a small drop.

Noise reduces the accuracy of objective function evaluation, and algorithms with fixed parameters such as MOEA/D and MOCOA are difficult to maintain accurate search under noise. The adaptive mechanism of the improved NSDE adjusts the mutation factor and crossover probability according to the real-time fitness noise level, automatically expands the search range in a high-noise environment, and focuses on fine optimization in a low-noise environment, maintaining a higher starting efficiency and a more stable increase. This shows that this method can effectively alleviate the negative impact of noise on the optimization process and achieve adaptive robust optimization in a noisy environment.

The comparison results of the original values and optimized values of the 9-dimensional decision variables are shown in Table 4.

Table 4

Table 4. Comparison results of the original values and optimized values of the 9-dimensional decision variables.

Table 4 shows the comparison results of the original and optimized values of the 9-dimensional decision variables. It can be seen that after optimization, all parameters are within a reasonable range and the overall performance is better. The gradients of the five quadrupole magnets were increased from the original Q1: 4.8, Q2: 5.0, Q3: 5.1, Q4: 4.9, and Q5: 5.2 T/m to the range of 6.0–6.4 T/m, respectively, which better enhanced the beam focusing capability; the deflection magnet current was increased from 150 A to 172 A, which is close to the upper limit but still within a reasonable range, which helps to improve the deflection accuracy; the lateral emittance of the incident beam was reduced from 1.5 mm·mrad to 0.9 mm·mrad, and the divergence angle was also optimized from 0.50 mrad to 0.42 mrad, indicating that the beam quality was significantly improved; the relative energy dispersion σ_p/p was compressed from 2.0 × 10⁻³ to 1.0 × 10⁻³, which is more conducive to energy stability; in terms of geometric and environmental conditions, the magnet spacing was reduced to 1.3 m, the adjustable beam limiter opening was narrowed to 9 mm, and the deflection angle was reduced from 1.5 mm·mrad to 0.9 mm·mrad. The angle was adjusted from 12° to 10°, both ensuring beam transmission while optimizing system compactness. Simultaneously, the pipeline pressure was reduced from 1.0 × 10⁻⁶ Pa to 8.0 × 10⁻⁷ Pa, approaching high vacuum conditions. Overall results show that the optimized parameter configuration significantly improves focusing intensity, beam quality, and system stability compared to the original solution.

5.6 Resource consumption and real-time analysis

The results of resource consumption and real-time analysis are shown in Figure 11. In Figure 11, the indicators used are CPU, GPU (Graphics Processing Unit) utilization, real-time feedback delay time and throughput to quantitatively analyze resource consumption and real-time performance. In the simulation results above, the real-time feedback delay time in the original medical heavy ion accelerator design reached 272.5 ms.

Figure 11

Figure 11. Resource consumption and real-time analysis results. (a) CPU, GPU utilization results. (b) Real-time feedback delay time and throughput results.

In Figure 11a, it can be seen that the improved NSDE algorithm has the most adequate performance in terms of resource utilization. The CPU utilization is 81.4% and the GPU utilization is 83.6%, which are higher than 68.3% and 41.7% of NSDE, 73.5% and 58.9% of NSGA-II, 65.7% and 37.2% of MOEA/D, 61.2% and 35.4% of MOPSO, and 76.8% and 67.1% of MOCOA. Among them, the CPU/GPU utilization of MOPSO is the lowest, indicating that its parallelism and computational intensity are both weak. The GPU resource call rate of NSDE is only 41.7%, while the improved NSDE almost fully utilizes the available computing resources, laying the foundation for subsequent real-time computing and high throughput.

In Figure 11b, the improved NSDE achieves a minimum value of 132.9 ms in real-time feedback delay, which is 139.6 ms shorter than the original medical heavy ion accelerator design, while other algorithms are all above 188 ms. MOCOA is 188.5 ms, NSGA-II is 197.4 ms, MOEA/D is 214.8 ms, NSDE is 224.6 ms, and MOPSO is the worst at 258.7 ms. In terms of throughput, the improved NSDE reaches 10.36 times/s, while MOCOA is only 8.87 times/s, NSGA-II is 8.43 times/s, MOEA/D is 7.92 times/s, NSDE is 7.81 times/s, and MOPSO is 6.47 times/s. The results show that the improved NSDE can make better use of hardware resources and complete more optimization iterations per unit time, significantly improving real-time performance and processing capabilities.

The improved NSDE integrates reverse learning initialization, adaptive parameter adjustment, and PSO local reinforcement. These modules can achieve efficient acceleration based on population vectorization operations in a parallel computing environment. The reverse learning and adaptive mechanism mainly involve vector operations and simple random number generation, which can be executed in parallel on the GPU; PSO local reinforcement uses matrix-level vector updates and speed calculations, which are also easy to assign to GPU threads. NSGA-II has multiple comparisons and dynamic congestion calculations in non-dominated sorting, MOEA/D frequently performs sub-problem decomposition and sub-population exchange, and MOPSO lacks efficient non-dominated sorting, which makes it difficult to fully utilize GPU resources. The improved NSDE also reduces redundant data copying and serial control flow on the CPU side, increasing CPU utilization to 81.4%.

Real-time feedback latency is subject to the efficiency of the computation-communication-scheduling pipeline. The improved NSDE reduces the number of GPU and CPU data transfers by merging key computation steps and adopts a batch simulation interpolation interface, so that most of the computations in each iteration are completed in the GPU memory, and the results are only sent back when necessary. Algorithms such as NSDE and MOEA/D lack this batching strategy and switch between the host and the device many times, resulting in high average latency. The improved NSDE uses multi-threaded asynchronous scheduling to start the next-generation of population evaluation in parallel, and interleaves the PSO local optimization with the main loop in parallel, further improving the throughput. The end-to-end scheduling and parallel optimization of this paper give it a significant performance advantage in real-time scenarios.

5.7 Ablation experiment

In this paper, the ablation experiment aims to systematically evaluate the actual contribution of each component in the improved NSDE algorithm, reverse learning initialization, adaptive mechanism, and PSO algorithm to the optimization effect of beam transmission efficiency. The experiment keeps the experimental data set, initial conditions and evaluation indicators consistent, and summarizes the marginal contribution of each module to the final performance improvement. The ablation experiment steps are as follows:

1. The complete improved algorithm is set as the baseline control group, including reverse learning initialization, adaptive mechanism, and PSO algorithm. This group is used to obtain the optimal beam transmission efficiency as a comparison standard.

2. Remove the reverse learning initialization separately, retain the adaptive mechanism and PSO algorithm, and further quantify the marginal benefit of reverse learning initialization.

3. The PSO algorithm auxiliary elite local enhancement module is removed, and the reverse learning initialization and adaptive mechanism are retained to evaluate the impact of global optimization ability on beam efficiency without local enhancement.

4. The adaptive mechanism and PSO algorithm modules can be removed, and only the reverse initialization is retained to observe the continuous optimization ability of the initial population diversity in the later stage when there is no dynamic adjustment of parameters and local reinforcement.

5. The reverse learning initialization, adaptive mechanism, and PSO algorithm are removed at the same time, and the most basic NSDE algorithm structure is returned to measure the joint improvement effect of all improved strategies and verify the improvement effect of reverse learning initialization.

The ablation experiment results are shown in Table 5.

Table 5

Table 5. Ablation experiment results.

In Table 5, as the key modules are gradually removed, the performance of the improved NSDE algorithm in various physical indicators has significantly decreased. The beam transmission efficiency achieved by the complete algorithm is 99.21%, the end beam spot size is 5.13 mm, the system power consumption is 153.68 kW, and the energy retention rate is 94.12%. When the PSO local enhancement is removed, the efficiency drops to 94.47%, the beam spot size increases to 5.69 mm, the power consumption increases to 161.21 kW, and the retention rate drops to 90.84%. When the adaptive mechanism is further removed, the efficiency drops to 91.56%, the beam spot increases to 5.97 mm, the power consumption increases to 166.44 kW, and the retention rate drops to 89.31%. Pure NSDE performs the worst, with an efficiency of only 89.02%, a beam spot size of 6.12 mm, a power consumption of 170.28 kW, and a retention rate of 88.45%.

Removing the reverse learning initialization (retaining the adaptive mechanism and PSO) resulted in a decrease in efficiency to 96.34%, an increase in the beam spot size to 5.44 mm, an increase in system power consumption to 157.86 kW, and a decrease in energy retention rate to 91.97%, indicating that the reverse learning initialization contributes significantly to the population diversity and the uniform distribution of the initial Pareto frontier. Compared with pure NSDE and the population without other modules, the marginal benefit of removing the reverse learning initialization alone on the optimization performance is obvious, proving that this module plays an irreplaceable role in improving the global coverage of the algorithm and avoiding early convergence.

The performance degradation caused by removing the PSO local enhancement module is the most significant, with an efficiency reduction of about 4.74%, an increase in the beam spot size of about 0.56 mm, an increase in system power consumption of about 7.53 kW, and a decrease in energy retention rate of about 3.28%. PSO can fine-tune the speed-position of non-dominated elite individuals, so that the algorithm has higher local search accuracy and target compromise ability in the late convergence stage. After removing this module, the algorithm lacks the ability to refine the elite solution, resulting in a relatively rough boundary of the optimal solution, making it difficult to further improve efficiency or reduce power consumption. The adaptive mechanism plays a secondary role, and reverse learning initialization also plays an important role in the algorithm.

When the adaptive mechanism and PSO are removed at the same time and only the reverse learning initialization is retained, the algorithm still has a certain initial diversity, but the efficiency further drops to 91.56%, indicating that it is difficult to maintain continuous optimization by relying solely on diversity. The removal of the adaptive mechanism means that the mutation factor and crossover probability are no longer dynamically adjusted according to the population state, resulting in the inability of the algorithm to effectively balance global exploration and local development in the middle and late stages, and the solution set is prone to localization, resulting in continuous deterioration of power consumption and beam spot size.

In Table 3, from the HV and SI indicators, it is found that after removing reverse learning, HV is reduced by 0.03 and SI is improved by 0.02, reflecting that it can fully improve the Pareto front coverage. Overall, reverse learning initialization, adaptive parameter control, and PSO local reinforcement work together to build an optimization path that is both “breadth” and “depth”. Z1 provides a starting point for global coverage, Z2 dynamically guides the search strategy, and Z3 enhances the accuracy of boundary solution. The ablation experiment clearly quantifies the marginal contribution of each module to beam performance and system energy consumption, and also verifies its necessity and complementarity in solving high-dimensional, multi-objective conflict problems.

5.8 Comparative experiment of different diversity enhancement initialization strategies

In order to further verify the role of the reverse learning initialization strategy in improving population diversity and optimizing performance, a comparative experiment was designed to compare the reverse learning initialization with two other common diversity enhancement initialization strategies, including: adversarial learning initialization (adversarial perturbation is used to enhance the coverage of the initial population); chaotic mapping-based initialization (introducing Logistic chaotic mapping to generate the initial population and improve the randomness of the distribution).

Experimental steps:

1. Under the improved NSDE algorithm framework, three initialization strategies are used to construct the population (only the initialization module is replaced, and other algorithm parameters are the same).

2. Each initialization strategy is run 30 times under the same decision variable dimension (9 dimensions) and the same experimental conditions, and the main performance indicators are counted.

3. Comparison of beam transmission efficiency, end beam spot size, system power consumption, energy retention rate and SPV (solution diversity index).

The comparison results of different diversity enhancement initialization strategies are shown in Table 6.

Table 6

Table 6. Comparison results of different diversity enhancement initialization strategies.

The results in Table 6 show that different initialization strategies have a significant impact on the multi-objective optimization performance of medical heavy ion accelerators. Reverse learning initialization achieves the best performance in beam transmission efficiency (99.21%), terminal spot size (5.13 mm), system power consumption (153.68 kW), energy retention rate (94.12%) and SPV (0.009), indicating that it has obvious advantages in improving population diversity and guiding the solution set to quickly approach the Pareto frontier. Adversarial initialization can enhance the boundary exploration ability, but the disturbance introduced in some areas leads to slightly inferior system power consumption and solution set uniformity (SPV 0.011). The initialization based on chaotic mapping is highly random, the initial diversity improvement effect is limited, and the convergence accuracy and boundary refinement ability are relatively weak (beam transmission efficiency 97.88%, SPV 0.013). In general, reverse learning initialization has more comprehensive advantages in balancing diversity, convergence speed and solution set refinement.

5.9 Analysis of the impact of local optimization frequency on computational complexity and performance cost ratio

In order to evaluate the trade-off between computational resource consumption and performance improvement of local optimization (PSO) in the multi-objective optimization process, experiments with different local optimization call frequencies were designed. The specific setting is: PSO local enhancement is performed once every T generations, and T is 5, 10, 20, 50 and 100. The average beam transmission efficiency, system running time and unit performance improvement computational cost ratio under each setting are statistically analyzed. The performance cost ratio is defined as: performance improvement (relative to the baseline without PSO local optimization) divided by the corresponding percentage of additional computing time. The experiment was conducted under the same data set and hardware environment, and the average was taken after 30 runs. The results are shown in Table 7.

Table 7

Table 7. Impact of local optimization frequency on performance and computational cost.

As shown in Table 7, with the increase in the frequency of local optimization calls (i.e., T decreases), the overall beam transmission efficiency shows an upward trend, but the computation time increases exponentially, leading to a gradual decrease in the performance-to-cost ratio. Notably, among the five frequency settings in this study, the T = 20 scheme achieved the optimal overall cost-effectiveness, with a significant improvement in beam transmission efficiency (from 94.47% to 99.21%), while the increase in computational cost remained within an acceptable range, resulting in a globally optimal “performance improvement/computational cost” ratio (0.052). In contrast, while higher frequencies like T = 10 or T = 5 slightly increased efficiency, the corresponding computational cost increased sharply, leading to a decrease in overall cost-effectiveness. Therefore, triggering a PSO local enhancement once every T = 20 cycles achieves the best balance between improving the quality of the multi-objective Pareto front and maintaining controllable overall operating costs. Based on this experimental result, this paper uses T = 20 as the final unified setting for the local optimization frequency and maintains consistency throughout all experiments.

5.10 Parameter robustness and sensitivity analysis experiment

In order to evaluate the robustness of the improved NSDE algorithm under different parameter configurations and its sensitivity to the main parameters, this paper designed a parameter scanning experiment. Three key parameters, the initial value of the inertia weight (ωmax), the learning factor (c), and the population size (Np), were selected based on experience and were taken in small, medium, and large ranges to analyze the impact on beam transmission efficiency, system power consumption, and hypervolume (HV). Each group of configurations was run independently 30 times and the average was taken. The experiment was carried out under the same data set, hardware, and multi-objective settings. The results are shown in Table 8.

Table 8

Table 8. Parameter robustness and sensitivity analysis experimental results.

The results in Table 8 show that the improved NSDE algorithm can maintain high beam transmission efficiency (≥97.7%) and low system power consumption (≤156.7 kW) within the range of parameter variation, and remains robust in comprehensive performance indicators such as HV. Parameter combination B (inertia weight 0.9, learning factor 1.8, population size 120) achieved the best transmission efficiency, minimum power consumption and highest HV, verifying the rationality of the parameter configuration. Combination C has a larger population size, but due to the increase in inertia weight, the convergence speed is slightly slower (average 105 generations). Combination A has a lower inertia weight, a smaller learning factor, and limited global exploration ability, resulting in a slightly lower HV than other groups. The experiment further verified the robustness of the improved NSDE to the main parameters, indicating that this method is suitable for multi-objective high-dimensional optimization problems under different optimization requirements, and the parameter settings can be flexibly adjusted according to resource and performance requirements.

5.11 More test case verification

In order to further verify the generalization ability and robustness of the improved NSDE algorithm in multi-objective optimization, this paper adds multiple sets of standard multi-objective optimization (MOO) benchmark test functions and experiments on real-world high-dimensional multi-objective problems. The selected test cases include the classic ZDT series (ZDT1, ZDT3), the DTLZ series (DTLZ2, DTLZ7) and a set of high-dimensional real industrial problems. All algorithms were run independently 30 times under the same parameter configuration and computing resources, and key indicators (the performance of beam transmission efficiency converted into the objective minimization or maximization indicators of the corresponding problem) were statistically analyzed and compared with existing algorithms. Table 9 shows the mean and standard deviation results of IGD and HV of the improved NSDE, NSDE, NSGA-II, and MOEA/D under each test case.

Table 9

Table 9. Test case verification results.

In Table 9, the results clearly demonstrate the superior performance of the improved NSDE algorithm in a variety of standard multi-objective optimization benchmarks and high-dimensional real-world industrial problems. The IGD index of the improved NSDE in all test cases is significantly lower than that of other algorithms. For example, the IGD in the ZDT1 problem is 0.027 ± 0.004, which is better than 0.034 ± 0.006 of NSDE and 0.038 ± 0.007 of NSGA-II, indicating that its solution set is closer to the true Pareto frontier. At the same time, the HV index also performs best. For example, the improved NSDE in the DTLZ2 problem reaches 0.480 ± 0.007, which is significantly better than 0.445 ± 0.010 of NSDE and 0.425 ± 0.012 of MOEA/D, indicating that its solution set covers a wider range and is more evenly distributed in the target space. In high-dimensional and complex real-world industrial problems (50 dimensions), the improved NSDE still maintains a low IGD (0.035 ± 0.005) and a high HV (0.462 ± 0.009), showing good generalization and robustness. The data fully proves that the improved strategy proposed in this paper has significant advantages in improving the convergence and diversity of the algorithm, and can effectively deal with the goal conflicts and high-dimensional challenges in multi-objective optimization.

5.12 PSO local optimization complexity changes with the number of targets and elite size

The curve of PSO local optimization complexity changing with the number of targets and elite size is shown in Figure 12.

Figure 12

Figure 12. PSO local optimization complexity curve with the number of targets and elite size.

In Figure 12, it can be observed that the complexity C increases linearly with the number of objective functions and the size of the elite, indicating that the two have a cumulative effect on the algorithm overhead. When the number of objectives is small (such as M = 1 or 2), the overall complexity is still controllable even if the elite size is large; in the case of multiple objectives (M = 4), the increase in the number of elites will significantly increase the computational burden, but it is still low. This figure provides a quantitative reference for optimizing the T value (PSO triggering frequency) and controlling the size of the elite subset, which helps to strike a balance between performance improvement and computational cost.

5.13 Comparison of the complexity of the baseline NSDE and the improved NSDE

In order to evaluate the impact of the improved NSDE algorithm proposed in this paper on the overall computational complexity, this paper designs a comparative experiment to statistically analyze the differences in the average number of operator calls and execution time of the two algorithms under unit iteration. The results are shown in Table 10.

Table 10

Table 10. Comparison of unit iteration complexity of different algorithms.

In Table 10, the improved algorithm increases the operator call and execution time by about 25% in unit iteration, which is mainly due to PSO local search and parameter adaptive update. However, compared with the significant improvement in optimization accuracy and solution quality, this complexity increase is acceptable, reflecting a good trade-off between performance and overhead.

It should be noted that in practical engineering optimization, such as the design of medical heavy ion accelerators involving FLUKA Monte Carlo simulations, each objective function evaluation can take days or even weeks. In such high-cost scenarios, although the improved NSDE algorithm increases computation time by approximately 25% per iteration, its significant improvements in convergence speed and solution accuracy through reverse learning initialization, adaptive parameter tuning, and PSO local enhancement can reduce the total number of iterations during the overall optimization process, potentially saving more computational resources over the entire optimization cycle. Furthermore, obtaining high-quality solutions is crucial for clinical and engineering safety, making the additional local computational overhead reasonable and necessary in practical applications. Of course, in specific engineering implementations, the performance improvement and computational cost can be balanced by adjusting the local optimization frequency, elite subset size, and parallel computing strategy to meet different optimization accuracy and time constraints.

5.14 Comparative experiment on the impact of different particle numbers (Nc) on algorithm performance and computational cost

To verify the rationality of the selection of the local reinforcement particle number Nc, this paper sets up five experimental groups with Nc = 10, 15, 20, 25, and 30, while keeping all other parameters completely consistent. Each group was run independently 20 times under the same hardware environment (Intel i7-12700, 32 GB RAM). The experimental procedure is as follows: First, the main NSDE framework and all adaptive mechanisms are fixed, and only the number of elite particles participating in the update in the PSO local reinforcement module, Nc, is changed; then, the average values ​of IGD and HV indices and the real-time feedback delay time of a single iteration are recorded for each Nc group; finally, the balance between convergence accuracy, Pareto front coverage, and computational cost of different Nc groups is compared to select the optimal Nc. The results are shown in Table 11.

Table 11

Table 11. Impact of different particle numbers nc on IGD, HV, and computational cost.

As shown in Table 11, when Nc increases from 10 to 20, both IGD and HV continuously improve, indicating that more sufficient local enhancement can improve the accuracy and coverage of the Pareto front. However, when Nc exceeds 20, the change in IGD tends to saturate, the improvement in HV is insufficient, and the computational latency increases significantly (from 132.9 ms to 204.5 ms). Therefore, Nc = 20 achieves the best balance between performance optimization and real-time computation cost, validating the rationality of the parameter setting in this paper.

6 Experimental discussion

The improved NSDE algorithm achieves near-optimal performance in all experimental scenarios, mainly due to the synergistic gain of the triple strategy. The reverse learning initialization provides the algorithm with a wider distribution of initial solutions, significantly enhancing the initial coverage of various regions of the Pareto frontier in the early stage. The adaptive mechanism dynamically adjusts the differential mutation factor and crossover probability according to the fitness of the population, realizes intelligent switching between the global exploration period and the local convergence period, and avoids premature convergence caused by fixed parameters. PSO local reinforcement fine-tunes the speed-position of the non-dominated elite solution, further compresses the precision error of the objective function value in the solution set, and improves the approximation of the solution set in terms of efficiency and energy retention. The above mechanisms complement each other in the high-dimensional coupled parameter space and multi-objective conflict environment, ensuring the stability, convergence speed and diversity coverage of the improved NSDE.

For the control rules of dynamically adjusting the mutation factor and crossover probability, the adaptive mechanism can effectively balance global exploration and local utilization, improve the convergence speed of the algorithm and the diversity of solutions, but the strict mathematical convergence proof is relatively complex and depends on the nature of the specific problem. This paper verifies the stability and convergence performance of the improved NSDE algorithm in multi-objective optimization tasks through a large number of numerical experiments. The results show that the algorithm can gradually approach the Pareto frontier during the iteration process and the solution set maintains good diversity. Moreover, the adjustment of the mutation factor and crossover probability is kept within the preset bounded interval, avoiding search oscillation or premature convergence caused by out-of-control parameters, and ensuring the stable operation of the algorithm in practical applications.

Regarding the universality of parameter settings, the parameter adjustment of the improved NSDE algorithm in this paper is mainly customized based on the specific characteristics of medical heavy ion accelerator beam transmission optimization and the characteristics of multi-objective conflicts. Parameters such as inertia weight, learning factor and population size are the core parameters of evolutionary and particle swarm optimization algorithms. Their reasonable range and adjustment mechanism also have certain reference value in other standard multi-objective test sets (such as DTLZ series and ZDT series). There are differences in the target dimensions, decision space complexity and target conflict degree of different problems, which will cause the optimal configuration of parameters to be offset. Therefore, when applied to other benchmark tests or actual problems, it is still necessary to perform targeted parameter tuning or adopt an adaptive parameter adjustment mechanism in combination with the specific problem characteristics to ensure the convergence and solution set diversity of the algorithm in different multi-objective optimization scenarios. Overall, the adaptive parameter adjustment framework proposed in this paper has good generalization potential and provides an effective idea for facing diversified optimization problems.

In the experiment of this paper, different decision variable dimensions significantly affect the effectiveness and synergy of each component module. As the dimension increases from 9 to 15, the search space expands sharply and the parameter coupling increases, resulting in a rapid performance decline of traditional algorithms (such as NSGA-II and MOEA/D), while the improved NSDE shows stronger robustness due to the combination advantages of its key modules. Reverse learning initialization effectively improves the initial diversity and global coverage of the population in high-dimensional space, especially significantly alleviates the early convergence problem when the dimension increases, and ensures the uniformity of the initial solution set. The adaptive mechanism can dynamically adjust the mutation factor and crossover probability under high-dimensional conditions, balance global exploration and local development, and avoid insufficient search capabilities in high-dimensional space. The contribution of PSO local enhancement is particularly prominent in high dimensions, which makes up for the shortcomings of local fine search of standard NSDE in the late convergence stage and improves the accuracy of high-dimensional Pareto frontier boundary solutions. Ablation experiments show that removing any module leads to a significant decline in performance when the dimension increases. The synergy of the three is particularly necessary in high-dimensional multi-objective optimization, and together constructs an optimization path in high-dimensional space that is both globally diverse and locally fine.

This study is the first to organically integrate reverse learning, adaptive control and PSO elite reinforcement into NSDE for the optimization of medical heavy ion accelerator beam transmission, which has important research significance. This paper constructs a multi-objective optimization framework that combines breadth exploration and depth refinement, which can be extended to other high-dimensional, multi-objective engineering systems. The research significantly improved the beam transmission efficiency and energy retention rate, reduced system power consumption and beam spot size fluctuations, and provided a feasible solution for the dose accuracy and operating cost control of heavy ion therapy systems. The algorithm module constructed by the experiment does not need to be manually adjusted for specific physical models, and can automatically adapt to different target scenarios through adaptive strategies, which improves the deployment efficiency of the algorithm in actual industrial/medical environments.

7 Conclusion

This paper adopts an improved NSDE algorithm that combines reverse learning initialization, adaptive parameter control and PSO local refinement to solve the high-dimensional coupling and multi-objective conflict problems in the beam transmission of medical heavy ion accelerators. (1) A diverse initial population can be constructed through symmetric mapping and non-dominated-crowding screening. (2) The differential mutation factor and crossover probability are adjusted in real time in the main loop to balance global jump out and local convergence. (3) PSO fine-tuning is applied to the Pareto elite solution every 20 generations to refine the target boundary. (1) A diverse initial population can be constructed through symmetric mapping and non-dominated-crowding screening. (2) The differential mutation factor and crossover probability are adjusted in real time in the main loop to balance global jump out and local convergence. (3) PSO fine-tuning is applied to the Pareto elite solution every 20 generations to refine the target boundary. Future work can introduce digital twins and real-time noise models to verify online robustness, develop lightweight parallel deployment solutions, and integrate decision maker preferences and closed-loop feedback to achieve online adaptive optimization of medical heavy ion accelerators.

This paper has achieved some small achievements. The following are the shortcomings of this paper and future research directions.

1. This paper is based on the Monte Carlo coupled interpolation model and does not fully consider the on-site time-varying disturbances, such as magnetic field drift and mechanical vibration. More real-time noise models or digital twin technologies can be introduced in the future to verify the real-time robustness.

2. The improved NSDE is heavily dependent on GPU and multi-core CPU. In the future, a lightweight version or multi-level parallel architecture can be studied to adapt to resource-constrained edge computing devices;

3. The current algorithm outputs a complete Pareto frontier and lacks a human-computer interactive decision support layer. In the future, it can be combined with decision maker preferences or a dynamic weighting mechanism based on utility functions to improve the efficiency of final solution selection;

4. The boundary limits of decision variables have not been systematically tested, such as stability under extremely high currents or extremely narrow apertures. Robustness and safety analysis under boundary conditions can be carried out in the future.

5. In the future, the algorithm can be integrated with the online data feedback closed loop to achieve adaptive real-time optimization of parameters under the accelerator operation state, and promote the transformation of the algorithm from offline design to online control.

Data availability statement

The original contributions presented in the study are included in the article/supplementary material, further inquiries can be directed to the corresponding author.

Author contributions

YY: Writing – review and editing, Writing – original draft. MZ: Writing – original draft, Data curation. KW: Writing – review and editing, Conceptualization.

Funding

The author(s) declared that financial support was not received for this work and/or its publication.

Conflict of interest

The author(s) declared that this work was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Generative AI statement

The author(s) declared that generative AI was not used in the creation of this manuscript.

Any alternative text (alt text) provided alongside figures in this article has been generated by Frontiers with the support of artificial intelligence and reasonable efforts have been made to ensure accuracy, including review by the authors wherever possible. If you identify any issues, please contact us.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

References

Ahmad, M. F., Isa, N. A. M., Lim, W. H., and Amg, K. M. (2022). Differential evolution with modified initialization scheme using chaotic oppositional based learning strategy. Alexandria Eng. J. 61 (12), 11835–11858. doi:10.1016/j.aej.2022.05.028

CrossRef Full Text | Google Scholar

Cao, P., and Huang, Q. (2024). Hybrid multi-strategy improved butterfly optimization algorithm. Appl. Sci. 14 (24), 1–14. doi:10.3390/app142411547

CrossRef Full Text | Google Scholar

Cao, D., Xu, Y., Yang, Z., Dong, H., and Li, X. (2023). An enhanced whale optimization algorithm with improved dynamic opposite learning and adaptive inertia weight strategy. Complex Intelligent Syst. 9 (1), 767–795. doi:10.1007/s40747-022-00827-1

CrossRef Full Text | Google Scholar

Cheng, L., Zhou, J. X., Hu, X., Mohamed, A. W., and Liu, Y. (2024). Adaptive differential evolution with fitness-based crossover rate for global numerical optimization. Complex Intelligent Syst. 10 (1), 551–576. doi:10.1007/s40747-023-01159-4

CrossRef Full Text | Google Scholar

Farda, I., and Thammano, A. (2023). An improved differential evolution algorithm for numerical optimization problems. HighTech Innovation J. 4 (2), 434–452. doi:10.28991/hij-2023-04-02-014

CrossRef Full Text | Google Scholar

Gad, A. G. (2022). Particle swarm optimization algorithm and its applications: a systematic review. Arch. Comput. Methods Eng. 29 (5), 2531–2561. doi:10.1007/s11831-021-09694-4

CrossRef Full Text | Google Scholar

Gao, X., He, F., Zhang, S., Luo, J., and Fan, B. (2024). A fast nondominated sorting-based MOEA with convergence and diversity adjusted adaptively. J. Supercomput. 80 (2), 1426–1463. doi:10.1007/s11227-023-05516-5

CrossRef Full Text | Google Scholar

Ge, Y., Zhong, Y., Murata, I., Tamaki, S., Yuan, N., Sun, Y., et al. (2024). Efficient optimization of an accelerator neutron source for neutron capture therapy using genetic algorithms. Med. Phys. 51 (9), 6445–6457. doi:10.1002/mp.17132

PubMed Abstract | CrossRef Full Text | Google Scholar

Han, M. C., Choi, S. H., Hong, C. S., Kim, Y. B., Koom, W. S., Kim, J. S., et al. (2024). The first Korean carbon-ion radiation therapy facility: current status of the heavy-ion therapy center at the onsei cancer center. Radiat. Oncol. J. 42 (4), 295–307. doi:10.3857/roj.2024.00206

PubMed Abstract | CrossRef Full Text | Google Scholar

Huang, C., Zhu, J., and Xu, M. (2025). DyS-MPADE: a novel multipopulation adaptive differential evolution methodology based on dynamic subpopulation. J. Comput. Des. Eng. 12 (3), 204–225. doi:10.1093/jcde/qwaf024

CrossRef Full Text | Google Scholar

Ji, F., Edelen, A., Roussel, R., Shen, X., Miskovich, S., Weathersby, S., et al. (2024). Multi-objective Bayesian active learning for MeV-ultrafast electron diffraction. Nat. Commun. 15 (1), 4726–4740. doi:10.1038/s41467-024-48923-9

PubMed Abstract | CrossRef Full Text | Google Scholar

Kuo, R. J., Song, P. F., Nguyen, T. P. Q., and Yang, T. J. (2023). An application of multi-objective simulation optimization to medical resource allocation for the emergency department in Taiwan. Ann. Operations Res. 326 (1), 199–221. doi:10.1007/s10479-023-05374-7

CrossRef Full Text | Google Scholar

Li, Z., and Tam, V. (2024). AdaGuiDE: an adaptive and guided differential evolution for continuous optimization problems. Appl. Intell. 54 (21), 10833–10911. doi:10.1007/s10489-024-05675-9

CrossRef Full Text | Google Scholar

Li, X., Wu, H., Yang, Q., Tan, S., Xue, P., and Yang, X. (2022). A multistrategy hybrid adaptive whale optimization algorithm. J. Comput. Des. Eng. 9 (5), 1952–1973. doi:10.1093/jcde/qwac092

CrossRef Full Text | Google Scholar

Liang, X., Mohammadi, H., Moreno, K. C., Beltran, C. J., and Holtzman, A. L. (2025). Heavy ion particle therapy in modern day radiation oncology. Hematology/Oncology Clin. 39 (2), 377–397. doi:10.1016/j.hoc.2024.11.007

PubMed Abstract | CrossRef Full Text | Google Scholar

Lin, A., Liu, D., Li, Z., Hasanien, H. M., and Shi, Y. (2023). Heterogeneous differential evolution particle swarm optimization with local search. Complex Intelligent Syst. 9 (6), 6905–6925. doi:10.1007/s40747-023-01082-8

CrossRef Full Text | Google Scholar

Liu, Y., Zhao, Y., He, Y., Zhang, Z., and Li, A. (2024a). Intelligent control ystem for the hard X-Ray nanoprobe beamline beam optimization based on automatic evolution algorithm and expert system. Sensors (Basel, Switz.) 24 (22), 7211–7227. doi:10.3390/s24227211

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, R., Fang, R., Zeng, T., Fei, H., Qi, Q., Zuo, P., et al. (2024b). A novel adaptive sand cat swarm optimization algorithm for feature selection and global optimization. Biomimetics 9 (11), 701–724. doi:10.3390/biomimetics9110701

PubMed Abstract | CrossRef Full Text | Google Scholar

Ma, H., Zhang, Y., Sun, S., Liu, T., and Shan, Y. (2023). A comprehensive survey on NSGA-II for multi-objective optimization and applications. Artif. Intell. Rev. 56 (12), 15217–15270. doi:10.1007/s10462-023-10526-z

CrossRef Full Text | Google Scholar

Pivi, M. T. F. (2024). The MedAustron particle therapy accelerator facility. Health Technol. 14 (5), 919–928. doi:10.1007/s12553-024-00840-z

CrossRef Full Text | Google Scholar

Qtaish, A., Braik, M., Albashish, D., Alshammari, M. T., Alreshidi, A., and Alreshidi, E. J. (2025). Enhanced coati optimization algorithm using elite opposition-based learning and adaptive search mechanism for feature selection. Int. J. Mach. Learn. Cybern. 16 (1), 361–394. doi:10.1007/s13042-024-02222-3

CrossRef Full Text | Google Scholar

Ramirez-Ochoa, D. D., Perez-Dominguez, L. A., Martínez-Gomez, E. A., and Luviano-Cruz, D. (2022). PSO, a swarm intelligence-based evolutionary algorithm as a decision-making strategy: a review. Symmetry 14 (3), 455–477. doi:10.3390/sym14030455

CrossRef Full Text | Google Scholar

Roussel, R., Hanuka, A., and Edelen, A. (2021). Multiobjective bayesian optimization for online accelerator tuning. Phys. Rev. Accel. Beams 24 (6), 1–14. doi:10.1103/physrevaccelbeams.24.062801

CrossRef Full Text | Google Scholar

Shami, T. M., El-Saleh, A. A., Alswaitti, M., Ai-Tashi, Q., Summakieh, M. A., and Mirjalili, S. (2022). Particle swarm optimization: a comprehensive survey. IEEE Access 10 (1), 10031–10061. doi:10.1109/access.2022.3142859

CrossRef Full Text | Google Scholar

Shayanmoghadam, A. A., Hosseinimotlagh, S. N., and Reisosadat, S. M. R. (2025). Heavy ion fusion simulation in key parameters optimization using multilayer fuel target and the effectiveness of ponderomotive force in ion acceleration. Indian J. Phys. 99 (5), 1909–1928. doi:10.1007/s12648-024-03397-9

CrossRef Full Text | Google Scholar

Sokol, O., and Durante, M. (2023). Carbon ions for hypoxic tumors: are we making the most of them?. Cancers 15 (18), 4494–4511. doi:10.3390/cancers15184494

PubMed Abstract | CrossRef Full Text | Google Scholar

Sun, L., Si, S., Ding, W., Xu, J., and Zhang, Y. (2023). BSSFS: binary sparrow search algorithm for feature selection. Int. J. Mach. Learn. Cybern. 14 (8), 2633–2657. doi:10.1007/s13042-023-01788-8

CrossRef Full Text | Google Scholar

Wang, M., Ma, Y., and Wang, P. (2022). Parameter and strategy adaptive differential evolution algorithm based on accompanying evolution. Inf. Sci. 607 (1), 1136–1157. doi:10.1016/j.ins.2022.06.040

CrossRef Full Text | Google Scholar

Wongsa, W., Puphasuk, P., and Wetweerapong, J. (2024). Differential evolution with adaptive mutation and crossover strategies for nonlinear regression problems. Bull. Electr. Eng. Inf. 13 (5), 3503–3514. doi:10.11591/eei.v13i5.6417

CrossRef Full Text | Google Scholar

Xu, C., Boltz, T., Mochihashi, A., Santamaria Garcia, A., Schuh, M., and Muller, A. S. (2023). Bayesian optimization of the beam injection process into a storage ring. Phys. Rev. Accel. Beams 26 (3), 1–10. doi:10.1103/physrevaccelbeams.26.034601

CrossRef Full Text | Google Scholar

Xu, H., Deng, Q., Zhang, Z., and Lin, S. (2025). A hybrid differential evolution particle swarm optimization algorithm based on dynamic strategies. Sci. Rep. 15 (1), 4518–4566. doi:10.1038/s41598-024-82648-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Yamada, S., Takiyama, H., Isozaki, Y., Shinoto, M., Ebner, D. K., Koto, M., et al. (2022). Carbon ion radiotherapy for locally recurrent rectal cancer of patients with prior pelvic irradiation. Ann. Surg. Oncol. 29 (1), 99–106. doi:10.1245/s10434-021-10876-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Yan, W., Bai, Y., Xu, R., and Neculaes, V. B. (2022). X-ray source design optimization using differential evolution algorithms—A case study. Rev. Sci. Instrum. 93 (5), 1–24. doi:10.1063/5.0079389

PubMed Abstract | CrossRef Full Text | Google Scholar

Yan, Y., Cheng, H., Li, Y., Fang, Y., Xia, Y., He, Q., et al. (2024). High transmission efficiency collection system for laser-accelerated proton beams based on permanent magnet quadrupoles prefocusing. Phys. Rev. Accel. Beams 27 (5), 1–9. doi:10.1103/physrevaccelbeams.27.052801

CrossRef Full Text | Google Scholar

Yang, X., and Qiu, Y. (2023). Research on improving gray wolf algorithm based on multi-strategy fusion. IEEE Access 11 (1), 66135–66149. doi:10.1109/access.2023.3289819

CrossRef Full Text | Google Scholar

Zhong, R., and Yu, J. (2024). DEA 2 H 2: differential evolution architecture based adaptive hyper-heuristic algorithm for continuous optimization. Clust. Comput. 27 (9), 12239–12266. doi:10.1007/s10586-024-04587-0

CrossRef Full Text | Google Scholar