A robust optimization method with successive linear programming for intensity-modulated radiation therapy

Intensity-modulated radiation therapy for cancer is considered to be effective when dealing with complicated tumour shapes because the dose distribution for each irradiation can be modulated. Fluence map optimization is often formulated as an optimization problem with dose volume constraints (DVCs). A linear programming (LP) method that approximated DVCs was proposed, and it was modified to the successive LP method (SLPM) to find a feasible treatment plan in a wider region. In the present paper, we propose a numerical method called SLPM-R (the SLPM with robustness) that enhances the SLPM using a robust optimization approach. We mathematically prove that the proposed method with extended LP problems has the favourable properties of the SLPM, even taking uncertainty in the influence matrix into consideration. In particular, when the optimal value of the LP problem is non-positive, the proposed SLPM-R guarantees that the output solution can satisfy all DVCs. Through numerical experiments, we observed that the proposed method found a feasible plan that the SLPM could not find. In addition, for a test case that even the SLPM-R failed, the largest deviations of 5.65 Gray in the SLPM was reduced to 3.15 Gray by the SLPM-R.


Introduction
Advances in technology have improved the accuracy and precision of radiation therapy, so that in addition to surgery and chemotherapy, it is now recognized as one of three major treatments.In particular, intensity-modulated radiation therapy (IMRT) is used in various cases, such as tumours with a concave shape, because the radiation intensity can be adjusted to the shapes by the use of a multileaf collimator (MLC), and IMRT is frequently used in oncology.In IMRT treatment planning, it is desirable that tumours (planning target volumes: PTVs) receive a reasonably high dose and healthy organs near tumours (organs at risk: OARs) receive a much lower dose [1], and such computation can be formulated as an optimization problem.
There are three optimization problems: beam angle optimization (BAO), fluence map optimization (FMO), and the leaf sequencing problem.An important goal of BAO is to find the best set of angles from different candidates based on an objective function that reflects the fitness of the treatment plan.(See, for example, [2,3].)Fluence map optimization is a problem of optimizing beam intensities.The leaf-sequencing problem is to determine the setting of MLC that conform to a given dose distribution.
Among many optimization aspects, the focus in the present paper is on an optimization problem related to FMO computation.Fluence map optimization sometimes contains specific constraints called dose volume constraints (DVCs).For example, an upper DVC of form U 0.1 Core = 25 can be roughly considered as a constraint in which the area that can receive a dose of 25 Gray or more is at most 10% of the organ Core.(Throughout the present paper, we use the Gray (Gy) as the unit of dose, and a precise definition of DVCs will be given in Section 2.) Therefore, it is effective in the FMO computation to identify the outliers (the areas that exceed the threshold).In 2002, Merritt et al. [4] proposed an iterative method that updates the outliers by solving linear programming (LP) problems.However, this method cannot directly take the fraction (like 10% of U 0.1 Core = 25) into consideration, hence, the areas identified by the method of Merritt et al. do not always reflect the fractions prescribed in DVCs.
If we express DVCs rigorously in optimization models, then DVCs require integer variables.Solving such optimization problems is known to be NP-hard [5].Therefore, a reduction in computation time is an important factor.In 2003, Romeijn et al. [6] introduced the concept of Conditional Value-at-Risk (C-VaR) to replace the DVCs with constraints that can be described in LP problems.Since LP problems can be solved in a polynomial time by interior-point methods, this C-VaR method can reduce the computation time compared to approaches that are strict to DVCs.However, the region that satisfies the C-VaR constraints is much narrower than that with the original DVCs, and this method sometimes fails to find a feasible plan.Kishimoto and Yamashita (2018) [7] relaxed the C-VaR constraints by detecting outliers so that the resultant LP problems always have a solution.Their successive LP method (SLPM) repeatedly solves LP problems updating the outliers and can find a feasible plan between the C-VaR constraints and the original DVCs.They reported based on numerical experiments on the C-Shape instance of TG119 [8] that the SLPM found beam intensities that satisfy all the DVCs, while the C-VaR approach by Romeijn et al. [9] failed.Compared with the computation time of 230 seconds by the iterative method of Merritt et al. [4], the SLPM completed its computation within 101 seconds.The number of iterations in the SLPM was less than that of Merritt et al., and this implies that the SLPM identified the outliers more efficiently.
On the other hand, in medical practice, the sequence of treatment planning includes various uncertainties, such as errors in the tissue densities in CT images, inaccuracy in delivering the correct dose, and positioning uncertainties due to patient movement during irradiation.These uncertainties should be addressed to develop treatment plans with more robustness.Chan et al. [10,11] developed an algorithm to reduce the total error by considering shape changes in organs and tumours with some uncertainties due to, for example, breathing or tumour motion.Stemkens et al. [12] proposed a framework to generate a subject-specific motion model on a voxel-by-voxel basis by performing a principal component analysis.However, the C-VaR method and the SLPM described earlier do not consider these uncertainties, so these two methods might be vulnerable to such uncertainties.
In the present paper, we propose a numerical method that combines the concept of robust optimization with the SLPM.It was proven in a previous paper [7] that the SLPM can find a treatment plan that satisfies all of the DVCs when the optimal values of its LP problems drop below zero.By extending the proof in the previous paper [7], we show that the proposed method still possesses this favourable property, even though the LP problems in the proposed method involve additional variables for robust optimization.
Through a numerical experiment with test instances of TG119 [8], the proposed method obtains solutions that satisfy all of the DVCs in more situations than the SLPM.Even when the proposed method cannot find a feasible solution, it is effective at reducing the largest deviations from DVCs.For the C-Shape instance with uncertainty, the largest deviations by the SLPM were 5.65 Gy and 4.41 Gy, while the deviation by proposed method was 3.15 Gy.In the Prostate instance, the deviations of 13.29 Gy and 9.43 Gy by the SLPM were reduced to 3.75 Gy in the proposed method.
In addition, when we extend the objective function in LP problems with penalty terms, we can give a higher priority to specific DVCs.For the Head and Neck dataset in the TG119 instances, all DVCs are satisfied by the use of the penalty terms, and for the MultiTarget dataset, we can reduce the deviations from the DVCs for the organ that receives the highest dose from 2.67 Gy to 1.48 Gy.
The remainder of the present paper is organized as follows.In Section 2, we introduce notation related to FMO and then explain existing methods.We describe the proposed method in Section 3 and discuss its mathematical properties.Section 4 shows the results of numerical experiments on the TG119 instances.In Section 5, we discuss an extension of the proposed method by incorporating the uncertainty in lung states considered by Chan et al. [10,11] and show that a corresponding problem in each iteration remains as an LP problem.Finally, the conclusions are given in Section 6.

Preliminaries and existing methods
In this section, we briefly introduce notation in FMO and then describe the C-VaR method proposed by Romeijn et al. [6,9], which approximates the difficult DVCs by linear constraints.In Section 2.2, we also briefly discuss the framework for the successive LP method (SLPM) proposed in a previous study [7].

Dose volume constraints and C-VaR type constraints
In IMRT optimization, in order to calculate the dose efficiently, the beams are discretized into small areas called beamlets.Similarly, the organs are also discretized into small volumes called voxels.Let S and J be the set of organs and the set of beamlets, respectively.We use I s to denote the set of voxels in s ∈ S. We can calculate the dose for the ith voxel in organ s ∈ S as z si = j∈J [D s ] ij x j , where x j is the intensity of beamlet j, and D s is an influence matrix for organ s, i.e. [D s ] ij represents the absorbed dose for voxel i in organ s from beamlet j at unit intensity.The influence matrices are also called dose description matrices.The size of D s is |I s | × |J |, where the notation |X| is used to denote the cardinality of a set X.
For the availability of a treatment plan, it is desirable to satisfy all of the DVCs.Dose volume constraints are classified into two types: lower DVCs and upper DVCs.Let A s ⊂ (0, 1) and A s ⊂ (0, 1) denote the sets of ratios used in the lower and upper DVCs, respectively, for organ s.A lower (upper) DVC on organ s with respect to a ratio α ∈ A s (α ∈ A s ) is a constraint that the fraction of voxels that receive at least L α s Gy should be no less than α (no more than α, respectively).More precisely, the lower DVC and an upper DVC can be formulated as respectively.In a mixed-integer linear programming formulation, the lower DVC is expressed as the following constraints: This formulation expresses a DVC rigorously but involves a binary variable for each DVC and each voxel.Therefore, the formulation often requires a long computation time.Fluence map optimization problems with DVCs have been proven to be NP-hard [5].
In order to reduce the computation cost, faster optimization approaches that do not involve binary variables are required.Based on the concept of C-VaR [13], Romeijn et al. [6,9] replaced the time-consuming lower and upper DVCs with cheaper linear C-VaR constraints of the following forms: Figure 1 illustrates the difference between a DVC and its C-VaR constraint.(This figure will also be used to illustrate the use of a hot spot in the approach of [7] described in the next subsection.)The horizontal and vertical axes represent the absorbed dose (Gy) and the percentage of voxels in a structure, respectively.The blue curve in Figure 1 is a dose volume histogram (DVH).For example, if the histogram passes a point (50 Gy, 95%), then 95% of voxels receive a dose of 50 Gy or higher.
For α ∈ A s , a DVC demands that the lowest dose received by the highest α fraction of voxels (the left end of the red and blue areas) be at most U α s .In contrast, the C-VaR constraint (1) requires the average dose received by the highest α fraction of voxels (the average of the red and blue areas, which is indicated by the dotted vertical line 'conventional CVaR α ') to be at most U α s .The average dose is larger than the lowest dose, thus any solution satisfying the C-VaR constraint of form (1) also satisfies the original DVC.(A mathematical proof of this property was given in a previous study [7].)The gap between the average and lowest doses implies that the feasible region of the C-VaR constraints is narrower than that of the DVCs.Therefore, it was pointed out in [7] that the C-VaR method [6,9] may discard a feasible solution of the original DVCs.

Successive linear programming method
In order to reduce the gap between a DVC and the corresponding C-VaR constraint, Kishimoto and Yamashita [7] proposed the successive linear programming method (SLPM) by introducing the concept of hot and cold spots.The average dose related to an upper DVC is affected strongly by a small number of voxels that receive extremely high doses.In the SLPM, such voxels are automatically detected as a hot spot.In Figure 1, a hot spot is illustrated as the blue area.The absorbed dose in the hot spot does not directly affect the satisfiability of the DVC, because the absorbed dose in the blue area in Figure 1 is not the lowest dose in the highest α fraction of voxels (the left end of the red area).Removing the hot spot from the computation of the average dose shifts the average to the left.Thus, the gap between the lowest dose and the average dose (in only the red area) will be tighter.Therefore, satisfying the DVC becomes easier.
A framework of the SLPM can be given as Algorithm 1.For the kth iteration, we use R α k,s ⊂ I s and R α k,s ⊂ I s to denote cold and hot spots, respectively, and we solve the following LP: , P α s , P α s , P s , P s are parameters that control the weights of DVCs.The constraint of (2c) is derived from bounds L s ≤ z si ≤ U s that each voxel should satisfy.The cold spot R α k,s and the hot spot R α k,s are removed from the computation of the averages in (2d) and (2e), respectively.These spots are updated by the rule given by (3) in Algorithm 1.
Kishimoto and Yamashita [7] proved the following proposition.

Proposition 2.1 ([7]):
The SLPM with (2) as the kth LP problem has the following three properties: (i) For each k ≥ 1, the kth LP (2) has an optimal solution.
(ii) If the optimal value t k in the kth LP (2) satisfies t k ≤ 0, then all DVCs are satisfied.(iii) The sequence {t k } is monotonically non-increasing, i.e.t k+1 ≤ t k for k ≥ 1.
Note that the optimal value t k in (2) can be regarded as the largest deviation from DVCs at the kth iteration adjusted by the parameters P α s , P α s , P s , and P s [7].
end for return x K Property (i) guarantees that Algorithm 1 can find a solution in each iteration, because the SLPM can detect the cold and hot spots adequately.This is different from the C-VaR method [6,9], which cannot output useful information if no solution can satisfy all of the C-VaR constraints, even when there is a solution that satisfies all DVCs.Furthermore, due to Property (iii), the SLPM reduces the deviation in each iteration, and this leads to a solution that satisfies all DVCs when t k ≤ 0 in Property (ii).Note that there is a possibility that the set of solutions that satisfy all of the DVCs is empty.For such a case, many approaches that relax some DVCs can be considered.As discussed in a previous study [7], Algorithm 1 can also partially give information for such a relaxation based on t k .

Robust optimization
For later discussions, we give a brief introduction to robust optimization.For more details, see Ben-tal et al. [14] and the references therein.Robust optimization has been applied to many situations.For example, Özmen et al. [15] employed robust optimization to analyze regulatory systems under polyhedral uncertainty.In [16], robust optimization was used in the context of a green agrifood supply chain.We can also find papers that utilize robustness in forestry researches, for example, by Orrego et al. [17] and Salmanmahiny et al. [18].
Roughly speaking, the concept of robust optimization is to optimize a given objective function over a feasible set that includes uncertainty.Suppose that we are solving an optimization problem where the input data are In global robust optimization, the objective function is optimized in a region in which the constraints are satisfied for any input data in Z.Thus, a robust counterpart of (4) can be formulated as min ( 5 ) and this can be equivalently rewritten as min For solving optimization problems that involve constraints like ( 5) or ( 6) in a short time, the perturbation set Z is usually a box (like −δ ≤ ξ l ≤ δ for each l = 1, . . ., L with a parameter δ > 0) [19] or an ellipsoid [14].It was shown in a previous study [14] that optimization problems of a linear objective function over the box or ellipsoid-shaped perturbation set can be converted into LP problems or second-order cone programming (SOCP) problems, respectively.Since computation time is an important factor in IMRT optimization, and the computation cost of LP problems is usually less than that of SOCP problems, our interest is the box-shaped perturbation set.

Proposed method
The SLPM [7] does not consider data that includes uncertainty of beam irradiations or a movement during treatment.Therefore, it is more practical to solve an optimization problem assuming that the data contains uncertainty.We propose a numerical method that combines the SLPM with the concept of robust optimization, and we will show that Proposition 2.1 can still hold in the proposed method.

Framework for proposed method
The framework for the proposed method is given in Algorithm 2 and is a modification of Algorithm 1 with a replacement of the kth LP problem (2) with its robust counterpart (8) below.We refer to this method as the SLPM with robustness (SLPM-R).
A treatment plan in clinical practice contains various uncertainties, such as measurement errors in CT, and numerical errors in the quadrant infinite beam (QIB) method [20].In the proposed method, we focus on the uncertainty in end for return x K the influence matrix and assume that other uncertainty elements are implicitly reflected in the uncertainty for the influence matrix.For example, Shan et al. [21] discussed a model in which uncertainties are considered in influence matrices.
Recall that the order of the influence matrix D s is |I s | × |J |.We assume that a perturbation set Z s is given as a box set, so that D s is described as where are parameter matrices and s is taken from the box-shaped perturbation set We can choose the same δ > 0 for all voxels by adjusting D 0 s and D s appropriately.For instance, if 3 In the numerical experiment described later herein, we will change δ to evaluate the effect of the uncertainty range.
In order to derive a robust counterpart of the LP problem (2) in the SLPM, we first focus on an upper bound for each voxel of form j∈J [D s ] ij x j ≤ U s + P s t.In similar steps from (4) to (6), we convert this constraint: where |D s | is the matrix that takes element-wise absolute values of D s .
In order to apply the same procedure to the inequalities in (2), we split the variable z si into two variables, z si and z si , for lower and upper DVCs, respectively.Consequently, we derive the robust counterpart of (2) as follows: Corresponding to the split of z si into z si and z si , the update rule (3) in the SLPM for the cold spots and hot spots is modified as (7) in the SLPM-R.
In (8), we use the simple objective function t in the same way as (2).We will consider a variant of this objective function in Section 3.3.

Properties of proposed method (SLPM-R)
The proposed method (SLPM-R) shares the basic framework with the SLPM of [7].However, it is not obvious as to whether the SLPM-R can maintain the three properties in Proposition 2.1, because we split z si into z si and z si .Therefore, we extend the proof in [7] along with the SLPM-R.As a result, the SLPM-R retains the properties, and this indicates that the SLPM-R can find a favourable treatment plan using robust optimization.
Proposition 3.1: The three properties (i), (ii), and (iii) in Proposition 2.1 hold in the SLPM-R.
To prove Proposition 3.1, we will use Lemma 3.1 below.Thus, we first give a proof for Lemma 3.1, and then prove Proposition 3.1.

Lemma 3.1:
For each s ∈ S and α ∈ A s , any feasible point in (8) Similarly, for each s ∈ S and α ∈ A s , any feasible point in (8) Proof: For each s ∈ S, α ∈ A s , we know that 0 < α < 1.Thus, it holds that Therefore, from (8f), we obtain where the last non-negativity is derived from an inequality p + (q − p) + ≥ 0, which holds for ∀p ∈ R, ∀q ≥ 0. Therefore, from (8g), we obtain U α s + P α s t ≥ 0.
We are now prepared to give a proof of Proposition 3.1.

Proof of Proposition 3.1.:
We start from Property (i).Since ( 8) is an LP problem, we can use the duality theorem [22].Therefore, it is sufficient to show two points: (a) a feasible solution exists in (8), and (b) there is a lower bound of the objective function t.We will prove these by induction.When k = 1, it holds for all s ∈ S that Thus, the denominators (8f) and (8g) are nonzeros.Let 8) has at least one feasible solution.Next, we verify a lower bound of the objective function t.From Lemma 3.1, the objective function t of LP ( 8) has a lower bound: We assume (a) and (b) for the kth LP and consider the k denote an optimal solution of the kth LP.For (8g), we temporarily assume that the number of voxels that are newly added to hot spots after the kth iteration is greater than or equal to but this inequality is inconsistent with (8g).Thus, we know that Regarding cold spots, we can also derive Thus, the denominators (8f) and (8g) are nonzero, and we can apply the same proof as in the first iteration to show (a) and (b) in the (k + 1)th LP.
Next, we consider Property (ii).Recall that t k is the optimal value t of the kth LP (8).When , and this indicates that the upper DVC for α ∈ A s holds for any s ∈ Z s .In regard to the lower DVCs, we can show this property in a similar manner.Finally, we discuss Property (iii).From Property (i), there exists an optimal solution for any k ≥ 1.Thus, it is sufficient to find a feasible solution in the (k + 1)th LP, the objective value of which is t k .We show that a feasible point in the (k + 1)th LP can be constructed with The objective value for this solution is t = t k , and (8b), (8c), (8d), and (8e) hold in the (k + 1)th LP, because these constraints are not affected by the updates of the hot and cold spots.
The constraint (8f) involves R α (k+1),s , but we can still show that Here, we used In the same way, the inequality (8g) holds.Therefore, we can find a feasible solution, the objective value of which is t k .This indicates t k+1 ≤ t k and completes the proof.

Proposed method with penalty terms
Although the SLPM can evaluate the deviation from DVCs by the optimal value t k , as discussed in [7], there remain voxels that receive much higher or lower doses.In order to reduce such voxel-wise deviations, we modify the objective function t in (8a) by adding penalty terms with thresholds θ α s , θ α s , θ s , θ s > 0 and weight parameters λ α s , λ α s , λ s , λ s ≥ 0 as follows: In addition, oncologists sometimes give a higher priority to PTVs than to healthy organs, and if there are several PTVs, then they prioritize a PTV that requires the highest absorbed dose.
As will be shown in the numerical results in Section 4, the penalty terms improve the solution quality for some cases by reaching a region where the simple objective function t cannot search.If we combine a framework of bilevel optimization [23], an automatic selection of the best weight parameters λ would be possible.However, it would complicate the analysis of the properties of the proposed method in Proposition 3.1, thus we would like to leave the automatic selection as a future work.

Numerical experiments
In this section, we discuss numerical experiments that compare the proposed method (SLPM-R) to the existing SLPM.The tests in the present study were performed on a Linux server with two Opteron 4386 (3.10 GHz) CPUs and 128 GB of memory.We used CPLEX 12.6.2 to solve the LP problems ( 2) and ( 8).We performed the experiments under the same conditions as [7], so the number of LP problems solved successively is fixed at five (K = 5).We also used the same irradiation settings (irradiations from five directions at 72 degrees each, 0 • , 72 • , 144 • , 216 • , and 288 • ).We then calculated the influence matrix D 0 s by the QIB method [20] with the default setting of CERR 3.0 [24].
We used test datasets named TG119 (Task Group 119) provided by the American Association of Physicists in Medicine (AAPM) [8], which contain four datasets: C-Shape, Head and Neck, Prostate, and MultiTarget.In Table 1, from left to right, the columns represent structure names, whether structures are PTVs, the number of voxels in each structure, information on DVCs, the index numbers of DVCs that will be used in the figures below, and the number of beamlets.
As mentioned in Section 3, the proposed method is focused on uncertainty in the influence matrix of the form D s = D 0 s + s • D s with s ∈ Z s .We set the matrix D 0 s as the influence matrix computed with CERR.The range of Z s is determined by the parameter δ, which we vary from 0.1 to 0.5.For preparing D s , let For each (i, j) ∈ P s , we take ε ij from the normal distribution N (0, 1), and we set [D s ] ij = min{max{ε ij , −1/δ}, 1/δ}D 0 s .On the other hand, for (i, j) / ∈ P s , we simply set [D s ] ij = 0. Thus, each element of D s = D 0 s + s • D s is nonnegative for any s ∈ Z s .In the numerical experiments, we vary γ from 0.1 to 1.We set the parameters P α s , P α s , P s , P s in the LP problems (2) and ( 8) as 1 in the manner described in a previous study [7].

Numerical results
We first report the numerical results for the SLPM-R with the simple objective function t of (8a).

C-Shape
Table 2 shows the deviations from the DVCs.The SLPM-R column represents the deviation for the proposed method, whereas the minus and plus columns represent the deviations for the existing method (SLPM) with the worst cases of uncertainty.More precisely, these two deviations correspond to the results of the SLPM with the influence matrices D 0 s − δ|D s | and D 0 s + δ|D s |.These two cases of the SLPM are hereinafter referred as the worst cases.In Table 2, a DVC is satisfied if the corresponding value is non-positive.For example, the value −0.31 for the column of the SLPM-R at (γ , δ) = (0.1, 0.1) with respect to L 0.95 OuterTarget = 50 indicates that more than 95% of the voxels in Outer Target receive 50.31 Gy or more.Thus, the DVC L 0.95 OuterTarget = 50 is satisfied.We use non-positive values to indicate the satisfaction of DVCs, since Property (ii) in Proposition 3.1 indicates that when t k ≤ 0, the obtained solution x k satisfies all DVCs.The bold numbers in Table 2 are used to highlight the best solution among the three columns, more precisely, the solution where the largest deviation from the DVCs is minimum among the SLPM-R, minus and plus columns.The highlighted solution also implies the smallest objective value t in (2a) and (8a).
For the smallest pair of parameters (γ , δ) = (0.1, 0.1), the SLPM-R can find a solution that satisfies all DVCs, but the minus case of the SLPM cannot.For   The right panel of Figure 2 is the result with respect to (γ , δ) = (0.5, 0.5).Due to the larger uncertainty, the ranges sandwiching the worst cases in the existing method turn out to be wider.In the SLPM-R, the deviation from the DVC U 0.1 OuterTarget = 55 becomes slightly worse, but the solid curve still passes near the DVC point, and the largest deviation in the three DVCs is small compared to the SLPM.
Through Table 2, we can also observe that the SLPM-R obtained the better solutions (the solution in the bold face in the table) than the SLPM for most cases.In particular, when δ is small, the SLPM-R is better than the SLPM through all γ .

Head and Neck
In the Head and Neck dataset, the PTV has 53,994 voxels, and this dataset is the largest in TG119.If we construct five successive LP problems with all voxels, then a heavy computation cost is necessary.We randomly selected 10,000 voxels from the PTV and OARs, because a previous study [7] reported that this size reduction does not remarkably affect the computed DVHs for the Head and Neck dataset.(If the reduction has a significant impact, then we may merge voxels by changing the discretization grid or consider re-calculation by focusing on only those areas with the largest differences.)Table 3 reports the deviations from DVCs.Both the SLPM and the SLPM-R find a solution that satisfies all DVCs when γ and δ are small.In contrast, when (γ , δ) = (0.6, 0.1) and (γ , δ) = (0.2, 0.2), the SLPM-R outputs a solution satisfying all DVCs, but the SLPM does not.
In Figure 3, we show two DVHs for parameter pairs (γ , δ) = (0.2, 0.2) and (γ , δ) = (0.8, 0.2).When (γ , δ) = (0.2, 0.2), the solution of the SLPM-R satisfies all DVCs, because the parameters regarding uncertainty are relatively small.On the other hand, when (γ , δ) = (0.8, 0.2), it is difficult for even the SLPM-R to find a feasible solution that satisfies all DVCs, because most components of the influence matrix contain errors.In particular, DVCs of Lt Parotid and Rt Parotid are more severe than those of the PTV and Cord.In the next subsection, we will show that the proposed method with the penalty terms described in Section 3.3 can find a solution.

Prostate
Table 4 shows the deviations from DVCs for the Prostate dataset.When the parameters γ and δ are small, the SLPM-R again outputs a solution satisfying all DVCs.In this dataset, the SLPM-R found better solutions than the SLPM through all γ and δ.
From the table, we can also observe that the lower DVC of Prostate and the DVCs of Rectum are in a trade-off relation.Therefore, it is difficult to find a solution that satisfies all DVCs.In contrast, DVCs of Bladder are easily satisfied.
This dataset requires the highest doses among the four datasets, whereas the number of beamlets is the least, Therefore, uncertainty in the influence matrix strongly affects the results.In fact, when (γ , δ) = (0.5, 0.3) in Figure 4, the deviations of the SLPM-R for the Prostate PTV are at most 3.75 Gy, whereas the largest deviations in the worst cases of the SLPM are 13.29 Gy and 9.43 Gy.

MultiTarget
As shown in Table 1, all three structures in this dataset are PTVs, and each PTV sets lower and upper DVCs.In Table 5, if we focus on the case of (γ , δ) = (0.5, 0.2), the largest deviations for the two worst cases by the SLPM are 7.54 Gy and 6.37 Gy, whereas the largest deviation by the SLPM-R is 6.17 Gy.Therefore, the SLPM-R can reduce the largest deviations for the worst cases.We observe that neither the SLPM nor the SLPM-R can satisfy all six DVCs.The difficulty with this dataset was discussed in a previous study [7].We confirm in Table 5 that Superior and Inferior are strongly affected by uncertainty.Since the irradiation is performed on the coordinate plane of z = 0 and only Center exists on that plane, Superior and Inferior are subject to large errors in uncertainty.Since we discuss FMO with given beam angles in the present paper, this argument is beyond the scope of the present paper, but one resolution would be to conduct irradiations from other angles, such as the coordinate plane of y = 0.
For the MultiTarget dataset, when δ is large, the advantage of the SLPM-R over the SLPM (plus) is not clear.The dose level of 12.5 Gy in the DVC L 0.99 Inferior = 12.5 is the lowest level in Table 1.Note that δ appears in the constraint (8b) in the SLPM-R, that is, j∈J [D 0 s − δ|D s |] ij x j = z si .Therefore, if is close to zero for some some organ i, the flexibility of x j to control z si in the low range around 12.5 Gy is limited, since the same intensity of x j needs to satisfy other DVCs simultaneously.However, the priority of DVCs with low dose can be considered lower than that with high dose, since the effect of beam irradiation to healthy tissues is stronger in the high-dose volumes.In Section 4.2, we will discuss such a priority by using the penalty term introduced in Section 3.3.
We compare the computation times for the SLPM and the SLPM-R.Table 6 shows the entire computation time, which includes the time to construct the input matrices for LP problems, to solve the successive five LP problems (K = 5) and to calculate the DVHs.In particular, solving LP problems occupies the most computation time (more than 90%).Note that the SLPM in [7] does not use γ .For the SLPM-R, the table reports the computation time for each γ .
From the results for C-Shape, adding uncertainty to the original influence matrix does not affect the computation time remarkably, because the computation times with γ = 0.1 and γ = 1.0 are almost the same, and a strong dependence of the computation time with increasing γ cannot be seen.The numbers of variables in the first LP problems in the SLPM and the SLPM-R are 38,549 and 76,680, respectively.The difference of 38,131 corresponds to s∈S {|I s | × (|A s | + |A s |)} due to splitting the variable z si into z si and z si in deriving the robust counterpart (8).Therefore, the number of variables is almost double in the SLPM-R, but Table 2 implies that the increase in computation time is only approximately 170−125 125 ∼ 36%.We can see similar tendencies for the Prostate and MultiTarget data.However, for the Head and Neck data, the SLPM-R is faster than the SLPM.Since the Head and Neck dataset involves the largest number of voxels and the SLPM solves the worst cases, the large deviations from D 0 s might affect the convergence of the interior-point method implemented in CPLEX.
The value of γ affects the number of variables in the LP problems ( 2) and (8) in the first iteration, while δ changes only the coefficients.Therefore, an impact on the computation time due to the change in δ is small compared to that in γ , so we include only δ = 0.2 in Table 6.
At the end of this subsection, we summarize the advantages of the SLPM-R.For the most cases in the three dataset (C-Shape, Head and Neck, and Prostate), the SLPM-R found better solutions than the SLPM as indicated in Tables 2-4.Actually, the largest deviations from DVCs were reduced by the SLPM-R.Though Table 6 indicates a tendency that the SLPM-R demanded a longer computation times than the SLPM, the increase was kept small compared to the increase in the number of variables.On the other hand, as we discussed in the MultiTarget dataset, when the coefficients were much close to zero due to uncertainty, the SLPM-R partially lost the flexibility.

Numerical results with penalty terms
In Table 7, we compare the deviations from DVCs in the results for the simple objective function t in (8a) and the results for the objective function with the penalty term (10) introduced in Section 3.3.Here, we fix the parameters at (γ , δ) = (0.5, 0.2).We use the thresholds of DVCs for parameters θ in (10), more precisely, θ α s = L α s , θ α s = U α s , θ s = L s , and θ s = U s , and we set the penalty weight λ in (10) as in Table 7.We should use large λ to give a high penalty for the deviations from DVCs.However, if λ is set too large, then the relative importance of minimizing t in the objective function becomes small.Thus, the effect of hot and cold spots will also be decreased.In order to examine the effectiveness of the penalty terms, we chose λ based on preliminary experiments in which we employed a grid search as λ ∈ {0.1, 0.2, 0.3, . . ., 0.9, 1.0} considering the deviations in Tables 2-5 and whether the organs are PTV or not.The results we discuss below are bold in Table 7 and they reflect the penalty terms.
We first focus on the results for Head and Neck.When we use the simple objective function (8a), the two worst cases of the SLPM, and even the SLPM-R, cannot satisfy all of the DVCs.In contrast, if we use the penalty term (10), then the SLPM-R can output a solution that satisfies all of the DVCs.In particular, the penalty terms are effective for reducing the deviations from U 0.5 LtParotid and U 0.5 RtParotid .
As for C-Shape, the parameters λ for Outer Target and Core were set to 0.4 and 0.1, respectively, in order to prioritize Outer Target, because the constraints of Outer Target are severe.As a result, the deviations for L 0.5 OuterTarget and U 0.5 OuterTarget were 0.52 Gy and 0.44 Gy in ( 8), but the penalty term in (10) enabled the SLPM-R to find a plan satisfies these DVCs.On the other hand, the deviation from the DVC on Core is worse.
In the Prostate dataset, we also set higher values λ = 0.6 for the Prostate PTV to give a higher priority.We can observe that all of the DVCs except only U 0.10 Rectum are satisfied.
In the MultiTarget dataset, it was difficult to find a solution satisfying all DVCs, as shown in Section 4.1.The three structures are PTVs in this dataset, and we should give a higher priority to a structure with a high threshold, because the change in high-dose areas can be considered to be more important.Thus, we set large parameters λ to the DVCs of Center.In accordance with such parameter settings, the SLPM-R with the penalty terms can reduce the deviations from DVCs on Center from 2.67 Gy and 2.28 Gy to 1.48 Gy and 1.11 Gy.
Compared to the worst cases by the existing method (SLPM), the proposed method reflects the parameters.For example, in the Prostate dataset, we set higher a λ to the PTV.Thus, the proposed method can satisfy the two DVCs on the PTV, whereas the worst cases of the SLPM can satisfy only one of the two DVCs.Similar results can be found in C-Shape and Head and Neck results.

Discussion
Since one advantage of IMRT is the capability of giving more flexibility to the control of beamlets, it is expected that IMRT can improve plans for moving structures, such as lungs.Here, we discuss an extension of our approach to deal with variations induced by breathing.
Chan et al. [10,11] proposed an algorithm to compute beamlet intensities based on an assumption that organs and tumours change their shapes with some uncertainties due to, for example, breathing or tumour motion.Let X be a finite set of lung states, and let P be a set of probability mass functions (PMFs): In other words, p(χ) is the probability that the lung takes a state χ ∈ X .Onak et al. [25] discussed a minimum relative entropy approach to infer a probability density function from noisy input data.Let p ∈ P be a nominal PMF in P. The set of lung states X can be divided into two sets U(⊂ X ) and X \U, where U specifies locations to which the realized PMF is allowed to deviate from the nominal one p.Chan et al. considered the following uncertainty set: P = p ∈ P : p(χ) − p(χ) ≤ p(χ) ≤ p(χ) + p(χ) ∀χ ∈ U, p(χ) = p(χ) ∀χ ∈ X \U , where p and p determine the interval in P, and they computed the absorbed dose of voxel i ∈ I s of an organ s as χ ∈X j∈J ( χ ,s ) ij p(χ)x j , where ( χ ,s ) ij is the (i, j)th element of the influence matrix χ,s for a state χ ∈ X .
We can extend the robust optimization model with the hot and cold spots (8) using the uncertainty of P. In (8), we considered D s as a perturbation to the influence matrix D s .In a similar way, we can introduce χ ,s as a perturbation to χ,s .By noting this, (8) Here, β si is the largest possible negative change caused by p ∈ P. Similarly, β si is the largest possible positive change.
The extended problem (11) is not a standard LP problem because β si and β si are determined by the lower-level LP problems (12) and (13).However, by following a procedure discussed by Bortfeld et al. [26] that exploits the duality theorem on LP, we can reformulate (11) as a standard LP problem.Therefore, (11) can be substantially solved with interior-point methods.
This approach can consider not only the uncertainty contained in the influence matrix, but also the uncertainty contained in the probability of states at the same time.In particular, as seen in Section 4, the advantage of the SLPM-R does not strongly depend on γ when δ is not large.Even if the perturbation χ ,s affects γ for each state s, stable results can be expected.Although this approach requires more variables than the proposed method and a longer computation time, it would derive more practical treatment plans.

Conclusions and future directions
In the present paper, we extended the SLPM with a framework of robust optimization.We mathematically showed that the proposed method maintains the three favourable properties of the SLPM.In particular, when the objective function in the LP problem is non-positive, the proposed method (SLPM-R) can satisfy all of the DVCs, even when the SLPM-R takes the uncertainty in the influence matrix into consideration.
Through numerical experiments, we observed that the SLPM-R provides a solution that reduces the deviation from DVCs compared to the SLPM, which can lead to a more suitable treatment plan.For the C-Shape dataset with parameters (γ , δ) = (0.1, 0.1), the SLPM-R found a solution that satisfies all DVCs, while the SLPM could not.For large parameters (γ , δ) = (0.5, 0.2), though neither the SLPM nor the SLPM-R found a solution, the largest deviation was reduced from 5.65 Gy and 4.41 Gy of the SLPM to 3.15 Gy of the SLPM-R.In the Prostate dataset, the SLPM-R found smaller deviations than the SLPM for all γ and δ.
In contrast, the SLPM-R requires more computation time than the SLPM due to the increased number of variables in LP problems.By introducing the penalty terms in the objective function, we can give priorities to DVCs, and the SLPM-R can find a solution that fulfills all of the DVCs for the Head and Neck dataset.In the MultiTarget dataset, a solution that satisfies all of the DVCs was not found; however, setting λ for the penalty terms enabled the SLPM-R to reduce the derivations from 2.67 Gy to 1.48 Gy.
Regarding future work, there are mainly two directions: development of more practical models and reduction of the computation time.With respect to the former, the proposed method can be applied to the lungs or to multiple influence matrices by assuming that each state changes stochastically, as discussed in Section 5. Furthermore, Chan et al. [11] developed an algorithm for dividing irradiation into small amounts with several steps so that the irradiation can be adjusted in later steps.We may combine this idea with the proposed method.
In regard to the computation time, one way is to accelerate an interior-point method by using the structure when formulating FMO as LP problems in a similar way to Enberg et al. [27].In particular, it may be possible to use the structure defined by the additional variables ζ and z.Arc-search type interior-point methods [28][29][30] can also be considered.

Figure 1 .
Figure 1.Dose volume histogram showing the dose volume constraint (DVCs) and C-VaR constraint.

Algorithm 2
Framework for proposed method: SLPM with robustness (SLPM-R)

Table 1 .
Detailed information on TG119 datasets.

Table 2 .
Deviations from dose volume constraints (DVCs) in C-Shape.

Table 3 .
Deviations from DVCs in Head and Neck.

Table 4 .
Derivations from DVCs in Prostate