- Email: [email protected]
Comparison of simulated annealing algorithms for conformal therapy treatment planning
Im. J. Radiation
Oncology
Biol.
Phys..
Vol.
Copyright Printed
Pergamon
3.1. No. 5. pp.
I091 - 1099.
1995
0 1995 Elscvier Science Inc. in the USA. All rtghts rebervcd OX@301ms $9.50 + .oo
0360-3016(95)02031-6
l
Treatment
Planning
and Plan Evaluation
COMPARISON OF SIMULATED CONFORMAL THERAPY ISAAC I. ROSEN, PH.D.,* MARK LANCER,
KAM
M.D.’
ANNEALING TREATMENT
S. LAM, PH.D.,+ AND
STEVEN
ALGORITHMS PLANNING
FOR
RICHARD G. LANE, PH.D.,+ M. MORRILL, PH.D.+
*Departmentof RadiationPhysics,The University of Texas M. D. AndersonCancerCenter, Houston, TX 77030, ‘Departmentof RadiationTherapy, The University of Texas Medical Branch, Galveston, TX 77550 Purpose: The efficiency of four fast simulated annealing algorithms for optimizing conformal radiation therapy treatment plans was studied and the resulting plans were compared with each other and to optimized conventional plans. Methods and Materials: Four algorithms were selected on the basis of their reported successes in solving other minimization problems: fast simulated annealing with a Cauchy generating function, fast simulated annealing with a Lorentzian generating function, variable step size generalized simulated annealing (VSGSA), and very fast simulated reannealing (VFSR). They were tested on six clinical cases using a multiple beam coplanar conformal treatment technique. Relative beam weights were computed that maximized the minimum tumor dose subject to dose-volume constraints on normal organ doses. Following some initial tuning of the annealing parameters, each algorithm was applied identically to each test case. Optimization tests were run using different random number sequences and different numbers of iterations. Results: The VSGSA algorithm consistently produced the best results. Using long run times, it generated -with the highest minimum tumor dose in five of the six cases. For the short run times, the VSGSA solutions averaged larger minimum tumor doses than those of the other algorithms for all six patients, with increases ranging from 0.4 to 5.9 Gy. For three of the patients, the conformal plan gave a clinically significant increase in the minimum tumor dose over the conventional plan, ranging from 8.2 to 13.0 Gy. In two other cases, there was little difference between the two treatment approaches. For one case, the optimized conventional plan was much better than the conformal plan because the conventional beam arrangement included wedges, which offset the multiple beam advantage of the conformal plans. Conclusions: For equal computing times of both long and short duration, the VSGSA algorithm consistently produced conformal plans that were superior to those produced by the other algorithms. The simple conformal technique used in this study showed a significant potential advantage in the treatment of abdominal tumors. In three of the cases, the conformal plans showed clinically important increases in tumor dose over optimized conventional plans. Treatment
planning,
Optimization,
Simulated
annealing,
Radiation
therapy,
Conformal
radiation
therapy.
Substantial work has been done over the past 20 years on the use of computers to optimize radiation therapy treatment plans. Because no single minimization algorithm is ideal for solving all problems, a variety of algorithms have been adapted based on the individual mathematical characteristics of different treatment strategies. For example, linear programming is ideally suited for objective functions and constraints that can be expressed in linear form; quadratic programming is most efficient for matching computed dose distributions to desired ideal ones when the number of variables is small; and inverse planning algorithms are specifically tailored to computing the nonuniform beam profiles that produce a desired dose
INTRODUCTION
Conformal and computer-controlled radiation therapy techniques offer the promise of improved local control and/or reduced treatment morbidity through the delivery of superior physical dose distributions. The variety of proposed conformal and computer-controlled techniques under study has increased as accelerator technology and computer technology have improved. Computer optimization of treatment plans will be essential for the full use of these sophisticated techniques. Indeed, virtually all of these techniques are too complex for the trial-and-error approach that is usedfor conventional treatment planning.
Reprint requeststo: IsaacI. Rosen,Ph.D., The University of
Acknowledgements-This
work was supported in part by Grant
CA 46634 from the National CancerInstitute, DHHS.
Texas M. D. Anderson Cancer Center, Department of Radiation Physics, Box 94, 1515 Holcombe Blvd., Houston, TX 77030.
Accepted 1091
for publication
7 July 1995.
1092
I. J. Radiation Oncology 0 Biology 0 Physics
pattern. Simulated annealing seemsmost useful for creating conformal therapy plans when the optimization goal and/or constraints are complicated nonanalytic functions. Simulated annealingis a mathematicalminimization processthat bearsa conceptualsimilarity to the physical annealing of solidsin which a substanceis slowly cooled from an initial high temperature-high energy state to a crystalline low temperature-minimumenergy state. Expanding on the ideas of Metropolis et al. (lo), Kirkpatrick (8) initially devised simulated annealing as a means of optimizing the placement of components in computer chips. Since then, simulated annealing methods have been applied to a wide variety of optimization problems in the fields of image processing(4), astronomy (5), financial modeling (7), nuclear medicine (20), and radiation therapy (2,9,11- 16,19,21-23). The term “temperature” is usedfor the simulatedannealing control parameter becausethe searchprocessstarts with a high degreeof randomness,correlating to a high temperature in a physical system, and then progressively lowers the degree of randomness,corresponding to the cooling of a physical system. Interest in simulated annealing minimization as a treatment plan optimization tool arose primarily because of its theoretical ability to find a global minimum in a complicated solution space. Unlike deterministic methods such as downhill simplex, simulated annealing can theoretically avoid being trapped in local minima. The statistical nature of its search allows the possibility of tunnelling through energy barriers, and use of the Metropolis criterion permits climbing out of local minima. Simulated annealing can easily incorporate cost functions that have nonanalytic forms, such astumor control probabilities and normal tissue complication probabilities. Similarly, these complicated functions can easily be used to constrain the solution space. Simulated annealing is also conceptually simple and relatively easy to implement. On the other hand, there are several important difficulties in using simulated annealing. First and foremost is that simulated annealing is not guaranteed to find the global optimum in a finite time. The ability to escape from a local minimum is a necessarycondition for finding the global minimum, but does not ensure it. Szu and Hartley (18) showed that a Cauchy generating function coupled with the proper annealing schedulewill converge to the global minimum in a finite, but not necessarily short, time. However, the treatment planning problem does not fit their mathematical framework becauseof the constraints on the solution space, the minimum constraint being that beam weights be nonnegative. Consequently, we have no way of knowing how close any given solution is to the global minimum. The efficacy of a simulated annealing algorithm can only be ascertained by testing it against other minimization algorithms for the class of problems under consideration. Second, because the simulated annealing method was originally designedfor use with discrete variables, adaptations for usewith continuous variables introduce the prob-
Volume 33, Number 5, 1995
lem of how to quantize the values of the variables. At each iteration, a new set of values for the independent variables is computed on the basis of a displacement vector whose direction and magnitude are selected according to the generating function probability distribution. Small displacement vectors are very inefficient for exploring large areas of the solution space, while large vectors are inefficient in searching the small area near the minimum. Ideally, the vector size should be large initially to allow for exploration of the overall topography of the function. Then it should shrink gradually as the search process proceeds so that the solution will converge to the global minimum. Finally, there are many possible generating functions and annealing schedules, each yielding a different result for a given number of iterations. The efficiency of a given generating function and annealing schedule will, in general, depend on the shape of the solution space to be explored. Ideally, the choice of generating function and annealing schedule should be governed by the characteristics of the problem to be solved. Algorithms are typically tested against mathematical expressions whose exact solutions are known. However, these results may not be indicative of performance with more complicated and less symmetric solution spaces.There is no way at present to predict which generating function and annealing schedule is most efficient for a given class of problems. Therefore, a variety of simulated annealing algorithms have been developed to deal with the various problems described above (1, 3, 6, 17, 18). The aims of this study were (a) to compare the efficiency of four different algorithms for producing optimized conformal treatment plans and (b) to compare optimized conformal treatment plans to optimized conventional plans for a seriesof abdominal tumors. The selected algorithms had been previously shown to be efficient at solving minimization problems in other fields. Two of them use global generating functions linked to global annealing schedules, and two use generating functions and annealing schedules that dynamically change as search results are evaluated. They were all implemented in our treatment planning system and then compared by optimizing conformal treatment plans for abdominal cancer cases selected from our clinic. It was assumed that the algorithms would all converge to the same or comparable solutions given enough time (number of iterations). We were interested in which algorithm would give the best results in about 30 min of computing, a time scale that is clinically useful. Because the shape of the solution space is usually unknown for conformal treatment techniques, one would expect algorithms that dynamically modify the generating function and/or annealing schedule to be the most efficient. The algorithms were also used to optimize conventional wedged three-field plans to evaluate whether the results of algorithm comparison were dependent on the number of optimization variables. Finally, the best conventional plan for each patient was
Simulatedannealingalgorithms for confotmal therapy treatmentplanningl I. I. Table 1. Organ dose-volume constraints for the six patients in this study Volume 5 dose limit
Organ
66% 66% 60% 75% 100%
Left kidney Right kidney Liver Small bowel Spinal cord
5 5 5 5 I
20 20 30 45 45
Gy Gy Gy Gy Gy
compared to the best conformal plan to see if the conforma1 technique using multiple static beams offered any dose distribution advantages over a well-optimized conventional beam arrangement. METHODS
AND
MATERIALS
Optimized treatment plans were generated for six patients with tumors in the abdomen selected from our clinic. Three patients had gastric tumors, two had pancreatic tumors, and one had a tumor in the colon. For each patient, a three-dimensional computer model of the normal anatomy and tumor volume was constructed from outlines drawn by a physician on the individual crosssections of the computed tomography study. Each normal organ was represented by dose constraint points randomly distributed throughout its volume. Because both dose computation time and optimization time increase as the number of dose constraint points increases, there is an incentive to minimize the number of points used. On the other hand, the number must be sufficiently large to accurately represent the distribution of dose in the organ. For these tests, the number of dose constraint points within each organ was determined by both the size of the organ and its proximity to the tumor volume. A minimum of 100 dose constraint points was required for each organ, and roughly the same number of points was used to represent the same organ in all cases. Point densities were typically about 0.5, 1, 1.5, and 5 points per cm’ for liver, small bowel, kidney, and spinal cord, respectively. For each patient the same dose constraint points were used with all of the annealing algorithms. The optimized plans maximized the minimum tumor dose subject to dose-volume constraints on the doses to the normal structures. For the treatment techniques under consideration, the minimum target dose occurs on the boundary of the target volume. Therefore, the target volume was represented in each CT slice by many dose points placed at 1 to 2 cm intervals along the perimeter of the target outline and a few randomly distributed within the target outline. Only solutions that met the normal tissue dose-volume constraints specified in Table 1 were considered feasible. These dose-volume constraints are equivalent to points lying on cumulative dose-volume histograms. For each set of beam weights, an approximate dose-volume histo-
ROSEN et al.
1093
gram was computed for each organ from the doses to the dose constraint points and then compared to the corresponding dose-volume constraint. Solutions that violated dose-volume constraints (nonfeasible solutions) were renormalized and used in place of the original nonfeasible solutions. Simulated annealing method Simulated annealing is an iterative process that searches for the value of the vector x that minimizes the value of the cost function F(x). The components of x are the independent variables x,, x2, . . . , x,,. For external beam treatment plan optimization, these independent variables are relative beam weights proportional to the doses delivered from beams or parts of beams. The cost function assigns a figure of merit to the plan resulting from a selected set of beam weights. Because simulated annealing is a minimization procedure, we maximized the tumor dose by choosing the negative of the minimum tumor dose as the cost function. Each iteration of the simulated annealing process has four basic steps and two relevant solutions. The process starts with an initial solution that can be obtained by any means, but must meet all constraints placed on the solution space. For treatment plan optimization problems, the initial solution is usually all beams with zero weight because this solution immediately meets all of the dose constraints on normal structures. This initial solution serves as the “current solution” for starting the simulated annealing search. The current solution is considered to be the point in the solution space closest to the global minimum. At each point in the process, the algorithm generates a “trial solution” and decides whether the trial solution should become the current solution. Ideally, the current solution will be the best solution when the simulated annealing process is terminated. However, because the Metropolis criterion occasionally moves the current solution to a higher cost function value and because the process is terminated after a finite number of iterations, the current solution at the end of a simulated annealing run may not have the lowest cost of any solution encountered during the run. Therefore, the best solution found in any iteration is usually retained and used as the final result. First, a trial solution x’ is generated from the current solution. A generating function is used to select the magnitude of a displacement vector Ax from the current solution x in the n-dimensional variable space. The direction of Ax is selected by randomly distributing its components among the independent variables. The components of x’ are then given by 4 ’ =x,
+ Ax;, (i = 1, . . . . n).
(Eq. 1)
Each of the four algorithms under consideration uses a different generating function. In general, the probability
1094
I. J. Radiation
Oncology
0 Biology
l Physics
of choosing a displacement vector of a given length decreases as the temperature decreases. The range of search is initially large with the hope that the well containing the global minimum will be discovered. Then, the search area decreases in an effort to find the minimum of the well containing the current solution. However, at every iteration there is still a finite probability of escaping the current well. Next, the trial solution is tested against constraint limits, and the value of the cost function is computed. If the trial solution meets all constraints, the algorithm continues on to the next step. If not, the trial solution is nonfeasible. It can be rejected immediately or it can be renormalized so that it does meet the constraints. Renormalization involves uniformly reducing all of the beam weights by the minimum amount needed to meet all of the constraint requirements. Renormalizing nonfeasible trial solutions is generally more efficient than discarding them. When nontrivial constraints are used, the minimum value of the cost function will fall at one of the constraint boundaries. Therefore, as the search approaches the minimum, the number of nonfeasible trial solutions increases rapidly. If these nonfeasible trial solutions are discarded, the efficiency of the algorithm decreases proportionally. Although renormalizing a nonfeasible solution is generally more computationally expensive than simply discarding it and selecting another, renormalization improves overall efficiency by providing a feasible trial solution at every iteration. Given a feasible trial solution, the algorithm must decide whether this trial solution represents an improvement over the current solution. If the cost function value for the trial solution is smaller than that for the current solution, then the trial solution is clearly superior and it becomes the new current solution. However, if the cost value for the trial solution is larger than that for the current solution, the trial solution is not immediately discarded. Rather, it is retained as the current solution with a probability given by the Metropolis criterion,
p = e -[m’k-w)l T
’
F(x’)
> F(x)
(Eq. 2)
where T is a “temperature” parameter. The probability of retaining a trial solution with higher cost decreases as the temperature of the system decreases and as the difference in cost values increases. Finally, the temperature of the system is reduced according to an annealing schedule that specifies a sequence of declining values for the temperature parameter as the optimization process proceeds. This annealing schedule requires an initial temperature value, a function describing the rate of temperature reduction, a final temperature, and the number of iterations to be performed at each temperature. A single temperature schedule may be used for all variables or each variable may have its own annealing schedule.
Volume
33, Number
5, 1995
Fast simulated annealing-Cauchy distribution generating function One of the earliest of the fast simulated annealing algorithms uses a generating function with a probability density given by a Cauchy distribution (18). The probability of obtaining a displacement vector of magnitude ) Ax) in an n-dimensional search space is given by
gOAx
=
W(k)
[(IA ()* + W(k)2]“‘+
“‘*
(Eq. 3)
where W(k) has been defined by Mageras and Mohan (9) for treatment plan optimization as W(k) =
W, (1 + WR)
(Eq. 4)
and k is the iteration number. The Cauchy distribution has a Gaussian-like peak and Lorentzian wings, allowing occasional long jumps during the search process. The parameter W,, defines the width of the Cauchy distribution, and the parameter R determines the rate at which the Cauchy distribution diminishes in size during the annealing process. Because there is no quick algorithm for calculating an n-dimensional Cauchy random number distribution, the magnitude of the step size for any generated random number was obtained by interpolation from a precalculated lookup table. The direction of the displacement vector was selected by randomly distributing its components among the independent variables. A sequence of n random numbers (u, 7 u2, . ..1 u,,) uniformly distributed in the interval [- 1, 1) was generated. These numbers define the components of a vector u in the n-dimensional space whose magnitude is
IllI = [(UT+ u; + * * - + Uf)]“2. 0% 5) The direction vector u, displacement vector Ax, and the current solution were then used to compute a new trial solution: X: = xi +
IAXI
.
0%. 6)
Any negative beam weights in the trial solution were set to zero. The temperature annealing schedule for this algorithm was given by T(k) =
” (1 + WI?)
@q. 7)
where T,, is the initial temperature and parameter R from Eq. 4 also determines the rate at which the temperature
Simulated
annealing
algorithms
for conformal
decreases. The temperature value was updated at every iteration. Based on tuning with a single case, we selected values for the control parameters of W, = 1, T,, = 0.2, and R = 500.
Fast simulated annealing-lorentzian generating function
distribution
A Lorentzian probability distribution in one-dimensional form is similar in shape to the Cauchy distribution, but is much simpler to compute. It also allows occasional large step sizes during the search process. The Lorentzian probability of obtaining a displacement vector of magnitude 1Ax 1 is given by
g(IAxl) =
W(k)* [( I Ax I >’ + W(k121 .
treatment
planning AXj
(W. 9)
where W(k) is the width of the Lorentzian distribution for the kth iteration. It has the same form as the width of the Cauchy distribution given by Eq. 4. The directions of displacement vectors and new trial solutions were generated in the same way as for the Cauchy distribution. The cooling schedule used here was the same as that of the fast annealing algorithm with the Cauchy distribution generating function (Eq. 7). Again, the values for the control parameters were determined by tuning the algorithm with one test case: W, = 10, T, = 0.2, and R = 500. The temperature value was also updated at every iteration.
Variable step size generalized simulated annealing (VSGSA) The variable step size generalized simulated annealing (VSGSA) algorithm is a modification by Sutter and Kalivas (17) of the generalized simulated annealing algorithm originally developed by Bohachevsky et al. (1). In this algorithm, the components of the displacement vector are not computed according to a theoretical probability distribution. Rather, they are calculated to maintain a relatively high efficiency of finding improved solutions as measured by the acceptance ratio, which is defined as the number of accepted trial solutions divided by the number of trial solutions evaluated within a specified block of iterations. An accepted trial solution is one that becomes the current solution by virtue of either a lower cost value or the Metropolis criterion. The change in beam weight for ith beam, AX, is calculated from a random number u, uniformly distributed in the interval [ - 1, 1) using the expression:
0 I. I. ROSEN et al. = hjUi(Ximax
-
1095 Ximin)-
(Eq.
10)
The minimum value for each beam weight, X,,inr was set to 0. The value of ximaxfor each beam was determined by finding the maximum weight that would generate a feasible solution when all other beam weights were set to zero. The values of Ximinand Ximaxwere determined in advance and not changed. The parameter hj controls the overall magnitude of step sizes for the jth block of iterations and it changes according to the acceptance ratio. If the acceptance ratio is greater than 0.9 for a block of iterations, then the current solution is considered to be far from the global minimum and the step size factor is increased according to
hj+, = mhj,
(Eq. 8)
Given a random number r uniformly distributed in the interval [0, I), the magnitude of a random displacement vector Ax is obtained from
1Ax ( = W(k) etan(O.%rr),
therapy
0% 11)
where h, is the old value for the jth block of iterations and m is a factor greater than 1, whose value is set by the user. If the acceptance ratio is less than 0.2, then the current solution is considered to be close to the global minimum and the step size factor is reduced according to
hj+, = h,lm.
0%. 12)
The initial value of h was chosen to be 0.20; the value of m was set at 1.05; and the acceptance ratio was calculated every 100 iterations. The probability of accepting a trial solution with a higher cost value is given by the function p = e~lcc,-c”)l~-tc,~c”,l 9
(Eq. 13)
where C, is the cost value of the trial solution, C,. is the cost value of the current solution, and C, is the global optimum value estimated by the user. The initial value for the temperature, T, was determined so that the probability P was between 0.2 and 0.7. Thereafter, the value of T was adjusted whenever a detrimental step occurred after the step size had been changed.
Very fast simulated
reannealing
(VFSR)
In the very fast simulated reannealing (VFSR) algorithm developed by Ingber (6), each independent variable has its own exponentially decreasing temperature, I;, given by
T;(k) = Ti(O)exp(-cik”“),
(Eq. 14)
where k is the iteration number, n is the number of independent variables, and ci is a user-defined parameter for the ith variable. The temperature of each variable is periodically resealed, depending on the sensitivity of the cost function to that variable. The purpose of resealing the temperature for each variable independently is so that the
IO96
I. J. Radiation Oncology 0 Biology 0 Physics
Volume 33, Number 5. 1995
generating function will give larger displacement vectors when searching relatively flat and insensitive dimensions and smaller displacement vectors in more rapidly varying dimensions. The sensitivity, s,, of the cost function to a variable is determined by sampling the pseudogradient of the steepestminimum found using numerically computed partial derivatives of the cost function with respect to that variable:
II
dF s; = . 8%
@I. 13 4
“Patient
Then, in each reannealing step, a new temperature T,‘(k) for each variable is computed, which is proportional to the old one according to the sensitivity of the variable,
where s,,, is the largest value of the sj computed during the reannealing step. This resealing of the temperature, or reannealing, can often improve the overall efficiency of the algorithm. However, any reduction in the number of iterations required for convergence comes at someexpense in increased computation time per iteration. The generating function for VFSR is temperature dependent. The change in weight for the ith beam is calculated from (Eq. 17) yi = sign(u, - 0.5)Ti[(L
+ -$‘“”
- 11
(Q.
18)
where ui is a random number in the interval [- 1,l) and the value of yi lies in the interval [ - 1, 1). Repeated values of Ax, are generated until a new value for the ith variable is obtained that is bounded by the interval [Ximinrx;~J. The values for X;,in and x,,,, are determined as for the VSGSA algorithm. Calculation of optimized treatment plans For the conformal plans, 36 coplanar beamsof 15 MV x-rays were deployed around the patient at 10” intervals, giving 36 variables for optimization. The conventional plans used anterior and lateral beams of 15 MV x-rays. Each beam direction included an open beam and two beams with 60” wedges oriented in opposite directions for a total of nine optimization variables. All beamswere shapedto conform to the projections of the tumor volume
’ VAXstation 4000 Model 60, Digital EquipmentCorporation, Maynard,
MA.
5
6
Fig. 1. Minimum target doses for the optimized conformal treatment plans generated using the different simulated annealing algorithms with long run times. The optimization objective was to maximize the minimum target dose subject to normal tissue dose-volume limits. The VSFR algorithm was run for 80,000 iterations and the others were run for 100,000 iterations each.
with a margin for beam penumbra. Doses were computed using a scatter-air ratio sector summation algorithm with a Batho power-law correction for tissue heterogeneities. Short and long optimization runs were performed for each patient. The short runs consisted of 8,000 iterations for the VFSR algorithm and 10,000 for each of the others. Preliminary tests showed that the VFSR algorithm was approximately 20% slower per iteration than the others. Therefore, the number of iterations for that algorithm was reduced so that elapsed times would be comparable for all four algorithms. The long runs consisted of 80,000 iterations for the VFSR algorithm and 100,000 iterations for each of the others. The long runs were performed to give an indication of how rapidly the short runs were converging to the best possible solution. Computation times’ for the short runs were approximately 20-30 CPU min for the conformal plans and 5-8 CPU min for the conventional plans. Becauseof the statistical nature of simulated annealing, running
the same treatment-planning
problem
with
two
different random number sequenceswill produce two different solutions. The longer the run time, the smaller will be the variation between the different solutions. Therefore, the short runs for the conformal plans were repeated five times for each patient using different random number sequences.The long runs were performed only once. Random numbers were obtained from the linear congruent random number generator in the VAXstation system software. RESULTS Results of the long optimization runs for the conformal plans are shown in Fig. 1. The minimum target dose is
Simulated
annealing
algorithms
for conformal
shown for each patient and each algorithm. Overall, the VSGSA algorithm gave the best results. For five of the six patients, the VSGSA plan had the highest minimum target dose, with improvements ranging from 0.1 to 3.8 Gy. Averaged over the six patients, the VSGSA plans were 0.8 Gy higher than the VFSR plans, 1.7 Gy higher than the Cauchy plans, and 2.1 Gy higher than the Lorentzian plans. For the conformal beam arrangements, the results of the five short runs using different random number sequences were averaged for each patient and for each simulated annealing algorithm. In Fig. 2, the average values and the corresponding standard deviations are shown. For each of the six patients, the VSGSA algorithm produced the highest average cost value: improvements over the VFSR plans ranged from 1.O to 3.3 Gy; improvements over the Cauchy plans ranged from 0.5 to 4.7 Gy; and improvements over the Lorentzian plans ranged from 0.4 to 5.9 Gy. The results of the short runs are compared to the results of the long run times for the VSGSA algorithm in Fig. 3. In five of the six patients, the VSGSA short run average was within 0.8 Gy of the long run solution. In the other case (Patient 5), the VSGSA short run average was 1.1 Gy less than the long run solution. For the conventional plans, the results produced by the four algorithms typically differed by less than 1% and the maximum difference was 1.5%. The variations produced by using different random number sequenceswere typically less than 0.5%. Although the differences were small, the VSGSA algorithm again gave the best average result for each case. Figure 4 shows the average computation times in CPU minutes and the corresponding standard deviations for the four algorithms for calculating optimized conformal plans using the short runs. The CPU times for the short runs
therapy
treatment
planning
l I. I. ROSEN
Bi ::: jjj i:: ii: I// g// ii: ii: ii_ :i: ::: ::. j/: ::: i:: ii: ::::: /Iii. -
1
2
2
3
Patient4
5
6
Fig. 2. Average minimum target doses for the optimized conforma1 treatment plans generated using short run times of the different simulated annealing algorithms. The vertical bars represent the averages and the error bars represent the standard deviations for five runs using different random number sequences. The optimization objective was to maximize the minimum target dose subject to normal tissue dose-volume limits. The VFSR algorithm was run for 8,000 iterations and the others were run for 10,000 iterations each.
3
Patient4
5
6
Fig. 3. Comparison of short run time results with long run time results for the VSGSA algorithm.
are quite similar for all algorithms and cases. With the exception of Patient 2, the time per iteration for the Cauthy, Lorentzian, and VSGSA algorithms were all within 4% of each other. In general, the compute times varied only slightly with different random number seeds. For Patient 6, the VFSR had an unusually large standard deviation. For this patient, the number of times the algorithm did reannealing varied vastly with different sequencesof random numbers. In three of the cases, the conformal plan was significantly better than the optimized conventional plan. The short VSGSA runs for the conformal and conventional techniques are compared in Fig. 5. For Patients 2, 3, and 5, the conformal plan increased the minimum tumor dose by 8.2 to 13.0 Gy over the conventional treatment. For Patients 1 and 4, there were only small differences (< 0.8 Gy) between the two treatment methods. For Patient 4, this was probably because the conventional plan dose was already very high (66.5 Gy). In one patient (Patient 6), the conventional plan was better than the conformal plan. This patient had a very large tumor, and the use of wedges in the conventional plan was a significant advantage that the conformal beam arrangement could not overcome, even by the useof many more beams. DISCUSSION
1
et a[.
AND CONCLUSIONS
One of the advantages of simulated annealing over other optimization methods is that it is easy to implement complex cost and constraint functions. However, the use of absolute constraints creates discontinuities in the solution space, which in turn, substantially increasesthe difficulty of solving the problem. For nontrivial constraints, the global minimum will always lie at the boundary of a constraint. As the solution approaches the boundary of the valid solution region, most increases in beam weights will generate solutions that violate the dose-volume con-
straints and, therefore, the probability of generating a feasible solution decreases. Renormalization of beam weights whenever organ dosesexceed dose-volume con-
I. J. RadiationOncology0 Biologyl Physics
1098
3s z .3 g
30
,g
20
!3
10
5
150 1
25
5
2
3
4
Patient
5
6
Fig. 4. Average computationtimes(CPU minutes)and standard deviationsfor calculating conformal plansusing short runs of the four algorithms-8,000 iterationsfor the VFSR algorithm and 10.000 for each of the others.
straints improves efficiency. Unfortunately, this modification violates the convergence criteria of simulated annealing algorithms and could causesolutions to be trapped in local minima. In the VSGSA algorithm, the calculated acceptance ratio gives a measure of the closenessof the current solution to the global minimum. With renormalization, however, the acceptance ratio is no longer related to the point at which it was evaluated. Nevertheless, we did find a valuable increase in overall efficiency when renormalization was implemented. An alternative method for constraining solutions is to add a penalty value to the cost function whenever solutions fall outside some defined limits (11). Such “soft” constraints may be used in place of or in conjunction with the “hard” limits described above. We have avoided using soft constraints because they represent compromisesbetween the absolute constraints and the objective, and the most appropriate compromise for any given problem will generally depend on the preferences of the physician and the details of the patient’s case. Rather than attempting to incorporate the physician’s preferences into the cost function, we present the best possible solution for the given hard constraints and require the physician to explicitly identify what constraints to relax if a higher tumor dose or other objective is desired. The use of soft constraints in place of hard constraints simplifies the simulated annealing optimization algorithm and probably gives better convergence becausethere are no discontinuities in the solution space. If soft constraints are used in addition to the hard ones, it is not clear that there is any additional benefit in terms of algorithmic efficiency. Hard constraints clearly separate solutions that are acceptable from those that are not. Soft constraints provide a smooth transition between those two types of solutions and in doing so become part of the definition of the optimum solution by the way in which the penalties are imposed. We found that the VSGSA algorithm consistently produced better results than the other algorithms. No simulated annealing technique can produce a guaranteed opti-
Volume33,Number5, 1995 ma1 solution with a single finite run. However, the VSGSA method performed best using both long and short run times for both conformal and conventional techniques. For the long runs, the VSGSA algorithm gave the highest minimum tumor dose for five of the six cases. For the short runs, the VSGSA solutions averaged over five runs gave minimum tumor doses that were larger by 1.8 to 5.9 Gy in 67% of the tests. An increase in tumor dose of 1.8 Gy is approximately equal to an increase of one treatment fraction and, therefore, could be clinically significant. Although only one cost function and constraint type were used in this study, we expect that the VSGSA algorithm will still be the most efficient with other cost and constraint functions. Because the performance of a given simulated annealing algorithm is dependent on the nature of the problem to be solved, it is appropriate to ask if the results presented would be applicable to other conformal therapy techniques. In other words, would the VSGSA algorithm still perform best for optimizing treatments with small numbers of noncoplanar or intensity-modulated beams?A definitive answer would require testing the various algorithms using those treatment techniques. However, for the objective function and constraints under consideration, the VSGSA algorithm performed best for 9 beam variables (conventional technique) and 36 beam variables (conformal technique) with the superiority of the VSGSA algorithm over the others increasing with the number of beams. Based on these observations, we expect that the VSGSA algorithm would still be most efficient when applied to the problem of planning with intensity-modulated beams. Conformal therapy using multiple static shapedcoplanar beamsappearsto offer clinical advantagesin the treatment of abdominal tumors. In half the patients studied the minimum tumor dose increased by more than 8 Gy for the samenormal tissue dose-volume constraints. It should be remembered that the conformal plans were compared to optimized conventional beam arrangements,not to standard plans. Although it was not evaluated as part of this study, the gain should be even greater when comparedto standard
ii: ::: ::: ii: /> /i: ::: .j/ .:: ii: :lj :i: i:: ::: ii: :;/ 1
::: ii: ::: ::: :// jj/ j// ii: ii: hsi -I
3
Patient4
5
6
Fig. 5. Comparisonof conformalplansand conventionalplans optimized using VSGSA.
Simulated
annealing
algorithms
for conformal
plans. In the one patient for whom the conventional plan was better, the large tumor size gave a substantial advantage to the use of wedges. Including wedges in the conformal technique is simple and should provide some additional increase in tumor dose at the expense of increased computer
therapy
treatment
planning
0 I. 1. ROSEN er al.
1099
time in planning and possibly increased treatment delivery time. The routine use of this conformal technique is no more difficult or complex than many conventional treatments but will require a computer-controlled accelerator and multileaf collimator for efficient automated delivery of treatments.
REFERENCES 1. Bohachevsky, I. 0.; Johnson, M. E.; Stein, M. L. Generalized simulated annealing for function optimization. Technometrics 28:209-217; 1986. 2. Bortfeld, T.; Schlegel, W. Optimization of beam orientations in radiation therapy: Some theoretical considerations. Phys. Med. Biol. 38:291-304; 1993. 3. Corana, A.; Marchesi, M.; Martini, C.; Ridella, S. Minimizing multimodal functions of continuous variables with the “simulated annealing” algorithm. ACM Trans. Math. Soft. 13:262-280; 1987. 4. Geman, S.; Geman, D. Stochastic relations, Gibbs distributions, and the Bayesian restoration of images. IEEE Trans. Patt. Anal. Mach. Intel. PAM1 6:721-741; 1984. 5. Groisman, G.; Parker, J. R. Computer-assisted photometry using simulated annealing. Comp. Phys. 7:87-96; 1993. 6. Ingber, L. Very fast simulated reannealing. Math. Comput. Modelling 12:967-973; 1989. 7. Ingber, L. Statistical-mechanical aids. to calculating termstructure models. Phys. Rev. [A] 42:7057-7064; 1990. 8. Kirkpatrick, S.; Gelatt, C. D.; Vecchi, M. P. Optimization by simulated annealing. Science 220:671-680; 1983. 9. Mageras, G. S.; Mohan, R. Application of fast simulated
annealing to optimization of conformal radiation treatments. Med. Phys. 20:639-647; 1993. 10. Metropolis, N.; Rosenbluth, A. W.; Rosenbluth,
M. N.;
14.
15. 16. 17.
18. 19. 20. 21.
Teller, A. H.; Teller, E. Equation of state calculationsby fast computing machines. J. Chem. Phys. 21:1087-1092; 1953. 11. Mohan, R.; Mageras, G. S.; Baldwin, B.; Brewster, L. J.; Kutcher, G. J.; Leibel, S.; Burman, C. M.; Ling, C. C.; Fuks, Z. Clinically relevant optimization of 3-D conformal treatments. Med. Phys. 19:933-944; 1992. 12. Merrill, S. M.; Lane, R. G.; Jacobson, G.; Rosen, I. I. Treatment planning optimization using constrained simulated annealing. Phys. Med. Biol. 36:1341-1361; 1991. 13. Morrill, S. M.; Lane, R. G.; Rosen, I. I. Constrained simu-
22. 23.
lated annealing for optimized radiation therapy treatment planning. Comp. Methods Prog. Biomed. 33: 135- 144; 1990. Niemierko, A. Random search algorithm (RONSC) for optimization of radiation therapy with both physical and biological end points and constraints. Int. J. Radiat. Oncol. Biol. Phys. 23:89-98; 1992. Sloboda, R. S. Optimization of brachytherapy dose distributions by simulated annealing. Med. Phys. 19:955-964; 1992. Sloboda, R. S.; Pearcey, R. G.; Gillan, S. J. Optimized low dose rate pellet configurations for intravaginal brachytherapy. Int. J. Radiat. Oncol. Biol. Phys. 26:499-511; 1993. Sutter, J. M.; Kalivas, J. H. Convergence of generalized simulated annealing with variable step size with application toward parameter estimations of linear and nonlinear models. Anal. Chem. 63:2383-2386; 1991. Szu, H.; Hartley, R. Fast simulated annealing. Phys. Len. [A] 122:157-162; 1987. Webb, S. Optimisation of conformal radiotherapy dose distributions by simulated annealing. Phys. Med. Biol. 34:1349-1370; 1989. Webb, S. SPECT reconstruction by simulated annealing. Phys. Med. Biol. 34:259-281; 1989. Webb, S. Optimization by simulated annealing of threedimensional conformal treatment planning for radiation fields defined by a multileaf collimator. Phys. Med. Biol. 36:1201-1226; 1991. Webb, S. Optimization of conformal radiotherapy dose distributions by simulated annealing: 2. Inclusion of scatter in the 2D technique. Phys. Med. Biol. 36: 1227-1237; 1991. Webb, S. Optimization by simulated annealing of threedimensional, conformal treatment planning for radiation fields defined by a multileaf collimator: II. Inclusion of two-dimensional modulation of the x-ray intensity. Phys. Med. Biol. 37:1689-1704; 1992.
Oncology
Biol.
Phys..
Vol.
Copyright Printed
Pergamon
3.1. No. 5. pp.
I091 - 1099.
1995
0 1995 Elscvier Science Inc. in the USA. All rtghts rebervcd OX@301ms $9.50 + .oo
0360-3016(95)02031-6
l
Treatment
Planning
and Plan Evaluation
COMPARISON OF SIMULATED CONFORMAL THERAPY ISAAC I. ROSEN, PH.D.,* MARK LANCER,
KAM
M.D.’
ANNEALING TREATMENT
S. LAM, PH.D.,+ AND
STEVEN
ALGORITHMS PLANNING
FOR
RICHARD G. LANE, PH.D.,+ M. MORRILL, PH.D.+
*Departmentof RadiationPhysics,The University of Texas M. D. AndersonCancerCenter, Houston, TX 77030, ‘Departmentof RadiationTherapy, The University of Texas Medical Branch, Galveston, TX 77550 Purpose: The efficiency of four fast simulated annealing algorithms for optimizing conformal radiation therapy treatment plans was studied and the resulting plans were compared with each other and to optimized conventional plans. Methods and Materials: Four algorithms were selected on the basis of their reported successes in solving other minimization problems: fast simulated annealing with a Cauchy generating function, fast simulated annealing with a Lorentzian generating function, variable step size generalized simulated annealing (VSGSA), and very fast simulated reannealing (VFSR). They were tested on six clinical cases using a multiple beam coplanar conformal treatment technique. Relative beam weights were computed that maximized the minimum tumor dose subject to dose-volume constraints on normal organ doses. Following some initial tuning of the annealing parameters, each algorithm was applied identically to each test case. Optimization tests were run using different random number sequences and different numbers of iterations. Results: The VSGSA algorithm consistently produced the best results. Using long run times, it generated -with the highest minimum tumor dose in five of the six cases. For the short run times, the VSGSA solutions averaged larger minimum tumor doses than those of the other algorithms for all six patients, with increases ranging from 0.4 to 5.9 Gy. For three of the patients, the conformal plan gave a clinically significant increase in the minimum tumor dose over the conventional plan, ranging from 8.2 to 13.0 Gy. In two other cases, there was little difference between the two treatment approaches. For one case, the optimized conventional plan was much better than the conformal plan because the conventional beam arrangement included wedges, which offset the multiple beam advantage of the conformal plans. Conclusions: For equal computing times of both long and short duration, the VSGSA algorithm consistently produced conformal plans that were superior to those produced by the other algorithms. The simple conformal technique used in this study showed a significant potential advantage in the treatment of abdominal tumors. In three of the cases, the conformal plans showed clinically important increases in tumor dose over optimized conventional plans. Treatment
planning,
Optimization,
Simulated
annealing,
Radiation
therapy,
Conformal
radiation
therapy.
Substantial work has been done over the past 20 years on the use of computers to optimize radiation therapy treatment plans. Because no single minimization algorithm is ideal for solving all problems, a variety of algorithms have been adapted based on the individual mathematical characteristics of different treatment strategies. For example, linear programming is ideally suited for objective functions and constraints that can be expressed in linear form; quadratic programming is most efficient for matching computed dose distributions to desired ideal ones when the number of variables is small; and inverse planning algorithms are specifically tailored to computing the nonuniform beam profiles that produce a desired dose
INTRODUCTION
Conformal and computer-controlled radiation therapy techniques offer the promise of improved local control and/or reduced treatment morbidity through the delivery of superior physical dose distributions. The variety of proposed conformal and computer-controlled techniques under study has increased as accelerator technology and computer technology have improved. Computer optimization of treatment plans will be essential for the full use of these sophisticated techniques. Indeed, virtually all of these techniques are too complex for the trial-and-error approach that is usedfor conventional treatment planning.
Reprint requeststo: IsaacI. Rosen,Ph.D., The University of
Acknowledgements-This
work was supported in part by Grant
CA 46634 from the National CancerInstitute, DHHS.
Texas M. D. Anderson Cancer Center, Department of Radiation Physics, Box 94, 1515 Holcombe Blvd., Houston, TX 77030.
Accepted 1091
for publication
7 July 1995.
1092
I. J. Radiation Oncology 0 Biology 0 Physics
pattern. Simulated annealing seemsmost useful for creating conformal therapy plans when the optimization goal and/or constraints are complicated nonanalytic functions. Simulated annealingis a mathematicalminimization processthat bearsa conceptualsimilarity to the physical annealing of solidsin which a substanceis slowly cooled from an initial high temperature-high energy state to a crystalline low temperature-minimumenergy state. Expanding on the ideas of Metropolis et al. (lo), Kirkpatrick (8) initially devised simulated annealing as a means of optimizing the placement of components in computer chips. Since then, simulated annealing methods have been applied to a wide variety of optimization problems in the fields of image processing(4), astronomy (5), financial modeling (7), nuclear medicine (20), and radiation therapy (2,9,11- 16,19,21-23). The term “temperature” is usedfor the simulatedannealing control parameter becausethe searchprocessstarts with a high degreeof randomness,correlating to a high temperature in a physical system, and then progressively lowers the degree of randomness,corresponding to the cooling of a physical system. Interest in simulated annealing minimization as a treatment plan optimization tool arose primarily because of its theoretical ability to find a global minimum in a complicated solution space. Unlike deterministic methods such as downhill simplex, simulated annealing can theoretically avoid being trapped in local minima. The statistical nature of its search allows the possibility of tunnelling through energy barriers, and use of the Metropolis criterion permits climbing out of local minima. Simulated annealing can easily incorporate cost functions that have nonanalytic forms, such astumor control probabilities and normal tissue complication probabilities. Similarly, these complicated functions can easily be used to constrain the solution space. Simulated annealing is also conceptually simple and relatively easy to implement. On the other hand, there are several important difficulties in using simulated annealing. First and foremost is that simulated annealing is not guaranteed to find the global optimum in a finite time. The ability to escape from a local minimum is a necessarycondition for finding the global minimum, but does not ensure it. Szu and Hartley (18) showed that a Cauchy generating function coupled with the proper annealing schedulewill converge to the global minimum in a finite, but not necessarily short, time. However, the treatment planning problem does not fit their mathematical framework becauseof the constraints on the solution space, the minimum constraint being that beam weights be nonnegative. Consequently, we have no way of knowing how close any given solution is to the global minimum. The efficacy of a simulated annealing algorithm can only be ascertained by testing it against other minimization algorithms for the class of problems under consideration. Second, because the simulated annealing method was originally designedfor use with discrete variables, adaptations for usewith continuous variables introduce the prob-
Volume 33, Number 5, 1995
lem of how to quantize the values of the variables. At each iteration, a new set of values for the independent variables is computed on the basis of a displacement vector whose direction and magnitude are selected according to the generating function probability distribution. Small displacement vectors are very inefficient for exploring large areas of the solution space, while large vectors are inefficient in searching the small area near the minimum. Ideally, the vector size should be large initially to allow for exploration of the overall topography of the function. Then it should shrink gradually as the search process proceeds so that the solution will converge to the global minimum. Finally, there are many possible generating functions and annealing schedules, each yielding a different result for a given number of iterations. The efficiency of a given generating function and annealing schedule will, in general, depend on the shape of the solution space to be explored. Ideally, the choice of generating function and annealing schedule should be governed by the characteristics of the problem to be solved. Algorithms are typically tested against mathematical expressions whose exact solutions are known. However, these results may not be indicative of performance with more complicated and less symmetric solution spaces.There is no way at present to predict which generating function and annealing schedule is most efficient for a given class of problems. Therefore, a variety of simulated annealing algorithms have been developed to deal with the various problems described above (1, 3, 6, 17, 18). The aims of this study were (a) to compare the efficiency of four different algorithms for producing optimized conformal treatment plans and (b) to compare optimized conformal treatment plans to optimized conventional plans for a seriesof abdominal tumors. The selected algorithms had been previously shown to be efficient at solving minimization problems in other fields. Two of them use global generating functions linked to global annealing schedules, and two use generating functions and annealing schedules that dynamically change as search results are evaluated. They were all implemented in our treatment planning system and then compared by optimizing conformal treatment plans for abdominal cancer cases selected from our clinic. It was assumed that the algorithms would all converge to the same or comparable solutions given enough time (number of iterations). We were interested in which algorithm would give the best results in about 30 min of computing, a time scale that is clinically useful. Because the shape of the solution space is usually unknown for conformal treatment techniques, one would expect algorithms that dynamically modify the generating function and/or annealing schedule to be the most efficient. The algorithms were also used to optimize conventional wedged three-field plans to evaluate whether the results of algorithm comparison were dependent on the number of optimization variables. Finally, the best conventional plan for each patient was
Simulatedannealingalgorithms for confotmal therapy treatmentplanningl I. I. Table 1. Organ dose-volume constraints for the six patients in this study Volume 5 dose limit
Organ
66% 66% 60% 75% 100%
Left kidney Right kidney Liver Small bowel Spinal cord
5 5 5 5 I
20 20 30 45 45
Gy Gy Gy Gy Gy
compared to the best conformal plan to see if the conforma1 technique using multiple static beams offered any dose distribution advantages over a well-optimized conventional beam arrangement. METHODS
AND
MATERIALS
Optimized treatment plans were generated for six patients with tumors in the abdomen selected from our clinic. Three patients had gastric tumors, two had pancreatic tumors, and one had a tumor in the colon. For each patient, a three-dimensional computer model of the normal anatomy and tumor volume was constructed from outlines drawn by a physician on the individual crosssections of the computed tomography study. Each normal organ was represented by dose constraint points randomly distributed throughout its volume. Because both dose computation time and optimization time increase as the number of dose constraint points increases, there is an incentive to minimize the number of points used. On the other hand, the number must be sufficiently large to accurately represent the distribution of dose in the organ. For these tests, the number of dose constraint points within each organ was determined by both the size of the organ and its proximity to the tumor volume. A minimum of 100 dose constraint points was required for each organ, and roughly the same number of points was used to represent the same organ in all cases. Point densities were typically about 0.5, 1, 1.5, and 5 points per cm’ for liver, small bowel, kidney, and spinal cord, respectively. For each patient the same dose constraint points were used with all of the annealing algorithms. The optimized plans maximized the minimum tumor dose subject to dose-volume constraints on the doses to the normal structures. For the treatment techniques under consideration, the minimum target dose occurs on the boundary of the target volume. Therefore, the target volume was represented in each CT slice by many dose points placed at 1 to 2 cm intervals along the perimeter of the target outline and a few randomly distributed within the target outline. Only solutions that met the normal tissue dose-volume constraints specified in Table 1 were considered feasible. These dose-volume constraints are equivalent to points lying on cumulative dose-volume histograms. For each set of beam weights, an approximate dose-volume histo-
ROSEN et al.
1093
gram was computed for each organ from the doses to the dose constraint points and then compared to the corresponding dose-volume constraint. Solutions that violated dose-volume constraints (nonfeasible solutions) were renormalized and used in place of the original nonfeasible solutions. Simulated annealing method Simulated annealing is an iterative process that searches for the value of the vector x that minimizes the value of the cost function F(x). The components of x are the independent variables x,, x2, . . . , x,,. For external beam treatment plan optimization, these independent variables are relative beam weights proportional to the doses delivered from beams or parts of beams. The cost function assigns a figure of merit to the plan resulting from a selected set of beam weights. Because simulated annealing is a minimization procedure, we maximized the tumor dose by choosing the negative of the minimum tumor dose as the cost function. Each iteration of the simulated annealing process has four basic steps and two relevant solutions. The process starts with an initial solution that can be obtained by any means, but must meet all constraints placed on the solution space. For treatment plan optimization problems, the initial solution is usually all beams with zero weight because this solution immediately meets all of the dose constraints on normal structures. This initial solution serves as the “current solution” for starting the simulated annealing search. The current solution is considered to be the point in the solution space closest to the global minimum. At each point in the process, the algorithm generates a “trial solution” and decides whether the trial solution should become the current solution. Ideally, the current solution will be the best solution when the simulated annealing process is terminated. However, because the Metropolis criterion occasionally moves the current solution to a higher cost function value and because the process is terminated after a finite number of iterations, the current solution at the end of a simulated annealing run may not have the lowest cost of any solution encountered during the run. Therefore, the best solution found in any iteration is usually retained and used as the final result. First, a trial solution x’ is generated from the current solution. A generating function is used to select the magnitude of a displacement vector Ax from the current solution x in the n-dimensional variable space. The direction of Ax is selected by randomly distributing its components among the independent variables. The components of x’ are then given by 4 ’ =x,
+ Ax;, (i = 1, . . . . n).
(Eq. 1)
Each of the four algorithms under consideration uses a different generating function. In general, the probability
1094
I. J. Radiation
Oncology
0 Biology
l Physics
of choosing a displacement vector of a given length decreases as the temperature decreases. The range of search is initially large with the hope that the well containing the global minimum will be discovered. Then, the search area decreases in an effort to find the minimum of the well containing the current solution. However, at every iteration there is still a finite probability of escaping the current well. Next, the trial solution is tested against constraint limits, and the value of the cost function is computed. If the trial solution meets all constraints, the algorithm continues on to the next step. If not, the trial solution is nonfeasible. It can be rejected immediately or it can be renormalized so that it does meet the constraints. Renormalization involves uniformly reducing all of the beam weights by the minimum amount needed to meet all of the constraint requirements. Renormalizing nonfeasible trial solutions is generally more efficient than discarding them. When nontrivial constraints are used, the minimum value of the cost function will fall at one of the constraint boundaries. Therefore, as the search approaches the minimum, the number of nonfeasible trial solutions increases rapidly. If these nonfeasible trial solutions are discarded, the efficiency of the algorithm decreases proportionally. Although renormalizing a nonfeasible solution is generally more computationally expensive than simply discarding it and selecting another, renormalization improves overall efficiency by providing a feasible trial solution at every iteration. Given a feasible trial solution, the algorithm must decide whether this trial solution represents an improvement over the current solution. If the cost function value for the trial solution is smaller than that for the current solution, then the trial solution is clearly superior and it becomes the new current solution. However, if the cost value for the trial solution is larger than that for the current solution, the trial solution is not immediately discarded. Rather, it is retained as the current solution with a probability given by the Metropolis criterion,
p = e -[m’k-w)l T
’
F(x’)
> F(x)
(Eq. 2)
where T is a “temperature” parameter. The probability of retaining a trial solution with higher cost decreases as the temperature of the system decreases and as the difference in cost values increases. Finally, the temperature of the system is reduced according to an annealing schedule that specifies a sequence of declining values for the temperature parameter as the optimization process proceeds. This annealing schedule requires an initial temperature value, a function describing the rate of temperature reduction, a final temperature, and the number of iterations to be performed at each temperature. A single temperature schedule may be used for all variables or each variable may have its own annealing schedule.
Volume
33, Number
5, 1995
Fast simulated annealing-Cauchy distribution generating function One of the earliest of the fast simulated annealing algorithms uses a generating function with a probability density given by a Cauchy distribution (18). The probability of obtaining a displacement vector of magnitude ) Ax) in an n-dimensional search space is given by
gOAx
=
W(k)
[(IA ()* + W(k)2]“‘+
“‘*
(Eq. 3)
where W(k) has been defined by Mageras and Mohan (9) for treatment plan optimization as W(k) =
W, (1 + WR)
(Eq. 4)
and k is the iteration number. The Cauchy distribution has a Gaussian-like peak and Lorentzian wings, allowing occasional long jumps during the search process. The parameter W,, defines the width of the Cauchy distribution, and the parameter R determines the rate at which the Cauchy distribution diminishes in size during the annealing process. Because there is no quick algorithm for calculating an n-dimensional Cauchy random number distribution, the magnitude of the step size for any generated random number was obtained by interpolation from a precalculated lookup table. The direction of the displacement vector was selected by randomly distributing its components among the independent variables. A sequence of n random numbers (u, 7 u2, . ..1 u,,) uniformly distributed in the interval [- 1, 1) was generated. These numbers define the components of a vector u in the n-dimensional space whose magnitude is
IllI = [(UT+ u; + * * - + Uf)]“2. 0% 5) The direction vector u, displacement vector Ax, and the current solution were then used to compute a new trial solution: X: = xi +
IAXI
.
0%. 6)
Any negative beam weights in the trial solution were set to zero. The temperature annealing schedule for this algorithm was given by T(k) =
” (1 + WI?)
@q. 7)
where T,, is the initial temperature and parameter R from Eq. 4 also determines the rate at which the temperature
Simulated
annealing
algorithms
for conformal
decreases. The temperature value was updated at every iteration. Based on tuning with a single case, we selected values for the control parameters of W, = 1, T,, = 0.2, and R = 500.
Fast simulated annealing-lorentzian generating function
distribution
A Lorentzian probability distribution in one-dimensional form is similar in shape to the Cauchy distribution, but is much simpler to compute. It also allows occasional large step sizes during the search process. The Lorentzian probability of obtaining a displacement vector of magnitude 1Ax 1 is given by
g(IAxl) =
W(k)* [( I Ax I >’ + W(k121 .
treatment
planning AXj
(W. 9)
where W(k) is the width of the Lorentzian distribution for the kth iteration. It has the same form as the width of the Cauchy distribution given by Eq. 4. The directions of displacement vectors and new trial solutions were generated in the same way as for the Cauchy distribution. The cooling schedule used here was the same as that of the fast annealing algorithm with the Cauchy distribution generating function (Eq. 7). Again, the values for the control parameters were determined by tuning the algorithm with one test case: W, = 10, T, = 0.2, and R = 500. The temperature value was also updated at every iteration.
Variable step size generalized simulated annealing (VSGSA) The variable step size generalized simulated annealing (VSGSA) algorithm is a modification by Sutter and Kalivas (17) of the generalized simulated annealing algorithm originally developed by Bohachevsky et al. (1). In this algorithm, the components of the displacement vector are not computed according to a theoretical probability distribution. Rather, they are calculated to maintain a relatively high efficiency of finding improved solutions as measured by the acceptance ratio, which is defined as the number of accepted trial solutions divided by the number of trial solutions evaluated within a specified block of iterations. An accepted trial solution is one that becomes the current solution by virtue of either a lower cost value or the Metropolis criterion. The change in beam weight for ith beam, AX, is calculated from a random number u, uniformly distributed in the interval [ - 1, 1) using the expression:
0 I. I. ROSEN et al. = hjUi(Ximax
-
1095 Ximin)-
(Eq.
10)
The minimum value for each beam weight, X,,inr was set to 0. The value of ximaxfor each beam was determined by finding the maximum weight that would generate a feasible solution when all other beam weights were set to zero. The values of Ximinand Ximaxwere determined in advance and not changed. The parameter hj controls the overall magnitude of step sizes for the jth block of iterations and it changes according to the acceptance ratio. If the acceptance ratio is greater than 0.9 for a block of iterations, then the current solution is considered to be far from the global minimum and the step size factor is increased according to
hj+, = mhj,
(Eq. 8)
Given a random number r uniformly distributed in the interval [0, I), the magnitude of a random displacement vector Ax is obtained from
1Ax ( = W(k) etan(O.%rr),
therapy
0% 11)
where h, is the old value for the jth block of iterations and m is a factor greater than 1, whose value is set by the user. If the acceptance ratio is less than 0.2, then the current solution is considered to be close to the global minimum and the step size factor is reduced according to
hj+, = h,lm.
0%. 12)
The initial value of h was chosen to be 0.20; the value of m was set at 1.05; and the acceptance ratio was calculated every 100 iterations. The probability of accepting a trial solution with a higher cost value is given by the function p = e~lcc,-c”)l~-tc,~c”,l 9
(Eq. 13)
where C, is the cost value of the trial solution, C,. is the cost value of the current solution, and C, is the global optimum value estimated by the user. The initial value for the temperature, T, was determined so that the probability P was between 0.2 and 0.7. Thereafter, the value of T was adjusted whenever a detrimental step occurred after the step size had been changed.
Very fast simulated
reannealing
(VFSR)
In the very fast simulated reannealing (VFSR) algorithm developed by Ingber (6), each independent variable has its own exponentially decreasing temperature, I;, given by
T;(k) = Ti(O)exp(-cik”“),
(Eq. 14)
where k is the iteration number, n is the number of independent variables, and ci is a user-defined parameter for the ith variable. The temperature of each variable is periodically resealed, depending on the sensitivity of the cost function to that variable. The purpose of resealing the temperature for each variable independently is so that the
IO96
I. J. Radiation Oncology 0 Biology 0 Physics
Volume 33, Number 5. 1995
generating function will give larger displacement vectors when searching relatively flat and insensitive dimensions and smaller displacement vectors in more rapidly varying dimensions. The sensitivity, s,, of the cost function to a variable is determined by sampling the pseudogradient of the steepestminimum found using numerically computed partial derivatives of the cost function with respect to that variable:
II
dF s; = . 8%
@I. 13 4
“Patient
Then, in each reannealing step, a new temperature T,‘(k) for each variable is computed, which is proportional to the old one according to the sensitivity of the variable,
where s,,, is the largest value of the sj computed during the reannealing step. This resealing of the temperature, or reannealing, can often improve the overall efficiency of the algorithm. However, any reduction in the number of iterations required for convergence comes at someexpense in increased computation time per iteration. The generating function for VFSR is temperature dependent. The change in weight for the ith beam is calculated from (Eq. 17) yi = sign(u, - 0.5)Ti[(L
+ -$‘“”
- 11
(Q.
18)
where ui is a random number in the interval [- 1,l) and the value of yi lies in the interval [ - 1, 1). Repeated values of Ax, are generated until a new value for the ith variable is obtained that is bounded by the interval [Ximinrx;~J. The values for X;,in and x,,,, are determined as for the VSGSA algorithm. Calculation of optimized treatment plans For the conformal plans, 36 coplanar beamsof 15 MV x-rays were deployed around the patient at 10” intervals, giving 36 variables for optimization. The conventional plans used anterior and lateral beams of 15 MV x-rays. Each beam direction included an open beam and two beams with 60” wedges oriented in opposite directions for a total of nine optimization variables. All beamswere shapedto conform to the projections of the tumor volume
’ VAXstation 4000 Model 60, Digital EquipmentCorporation, Maynard,
MA.
5
6
Fig. 1. Minimum target doses for the optimized conformal treatment plans generated using the different simulated annealing algorithms with long run times. The optimization objective was to maximize the minimum target dose subject to normal tissue dose-volume limits. The VSFR algorithm was run for 80,000 iterations and the others were run for 100,000 iterations each.
with a margin for beam penumbra. Doses were computed using a scatter-air ratio sector summation algorithm with a Batho power-law correction for tissue heterogeneities. Short and long optimization runs were performed for each patient. The short runs consisted of 8,000 iterations for the VFSR algorithm and 10,000 for each of the others. Preliminary tests showed that the VFSR algorithm was approximately 20% slower per iteration than the others. Therefore, the number of iterations for that algorithm was reduced so that elapsed times would be comparable for all four algorithms. The long runs consisted of 80,000 iterations for the VFSR algorithm and 100,000 iterations for each of the others. The long runs were performed to give an indication of how rapidly the short runs were converging to the best possible solution. Computation times’ for the short runs were approximately 20-30 CPU min for the conformal plans and 5-8 CPU min for the conventional plans. Becauseof the statistical nature of simulated annealing, running
the same treatment-planning
problem
with
two
different random number sequenceswill produce two different solutions. The longer the run time, the smaller will be the variation between the different solutions. Therefore, the short runs for the conformal plans were repeated five times for each patient using different random number sequences.The long runs were performed only once. Random numbers were obtained from the linear congruent random number generator in the VAXstation system software. RESULTS Results of the long optimization runs for the conformal plans are shown in Fig. 1. The minimum target dose is
Simulated
annealing
algorithms
for conformal
shown for each patient and each algorithm. Overall, the VSGSA algorithm gave the best results. For five of the six patients, the VSGSA plan had the highest minimum target dose, with improvements ranging from 0.1 to 3.8 Gy. Averaged over the six patients, the VSGSA plans were 0.8 Gy higher than the VFSR plans, 1.7 Gy higher than the Cauchy plans, and 2.1 Gy higher than the Lorentzian plans. For the conformal beam arrangements, the results of the five short runs using different random number sequences were averaged for each patient and for each simulated annealing algorithm. In Fig. 2, the average values and the corresponding standard deviations are shown. For each of the six patients, the VSGSA algorithm produced the highest average cost value: improvements over the VFSR plans ranged from 1.O to 3.3 Gy; improvements over the Cauchy plans ranged from 0.5 to 4.7 Gy; and improvements over the Lorentzian plans ranged from 0.4 to 5.9 Gy. The results of the short runs are compared to the results of the long run times for the VSGSA algorithm in Fig. 3. In five of the six patients, the VSGSA short run average was within 0.8 Gy of the long run solution. In the other case (Patient 5), the VSGSA short run average was 1.1 Gy less than the long run solution. For the conventional plans, the results produced by the four algorithms typically differed by less than 1% and the maximum difference was 1.5%. The variations produced by using different random number sequenceswere typically less than 0.5%. Although the differences were small, the VSGSA algorithm again gave the best average result for each case. Figure 4 shows the average computation times in CPU minutes and the corresponding standard deviations for the four algorithms for calculating optimized conformal plans using the short runs. The CPU times for the short runs
therapy
treatment
planning
l I. I. ROSEN
Bi ::: jjj i:: ii: I// g// ii: ii: ii_ :i: ::: ::. j/: ::: i:: ii: ::::: /Iii. -
1
2
2
3
Patient4
5
6
Fig. 2. Average minimum target doses for the optimized conforma1 treatment plans generated using short run times of the different simulated annealing algorithms. The vertical bars represent the averages and the error bars represent the standard deviations for five runs using different random number sequences. The optimization objective was to maximize the minimum target dose subject to normal tissue dose-volume limits. The VFSR algorithm was run for 8,000 iterations and the others were run for 10,000 iterations each.
3
Patient4
5
6
Fig. 3. Comparison of short run time results with long run time results for the VSGSA algorithm.
are quite similar for all algorithms and cases. With the exception of Patient 2, the time per iteration for the Cauthy, Lorentzian, and VSGSA algorithms were all within 4% of each other. In general, the compute times varied only slightly with different random number seeds. For Patient 6, the VFSR had an unusually large standard deviation. For this patient, the number of times the algorithm did reannealing varied vastly with different sequencesof random numbers. In three of the cases, the conformal plan was significantly better than the optimized conventional plan. The short VSGSA runs for the conformal and conventional techniques are compared in Fig. 5. For Patients 2, 3, and 5, the conformal plan increased the minimum tumor dose by 8.2 to 13.0 Gy over the conventional treatment. For Patients 1 and 4, there were only small differences (< 0.8 Gy) between the two treatment methods. For Patient 4, this was probably because the conventional plan dose was already very high (66.5 Gy). In one patient (Patient 6), the conventional plan was better than the conformal plan. This patient had a very large tumor, and the use of wedges in the conventional plan was a significant advantage that the conformal beam arrangement could not overcome, even by the useof many more beams. DISCUSSION
1
et a[.
AND CONCLUSIONS
One of the advantages of simulated annealing over other optimization methods is that it is easy to implement complex cost and constraint functions. However, the use of absolute constraints creates discontinuities in the solution space, which in turn, substantially increasesthe difficulty of solving the problem. For nontrivial constraints, the global minimum will always lie at the boundary of a constraint. As the solution approaches the boundary of the valid solution region, most increases in beam weights will generate solutions that violate the dose-volume con-
straints and, therefore, the probability of generating a feasible solution decreases. Renormalization of beam weights whenever organ dosesexceed dose-volume con-
I. J. RadiationOncology0 Biologyl Physics
1098
3s z .3 g
30
,g
20
!3
10
5
150 1
25
5
2
3
4
Patient
5
6
Fig. 4. Average computationtimes(CPU minutes)and standard deviationsfor calculating conformal plansusing short runs of the four algorithms-8,000 iterationsfor the VFSR algorithm and 10.000 for each of the others.
straints improves efficiency. Unfortunately, this modification violates the convergence criteria of simulated annealing algorithms and could causesolutions to be trapped in local minima. In the VSGSA algorithm, the calculated acceptance ratio gives a measure of the closenessof the current solution to the global minimum. With renormalization, however, the acceptance ratio is no longer related to the point at which it was evaluated. Nevertheless, we did find a valuable increase in overall efficiency when renormalization was implemented. An alternative method for constraining solutions is to add a penalty value to the cost function whenever solutions fall outside some defined limits (11). Such “soft” constraints may be used in place of or in conjunction with the “hard” limits described above. We have avoided using soft constraints because they represent compromisesbetween the absolute constraints and the objective, and the most appropriate compromise for any given problem will generally depend on the preferences of the physician and the details of the patient’s case. Rather than attempting to incorporate the physician’s preferences into the cost function, we present the best possible solution for the given hard constraints and require the physician to explicitly identify what constraints to relax if a higher tumor dose or other objective is desired. The use of soft constraints in place of hard constraints simplifies the simulated annealing optimization algorithm and probably gives better convergence becausethere are no discontinuities in the solution space. If soft constraints are used in addition to the hard ones, it is not clear that there is any additional benefit in terms of algorithmic efficiency. Hard constraints clearly separate solutions that are acceptable from those that are not. Soft constraints provide a smooth transition between those two types of solutions and in doing so become part of the definition of the optimum solution by the way in which the penalties are imposed. We found that the VSGSA algorithm consistently produced better results than the other algorithms. No simulated annealing technique can produce a guaranteed opti-
Volume33,Number5, 1995 ma1 solution with a single finite run. However, the VSGSA method performed best using both long and short run times for both conformal and conventional techniques. For the long runs, the VSGSA algorithm gave the highest minimum tumor dose for five of the six cases. For the short runs, the VSGSA solutions averaged over five runs gave minimum tumor doses that were larger by 1.8 to 5.9 Gy in 67% of the tests. An increase in tumor dose of 1.8 Gy is approximately equal to an increase of one treatment fraction and, therefore, could be clinically significant. Although only one cost function and constraint type were used in this study, we expect that the VSGSA algorithm will still be the most efficient with other cost and constraint functions. Because the performance of a given simulated annealing algorithm is dependent on the nature of the problem to be solved, it is appropriate to ask if the results presented would be applicable to other conformal therapy techniques. In other words, would the VSGSA algorithm still perform best for optimizing treatments with small numbers of noncoplanar or intensity-modulated beams?A definitive answer would require testing the various algorithms using those treatment techniques. However, for the objective function and constraints under consideration, the VSGSA algorithm performed best for 9 beam variables (conventional technique) and 36 beam variables (conformal technique) with the superiority of the VSGSA algorithm over the others increasing with the number of beams. Based on these observations, we expect that the VSGSA algorithm would still be most efficient when applied to the problem of planning with intensity-modulated beams. Conformal therapy using multiple static shapedcoplanar beamsappearsto offer clinical advantagesin the treatment of abdominal tumors. In half the patients studied the minimum tumor dose increased by more than 8 Gy for the samenormal tissue dose-volume constraints. It should be remembered that the conformal plans were compared to optimized conventional beam arrangements,not to standard plans. Although it was not evaluated as part of this study, the gain should be even greater when comparedto standard
ii: ::: ::: ii: /> /i: ::: .j/ .:: ii: :lj :i: i:: ::: ii: :;/ 1
::: ii: ::: ::: :// jj/ j// ii: ii: hsi -I
3
Patient4
5
6
Fig. 5. Comparisonof conformalplansand conventionalplans optimized using VSGSA.
Simulated
annealing
algorithms
for conformal
plans. In the one patient for whom the conventional plan was better, the large tumor size gave a substantial advantage to the use of wedges. Including wedges in the conformal technique is simple and should provide some additional increase in tumor dose at the expense of increased computer
therapy
treatment
planning
0 I. 1. ROSEN er al.
1099
time in planning and possibly increased treatment delivery time. The routine use of this conformal technique is no more difficult or complex than many conventional treatments but will require a computer-controlled accelerator and multileaf collimator for efficient automated delivery of treatments.
REFERENCES 1. Bohachevsky, I. 0.; Johnson, M. E.; Stein, M. L. Generalized simulated annealing for function optimization. Technometrics 28:209-217; 1986. 2. Bortfeld, T.; Schlegel, W. Optimization of beam orientations in radiation therapy: Some theoretical considerations. Phys. Med. Biol. 38:291-304; 1993. 3. Corana, A.; Marchesi, M.; Martini, C.; Ridella, S. Minimizing multimodal functions of continuous variables with the “simulated annealing” algorithm. ACM Trans. Math. Soft. 13:262-280; 1987. 4. Geman, S.; Geman, D. Stochastic relations, Gibbs distributions, and the Bayesian restoration of images. IEEE Trans. Patt. Anal. Mach. Intel. PAM1 6:721-741; 1984. 5. Groisman, G.; Parker, J. R. Computer-assisted photometry using simulated annealing. Comp. Phys. 7:87-96; 1993. 6. Ingber, L. Very fast simulated reannealing. Math. Comput. Modelling 12:967-973; 1989. 7. Ingber, L. Statistical-mechanical aids. to calculating termstructure models. Phys. Rev. [A] 42:7057-7064; 1990. 8. Kirkpatrick, S.; Gelatt, C. D.; Vecchi, M. P. Optimization by simulated annealing. Science 220:671-680; 1983. 9. Mageras, G. S.; Mohan, R. Application of fast simulated
annealing to optimization of conformal radiation treatments. Med. Phys. 20:639-647; 1993. 10. Metropolis, N.; Rosenbluth, A. W.; Rosenbluth,
M. N.;
14.
15. 16. 17.
18. 19. 20. 21.
Teller, A. H.; Teller, E. Equation of state calculationsby fast computing machines. J. Chem. Phys. 21:1087-1092; 1953. 11. Mohan, R.; Mageras, G. S.; Baldwin, B.; Brewster, L. J.; Kutcher, G. J.; Leibel, S.; Burman, C. M.; Ling, C. C.; Fuks, Z. Clinically relevant optimization of 3-D conformal treatments. Med. Phys. 19:933-944; 1992. 12. Merrill, S. M.; Lane, R. G.; Jacobson, G.; Rosen, I. I. Treatment planning optimization using constrained simulated annealing. Phys. Med. Biol. 36:1341-1361; 1991. 13. Morrill, S. M.; Lane, R. G.; Rosen, I. I. Constrained simu-
22. 23.
lated annealing for optimized radiation therapy treatment planning. Comp. Methods Prog. Biomed. 33: 135- 144; 1990. Niemierko, A. Random search algorithm (RONSC) for optimization of radiation therapy with both physical and biological end points and constraints. Int. J. Radiat. Oncol. Biol. Phys. 23:89-98; 1992. Sloboda, R. S. Optimization of brachytherapy dose distributions by simulated annealing. Med. Phys. 19:955-964; 1992. Sloboda, R. S.; Pearcey, R. G.; Gillan, S. J. Optimized low dose rate pellet configurations for intravaginal brachytherapy. Int. J. Radiat. Oncol. Biol. Phys. 26:499-511; 1993. Sutter, J. M.; Kalivas, J. H. Convergence of generalized simulated annealing with variable step size with application toward parameter estimations of linear and nonlinear models. Anal. Chem. 63:2383-2386; 1991. Szu, H.; Hartley, R. Fast simulated annealing. Phys. Len. [A] 122:157-162; 1987. Webb, S. Optimisation of conformal radiotherapy dose distributions by simulated annealing. Phys. Med. Biol. 34:1349-1370; 1989. Webb, S. SPECT reconstruction by simulated annealing. Phys. Med. Biol. 34:259-281; 1989. Webb, S. Optimization by simulated annealing of threedimensional conformal treatment planning for radiation fields defined by a multileaf collimator. Phys. Med. Biol. 36:1201-1226; 1991. Webb, S. Optimization of conformal radiotherapy dose distributions by simulated annealing: 2. Inclusion of scatter in the 2D technique. Phys. Med. Biol. 36: 1227-1237; 1991. Webb, S. Optimization by simulated annealing of threedimensional, conformal treatment planning for radiation fields defined by a multileaf collimator: II. Inclusion of two-dimensional modulation of the x-ray intensity. Phys. Med. Biol. 37:1689-1704; 1992.