A conjugate gradient-assisted multi-objective evolutionary algorithm for fluence map optimization in radiotherapy treatment

Intensity-modulated radiotherapy (IMRT) is one of the most applied techniques for cancer radiotherapy treatment. The fluence map optimization is an essential part of IMRT plan designing, which has a significant impact on the radiotherapy treatment effect. In fact, the treatment planing of IMRT is an inverse multi-objective optimization problem. Existing approaches of solving the fluence map optimization problem (FMOP) obtain a satisfied treatment plan via trying different coupling weights, the optimization process needs to be conducted many times and the coupling weight setting is completely based on the experience of a radiation physicist. For fast obtaining diverse high-quality radiotherapy plans, this paper formulates the FMOP into a three-objective optimization problem, and proposes a conjugate gradient-assisted multi-objective evolutionary algorithm (CG-MOEA) to solve it. The proposed algorithm does not need to set the coupling weights and can produce the diverse radiotherapy plans within a single run. Moreover, the convergence speed is further accelerated by an adaptive local search strategy based on the conjugate-gradient method. Compared with five state-of-the-art multi-objective evolutionary algorithms (MOEAs), the proposed CG-MOEA can obtain the best hypervolume (HV) values and dose–volume histogram (DVH) performance on five clinical cases in cancer radiotherapy. Moreover, the proposed algorithm not only obtains the more optimal solution than traditional method used to solve the FMOP, but also can find diverse Pareto solution set which can be provided to radiation physicist to select the best treatment plan. The proposed algorithm outperforms dose-volume histogram state-of-the-art multi-objective evolutionary algorithms and traditional method for FMOP on five clinical cases in cancer radiotherapy.


Introduction
Malignant tumor (also called cancer) is the leading cause of death worldwide. The global cancer burden has risen to 9.3 million new cases and 10.0 million deaths in 2020 [1]. One in five people suffer from the disease in lifetime, one in eight men and one in eleven women die from cancer [2].
Radiotherapy is one of the most three treatment practices, especially when the tumor is localized in a one part of the body. Nearly seventy percent of cancer patients receive radiotherapy during their illness, sometimes combined with surgery or chemotherapy. Radiotherapy damage the malignant tumor cells using ionizing radiation. However, it is inescapable that the radiation may damage healthy tissue, even using the most advanced treatment technology. Radiation to healthy tissue may lead to radiotherapy complications, which may have a profound and lasting negative impact on the patient's well-being. Hence, the purpose of the radiotherapy treatment is to deliver precise high radiation dose to tumor target volume, while the absorbed dose in surrounding organ at risk (OAR) should simultaneously be minimized to reduce damage to healthy organs. It can be found that the planning of radiotherapy treatment needs to balance the high dose of the tumor and different dose constraints of many healthy organs with the purpose to achieve high quality life for the patient. Moreover, it requires balancing 10-30 highly correlated objectives or constraints, because each patient is anatomically unique and needs a unique individualized plan.
IMRT is one of the most commonly used delivery techniques in clinic. By adjusting the intensity of each radiation field, IMRT makes the shape of the high dose area conformal with the tumor target volume, while the dose of other parts is as low as possible. Figure 1 shows the treatment machine of medical linear accelerator, where the treatment planning for IMRT is to determine the machine parameters, including orientation of radiation field (also called beam angle) and the fluence map of each field. These parameters are commonly determined using forward planning and inverse planning. Forward planning based on the experience of the treatment planner is widely used for the selection of the orientation of radiation field, while inverse planning is the only feasible mean to determine the fluence map of each field. Inverse planning transforms the determination of the treatment parameters to an optimization problem where the expected dose distribution (the prescription to all targets and dose constraints to elsewhere) is given and the fluence map that realizes this dose distribution is unknown [3]. The optimization method was used to determine the fluence map that best produces the desired dose. In fact, FMOP is a multiobjective optimization problem that should balance the high radiation dose to tumor targets, while the dose in surrounding OARs should be minimized simultaneously. Clinically, all the OARs and targets collectively called region of interest (ROI). However, most of existing methods translate the FMOP into a single-objective optimization problem by calculating the weighted summation of all objectives and the constraints [4], and then optimize it using the analytical optimization method such as gradient-based method [5,6]. The coupling weights used in these approaches should be determined before the optimization. However, it is difficult to determine these weights in advance. Hence, the optimization procedure is a trial and error process, resulting in a timeconsuming and laborious planning.
To overcome the above issue, we propose a conjugate gradient-assisted multi-objective evolutionary algorithm (CG-MOEA) for the FMOP. The proposed algorithm can provide the Pareto set of fluence map set to radiation physicist within a single run. Moreover, the conjugategradient method (CG) as a local search strategy is adaptively embedded into the evolutionary algorithm for accelerating the convergence speed. Specifically, the contributions of this paper are summarized as follows: • The FMOP is formulated into a multi-objective optimization problem, where many dose objectives and constraints to different ROI volumes are modeled to three objectives optimization problem by analysing the relation of all the dose requirements. The problem formulation is suitable for different clinical cases and can be solved by multiobjective evolutionary algorithm. • An effective multi-objective evolutionary algorithm is proposed to solve the constructed multi-objective FMOP.
In the proposed algorithm, the CG is used as a local search to speed up the convergence, and an adaptive strategy for determining the frequency of using CG is proposed to keep population diversity. The proposed algorithm can guarantee good convergence and diversity of the obtained Pareto solutions. • The proposed algorithm is compared with several stateof-the-art large-scale optimization algorithms and one traditional analytic method on five clinical cases. Experimental results demonstrate that the proposed algorithm can provide a set of radiotherapy plans with good diversity and convergence within a single run. The radiation physicist only need select the satisfy plan in the Pareto set, without specifying weights or the priority list of dose requirements.
The rest of the paper is organized as follows. In the next section, the preliminaries about the FMOP are first given and then some related work for solving the FMOP is presented. In the following section, the multi-objective problem formulation for FMOP and the details of the proposed algorithm CG-MOEA are presented. In the next section, the empirical results and discussions are reported. Finally, the last section makes a conclusion about this paper and points out some future research issues.

Preliminaries of FMOP
In the optimization of FMOP, the radiation field is discretized into beamlets as shown in Fig. 2, which are the decision variables of the FMOP. The numerical value of a single beamlet represents the intensity fluence that passes through the grid element of the beamlet, and delivers the dose to the patient body. The body of a patient is discretized into voxels. Rather considering the dose to the entire patient body, the desired dose distribution is given per ROI including all the OARs and targets, and optimized by one or more costfunctions. The voxel sampling resolution in the ROI volumes is usually lower than the original resolution of the computed tomography (CT). To keep the optimization problem manageable, different ROI volumes are sampled with variable resolutions, so that important volumes can be sampled more accurately. The relation between the beamlet x and the voxel dose d is expressed as follows: where w is the dose influence matrix representing the dose influence of each beamlet to every voxel, which is computed by dose computation algorithm based on the CT scans and radiation field parameters, such as beam energy and angle.
For better understanding the parameters in the FMOP, the fluence map of the first beam is presented in Fig. 2.
The desired dose distribution is often given by radiation oncologists based on dose-based functions or functions of DVH. The usually used function type is shown in Table 1. Also the symbol descriptions in this study is presented in Table 2. For a target, one MaxDose, one MinDose and one UniformDose are often given in clinic. For an OAR, one Max-Dose is usually used. A dose-based function gives a penalty for deviations between the calculated dose d and a reference dose leveld according to [7], which is expressed as Eq. (2).
where d calculated according to Eq. (1),d is given by radiation oncologist and different ROI volumes are given with various dose. The MinDose function only penalizes the dose below the reference level and the MaxDose function penalizes the dose above the reference level. Where the δ(y) = min{y, 0} for MinDose, δ = max{y, 0} for Max-Dose.
Functions of dose-volume histogram involve two kinds of MinDVH and MaxDVH. MaxDVH functions are commonly used for OARs that exhibit a maximum sub-volume where the dose is equal to or higher than the reference dose. Such as lung or kidney, it is allowable to deliver a high dose to a subvolume as long as this sub-volume is less than the permitted volume. MinDVH function are used when the target volume overlaps with a dose-limiting OAR, where it can accept a controlled underdosage of a sub-volume of the target. The DVH function penalizes the deviation between the DVH curve associated with a given sub-volume r and the reference dose leveld according to Eq. (3) which is defined as [8]:

Related work for the FMOP
FMOP of radiotherapy has always been a multi-objective optimization problem with many constraints. It is a major challenge to directly solve such multi-objective problem from a mathematical perspective. Therefore, the single objective optimization using unconstrained quadratic programming by weighted summation each objectives has been widely used in commercial treatment planning system (TPS) [9,10]. However, the commonly used TPS requires radiation physicist to input not only desired dosimetric goals but also the weight factor of each goal, which represent the importance of goals and are difficult to assign before optimizing. The radiation physicist must manually adjust these parameters through a trial-and-error process [11,12], which results in a time-consuming and laborious planning proce-dure. Moreover, the quality of the final plan highly depends on the radiation physicist experience [13][14][15]. Some studies about automatic multi-criteria plan optimization have been attempted in recent years as shown in Table 3. The approaches using predefined constraint priority list and specific constraint adjustment mechanism have been studied for solving the FMOP. Wilkens et al. [16] proposed a four-step procedure and evaluated head and neck cases with the proposed procedure. Jee et al. [17] implemented a lexicographic ordering (LO) method, which was applied to prostate, head and neck cases. The optimization objectives and constraints are pre-categorized into four priority levels based on the clinical importance. Breedveld et al. [18] used an automatic constraint adjustment mechanism based on a wish list to produce a plan with all constraints required as well as possible. Based on the lexicographic ordering, Breedveld et al. [19] later proposed an improved automatic multi-criteria plan optimization model and showed the higher quality than that of the original plans from the dosimetrists. However, most  of those studies commonly generate the sub-optimal plan because of an incomplete mathematical solution space, and need to determine the importance parameters or priority or wish list of all the objective and constraint in advance. Therefore, we propose a multi-objective evolutionary optimization framework for the FMOP without the specification of weighting factors or priority list. First, we transform the FMOP into a three-objective optimization problem, and then propose an effective multi-objective evolutionary algorithm to obtain the Pareto solution set. Multi-objective optimization is one of the most popular fields in evolutionary computation, and many effective algorithms have been proposed [20,21]. Based on the evolutionary optimizer in [22], the beam angle optimization as demonstrated in Fig. 3 is presented in a monolithic multi-criteria setting. However, it is necessary to point out that there are fewer MOEAs proposed for solving FMOP. To the best of our knowledge, this work is the first study to introduce the multi-objective evolutionary optimization approach for the fluence map optimization in radiotherapy.

A conjugate gradient-assisted multi-objective evolutionary algorithm for FMOP
In this section, we will propose a conjugate gradient assisted multi-objective evolutionary algorithm for FMOP and the framework of the proposed approach is depicted in Fig. 4. As can be seen from the figure, the FMOP is first transformed into a three-objective optimization problem based on the dose requirements given by radiation oncologist. Then, the CG-MOEA is proposed to find the fluence map optimal set of treatment plan.

Multi-objective optimization problem formulation for FMOP
For the FMOP, it contains five types of constraints or objectives functions mentioned in "Preliminaries of FMOP", including the Mindose, Maxdose, MinDVH, MaxDVH, and UniformDose. An OAR often includes two or more types of constraints as shown in Table 4. A treatment target often Fig. 4 The framework of the proposed multi-objective evolutionary for FMOP contains more than two objectives or constraints, where treatment targets usually include PTV, CTV and multiple GTVs. It can be seen that the FMOP have many dose objectives or constraints, so how to transform these clinical dose requirements into a multi-objective problem is very difficult and a key issue. In the paper, all the dose objectives and dose-volume histogram constraints is formulated into a three-objective optimization problem, where all the Maxdose and MaxDVH functions are aggregated into one objective, all the Mindose and MinDVH functions are aggregated into the second objective, and all the UniformDose functions of targets are summed to the third objective. The minimum required constraint function will become easy to be satisfied when the dose is increased. On the contrary, the maximum required constraint function is more likely to be optimized when the dose is reduced. Besides, the uniform dose objectives is determined by the overall dose distribution. For this reason, these constraints are aggregated into three objectives to be optimized by the multi-objective evolutionary algorithm, and the influence of constraints aggregation will be described in the experimental section. Finally, the multi-objective optimization problem of FMOP can be formulated as Eq. (4).
where x is a solution containing all the fluence map of every field, d v (x) is the computed dose of all the voxels in the ROI volume v, when v ∈ Targets, v may be PTV, CTV or GTV, when v ∈ OARs, v represents each organ in risk that is be protected during radiotherapy treatment.d v is the reference dose level and ratio v is the specified volume ratio of sampling voxels dose exceeding the reference dose level for volume v. Different reference dose level and volume ratio are given to each volume v. MinDose, MaxDose and UniformDose function are formulated as Eq.
(2). The MinDVH and MaxDVH functions are shown as Eq. (3). Many multi-objective optimization method have been proposed [23], and MOEAs have been widely employed to solving multi-objective optimization problems through the population-based stochastic search process [24][25][26]. In the last 20 years, various multi-objective evolutionary algorithms have been proposed for complex optimization problems. In general, two important components in MOEAs have a strong Table 4 Dose requirements and radiation field parameters of five clinical cancer cases impact on the solving effect of the algorithm for different kinds of problems, the search strategy and the environmental selection strategy. The former involves how to produce promising solutions by designing effective search operators, including MOEA/DVA [27], MOEA/PSL [26] and SparseEA [28], especially MOEA/PSL can solve multi-objective optimization problems with sparse optimal solutions. The latter considers how to survive the high-quality solutions into the next generation. For example, in [29], the dominance relation is applied to select promising solutions in the environmental selection stage. Besides, the association relation between reference vectors or points is also used for environmental selection, such as MOEA/D [30], AR-MOEA [25], NSGA-III [31] and so on. Therefore, in the proposed algorithm, we elaborate these two components for solving the specific FMOP.

Framework of the proposed CG-MOEA
As mentioned above, the advantage of the multi-objective evolutionary algorithm is that the population-based stochastic search strategy can provide diverse solutions at a time.
And the conjugate gradient method can be used to accelerate the search process due to the gradient of the fluence map optimization. Thus, a conjugate gradient-assisted multi-objective evolutionary algorithm is designed for solving the FMOP. To be specific, the conjugate gradient method as a local search operator is adaptively embedded into the evolutionary process to generate promising solutions and thus speed up the convergence. In the meantime, a new environmental selection is designed for the FMOP, which can adaptively balance the convergence and the diversity during the evolutionary process.
The procedure of the CG-MOEA is detailed in the Algorithm 1. First, the population, Population, is randomly initialized. Then the sparsity solution generator suggested in [26] is employed to generate new solutions, O f f spring. In the third step, the conjugate gradient method is applied to part of solutions, which includes two parts: determinate the number of solutions using the conjugate gradient optimization approach and generate accelerator solutions by the conjugate gradient. Finally, the population, Population, is updated by all new solutions. And the details of each step will be described in the following sections.

Adaptive conjugate gradient approach
The conjugate gradient approach is characterized by the properties of fast convergence and low memory requirements [32]. The effectiveness of conjugate gradient approach has been widely validated in many treatment planning systems for the FMOP [33]. In this paper, we use the conjugate gradient approach as a local search strategy for accelerating the evolutionary search and generating the promising solutions with good convergence. According to the conjugate gradient approach mentioned in [34], the new solution is generated by the following rule listed as Eq. (5).
where x 0 ∈ R n is initial value of decision variable that is a solution need to be optimized by CG in this study and k is the number of update times. Besides, α k is a search step at x k calculated according to Golden section algorithm, d k is the search direction of kth step and its updating rule is listed as Eq. (6).
where the d 0 is the negative gradient at x 0 , g k as the gradient ∇ f (x k ) of f at x k that is a decision vector and β k is the updating parameter updating by Eq. (7), which is suggested by Fletcher and Reeves [32].
Since the convergence speed of the gradient-based search is faster than the stochastic search method suggested in [26], the frequency of using conjugate gradient method has a significant impact on the solution diversity of the population. Therefore, the conjugate gradient process is conducted every Max Eval * β evaluation, where the value of parameter β will be determined in the experimental section. Considering the influence of the number of accelerator solutions generated by conjugate gradient approach, an adaptive strategy is proposed based on the change of the population diversity. The details of the proposed adaptive strategy can be found in Algorithm 2.

Algorithm 2 AdaptiveConjugateGradient
Require: Current population Population, the diversity degree Div and the number of accelerator solution in last generation η. To be specific, η as the number of the accelerator solutions is reduced when the diversity of the current population is worse than the last updated population; otherwise, it will be increased. And the details of the diversity indicator suggested in [35] is described in Eq. (8).

S P(P)
where d i is the minimum distance between solution i and the other solution in population P, and d is mean value of all d i . Once the number of the accelerator solutions η is determined, solutions of η will be randomly selected from the current population and optimized by the conjugate gradient method.

Environmental selection
In the environmental selection procedure, the solution set O f f spring generated by the evolutionary search operator and Con Population optimized by the conjugate gradient method is used to update Population. Moreover, a new fitness function is designed to balance the diversity and convergence of the population, which is detailed in Eq. (9).
where Div(x i ) is the diversity of solution x i , Cov(x i ) calculated according to Eq. (10) represents the convergence of ith solution. Clearly, the selection mechanism prefers to select solutions with good diversity in the early stage of search process. Considering the convergence of solutions, the solutions with better convergence will be selected when the evolution times is meeting Max Eval. And the convergence indicator inspired by MOEA/D in [30] is employed to calculate the convergence of each solution, which is described as follows. Cov where the maximum kth objective value is set as convergence fitness. After obtaining the convergence of the ith solution, its diversity is calculated by Eq. (11).
where dist(x i , x j ) i = j represents the Euclidean distance between x i and its closet solution in Population. In the selection process, if the number of current non-dominated solutions is larger than the population size N , the solution with worst diversity calculated by Eq. (11) will be removed until the number of non-dominated solutions is equal to N ; otherwise, the solutions will be sorted according to Eq. (9). The details of the selection procedure are given in Algorithm 3. Remove the solution with the worst diversity until the solution number is equal to N , and update fitness. 8: else 9:

Algorithm 3 EnvironmentalSelection
All non-dominated solutions is set as Population. 10: end if 11: Return Population.
The values with bold mean the best value among the compared algorithms

Computational complexity analysis
The time complexity analysis of the proposed CG-MOEA is presented in this section. The most time-consuming operations in CG-MOEA include the solution generator in [26], the adaptive conjugate gradient strategy, and environmental selection. The time complexity of generating solutions is O(N · D · K ), where N , D, and K represent the population size, the number of decision variables, and the hidden layer size, respectively [26]. For the adaptive conjugate gradient strategy, the time complexity for the adaptive process is β · Max Eval N , and the conjugate gradient method requires I · η · D, where I is the iteration number of conjugate gradient process. Hence, the time complexity of the adaptive conjugate strategy is O( β·Max Eval N +I ·η· D), which is approximately equal to O(N + D) due to the small value of β · Max Eval , I and η. In the environmental selection stage, the time complexity of non-dominated sorting is M N 2 [29] and the time complexity of updating solutions is N 3 . Thus, the overall time complexity is N (D · K + N 2 ) for executing the proposed CG-MOEA once.

Experimental results and analysis
In this section, we verify the effectiveness of the proposed multi-objective FMOP model and analyze the performance of the proposed CG-MOEA. Specifically, the experiments are set as follows.

Datasets
In our experiments, we adopt five clinical IMRT cases with different sites. The characteristics of five clinical cases and the prescription requirements of a clinical treatment plan are shown in Table 4. The esophageal cancer case contains nine constraints, including four constraints for PTV-Plan, one constraint for the spinal cord, two constraints for two lungs, and two constraints for the heart. And five beams are generated

Comparison algorithms
To verify the performance of the proposed CG-MOEA, we first select five well-known multi-objective evolutionary algorithms, including NSGA-II [29], MOEA/D [30], RVEA [41], MOEA/DVA [27] and MOEA/PSL [26] to analyze the algorithm diversity and convergence. NSGA-II selects solutions based on the Pareto-dominance relation. MOEA/D transforms the multi-objective optimization problem into single-objective optimization problem by conducting the decomposition strategy in objective space. RVEA develops a set of reference vectors to maintain the diversity of solutions. Different from the general evolutionary algorithms mentioned above, MOEA/DVA is designed for multi-objective problems with a large number of decision variables by decision variable analysis. Besides, MOEA/PSL adopts the learning-based solution generator and traditional operator to search optimal solutions with the sparse property. Besides, the parameters of all comparisons are set to defaults according to [26,27,29,30,41], respectively. Finally, to evaluate the clinical effectiveness, we adopt the conjugate gradient that is commonly used for FMOP in clinic, in which FMOP is translated to a unconstrained quadratic programming by weighted summation each objective.

Genetic operators
In NSGA-II, MOEA/D, RVEA and MOEA/DVA, the simulated binary crossover [36] and the polynomial mutation [37] are used to generate new solutions, where the crossover probability is set to 1, the mutation probability of each dimension is set to 1 D (D is the number of decision variables), and the distribution index is set to 20. In the MOEA/DVA, a differential operator is employed to generate solutions [38], where    The values with bold mean the best value among the compared algorithms the probability of crossover C R is set to 1, and the scale factor F is set to 0.5. MOEA/PSL adopts two types of solution generators, including the simulated binary crossover and polynomial mutation for real variation, and the one point crossover and bitewise mutation for binary variation, the probability of crossover and mutation are set 1 and 1 D , respectively.

Parameter settings
For the fairness of comparison, the population size N and the maximum number of evolutions Max Eval are set to the same values in all comparison multi-objective evolutionary algorithms, where N is set to 100 and Max Eval is set to 1,000,000. For the single objective optimization algorithm CG compared with the CG-MOEA, all the weights of all the dose objectives and constraints are selected by tryerror-try until all the constraints and objectives meet clinical requirements. The iteration of the conjugate gradient method embedded in the proposed CG-MOEA is set to 2 for the fast convergence.

Performance metrics
To measure the quality of solutions obtained by the six multiobjective evolutionary algorithms, the hypervolume (HV) that is widely used performance indicator [39] is adopted in this paper. The HV indicator represents the volume of the region dominated by the obtained optimal solutions, which can be able to account for both the convergence and the diversity of the achieved non-dominated solutions. The larger the HV value, the better the quality of the solution set is. Let P be the set of final non-dominated solutions obtained in the objective space by an algorithm, and z = (z 1 , z 2 , . . . , z n ) T be a reference point in the objective space which is dominated by all Pareto optimal solutions. Then the hypervolume indicator shown as Eq. (12) measures the volume of the region dominated by P and bounded by z.
where n is the number of objectives and f is the objective vector.
The second metric is DVH of all the ROI volumes that is often used in radiotherapy to measure the quality of a radiotherapy plan [13]. A DVH curve represents the dose statistics of a ROI, where the horizontal axis represents the dose value and the vertical axis represents the volume. Therefore, it depicts the portion of the volume of the ROI that receives the amount of dose or higher. For OARs, the curves are ide- ally close to the origin, where the part of the volume that receives a high dose is minimal. For targets, on which often a minimum and a maximum dose are imposed, the curve is expected to be at the right, with a steep slope downward at the end. Furthermore, the Wilcoxon rank-sum test with a significance level is employed to perform statistical analysis, where "-" indicates a significantly better performance of CG-MOEA against another comparison MOEA, and "+" indicates CG-MOEA performs significantly worse against another MOEA, and the symbol "≈" indicates a significantly similar performance of CA-EMOA to another MOEA.

Experimental results
The construction of the optimization model is a key step in multi-objective optimization, hence we first verify the effectiveness of the proposed multi-objective FMOP model. Moreover, to verify the performance of the proposed CG-MOEA, the proposed algorithm is compared with several existing multi-objective evolutionary algorithms and single objective algorithm on some clinical cases, the effectiveness of the proposed adaptive conjugate gradient approach is demonstrated, and the sensitivity analysis on the parameter β is performed. All of the experiments are conducted on the evolutionary optimization platform PlatEMO [40].

Effectiveness of the proposed FMOP model
In this study, all constraints are regarded as objective functions and all objectives are transformed into three objectives to be optimized, including the minimum objective function, the maximize the objective function, and the uniform dose function. To investigate the influence of the problem transform strategy, two variants of the problem transformation, CG-MOEA α and CG-MOEA β , are compared with the proposed method. CG-MOEA α indicates five objective functions determined by the five types of constraints and objectives while CG-MOEA β considers all constraints as objective functions independently. For the fairness of comparison, the way of calculating HV value is the same as CG-MOEA β .
As shown in Table 5, CG-MOEA with the proposed problem transformation performs better than the other two variants on the optimization problem with a large number of objectives, such as nasopharyngeal carcinoma and cervical cancer. For lung cancer and kidney cancer, the performance of CG-MOEA β and CG-MOEA does not have significant differences due to the small number of objectives. Besides, three problem transformation strategies obtain almost the same optimal solutions, even CG-MOEA β perform best than the other two algorithms. A large number of objective functions can significantly increase the difficulty of the optimization problem, and the proposed problem transformation method can handle this issue and also perform well on the problem with a small number of objectives. Table 6 presents median HV and standard deviation values obtained by NSGA-II, MOEA/D, RVEA, MOEA/DVA, MOEA/PSL and the proposed CG-MOEA on five IMRT cancer cases, including nasopharyngeal carcinoma, lung cancer, cervical cancer, esophageal cancer, and kidney cancer, where thirty independent runs are performed for each algorithm on each test instance. The decision variables of the real-world cancer problem cannot be easily divided into two types of convergence-related and diversity-related, hence MOEA/DVA cannot obtain satisfactory results on all cancer cases. MOEA/PSL can effectively generate sparsity optimal solutions and it can provide good solutions on a large-scale multi-objective optimization problem, and MOEA/PSL performs slightly worse than the proposed CG-MOEA. The reference vector-based MOEA/D and RVEA obtain similar experimental results on all cases due to the highly irregular Pareto shape of all cancer cases and the similar environmental mechanism. Besides, NSGA-II achieves the medium results among six algorithms. Furthermore,the proposed CG-MOEA obtains the best HV values on all cases, which shows the superiority of the proposed algorithm comparing with other algorithms. Figures 5, 6, 7, 8 and 9 demonstrate the non-dominated solution sets with median HV value among 30 runs obtained by six multi-objective evolutionary algorithms on nasopharyngeal carcinoma, lung cancer, cervical cancer, esophageal cancer, and kidney cancer. The proposed CG-MOEA obtains the approximate solution sets with good diversity and convergence. MOEA/PSL performs slightly worse than CG-MOEA, especially for lung cancer. Although the uniform distribution solutions have been obtained by the adaptive reference vector-based RVEA on five cases with highly irregular PF shape, the convergence of RVEA is worse than other algorithms, except MOEA/DVA. The diversity of solutions Furthermore, for analyzing the quality of radiotherapy treatment plans from the clinical perspective, some representative compromising radiotherapy treatment plans are selected from each optimal solution set and demonstrated by DVH in Figs. 10, 11, 12, 13 and 14. There are the target preferred plan, the OAR preferred plan and the compromising treatment plan. Obviously, the diversity of MOEA/D and RVEA is significantly worse than other four algorithms, and the treatment plans obtained by these two algorithms cannot be accepted by the radiation physicist due to the worse convergence on the objectives. Although NSGA-II and MOEA/DVA obtain the optimal solutions with good diversity, the convergence of both two algorithms is also unacceptable to treat patients. MOEA/PSL performs best among all four comparing algorithms, but the performance on the most of cases (except kidney cancer case) is far less than that of CG-MOEA. Figures 15 and 16 show the more detailed DVH with all ROIs of esophageal cancer and lung cancer, the dose drops sharply when the minimum dose target is exceeded. Clearly, the performance of the proposed CG-MOEA is significantly better than comparing algorithms on all cases, and its optimal solutions also meet these requirements of radiotherapy treatment.

Comparing the CG-MOEA with commonly used method
To further analyse the effectiveness of CG-MOEA, Figs. 17 and 18 show the DVH of two cancer cases compared against the traditional single objective algorithm that is pure conjugate gradient method. It can be seen that the proposed multi-objective algorithm can find multiple Pareto solutions, moreover, for lung cancer case, the proposed CG-MOEA can find the better treatment plans than traditional CG. Although the optimized dose distribution for esophageal cancer is similar to the traditional single conjugate gradient method. It should be noted that the solution of CG method optimization is the clinical treatment plan obtained by trial-error-trial adjusting the weights of different objectives and constraints based on experience.

Effectiveness of the proposed adaptive conjugate gradient approach
The superiority of the proposed CG-MOEA is mainly owing to the assistance of the gradient information using the conjugate gradient method. For a better combination of the conjugate gradient method and multi-objective evolutionary algorithm, a novel adaptive conjugate gradient method is proposed to accelerate the convergence of the evolutionary multi-objective algorithm. To demonstrate the effectiveness of the proposed adaptive conjugate method, CG-MOEA is compared with CG-MOEA-WO (CG-MOEA without adaptive conjugate gradient) under the same experimental setting. Table 7 gives the mean HV values of the CG-MOEA and its variant version CG-MOEA-WO on nasopharyngeal carcinoma, lung cancer, cervical cancer, esophageal cancer, and kidney cancer with 30 run times. The proposed CG-MOEA performs better than its variant, especially for the slightly complex problem, including lung cancer, esophageal cancer, and kidney cancer. Besides, Fig. 19 shows the number of accelerators determined by the adaptive strategy in each round, it can be seen that the frequency of using the conjugate In the later stage, the significance of population convergence is more important than the diversity and the higher frequency of using the conjugate gradient is needed as shown in Fig. 19.

Parameter sensitivity analysis
The influence of the frequency of employing the adaptive conjugate gradient method is investigated in this section.
As shown in Fig. 20, the performance of the proposed CG-MOEA on cervical cancer and esophageal cancer do not have significant change while the frequency β varies between 0.1 and 0.5. The HV value decreases quickly when the gradient operation time is reduced for nasopharyngeal carcinoma and kidney cancer cases. Although its performance becomes better when the frequency β is set to 0.4, the whole trend of the HV curve is reduced with the increasing of β. Over the influence of β on all cancer cases, the frequency is finally set to 0.1.

Conclusions and future work
In this paper, we have proposed to transform the fluence map optimization problem in IMRT into a multi-objective optimization problem, including the minimization objective, the maximization objective and the uniform dose objective. Moreover, we have proposed a conjugate gradient-assisted multi-objective evolutionary algorithm for solving it, in which the adaptive conjugate gradient method is embedded into the solution generating process for the fast convergence In the future research, more clinical objective or constraint functions will be considered, such as the functions based on biological effectiveness. It is interesting to design an automatical objective/constraint aggregation approach when solving other problems with unknown objectives or constraints. More advanced evolutionary algorithms for complex multi-objective optimization problems in IMRT will be developed.    The values with bold mean the best value among the compared algorithms