Affinely adjustable robust optimization for radiation therapy under evolving data uncertainty via semi-definite programming

Abstract

Static robust optimization has played an important role in radiotherapy, where the decisions aim to safeguard against all possible realizations of uncertainty. However, it may lead to overly conservative decisions or too expensive treatment plans, such as delivering significantly more dose than necessary. Motivated by the success of adjustable robust optimization in reducing highly conservative decision-making of static robust optimization in applications, in this paper, we present an affinely adjustable robust optimization (AARO) model for hypoxia-based radiation treatment planning in the face of evolving data uncertainty. We establish an exact semi-definite program reformulation of the model under a so-called affine decision rule and evaluate our model and approach on a liver cancer case as a proof-of-concept. Our AARO model incorporates uncertainties both in dose influence matrix and re-oxygenation data as well as inexactness of the revealed (re-oxygenation) data. Our numerical experiments demonstrate that the adjustable model successfully handles uncertainty in both re-oxygenation and the dose matrix. They also show that, by utilizing information halfway through the treatment plan, the adjustable solutions of the AARO model outperform a static method while maintaining a similar total dose.

Keywords:

• Affinely adjustable robust optimization
• semi-definite programming
• radiation therapy optimization
• evolving uncertainty

1. Introduction

Radiation therapy is one of the most widely used treatment options for cancer. Long-term side effects can occur if treatments are not carefully designed due to the exposure of healthy tissues surrounding the tumour. The aim of the radiation treatment planning is to deliver curative doses of radiation to tumours while minimizing the risk of side effects [1–6]. The changes in patient's conditions, such as the changes in the patient's cell responses to external radiation when, for instance, the oxygenation is reduced, often result in uncertainty in dose prescription values and inexactness in oxygenation data because of measurement errors. The changes in cell oxygenation directly impact the response to radiation over time, known as the hypoxia effect [7,8].

Robust optimization (RO) is a deterministic approach to decision-making in the face of data uncertainty and is well-suited to handle static decision-making problems in the face of fixed uncertainty at hand [9]. In the (static) RO approach, all decision variables represent ‘here-and-now’ decisions, which means they should obtain specific numerical values as a result of the problem being solved before the actual uncertain parameter values ‘reveal themselves’ [10–12]. A series of works in the late 1990s has shown that robust optimization is a highly successful computationally tractable methodology for solving many classes of optimization problems, such as linear as well as some non-linear problems, under data uncertainty [13]. Since then there have been many further developments to the point where robust optimization has now become a high-potential technique to solve wide-ranging decision-making problems in the face of data uncertainty [14,15].

A recent advance in RO was the development of duality-based conic programming relaxation schemes for solving robust optimization problems, employing the recent numerically tractable convex polynomials [11]. More recently, highly successful first-order methods have been proposed for robust optimization problems that allow larger-scale problems to be solved [16]. For a review of RO and its applications, we refer the reader to [13,17].

Static RO models with fixed uncertainty sets have played important roles in designing radiation treatment plans where a tumour receives a sufficiently high dose to sterilize the tumourous cells while limiting dose exposure to surrounding organs at risk as much as possible to minimize the risk of side effects (see [18] and other references therein). These models aim to safeguard against all possible realizations of uncertainty and may lead to overly conservative decisions or too expensive treatment plans, such as delivering significantly more dose than necessary. Recently, Nohadani and Roy [4] have presented an elegant static robust optimization (RO) approach for a hypoxia-based radiation therapy model with evolving (time-dependent) uncertainty sets and have successfully evaluated their approach and the model on a prostate cancer case. We refer the reader to [19] for a review on robust treatment planning and [20] for the implementation of RO in commercial treatment planning systems for clinical use.

Adjustable robust optimization (ARO), which is an extension of static RO [9,11] to reduce overly conservative decision-making in application, is a deterministic methodology to handle multi-stage optimization problems. In ARO, some of the decision variables are ‘here-and-now ’, while others represent ‘wait-and-see’ decisions, which are assigned numerical values once some of the uncertain parameters have become known [9,21,22]. The significance of ARO is that its worst-case objective value is no worse than the corresponding value of the static RO [10,12,23]. Unfortunately, it is known that solving ARO problems is, in general, computationally very hard. Nevertheless, affinely adjustable robust optimization (AARO) models with fixed recourse that employ the so-called affine decision rules for wait-and-see variables have resulted in computationally tractable reformulations for many classes of functions and uncertainty sets [9,10,13,24–26]. A recent survey on adjustable robust optimization and its applications is available in [27] and extensions of AARO to quadratically adjustable robust optimization may be found in [12,23].

Motivated by the success of adjustable robust optimization over static robust optimization in applications, in this paper, we present an AARO model for hypoxia-based radiation treatment planning in the face of evolving data uncertainty, establish an exact semi-definite program reformulation of the model under a so-called affine decision rule, and evaluate our model and approach on a clinical liver cancer case as an illustration of potential clinical impact.

The main novelty in relation to recent static robust optimization approaches [4,18] is that (a) our AARO model incorporates uncertainties both in the dose-influence matrix and re-oxygenation data and (b) the decision variables on beamlet densities of the model are allowed to depend on the uncertain data. This results in an equivalent conic linear optimization problem with semi-definite as well as second-order cone constraints, rather than a second-order cone program [4] or a linear program [28].¹ Also, our model integrates inexactness of the revealed (re-oxygenation) data in the construction of the uncertainty sets. For recent research on adjustable robust optimization with inexactly revealed data, see [12,24]. Our approach employs semi-definite programming to solve the reformulation of the model under an affine decision rule which is otherwise known to be intractable computationally.

Our contributions to radiation therapy optimization are itemized as follows:

• We present, to the best of our knowledge for the first time, an affinely adjustable robust optimization model for a hypoxia-based radiation therapy planning problems with uncertainties both in dose influence matrix and re-oxygenation data. Our model is a variant of the static RO model by Nohadani and Roy [4] in the sense that the decision variables on beamlet densities are allowed to be affinely adjustable with respect to the uncertain parameter. Also, our model employs uncertainty sets which not only capture the evolving uncertain conditions but also incorporates inexactness of the revealed (re-oxygenation) data (see Section 3 for further details). Interested readers are directed to [3] for a more recent AARO approach to treatment length-optimization in radiation therapy, accounting for inexactness of biomarker information.

• We establish an equivalence between the AARO model and its computationally tractable conic linear programming problem reformulation of the AARO model under an affine decision rule. This is done by invoking the celebrated S-lemma [29,30] to convert the numerically intractable robust constraints into an equivalent system of semi-definite constraints and second-order cone constraints, following the standard approach pioneered by Ben-Tal et al. [9]. The reformulation may permit one to solve computationally expensive treatment planning models of clinical cancer cases because these conic models can be efficiently solved via commonly used semi-definite program software such as MOSEK [31].

• Finally, we implement a simplified form of our AARO model with dose influence matrix and re-oxygenation uncertainties on a liver cancer case as proof-of-concept. Our numerical experiments demonstrate that the adjustable model successfully handles uncertainty in both re-oxygenation and the dose matrix. They also show that, by utilizing information halfway through the treatment plan, the adjustable solutions of the AARO model outperform a static method while maintaining a similar total dose. The results are also within the comparable performance of the hypoxia linear program model with nominal data.

The paper is organized as follows. Section 2 presents preliminaries associated with the hypoxia-based radiation treatment planning optimization. Section 3 derives both the affinely adjustable robust optimization model and the corresponding static model incorporating both the dose influence matrix and re-oxygenation data uncertainties. Section 4 establishes the equivalence between the AARO model and its computationally tractable conic linear program reformulation of the AARO model under an affine decision rule. Section 5 describes the implementation of our AARO and static models, and their results in application to a liver cancer case. Appendix 1 provides the dose-volume histograms that are commonly used in treatment plan evaluations; Appendix 2 gives the simplified conic formulation for AARO model that is employed in the liver cancer case; Appendix 3 gives the equivalent second-order cone program reformulation for the static robust model.

2. Preliminaries on radiation treatment planning optimization

We begin this section by explaining the process and principle of radiotherapy, fixing the notation and preliminaries, and by recalling some recent successful developments of radiation therapy modelling such as the ones examined in [4,18,32]. The radiotherapy planning process involves finding the minimum effective dosage to deliver a prescribed amount of radiation to the target tumour region by solving an optimization problem to determine the dose intensity of every beamlet.

The patient interior organ structure, such as organs at risk, tumour targets and normal tissue, is mapped by voxels of a fixed size. These voxels then have constraints assigned to them in the optimization problem, creating upper and lower bounds on radiation dosage. Note that a lack of oxygen reduces the responsiveness of the tumour to the radiation, and this condition is known as hypoxia, referred sometimes as radiosensitivity. It directly affects the outcome of radiation treatment. Based on these principles of radiotherapy, we first present a basic hypoxia-based robust optimization model for radiation therapy planning.

We assume that the set of all voxels of interest, denoted by $\mathcal{S}$, can be logically partitioned into three mutually exclusive sets, $\mathcal{T}\cup \mathcal{C}\cup \mathcal{N}$, where $\mathcal{T}$ represents the subset of voxels that are within the target (or tumour), $\mathcal{C}$ represents the set of critical organ voxels that are not in $\mathcal{T}$ and $\mathcal{N}$ represents the set of normal tissue voxels, which is the subset of out-of-target voxels and also not in the set of critical organ voxels. Note that the set of hypoxia tumour voxels is denoted by $\mathcal{H}$, which is a subset of $\mathcal{T}$.

The basic optimization model in radiation therapy often takes the following form where the weighted sum of doses delivered to voxels is minimized subject to constraints (e.g. see [18,32,33]): $\begin{array}{rl}\underset{w=(w_p{)}_{p=1}^n\in {\mathbb{R}}^n}{min}& \sum _{p=1}^n\sum _{(\alpha ,\beta )\in \mathcal{S}}{w}_p{D}_{\mathrm{\alpha \beta }}^p\\ \text{subject to}& w_p\ge 0,\\ \text{(Tumour constraint)}& \sum _{p=1}^n{w}_p{D}_{\mathrm{\alpha \beta }}^p\ge {\gamma }_{\mathrm{lb}},\text{for all }(\alpha ,\beta )\in \mathcal{T},\\ \text{(Critical organ constraint)}& \sum _{p=1}^n{w}_p{D}_{\mathrm{\alpha \beta }}^p\le {\gamma }_{\mathrm{ub}}^c,\text{for all }(\alpha ,\beta )\in \mathcal{C},\\ \text{(Normal tissue constraint)}& \sum _{p=1}^n{w}_p{D}_{\mathrm{\alpha \beta }}^p\le {\gamma }_{\mathrm{ub}}^n,\text{for all }(\alpha ,\beta )\in \mathcal{N},\\ \text{(Hypoxic region constraint)}& \sum _{p=1}^n{w}_p{D}_{\mathrm{\alpha \beta }}^p\ge {\rho }_{\mathrm{\alpha \beta }}{\gamma }_{\mathrm{lb}},\text{for all }(\alpha ,\beta )\in \mathcal{H},\end{array}$where the decision variable $w=(w_p{)}_{p=1}^n$ is the vector of beamlet density and $D_{\mathrm{\alpha \beta }}^p$ represents the dose delivered to voxel $(\alpha ,\beta )$ by a unit weight of beamlet p. Note that the constraints describe lower (lb) or upper (ub) prescribed dose limits ${\gamma }_{\mathrm{lb}}$, ${\gamma }_{\mathrm{ub}}^c$ and ${\gamma }_{\mathrm{ub}}^n$ to target voxels, critical organ voxels and normal tissue voxels, respectively. The hypoxia constraint is modified by a factor to account for the cell response to the dose. So, a factor ${\rho }_{\mathrm{\alpha \beta }}$ is imposed on the clinically prescribed dose for any hypoxia tumour voxels, that is, for all $(\alpha ,\beta )\in \mathcal{H}$.

The radiation treatments are often conducted over an extended time period (typically several days or weeks) in order to increase the recovery probability of healthy tissue, and the clinical observations are generally measured in increments. We consider an optimization model which takes into account of reoxygenation over a given time period. Also, note that the effect of hypoxia changes over time also due to reoxygenation.

In the interest of simplifying the descriptions of the model, we will assume throughout that the pre-fixed length between each observation is one unit where, for instance, one unit may correspond to one day or one week. Denote the total length of the treatment time by T units.

Now, we consider the following hypoxia-based model, which is a slight variant of the model examined recently in [4]: $\begin{array}{rl}\underset{w_p^t\in \mathbb{R},z^t\in \mathbb{R}}{min}& \sum _{t=0}^T\sum _{(\alpha ,\beta )\in \mathcal{S}}\sum _{p=1}^n{w}_p^t{D}_{\mathrm{\alpha \beta }}^{p,t}\\ \text{subject to}& \sum _{t=0}^T\sum _{p=1}^n{w}_p^t{D}_{\mathrm{\alpha \beta }}^{p,t}\ge {\gamma }_{\mathrm{lb}},\text{for all }(\alpha ,\beta )\in \mathcal{T},\\ & \sum _{t=0}^T\sum _{p=1}^n{w}_p^t{D}_{\mathrm{\alpha \beta }}^{p,t}\le {\gamma }_{\mathrm{ub}}^c,\text{for all }(\alpha ,\beta )\in \mathcal{C},\\ & \sum _{t=0}^T\sum _{p=1}^n{w}_p^t{D}_{\mathrm{\alpha \beta }}^{p,t}\le {\gamma }_{\mathrm{ub}}^n,\text{for all }(\alpha ,\beta )\in \mathcal{N},\\ & \sum _{p=1}^n{w}_p^t{D}_{\mathrm{\alpha \beta }}^{p,t}\ge {\rho }_{\mathrm{\alpha \beta }}^t{z}^t,\text{for all }(\alpha ,\beta )\in \mathcal{H},t=0,1,\dots ,T,\\ & \sum _{t=0}^T{z}^t\ge {\gamma }_{\mathrm{lb}},\\ & w_p^t\ge 0,t=0,1,\dots ,T,p=1,\dots ,n,\\ & z^t\ge 0,t=0,1,\dots ,T,\end{array}$where $w_p^t$ denotes the density of the beamlet p at the time t, and the real number $D_{\mathrm{\alpha \beta }}^{p,t}\in \mathbb{R}$ represents the given dose matrix entry at voxel $(\alpha ,\beta )$ for beamlet p at the time t, where $p=1,\dots ,n$ and $t=0,\dots ,T$. Recall that t = 0 means the starting time and t = 1 means one unit since the starting time.

Here, the variables $z^t\in \mathbb{R}$ represent the fraction of the total dosage delivered to the hypoxia tumour voxels at the time t, $t=0,\dots ,T$. This is necessary since each fraction is affected by its own hypoxia coefficient ${\rho }_{\mathrm{\alpha \beta }}^t$, and so cannot be optimally split into, say, equal dosages, and instead needs to be determined as part of the model.

Our model is a variant to the one considered in [4] in the sense that we introduce new variables $z^t$ to allow us to determine the optimal fraction delivered at the time t with $t=0,\dots ,T$ (instead of fixing them as known equal percentages). Notice that the pre-fixed equal percentages choice of $z^t$ used in [4] is a feasible solution for the above model. Thus, our model produces a solution with a smaller amount of total dose of radiation emitted to the whole tissues (that is, a better objective value).

3. Adjustable and static robust optimization models

In this section, we present adjustable as well as static hypoxia-based robust optimization models. Recall that the elements of dose influence matrix and the hypoxia factor ${\rho }_{\mathrm{\alpha \beta }}^t$ are pre-determined by clinicians or physicists. They are often uncertain due to prediction error and patient's condition. To model this, we assume that the dose matrix data is uncertain and satisfies $D_{\mathrm{\alpha \beta }}^{p,t}={\overline{d}}_{\mathrm{\alpha \beta }}^{p,t}+{\overline{v}}_{\mathrm{\alpha \beta }}^{T}u,t=0,\dots ,T,p=1,\dots ,n,(\alpha ,\beta )\in \mathcal{S},$where ${\overline{d}}_{\mathrm{\alpha \beta }}^{p,t}\in \mathbb{R}$ is the nominal dose matrix entry at voxel $(\alpha ,\beta )$ for beamlet p at the time t with $p=1,\dots ,n$ and $t=0,\dots ,T$; ${\overline{v}}_{\mathrm{\alpha \beta }}\in {\mathbb{R}}^q$ and q is a positive integer; $u\in {\mathbb{R}}^q$ is the uncertain parameter and $\Vert u\Vert \le r$ with r>0. Note that $\Vert u\Vert$ denotes the Euclidean norm of u, that is, $\Vert u\Vert =\sqrt{u^{T}u}$. Recall also that, for two convex sets A and B, the Minkowski sum of A and B is denoted as A + B and is defined as $A+B=\{a+b:a\in A,b\in B\}$.

In this paper, we assume that the hypoxia discount parameter ${\rho }_{\mathrm{\alpha \beta }}^t$ in the re-oxygenation process is also uncertain with ${\rho }_{\mathrm{\alpha \beta }}^t$ belonging to the evolving uncertainty set ${\mathcal{U}}_{\mathrm{\alpha \beta }}^t$ with inexactly revealed data, defined as (1) $\begin{array}{rl}& {\mathcal{U}}_{\mathrm{\alpha \beta }}^t\\ & =\{\begin{array}{ll}\{{\rho }_{\mathrm{\alpha \beta }}^t:|{\rho }_{\mathrm{\alpha \beta }}^t-(({\rho }_0{)}_{\mathrm{\alpha \beta }}+\eta t)|\le \mathrm{\Gamma t}\}& \text{if }t<t_k,\\ \{{\rho }_{\mathrm{\alpha \beta }}^t:|{\rho }_{\mathrm{\alpha \beta }}^t-({\rho }_{\mathrm{\alpha \beta }}^{\mathrm{observed}}+\eta (t-t_k))|\le \mathrm{\Gamma }(t-t_k)\}+[-\nu ,\nu ]& \text{if }t\ge t_k,\end{array}\end{array}$(1) where $({\rho }_0{)}_{\mathrm{\alpha \beta }}\in \mathbb{R}$, ${\rho }_{\mathrm{\alpha \beta }}^{\mathrm{observed}}\in \mathbb{R}$ and $\eta ,\mathrm{\Gamma },\nu \ge 0$ are given, and $t_k$ is the time where the intermediate observation is made. Note that, when $t<t_k$, the uncertainty set takes the form ${\mathcal{U}}_{\mathrm{\alpha \beta }}^t=\{{\rho }_{\mathrm{\alpha \beta }}^t:|{\rho }_{\mathrm{\alpha \beta }}^t-(({\rho }_0{)}_{\mathrm{\alpha \beta }}+\eta t)|\le \mathrm{\Gamma t}\}.$At the observation time $t_k$, the available inexactly revealed data, ${\rho }_{\mathrm{\alpha \beta }}^{\mathrm{observed}}$ for the hypoxia parameter, is used as estimate for the true value of the unknown hypoxia parameter. The uncertainty set after the intermediate observation (that is, when $t\ge t_k$) is defined as ${\mathcal{U}}_{\mathrm{\alpha \beta }}^t=\{{\rho }_{\mathrm{\alpha \beta }}^t:|{\rho }_{\mathrm{\alpha \beta }}^t-({\rho }_{\mathrm{\alpha \beta }}^{\mathrm{observed}}+\eta (t-t_k))|\le \mathrm{\Gamma }(t-t_k)\}+[-\nu ,\nu ],$where $\nu \ge 0$ controls the level of inexactness of the observation at the time $t_k$ (see Figure 1 below for an illustration). We note that, in the case where the observation is exact (that is, $\nu =0$) and the dose matrix is free of uncertainty, such an evolving data uncertainty set with exactly revealed observation data has been introduced and used in [4] for static robust optimization models.

Figure 1. Uncertainty sets with inexactness in hypoxia observation data.

The construction of our uncertainty set ${\mathcal{U}}_{\mathrm{\alpha \beta }}^t$ illustrated in Figure 1 was inspired by [4]. The right-hand side of the figure highlights the inexactness of the observed data, whereas in [4] observation data were assumed to be exact.

Affine decision rule. Note that, in contrast to [4], we allow the decision variable $w_p^t$ on beamlet density to be a function of the uncertain parameter $u\in {\mathbb{R}}^q$, and assume that $w_p^t$ is adjustable and satisfies the following affine decision rule: $w_p^t(u)={\eta }_p^t+(a_p^t{)}^{T}u,t=0,1,\dots ,T,p=1,\dots ,n,$where ${\eta }_p^t\in \mathbb{R}$ and $a_p^t\in {\mathbb{R}}^q$; while the decision variables $z^t$, $t=0,\dots ,T$, are non-adjustable.

Then, the corresponding affinely adjustable robust optimization model under dose and re-oxygenation uncertainties can be written as

(2) $\begin{array}{rl}(P_{\mathrm{AARO}})\underset{{\eta }_p^t\in \mathbb{R},a_p^t\in {\mathbb{R}}^q,z^t\in \mathbb{R}}{min}& \underset{\Vert u\Vert \le r}{max}\{\sum _{t=0}^T\sum _{(\alpha ,\beta )\in \mathcal{S}}\sum _{p=1}^n{w}_p^t(u)[{\overline{d}}_{\mathrm{\alpha \beta }}^{p,t}+{\overline{v}}_{\mathrm{\alpha \beta }}^{T}u]\}\\ \text{subject to}& \sum _{t=0}^T\sum _{p=1}^n{w}_p^t(u)[{\overline{d}}_{\mathrm{\alpha \beta }}^{p,t}+{\overline{v}}_{\mathrm{\alpha \beta }}^{T}u]\ge {\gamma }_{\mathrm{lb}},\mathrm{\forall }\Vert u\Vert \le r,(\alpha ,\beta )\in \mathcal{T},\\ & \sum _{t=0}^T\sum _{p=1}^n{w}_p^t(u)[{\overline{d}}_{\mathrm{\alpha \beta }}^{p,t}+{\overline{v}}_{\mathrm{\alpha \beta }}^{T}u]\le {\gamma }_{\mathrm{ub}}^c,\mathrm{\forall }\Vert u\Vert \le r,(\alpha ,\beta )\in \mathcal{C},\\ & \sum _{t=0}^T\sum _{p=1}^n{w}_p^t(u)[{\overline{d}}_{\mathrm{\alpha \beta }}^{p,t}+{\overline{v}}_{\mathrm{\alpha \beta }}^{T}u]\le {\gamma }_{\mathrm{ub}}^n,\mathrm{\forall }\Vert u\Vert \le r,(\alpha ,\beta )\in \mathcal{N},\\ & \sum _{p=1}^n{w}_p^t(u)[{\overline{d}}_{\mathrm{\alpha \beta }}^{p,t}+{\overline{v}}_{\mathrm{\alpha \beta }}^{T}u]\ge {\rho }_{\mathrm{\alpha \beta }}^t{z}^t,\\ & \mathrm{\forall }{\rho }_{\mathrm{\alpha \beta }}^t\in {\mathcal{U}}_{\mathrm{\alpha \beta }}^t,\mathrm{\forall }\Vert u\Vert \le r,(\alpha ,\beta )\in \mathcal{H},t=0,\dots ,T,\\ & \sum _{t=0}^T{z}^t\ge {\gamma }_{\mathrm{lb}},z^t\ge 0,t=0,1,\dots ,T,\\ & w_p^t(u)={\eta }_p^t+(a_p^t{)}^{T}u,t=0,1,\dots ,T,\text{(Affine Decision Rule)}\\ & w_p^t(u)\ge 0,\mathrm{\forall }\Vert u\Vert \le r.\end{array}$(2) As we will see later (see Section 4), this affinely adjustable robust model problem $(P_{\mathrm{AARO}})$ can be equivalently reformulated as a conic linear program with semi-definite constraints and second order cone constraints, and so, can be efficiently solved via commonly used semi-definite program solver such as MOSEK [31].

Static robust model under dose and re-oxygenation uncertainties. In the static case, that is, the decision variable $w_p^t$ is uncertainty free, one has $w_p^t(u)={\eta }_p^t$. Then, the corresponding affinely adjustable robust optimization simplified as

$\begin{array}{rl}(P_{\mathrm{Static}})\underset{{\eta }_p^t\in \mathbb{R},z^t\in \mathbb{R}}{min}& \underset{\Vert u\Vert \le r}{max}\{\sum _{t=0}^T\sum _{(\alpha ,\beta )\in \mathcal{S}}\sum _{p=1}^n{\eta }_p^t[{\overline{d}}_{\mathrm{\alpha \beta }}^{p,t}+{\overline{v}}_{\mathrm{\alpha \beta }}^{T}u]\}\\ \text{subject to }& {\displaystyle \sum _{t=0}^T\sum _{p=1}^n{\eta }_p^t[{\overline{d}}_{\mathrm{\alpha \beta }}^{p,t}+{\overline{v}}_{\mathrm{\alpha \beta }}^{T}u]\ge {\gamma }_{\mathrm{lb}},\mathrm{\forall }\Vert u\Vert \le r,(\alpha ,\beta )\in \mathcal{T},}\\ & {\displaystyle \sum _{t=0}^T\sum _{p=1}^n{\eta }_p^t[{\overline{d}}_{\mathrm{\alpha \beta }}^{p,t}+{\overline{v}}_{\mathrm{\alpha \beta }}^{T}u]\le {\gamma }_{\mathrm{ub}}^c,\mathrm{\forall }\Vert u\Vert \le r,(\alpha ,\beta )\in \mathcal{C},}\\ & \sum _{t=0}^T\sum _{p=1}^n{\eta }_p^t[{\overline{d}}_{\mathrm{\alpha \beta }}^{p,t}+{\overline{v}}_{\mathrm{\alpha \beta }}^{T}u]\le {\gamma }_{\mathrm{ub}}^n,\mathrm{\forall }\Vert u\Vert \le r,(\alpha ,\beta )\in \mathcal{N},\\ & {\displaystyle \sum _{p=1}^n{\eta }_p^t[{\overline{d}}_{\mathrm{\alpha \beta }}^{p,t}+{\overline{v}}_{\mathrm{\alpha \beta }}^{T}u]\ge {\rho }_{\mathrm{\alpha \beta }}^t{z}^t,}\\ & \mathrm{\forall }{\rho }_{\mathrm{\alpha \beta }}^t\in {\mathcal{U}}_{\mathrm{\alpha \beta }}^t,\mathrm{\forall }\Vert u\Vert \le r,(\alpha ,\beta )\in \mathcal{H},\\ & {\displaystyle \sum _{t=0}^T{z}^t\ge {\gamma }_{\mathrm{lb}},{\eta }_p^t\ge 0,z^t\ge 0,t=0,\dots ,T.}\end{array}$

The static model $(P_{\mathrm{Static}})$ is a semi-infinite optimization problem which is, in general, not numerically tractable. Using a standard approach as in [9], this model can be equivalently reformulated as a second-order cone program (SOCP), which can efficiently be solved by interior point methods. This SOCP reformulation will also be used in Section 5 to study the liver cancer case. For completeness, the equivalent reformulation as an SOCP for the static problem $(P_{\mathrm{Static}})$ is given in Appendix 3.

For the corresponding static robust model with only re-oxygenation uncertainty (that is, there is no dose matrix uncertainty), ${\overline{v}}_{\mathrm{\alpha \beta }}=0$ for all $(\alpha ,\beta )$. So, the model can be further simplified as the following linear program:

$\begin{array}{rl}(P_{\mathrm{hypoxia}\mathrm{LP}})\underset{{\eta }_p^t\in \mathbb{R},z^t\in \mathbb{R}}{min}& \sum _{t=0}^T\sum _{(\alpha ,\beta )\in \mathcal{S}}\sum _{p=1}^n{\eta }_p^t{\overline{d}}_{\mathrm{\alpha \beta }}^{p,t}\\ \text{subject to}& {\displaystyle \sum _{t=0}^T\sum _{p=1}^n{\eta }_p^t{\overline{d}}_{\mathrm{\alpha \beta }}^{p,t}\ge {\gamma }_{\mathrm{lb}},\text{for all }(\alpha ,\beta )\in \mathcal{T},}\\ & {\displaystyle \sum _{t=0}^T\sum _{p=1}^n{\eta }_p^t{\overline{d}}_{\mathrm{\alpha \beta }}^{p,t}\le {\gamma }_{\mathrm{ub}}^c,\text{for all }(\alpha ,\beta )\in \mathcal{C},}\\ & {\displaystyle \sum _{t=0}^T\sum _{p=1}^n{\eta }_p^t{\overline{d}}_{\mathrm{\alpha \beta }}^{p,t}\le {\gamma }_{\mathrm{ub}}^n,\text{for all }(\alpha ,\beta )\in \mathcal{N},}\\ & {\displaystyle \sum _{p=1}^n{\eta }_p^t{\overline{d}}_{\mathrm{\alpha \beta }}^{p,t}\ge [({\rho }_0{)}_{\mathrm{\alpha \beta }}+\mathrm{\eta t}+\mathrm{\Gamma t}]z^t,}\\ & \text{for all }(\alpha ,\beta )\in \mathcal{H},t=0,\dots ,t_k-1,\\ & {\displaystyle \sum _{p=1}^n{\eta }_p^t{\overline{d}}_{\mathrm{\alpha \beta }}^{p,t}\ge [{\rho }_{\mathrm{\alpha \beta }}^{\mathrm{observed}}+\eta (t-t_k)+\mathrm{\Gamma }(t-t_k)+\nu ]z^t,}\\ & \text{for all }(\alpha ,\beta )\in \mathcal{H},t=t_k,\dots ,T,\\ & {\displaystyle \sum _{t=0}^T{z}^t\ge {\gamma }_{\mathrm{lb}},{\eta }_p^t\ge 0,z^t\ge 0,t=0,\dots ,T.}\end{array}$

4. Conic LP reformulation of radiation therapy AARO model

In this section, we establish an equivalent conic linear program reformulation for the affinely adjustable robust optimization model $(P_{\mathrm{AARO}})$. The model $(P_{\mathrm{AARO}})$ is a semi-infinite optimization problem which is, in general, not numerically tractable, whereas the conic linear program reformulation can efficiently be solved by commonly available software.

To do this, consider the problem $(P_{\mathrm{AARO}})$ given in (2(2) $\begin{array}{rl}(P_{\mathrm{AARO}})\underset{{\eta }_p^t\in \mathbb{R},a_p^t\in {\mathbb{R}}^q,z^t\in \mathbb{R}}{min}& \underset{\Vert u\Vert \le r}{max}\{\sum _{t=0}^T\sum _{(\alpha ,\beta )\in \mathcal{S}}\sum _{p=1}^n{w}_p^t(u)[{\overline{d}}_{\mathrm{\alpha \beta }}^{p,t}+{\overline{v}}_{\mathrm{\alpha \beta }}^{T}u]\}\\ \text{subject to}& \sum _{t=0}^T\sum _{p=1}^n{w}_p^t(u)[{\overline{d}}_{\mathrm{\alpha \beta }}^{p,t}+{\overline{v}}_{\mathrm{\alpha \beta }}^{T}u]\ge {\gamma }_{\mathrm{lb}},\mathrm{\forall }\Vert u\Vert \le r,(\alpha ,\beta )\in \mathcal{T},\\ & \sum _{t=0}^T\sum _{p=1}^n{w}_p^t(u)[{\overline{d}}_{\mathrm{\alpha \beta }}^{p,t}+{\overline{v}}_{\mathrm{\alpha \beta }}^{T}u]\le {\gamma }_{\mathrm{ub}}^c,\mathrm{\forall }\Vert u\Vert \le r,(\alpha ,\beta )\in \mathcal{C},\\ & \sum _{t=0}^T\sum _{p=1}^n{w}_p^t(u)[{\overline{d}}_{\mathrm{\alpha \beta }}^{p,t}+{\overline{v}}_{\mathrm{\alpha \beta }}^{T}u]\le {\gamma }_{\mathrm{ub}}^n,\mathrm{\forall }\Vert u\Vert \le r,(\alpha ,\beta )\in \mathcal{N},\\ & \sum _{p=1}^n{w}_p^t(u)[{\overline{d}}_{\mathrm{\alpha \beta }}^{p,t}+{\overline{v}}_{\mathrm{\alpha \beta }}^{T}u]\ge {\rho }_{\mathrm{\alpha \beta }}^t{z}^t,\\ & \mathrm{\forall }{\rho }_{\mathrm{\alpha \beta }}^t\in {\mathcal{U}}_{\mathrm{\alpha \beta }}^t,\mathrm{\forall }\Vert u\Vert \le r,(\alpha ,\beta )\in \mathcal{H},t=0,\dots ,T,\\ & \sum _{t=0}^T{z}^t\ge {\gamma }_{\mathrm{lb}},z^t\ge 0,t=0,1,\dots ,T,\\ & w_p^t(u)={\eta }_p^t+(a_p^t{)}^{T}u,t=0,1,\dots ,T,\text{(Affine Decision Rule)}\\ & w_p^t(u)\ge 0,\mathrm{\forall }\Vert u\Vert \le r.\end{array}$(2) ). We associate it with the following conic linear optimization problem: $\begin{array}{rl}(P_{\mathrm{AARO\_SDP}})& \underset{\begin{array}{c}s\in \mathbb{R},{\eta }_p^t\in \mathbb{R},a_p^t\in {\mathbb{R}}^q,z^t\in \mathbb{R}\\ {\lambda }_1,{\lambda }_2^{\alpha ,\beta },{\lambda }_3^{\alpha ,\beta },{\lambda }_4^{\alpha ,\beta }\ge 0\\ {\lambda }_{5,t}^{\alpha ,\beta }\ge 0\end{array}}{min}s\\ \text{subject to}& {\eta }_p^t-r\Vert a_p^t\Vert \ge 0,t=0,1,\dots ,T\\ & \left(\begin{array}{cc}{\displaystyle -\sum _{(\alpha ,\beta )\in \mathcal{S}}{W}_{\mathrm{\alpha \beta }}+{\lambda }_1{I}_q}& {\displaystyle -\sum _{(\alpha ,\beta )\in \mathcal{S}}{w}_{\mathrm{\alpha \beta }}/2}\\ {\displaystyle -\sum _{(\alpha ,\beta )\in \mathcal{S}}{w}_{\mathrm{\alpha \beta }}^T/2}& {\displaystyle -\sum _{(\alpha ,\beta )\in \mathcal{S}}{\theta }_{\mathrm{\alpha \beta }}+s-{\lambda }_1{r}^2}\end{array}\right)⪰0\\ & \left(\begin{array}{cc}{W}_{\mathrm{\alpha \beta }}+{\lambda }_2^{\alpha ,\beta }{I}_q& w_{\mathrm{\alpha \beta }}/2\\ w_{\mathrm{\alpha \beta }}^T/2& {\theta }_{\mathrm{\alpha \beta }}-{\gamma }_{\mathrm{lb}}-{\lambda }_2^{\alpha ,\beta }{r}^2\end{array}\right)⪰0,\\ & \text{for all }(\alpha ,\beta )\in \mathcal{T},\\ & \left(\begin{array}{cc}-W_{\mathrm{\alpha \beta }}+{\lambda }_3^{\alpha ,\beta }{I}_q& -w_{\mathrm{\alpha \beta }}/2\\ -w_{\mathrm{\alpha \beta }}^T/2& -{\theta }_{\mathrm{\alpha \beta }}+{\gamma }_{\mathrm{ub}}^c-{\lambda }_3^{\alpha ,\beta }{r}^2\end{array}\right)⪰0,\\ & \text{for all }(\alpha ,\beta )\in \mathcal{C},\\ & \left(\begin{array}{cc}-W_{\mathrm{\alpha \beta }}+{\lambda }_4^{\alpha ,\beta }{I}_q& -w_{\mathrm{\alpha \beta }}/2\\ -w_{\mathrm{\alpha \beta }}^T/2& -{\theta }_{\mathrm{\alpha \beta }}+{\gamma }_{\mathrm{ub}}^n-{\lambda }_4^{\alpha ,\beta }{r}^2\end{array}\right)⪰0,\\ & \text{for all }(\alpha ,\beta )\in \mathcal{N},\\ & \left(\begin{array}{cc}{W}_{\mathrm{\alpha \beta }}^t+{\lambda }_{5,t}^{\alpha ,\beta }{I}_q& w_{\mathrm{\alpha \beta }}^t/2\\ (w_{\mathrm{\alpha \beta }}^t{)}^T/2& {\theta }_{\mathrm{\alpha \beta }}^t-[({\rho }_0{)}_{\mathrm{\alpha \beta }}+\mathrm{\eta t}+\mathrm{\Gamma t}]z^t-{\lambda }_{5,t}^{\alpha ,\beta }{r}^2\end{array}\right)⪰0,\\ & \text{for all }(\alpha ,\beta )\in \mathcal{H},t=0,1,\dots ,t_k-1\\ & \left(\begin{array}{cc}{W}_{\mathrm{\alpha \beta }}^t+{\lambda }_{5,t}^{\alpha ,\beta }{I}_q& w_{\mathrm{\alpha \beta }}^t/2\\ (w_{\mathrm{\alpha \beta }}^t{)}^T/2& {\theta }_{\mathrm{\alpha \beta }}^t-[{\rho }_{\mathrm{\alpha \beta }}^{\mathrm{observed}}+\eta (t-t_k)\\ & +\mathrm{\Gamma }(t-t_k)+\nu ]z^t-{\lambda }_{5,t}^{\alpha ,\beta }{r}^2\end{array}\right)⪰0,\\ & \text{for all }(\alpha ,\beta )\in \mathcal{H},t=t_k,\dots ,T\\ & {\displaystyle {\theta }_{\mathrm{\alpha \beta }}=\sum _{t=0}^T\sum _{p=1}^n{\eta }_p^t{\overline{d}}_{\mathrm{\alpha \beta }}^{p,t},w_{\mathrm{\alpha \beta }}=\sum _{t=0}^T\sum _{p=1}^n({\overline{d}}_{\mathrm{\alpha \beta }}^{p,t}{a}_p^t+{\eta }_p^t{\overline{v}}_{\mathrm{\alpha \beta }}),}\\ & {\displaystyle W_{\mathrm{\alpha \beta }}=\sum _{t=0}^T\sum _{p=1}^n\frac{a_p^t{\overline{v}}_{\mathrm{\alpha \beta }}^T+{\overline{v}}_{\mathrm{\alpha \beta }}(a_p^t{)}^T}{2},}\\ & {\displaystyle {\theta }_{\mathrm{\alpha \beta }}^t=\sum _{p=1}^n{\eta }_p^t{\overline{d}}_{\mathrm{\alpha \beta }}^{p,t},w_{\mathrm{\alpha \beta }}^t=\sum _{p=1}^n({\overline{d}}_{\mathrm{\alpha \beta }}^{p,t}{a}_p^t+{\eta }_p^t{\overline{v}}_{\mathrm{\alpha \beta }}),}\\ & W_{\mathrm{\alpha \beta }}^t=\sum _{p=1}^n\frac{a_p^t{\overline{v}}_{\mathrm{\alpha \beta }}^T+{\overline{v}}_{\mathrm{\alpha \beta }}(a_p^t{)}^T}{2},\\ & {\displaystyle \sum _{t=0}^T{z}^t\ge {\gamma }_{\mathrm{lb}},z^t\ge 0.}\end{array}$ Note that, for a symmetric matrix A, $A⪰0$ means that A is positive semi-definite. Also, $I_q$ denotes the $(q\times q)$ identity matrix.

Now, we show that the adjustable robust optimization problem $(P_{\mathrm{AARO}})$ can be equivalently reformulated as $(P_{\mathrm{AARO\_SDP}})$, and one can find a solution of $(P_{\mathrm{AARO}})$ by solving $(P_{\mathrm{AARO\_SDP}})$. It is worth noting that $(P_{\mathrm{AARO\_SDP}})$ is a conic linear optimization problem with second-order cone and linear matrix constraints, and so, can be efficiently solved by commonly used semi-definite program solver such as MOSEK [31]. Note also that the AARO models with the affine decision rule will lead to reformulations involving semi-definite constraints rather than second-order cone constraints or linear constraints as in [4,28].

Let us recall the following two useful lemmas which play key roles in reformulating $(P_{\mathrm{AARO}})$ as a conic linear program.

Lemma 4.1

S-lemma [29,30]

Let f and g be quadratic functions on ${\mathbb{R}}^q$. Suppose that $\{u\in {\mathbb{R}}^q:g(u)<0\}\ne \mathrm{\varnothing }$. Then, the following statements (i) and (ii) are equivalent:

1. $g(u)\le 0\Rightarrow f(u)\ge 0;$

2. there exists ${\lambda }_1\ge 0$ such that $f(u)+{\lambda }_{1}g(u)\ge 0\text{ for all }u\in {\mathbb{R}}^q.$

Lemma 4.2

Representation of non-negativity for quadratic functions [9]

Let h be a quadratic function on ${\mathbb{R}}^q$ such that $h(u)=u^T\mathrm{Mu}+b^{T}u+\gamma$, where M is a symmetric $(q\times q)$ matrix, $b\in {\mathbb{R}}^q$ and $\gamma \in \mathbb{R}$. Then $h(u)\ge 0$ for all $u\in {\mathbb{R}}^q$ if and only if $\left(\begin{array}{cc}M& b/2\\ b^T/2& \gamma \end{array}\right)⪰0.$

Theorem 4.1

Exact conic formulation for adjustable robust optimization problem

Consider the affinely adjustable robust optimization problem $(P_{\mathrm{AARO}})$ and the conic linear program $(P_{\mathrm{AARO\_SDP}})$. Then, $min(P_{\mathrm{AARO}})=min(P_{\mathrm{AARO\_SDP}})$. Moreover, $({\eta }_p^t,a_p^t,z^t)\in \mathbb{R}\times {\mathbb{R}}^q\times \mathbb{R}$, $t=0,\dots ,T$ and $p=1,\dots ,n$, is a solution for $(P_{\mathrm{AARO}})$ if and only if there exist $s\in \mathbb{R}$, ${\lambda }_1\ge 0$, ${\lambda }_2^{\alpha ,\beta }\ge 0$ with $(\alpha ,\beta )\in \mathcal{T}$, ${\lambda }_3^{\alpha ,\beta }\ge 0$ with $(\alpha ,\beta )\in \mathcal{C}$, ${\lambda }_4^{\alpha ,\beta }\ge 0$ with $(\alpha ,\beta )\in \mathcal{N}$ and ${\lambda }_{5,t}^{\alpha ,\beta }\ge 0$ with $(\alpha ,\beta )\in \mathcal{H}$ and $t=0,\dots ,T$, such that $(s,{\eta }_p^t,a_p^t,z^t,{\lambda }_1,{\lambda }_2^{\alpha ,\beta },{\lambda }_3^{\alpha ,\beta },{\lambda }_4^{\alpha ,\beta },{\lambda }_{5,t}^{\alpha ,\beta })$ is a solution for $(P_{\mathrm{AARO\_SDP}})$.

Proof.

Firstly, we note that the problem $(P_{\mathrm{AARO}})$ is equivalent to the following adjustable robust problem where the objective function is uncertainty free: $\begin{array}{ll}{\displaystyle \underset{\begin{array}{c}{\eta }_p^t\in \mathbb{R},a_p^t\in {\mathbb{R}}^q\\ z^t\ge 0,s\in \mathbb{R}\end{array}}{min}}& s\\ \text{subject to}& {\displaystyle \sum _{t=0}^T\sum _{(\alpha ,\beta )\in \mathcal{S}}\sum _{p=1}^n{w}_p^t(u)[{\overline{d}}_{\mathrm{\alpha \beta }}^{p,t}+{\overline{v}}_{\mathrm{\alpha \beta }}^{T}u]\le s,\mathrm{\forall }\Vert u\Vert \le r,}\\ & {\displaystyle \sum _{t=0}^T\sum _{p=1}^n{w}_p^t(u)[{\overline{d}}_{\mathrm{\alpha \beta }}^{p,t}+{\overline{v}}_{\mathrm{\alpha \beta }}^{T}u]\ge {\gamma }_{\mathrm{lb}},\mathrm{\forall }\Vert u\Vert \le r,(\alpha ,\beta )\in \mathcal{T},}\end{array}$ $\begin{array}{ll}& {\displaystyle \sum _{t=0}^T\sum _{p=1}^n{w}_p^t(u)[{\overline{d}}_{\mathrm{\alpha \beta }}^{p,t}+{\overline{v}}_{\mathrm{\alpha \beta }}^{T}u]\le {\gamma }_{\mathrm{ub}}^c,\mathrm{\forall }\Vert u\Vert \le r,(\alpha ,\beta )\in \mathcal{C},}\\ & {\displaystyle \sum _{t=0}^T\sum _{p=1}^n{w}_p^t(u)[{\overline{d}}_{\mathrm{\alpha \beta }}^{p,t}+{\overline{v}}_{\mathrm{\alpha \beta }}^{T}u]\le {\gamma }_{\mathrm{ub}}^n,\mathrm{\forall }\Vert u\Vert \le r,(\alpha ,\beta )\in \mathcal{N},}\\ & {\displaystyle \sum _{p=1}^n{w}_p^t(u)[{\overline{d}}_{\mathrm{\alpha \beta }}^{p,t}+{\overline{v}}_{\mathrm{\alpha \beta }}^{T}u]\ge {\rho }_{\mathrm{\alpha \beta }}^t{z}^t,}\\ & \mathrm{\forall }{\rho }_{\mathrm{\alpha \beta }}^t\in {\mathcal{U}}_{\mathrm{\alpha \beta }}^t,\mathrm{\forall }\Vert u\Vert \le r,(\alpha ,\beta )\in \mathcal{H},\\ & {\displaystyle \sum _{t=0}^T{z}^t\ge {\gamma }_{\mathrm{lb}},z^t\ge 0,t=0,\dots ,T,}\\ & w_p^t(u)={\eta }_p^t+(a_p^t{)}^{T}u,t=0,1,\dots ,T,\\ & w_p^t(u)\ge 0,\mathrm{\forall }\Vert u\Vert \le r.\end{array}$For each fixed $(\alpha ,\beta )\in \mathcal{S}$, observe that $\begin{array}{rl}& \sum _{t=0}^T\sum _{p=1}^n{w}_p^t(u)[{\overline{d}}_{\mathrm{\alpha \beta }}^{p,t}+{\overline{v}}_{\mathrm{\alpha \beta }}^{T}u]\\ & =\sum _{t=0}^T\sum _{p=1}^n({\eta }_p^t+(a_p^t{)}^{T}u)({\overline{d}}_{\mathrm{\alpha \beta }}^{p,t}+{\overline{v}}_{\mathrm{\alpha \beta }}^{T}u)\\ & =\left(\sum _{t=0}^T\sum _{p=1}^n{\eta }_p^t{\overline{d}}_{\mathrm{\alpha \beta }}^{p,t}\right)+{[\sum _{t=0}^T\sum _{p=1}^n({\overline{d}}_{\mathrm{\alpha \beta }}^{p,t}{a}_p^t+{\eta }_p^t{\overline{v}}_{\mathrm{\alpha \beta }})]}^{T}u\\ & +u^T\left(\sum _{t=0}^T\sum _{p=1}^n\frac{a_p^t{\overline{v}}_{\mathrm{\alpha \beta }}^T+{\overline{v}}_{\mathrm{\alpha \beta }}(a_p^t{)}^T}{2}\right)u\\ & ={\theta }_{\mathrm{\alpha \beta }}+w_{\mathrm{\alpha \beta }}^{T}u+u^T{W}_{\mathrm{\alpha \beta }}u,\end{array}$where ${\theta }_{\mathrm{\alpha \beta }}=\sum _{t=0}^T\sum _{p=1}^n{\eta }_p^t{\overline{d}}_{\mathrm{\alpha \beta }}^{p,t},w_{\mathrm{\alpha \beta }}=\sum _{t=0}^T\sum _{p=1}^n({\overline{d}}_{\mathrm{\alpha \beta }}^{p,t}{a}_p^t+{\eta }_p^t{\overline{v}}_{\mathrm{\alpha \beta }})$and $W_{\mathrm{\alpha \beta }}=\sum _{t=0}^T\sum _{p=1}^n\frac{a_p^t{\overline{v}}_{\mathrm{\alpha \beta }}^T+{\overline{v}}_{\mathrm{\alpha \beta }}(a_p^t{)}^T}{2}.$So, $\begin{array}{rl}& \sum _{t=0}^T\sum _{(\alpha ,\beta )\in \mathcal{S}}\sum _{p=1}^n{w}_p^t(u)[{\overline{d}}_{\mathrm{\alpha \beta }}^{p,t}+{\overline{v}}_{\mathrm{\alpha \beta }}^{T}u]\\ & =\sum _{(\alpha ,\beta )\in \mathcal{S}}\sum _{t=0}^T\sum _{p=1}^n{w}_p^t(u)[{\overline{d}}_{\mathrm{\alpha \beta }}^{p,t}+{\overline{v}}_{\mathrm{\alpha \beta }}^{T}u]\\ & =\left(\sum _{(\alpha ,\beta )\in \mathcal{S}}{\theta }_{\mathrm{\alpha \beta }}\right)+{\left(\sum _{(\alpha ,\beta )\in \mathcal{S}}{w}_{\mathrm{\alpha \beta }}\right)}^{T}u+u^T\left(\sum _{(\alpha ,\beta )\in \mathcal{S}}{W}_{\mathrm{\alpha \beta }}\right)u.\end{array}$Now, the constraint $\sum _{t=0}^T\sum _{(\alpha ,\beta )\in \mathcal{S}}\sum _{p=1}^n{w}_p^t(u)[{\overline{d}}_{\mathrm{\alpha \beta }}^{p,t}+{\overline{v}}_{\mathrm{\alpha \beta }}^{T}u]\le s,\mathrm{\forall }\Vert u\Vert \le r,$can be equivalently rewritten as $\begin{array}{rl}& \Vert u{\Vert }^2\le r^2\Rightarrow (s-\sum _{(\alpha ,\beta )\in \mathcal{S}}{\theta }_{\mathrm{\alpha \beta }})-{\left(\sum _{(\alpha ,\beta )\in \mathcal{S}}{w}_{\mathrm{\alpha \beta }}\right)}^{T}u-u^T\left(\sum _{(\alpha ,\beta )\in \mathcal{S}}{W}_{\mathrm{\alpha \beta }}\right)\\ & u\ge 0.\end{array}$By Lemma 4.1, this implication is further equivalent to the condition that there exists ${\lambda }_1\ge 0$ such that $\begin{array}{rl}& (s-\sum _{(\alpha ,\beta )\in \mathcal{S}}{\theta }_{\mathrm{\alpha \beta }})-{\left(\sum _{(\alpha ,\beta )\in \mathcal{S}}{w}_{\mathrm{\alpha \beta }}\right)}^{T}u-u^T\left(\sum _{(\alpha ,\beta )\in \mathcal{S}}{W}_{\mathrm{\alpha \beta }}\right)u\\ & +{\lambda }_1(\Vert u{\Vert }^2-r^2)\ge 0\text{for all }u\in {\mathbb{R}}^q.\end{array}$Using Lemma 4.2, this inequality can be equivalently rewritten as $\left(\begin{array}{cc}{\displaystyle -\sum _{(\alpha ,\beta )\in \mathcal{S}}{W}_{\mathrm{\alpha \beta }}+{\lambda }_1{I}_q}& {\displaystyle -\sum _{(\alpha ,\beta )\in \mathcal{S}}{w}_{\mathrm{\alpha \beta }}/2}\\ {\displaystyle -\sum _{(\alpha ,\beta )\in \mathcal{S}}{w}_{\mathrm{\alpha \beta }}^T/2}& {\displaystyle -\sum _{(\alpha ,\beta )\in \mathcal{S}}{\theta }_{\mathrm{\alpha \beta }}+s-{\lambda }_1{r}^2}\end{array}\right)⪰0.$Similarly, we have the following equivalences:

• For each $(\alpha ,\beta )\in \mathcal{T}$, the constraint $\sum _{t=0}^T\sum _{p=1}^n{w}_p^t(u)[{\overline{d}}_{\mathrm{\alpha \beta }}^{p,t}+{\overline{v}}_{\mathrm{\alpha \beta }}^{T}u]\ge {\gamma }_{\mathrm{lb}},\mathrm{\forall }\Vert u\Vert \le r,$is equivalent to the matrix inequality condition that, there exists ${\lambda }_2^{\alpha ,\beta }\ge 0$ such that $\left(\begin{array}{cc}{W}_{\mathrm{\alpha \beta }}+{\lambda }_2^{\alpha ,\beta }{I}_q& w_{\mathrm{\alpha \beta }}/2\\ w_{\mathrm{\alpha \beta }}^T/2& {\theta }_{\mathrm{\alpha \beta }}-{\gamma }_{\mathrm{lb}}-{\lambda }_2^{\alpha ,\beta }{r}^2\end{array}\right)⪰0.$

• For each $(\alpha ,\beta )\in \mathcal{C}$, the constraint $\sum _{t=0}^T\sum _{p=1}^n{w}_p^t(u)[{\overline{d}}_{\mathrm{\alpha \beta }}^{p,t}+{\overline{v}}_{\mathrm{\alpha \beta }}^{T}u]\le {\gamma }_{\mathrm{ub}}^c,\mathrm{\forall }\Vert u\Vert \le r,$is equivalent to the condition that there exists ${\lambda }_3^{\alpha ,\beta }\ge 0$ such that $\left(\begin{array}{cc}-W_{\mathrm{\alpha \beta }}+{\lambda }_3^{\alpha ,\beta }{I}_q& -w_{\mathrm{\alpha \beta }}/2\\ -w_{\mathrm{\alpha \beta }}^T/2& -{\theta }_{\mathrm{\alpha \beta }}+{\gamma }_{\mathrm{ub}}^c-{\lambda }_3^{\alpha ,\beta }{r}^2\end{array}\right)⪰0.$

• For each $(\alpha ,\beta )\in \mathcal{N}$, the constraint $\sum _{t=0}^T\sum _{p=1}^n{w}_p^t(u)[{\overline{d}}_{\mathrm{\alpha \beta }}^{p,t}+{\overline{v}}_{\mathrm{\alpha \beta }}^{T}u]\le {\gamma }_{\mathrm{ub}}^n,\mathrm{\forall }\Vert u\Vert \le r,$is equivalent to the condition that there exists ${\lambda }_4^{\alpha ,\beta }\ge 0$ such that $\left(\begin{array}{cc}-W_{\mathrm{\alpha \beta }}+{\lambda }_4^{\alpha ,\beta }{I}_q& -w_{\mathrm{\alpha \beta }}/2\\ -w_{\mathrm{\alpha \beta }}^T/2& -{\theta }_{\mathrm{\alpha \beta }}+{\gamma }_{\mathrm{ub}}^n-{\lambda }_4^{\alpha ,\beta }{r}^2\end{array}\right)⪰0.$

• Since $z^t\ge 0$, for each $(\alpha ,\beta )\in \mathcal{H}$, the constraint $\sum _{p=1}^n{w}_p^t(u)[{\overline{d}}_{\mathrm{\alpha \beta }}^{p,t}+{\overline{v}}_{\mathrm{\alpha \beta }}^{T}u]\ge {\rho }_{\mathrm{\alpha \beta }}^t{z}^t,\mathrm{\forall }{\rho }_{\mathrm{\alpha \beta }}^t\in {\mathcal{U}}_{\mathrm{\alpha \beta }}^t,\mathrm{\forall }\Vert u\Vert \le r,$is equivalent to the condition that, for $t=0,1,\dots ,t_k-1$, $\sum _{p=1}^n{w}_p^t(u)[{\overline{d}}_{\mathrm{\alpha \beta }}^{p,t}+{\overline{v}}_{\mathrm{\alpha \beta }}^{T}u]\ge [({\rho }_0{)}_{\mathrm{\alpha \beta }}+\mathrm{\eta t}+\mathrm{\Gamma t}]z^t,\mathrm{\forall }\Vert u\Vert \le r,$and, for $t=t_k,\dots ,T$, $\begin{array}{rl}& \sum _{p=1}^n{w}_p^t(u)[{\overline{d}}_{\mathrm{\alpha \beta }}^{p,t}+{\overline{v}}_{\mathrm{\alpha \beta }}^{T}u]\ge [{\rho }_{\mathrm{\alpha \beta }}^{\mathrm{observed}}+\eta (t-t_k)+\mathrm{\Gamma }(t-t_k)+\nu ]z^t,\\ & \mathrm{\forall }\Vert u\Vert \le r.\end{array}$Similarly, as before, the latter constraints can be further equivalently reformulated as the matrix inequality constraints that there exist ${\lambda }_{5,t}^{\alpha ,\beta }\ge 0$ such that, for $t=0,1,\dots ,t_k-1$, $\left(\begin{array}{cc}{W}_{\mathrm{\alpha \beta }}^t+{\lambda }_{5,t}^{\alpha ,\beta }{I}_q& w_{\mathrm{\alpha \beta }}^t/2\\ (w_{\mathrm{\alpha \beta }}^t{)}^T/2& {\theta }_{\mathrm{\alpha \beta }}^t-[({\rho }_0{)}_{\mathrm{\alpha \beta }}+\mathrm{\eta t}+\mathrm{\Gamma t}]z^t-{\lambda }_{5,t}^{\alpha ,\beta }{r}^2\end{array}\right)⪰0,$

  and, for $t=t_k,\dots ,T$,

  $\begin{array}{rl}& \left(\begin{array}{cc}{W}_{\mathrm{\alpha \beta }}^t+{\lambda }_{5,t}^{\alpha ,\beta }{I}_q& w_{\mathrm{\alpha \beta }}^t/2\\ (w_{\mathrm{\alpha \beta }}^t{)}^T/2& {\theta }_{\mathrm{\alpha \beta }}^t-[{\rho }_{\mathrm{\alpha \beta }}^{\mathrm{observed}}+\eta (t-t_k)+\mathrm{\Gamma }(t-t_k)+\nu ]z^t-{\lambda }_{5,t}^{\alpha ,\beta }{r}^2\end{array}\right)⪰0.\end{array}$

Finally, note that the constraint $w_p^t(u)\ge 0,\mathrm{\forall }\Vert u\Vert \le r$ with $w_p^t(u)={\eta }_p^t+(a_p^t{)}^{T}u$ is equivalent to ${\eta }_p^t-r\Vert a_p^t\Vert \ge 0$. Therefore, the conclusion follows by combining the above equivalent reformulations of the constraints.

Remark 4.1

Solving large-scale conic reformulation problems

Note that our exact conic formulation $(P_{\mathrm{AARO\_SDP}})$ involves the variables ${\lambda }_2^{\alpha ,\beta },\dots ,{\lambda }_4^{\alpha ,\beta },{\lambda }_{5,t}^{\alpha ,\beta }$. In the case where the number of voxel is huge, that is, the range of $(\alpha ,\beta )$ is large, the exact conic formulation may lead to a large size of a semi-definite program. The resulting large-size semi-definite program (SDP) can be expansive or impractical to solve due to the limitation of the state-of-art of SDP software. This motivates us to examine, in our implementation, a modified conic formulation $(P_{\mathrm{SDP}})$ where we require the variables ${\lambda }_2,{\lambda }_3,{\lambda }_4,{\lambda }_{5,t}$ are independent of the choice of voxel $(\alpha ,\beta )$. See Appendix 2 for the explicit formulation of the modified conic formulation $(P_{\mathrm{SDP}})$.

From the constructions of $(P_{\mathrm{AARO\_SDP}})$ and $(P_{\mathrm{SDP}})$, one sees that the feasible region of $(P_{\mathrm{AARO\_SDP}})$ is larger than $(P_{\mathrm{SDP}})$ and both problem shares the same objective function. So, $min(P_{\mathrm{AARO\_SDP}})\le min(P_{\mathrm{SDP}})$. This together with the preceding theorem shows that the objective value of the modified conic formulation $(P_{\mathrm{SDP}})$ serves as an upper bound for the objective value of the adjustable robust optimization problem $(P_{\mathrm{AARO}})$, that is, $min(P_{\mathrm{AARO}})\le min(P_{\mathrm{SDP}})$. In other words, the modified conic formulation $(P_{\mathrm{SDP}})$ provides an upper bound (upper estimate) to the true minimum total dose emitted to the whole tissues in the AARO model $(P_{\mathrm{AARO}})$.

We also note that the size of the modified problem is, in general, much smaller than the exact conic formulation $(P_{\mathrm{SDP}})$. In our implementation, the modified conic problem $(P_{\mathrm{SDP}})$ can be solved efficiently in a reasonable time with improved performance comparing with the static model.

5. Evaluating AARO approach on liver cancer case

In this section, we evaluate three models: the hypoxia linear program model with nominal data, $(P_{\mathrm{hypoxia}\mathrm{LP}})$, and two models for handling both dose matrix uncertainty and re-oxygenation uncertainty, namely the static robust model $(P_{\mathrm{static}})$, and the modified conic formulation of the adjustable robust model $(P_{\mathrm{SDP}})$, see Appendix 2. The linear program model and the static robust model have been considered in the literature [4,32], and so, they serve as a benchmark for our proposed adjustable robust model $(P_{\mathrm{SDP}})$. We apply each of these to a simulated liver cancer case using freely available data to simulate both hypoxia and dose influence matrix uncertainty over the treatment period. We evaluate our models' performances using the common metrics of dose criteria [5], volume criteria, equivalent uniform dose (EUD) and dose-volume histograms (DVHs) across the total volume, the planning target volume (PTV), and the organs-at-risk (OARs).

Patient and Source Data. For our study, we use a dataset provided by [6] for a liver cancer case. We set a lower threshold on the tumour of 55 Gy based on a prescribed dose of 60 Gy to the tumour. Note that Gray (Gy) is the unit used to measure the amount of radiation. We set the sole OAR to be the heart, with a limiting dose of 70 Gy in line with the recommendations in [34]. The number of voxels in the tumour and heart are 6954 and 33021, respectively, with a total voxel count of 7910952. All three models follow a treatment plan, where T = 3 and the observation of the hypoxia value is taken at $t_k=2$. To reduce computational complexity, following [35], we eliminate the set of voxels in the liver cancer data set which, by prior knowledge, will never receive a dose greater than zero. Rather than enforcing a fractionated approach with dose limits ${\gamma }_{\mathrm{lb}}/(T+1)$, we utilize the parameter $z^t$ with $t=0,\dots ,T$, to allow the optimization program to instead select the optimal dosage distribution over the treatment period.

For our conic programming models (SDP and SOCP), we used the commercial solver MOSEK [31] (Mosek ApS, Copenhagen) with the optimization toolbox YALMIP [36], implemented in Matlab 2020 (MathWorks, Inc., Natick, Massachusetts). For linear programming models, we used the inbuilt solver linprog from Matlab 2020. All calculations were done on a machine with an AMD Ryzen 7 3700X (4.00 Ghz) processor and 32 GB RAM.

The choice of couch-gantry angle combinations were determined via preprocessing to reduce model complexity. The hypoxia linear program model is first solved using $\mathtt{\text{linprog}}$ in Matlab and the total dosage for the beamlets of each couch-gantry angle is summed up. The five largest couch-gantry angle contributors were selected and the linear program re-solved to ensure feasibility. This resulted in 5 beams from gantry angles ${16}^{\circ }$, ${30}^{\circ }$, ${270}^{\circ }$, ${279}^{\circ }$, and ${306}^{\circ }$ , and couch angles $-{90}^{\circ }$, $-{59}^{\circ }$, $0^{\circ }$, ${15}^{\circ }$, and ${32}^{\circ },$ respectively.

Model setups and uncertainty simulation. Following the setup of existing literature [12,37], we simulated 100 instances of realized uncertainty in the dose influence matrix and hypoxia at the starting time 0 and the time $t_k=2$ in order to compare the affinely adjustable robust model, the static robust model and the hypoxia linear program model.

To simulate re-oxygenation, we followed the approach of [4] (outlined in Section 3) to construct uncertainty sets with the following parameters: re-oxygenation rate $\eta =2.8\mathrm{\%}$, uncertainty budget rate $\mathrm{\Gamma }=1.6\mathrm{\%}$, initial observed hypoxia parameter $({\rho }_0{)}_{\mathrm{\alpha \beta }}=1.2$ for all $(\alpha ,\beta )$, and inexactness parameter $\nu =0.05$. We used the above constructed uncertainty sets to form the static robust model and the affinely adjustable robust model.

Our realized hypoxia is an ‘ideal value’ of the hypoxia parameter assuming we have perfect knowledge, representing the best-case scenario. To do this, we first generated the realized hypoxia for the hypoxia linear program model by simulating an observed hypoxia ${\rho }_{\mathrm{\alpha \beta }}^{\mathrm{observed}}$ at $t=t_k=2$ with a uniformly random selection from the outer cone (see Figure 1). Then, we obtained the realized hypoxia parameter by uniformly randomly selecting a sample value from the inexact uncertainty set which is an interval centred at ${\rho }_{\mathrm{\alpha \beta }}^{\mathrm{observed}}$ with radius being the inexactness parameter ν. The hypoxic region for all three models is set to 54.4% of the tumour region using uniform random selection, in accordance to previous studies [38].

To simulate dose influence matrix uncertainty, we set the dimension q = 1, set r as a constant equal to 1% of the median dosage value of our dose matrix, and set ${\overline{v}}_{\mathrm{\alpha \beta }}$ to a constant value of 1 for all $(\alpha ,\beta )$. We set ${\overline{d}}_{\mathrm{\alpha \beta }}^{p,t}$ to the values taken from the dose matrix values provided in the dataset for all three models. The values for the uncertainty parameter u were simulated by uniformly sampling from the 1-dimensional ball $\mathcal{B}=\{u:\Vert u\Vert \le r\}$ which is just the interval $[-r,r]$.

We emphasize that the purpose of our experiments is to evaluate the potential for adjustable robust methods in radiation treatment planning, and so our choice of r is merely a simple heuristic and is likely far from optimal. Estimation of the size of the uncertainty set is a common problem in robust optimization in general and may often rely on expert opinion.²However, some other simple heuristics may be employed. The choice of median value of the dose influence matrix is a result of the distribution and scale of values in the dose influence matrix, with a large number of small values and a small number of values multiple magnitudes larger.

Results. The results reflect cumulative dose over the treatment period. The total dose, dose to the PTV region and dose to the heart across the simulations are shown in Figures 2–4. The mean dose values across the simulations are summarized in Table 1 which shows that the adjustable robust model outperforms the static robust model in all the measures (total dosage, PTV dosage and heart dosage) with a 3.8% reduction in dosage to the heart and a 0.03% reduction in total dosage. Also, the adjustable robust model only has an 8% increase in total dosage compared to the ideal hypoxia model. This is also illustrated in Figures 2–4.

Figure 2. Overall total dose.

Figure 3. PTV dose.

Figure 4. Heart dose.

Table 1. Dose summary.

Dose (Gy)	Hypoxia	Static robust	Adjustable robust
Mean total dose	3.183e+06	3.438e+06	3.437e+06
Mean PTV dose	4.711e+05	5.069e+05	5.062e+05
Mean heart dose	8.827e+04	1.014e+05	9.758e+04

Our results demonstrate that the adjustable model successfully handles uncertainty in both re-oxygenation and the dose matrix. It also shows that, for a comparable total dose, it provides lower doses to critical organs, outperforming a static method by utilizing information halfway through the treatment plan (i.e. the hypoxia observation taken at time $t_k$).

This success comes with a small cost with reference to the ideal (nominal) solution obtained by solving the linear program $(P_{\mathrm{hypoxia}\mathrm{LP}})$ with exact knowledge on the data. We emphasise that in practice, the hypoxia solution cannot be achieved, as it requires perfect knowledge on the (uncertain) dose matrix before any treatment has begun.

Our two competing objectives are to achieve the minimum dosage to the tumour region while minimizing dosage to OAR. To measure the clinical quality of the results, we use typical criteria such as dose-volume histograms (DVHs), dose-criteria $D_x$ and volume-criteria $V_x$ specifications. Here, $D_{x\mathrm{\%}}$ is the minimum dose received by $\ge x\mathrm{\%}$ of the structure and $V_x$ is the minimum volume percentage receiving at least x Gy. For the tumour, the median dose $D_{50\mathrm{\%}}$ is recommended to be the prescribed dose and $D_{95\mathrm{\%}}$ is suggested at 95% of the prescribed dose [2]. Table 2 demonstrates the results of our three methods over the PTV region in terms of these measures. Table 3 presents volume-criteria metrics for dosage values $5,10,20,30,40,$ and 50 over the sole critical organ of our investigation, the heart. We are also interested in the equivalent uniform dose (EUD) to the tumour, a measure of dose homogeneity. Given our tumour region $\mathcal{T}$, $\mathrm{EUD}={(\frac{1}{|\mathcal{T}|}\sum _{(\alpha ,\beta )\in \mathcal{T}}\sum _{p,t}({\overline{d}}_{\mathrm{\alpha \beta }}^{p,t}{w}_p^t{)}^{\zeta })}^{\frac{1}{\zeta }},$where ${\overline{d}}_{\mathrm{\alpha \beta }}^{p,t}$ is the nominal dose matrix entries obtained from the dataset and $w_p^t$ are the beamlet density with $p=1,\dots ,n$ and $t=0,\dots ,T$. We have used $\zeta =-10$ for the EUD of the PTV region, as is common practice [5].

Table 2. PTV plan evaluation.

Dose criteria (Gy)	Hypoxia	Static robust	Adjustable robust
$D_{50\mathrm{\%}}\ge 60$	66.44	71.87	71.71
$D_{95\mathrm{\%}}\ge 57$	61.01	65.82	65.74
$\mathrm{EUD}=60$	65.73	70.94	70.85

Table 3. Heart plan evaluation.

Heart dose criteria (Gy)	Hypoxia	Static robust	Adjustable robust
$V_5$	13.03%	14.96%	13.90%
$V_{10}$	7.03%	7.99%	7.72%
$V_{20}$	2.96%	3.77%	3.62%
$V_{30}$	1.60%	2.25%	2.19%
$V_{40}$	0.95%	1.33%	1.34%
$V_{50}$	0.48%	0.75%	0.76%

It is clear from Tables 2 and 3 as well as Figures 3 and 4 that the adjustable model performs better than the static model, especially in the case of dosage to the heart. Note that our models do not directly optimize the dosage to critical organs (our objective function uses the total dosage as opposed to critical organ dosage) and so there is no guarantee that an adjustable model will prescribe less dosage to the critical organs than a static model (see Table 3), static model performs slightly better over larger volumes of the heart. However, our results have shown that the adjustable model still successfully prescribes dosage over smaller volumes of the critical organ, with $V_5$ represented over 7% less volume for the adjustable model versus the static model. In addition, we have also included dose-volume histograms for these results in Appendix 1.

6. Conclusion and future work

We presented an AARO model that incorporates uncertainties both in dose influence matrix and re-oxygenation data as well as inexactness of the revealed (re-oxygenation) data. We allowed, more generally, the beamlet densities to be affinely adjustable with respect to the uncertainty parameter. We evaluated our model on a simulated liver cancer case as a proof-of-concept.

Our experiments have demonstrated that the adjustable model successfully handles uncertainty in both re-oxygenation and the dose matrix. By utilizing information halfway through the treatment plan, we have shown that the adjustable solutions of the AARO model outperformed a static method on the two counts – PTV dosage and OAR dosage, while maintaining a similar total dose. They are within comparable performance of the hypoxia linear program model with nominal data.

It would be of great interest to extend our study to models that incorporate intermediate information in treatment planning with accrued additional costs, such as the intermediate imaging costs. This will be an interesting topic for further research.

Acknowledgments

The authors are grateful to the referees for their comments and suggestions on the paper.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Data availability statement

We used a dataset provided by [6] for our study on a liver cancer case. The data generated during the current study are available from the fourth author on reasonable request.

Additional information

Funding

Research was partially supported by a grant from the Australian Research Council.

Notes

1 This paper considered a multi-stage robust optimization model, incorporating intermediate information in treatment planning with accrued additional cost. This multi-stage model results in an equivalent linear optimization problem.

2 A better estimate of the size of the uncertainty set r for our problem could be found by applying max-pooling over the 3-dimensional neighbourhood surrounding each voxel in the treatment volume, which can be interpreted as describing the maximal possible change in dosage throughout the treatment due to patient movement (under certain assumptions). We leave the evaluation of this estimate and other such choices for future work.

Appendices

Appendix 1. Dose-volume histograms

We also provide dose-volume histograms (DVHs) for the interest of the medical community. The DVH bands in the shaded region represent the maximum and minimum dose of the model solutions across the simulations, while the centre line represents the median dose of the model solutions. For the hypoxia model, the minimum and maximum dose has a larger variation and hence more visible than in the adjustable model, especially in the PTV DVH (Figure A1). In the heart DVH (Figure A2), the adjustable model and the hypoxia model are consistent.

Appendix 2. Simplified conic model

In this appendix, we give the explicit form of the modified conic formulation for the adjustable robust optimization model which we denoted as $(P_{\mathrm{SDP}})$. We note that, by comparing with the exact conic reformulation $(P_{\mathrm{AARO\_SDP}})$, the only difference is that the variables ${\lambda }_2,{\lambda }_3,{\lambda }_4,{\lambda }_{5,t}$ in $(P_{\mathrm{SDP}})$ are independent of the choice of voxel $(\alpha ,\beta )$. $\begin{array}{r}\begin{array}{ll}& (P_{\mathrm{SDP}})\\ & {\displaystyle \underset{\begin{array}{c}s\in \mathbb{R},{\eta }_p^t\in \mathbb{R},a_p^t\in {\mathbb{R}}^q,z^t\in \mathbb{R}\\ {\lambda }_1,{\lambda }_2,{\lambda }_3,{\lambda }_4,{\lambda }_{5,t}\ge 0\end{array}}{min}}& s\\ & \text{subject to }& {\eta }_p^t-r\Vert a_p^t\Vert \ge 0,t=0,1,\dots ,T\\ & & \left(\begin{array}{cc}{\displaystyle -\sum _{(\alpha ,\beta )\in \mathcal{S}}{W}_{\mathrm{\alpha \beta }}+{\lambda }_1{I}_q}& {\displaystyle -\sum _{(\alpha ,\beta )\in \mathcal{S}}{w}_{\mathrm{\alpha \beta }}/2}\\ {\displaystyle -\sum _{(\alpha ,\beta )\in \mathcal{S}}{w}_{\mathrm{\alpha \beta }}^T/2}& {\displaystyle -\sum _{(\alpha ,\beta )\in \mathcal{S}}{\theta }_{\mathrm{\alpha \beta }}+s-{\lambda }_1{r}^2}\end{array}\right)⪰0,\\ & & \left(\begin{array}{cc}{W}_{\mathrm{\alpha \beta }}+{\lambda }_2{I}_q& w_{\mathrm{\alpha \beta }}/2\\ w_{\mathrm{\alpha \beta }}^T/2& {\theta }_{\mathrm{\alpha \beta }}-{\gamma }_{\mathrm{lb}}-{\lambda }_2{r}^2\end{array}\right)⪰0,\text{for all }(\alpha ,\beta )\in \mathcal{T},\\ & & \left(\begin{array}{cc}-W_{\mathrm{\alpha \beta }}+{\lambda }_3{I}_q& -w_{\mathrm{\alpha \beta }}/2\\ -w_{\mathrm{\alpha \beta }}^T/2& -{\theta }_{\mathrm{\alpha \beta }}+{\gamma }_{\mathrm{ub}}^c-{\lambda }_3{r}^2\end{array}\right)⪰0,\text{for all }(\alpha ,\beta )\in \mathcal{C},\\ & & \left(\begin{array}{cc}-W_{\mathrm{\alpha \beta }}+{\lambda }_4{I}_q& -w_{\mathrm{\alpha \beta }}/2\\ -w_{\mathrm{\alpha \beta }}^T/2& -{\theta }_{\mathrm{\alpha \beta }}+{\gamma }_{\mathrm{ub}}^n-{\lambda }_4{r}^2\end{array}\right)⪰0,\text{for all }(\alpha ,\beta )\in \mathcal{N},\\ & & \left(\begin{array}{cc}{W}_{\mathrm{\alpha \beta }}^t+{\lambda }_{5,t}{I}_q& w_{\mathrm{\alpha \beta }}^t/2\\ (w_{\mathrm{\alpha \beta }}^t{)}^T/2& {\theta }_{\mathrm{\alpha \beta }}^t-[({\rho }_0{)}_{\mathrm{\alpha \beta }}+\mathrm{\eta t}+\mathrm{\Gamma t}]z^t-{\lambda }_{5,t}{r}^2\end{array}\right)⪰0,\\ & & \text{for all }(\alpha ,\beta )\in \mathcal{H},t=0,1,\dots ,t_k-1\\ & & \left(\begin{array}{cc}{W}_{\mathrm{\alpha \beta }}^t+{\lambda }_{5,t}{I}_q& w_{\mathrm{\alpha \beta }}^t/2\\ (w_{\mathrm{\alpha \beta }}^t{)}^T/2& {\theta }_{\mathrm{\alpha \beta }}^t-[({\rho }_{t_k}^{\mathrm{observe}}{)}_{\mathrm{\alpha \beta }}+\eta (t-t_k)\\ & +\mathrm{\Gamma }(t-t_k)+\nu ]z^t-{\lambda }_{5,t}{r}^2\end{array}\right)⪰0,\\ & & \text{for all }(\alpha ,\beta )\in \mathcal{H},t=t_k,\dots ,T\\ & & {\displaystyle {\theta }_{\mathrm{\alpha \beta }}=\sum _{t=0}^T\sum _{p=1}^n{\eta }_p^t{\overline{d}}_{\mathrm{\alpha \beta }}^{p,t},w_{\mathrm{\alpha \beta }}=\sum _{t=0}^T\sum _{p=1}^n({\overline{d}}_{\mathrm{\alpha \beta }}^{p,t}{a}_p^t+{\eta }_p^t{\overline{v}}_{\mathrm{\alpha \beta }}),}\\ & & {\displaystyle W_{\mathrm{\alpha \beta }}=\sum _{t=0}^T\sum _{p=1}^n\frac{a_p^t{\overline{v}}_{\mathrm{\alpha \beta }}^T+{\overline{v}}_{\mathrm{\alpha \beta }}(a_p^t{)}^T}{2},}\\ & & {\displaystyle {\theta }_{\mathrm{\alpha \beta }}^t=\sum _{p=1}^n{\eta }_p^t{\overline{d}}_{\mathrm{\alpha \beta }}^{p,t},w_{\mathrm{\alpha \beta }}^t=\sum _{p=1}^n({\overline{d}}_{\mathrm{\alpha \beta }}^{p,t}{a}_p^t+{\eta }_p^t{\overline{v}}_{\mathrm{\alpha \beta }}),}\end{array}\end{array}$ $\begin{array}{r}\begin{array}{lll}& & W_{\mathrm{\alpha \beta }}^t=\sum _{p=1}^n\frac{a_p^t{\overline{v}}_{\mathrm{\alpha \beta }}^T+{\overline{v}}_{\mathrm{\alpha \beta }}(a_p^t{)}^T}{2},\\ & & {\displaystyle \sum _{t=0}^T{z}^t\ge {\gamma }_{\mathrm{lb}},z^t\ge 0.}\end{array}\end{array}$

Appendix 3. SOCP reformulation for the static robust model

In this appendix, we give the explicit form of the second-order cone program reformulation for the static robust model problem, which we denote it as $(P_{\mathrm{Static\_SOCP}})$. $\begin{array}{r}\begin{array}{lll}(P_{\mathrm{Static\_SOCP}})& \underset{s\in \mathbb{R},{\eta }_p^t\in \mathbb{R},z^t\in \mathbb{R}}{min}& s\\ & \text{subject to }& {\displaystyle \sum _{(\alpha ,\beta )\in \mathcal{S}}{\overline{\alpha }}_{\mathrm{\alpha \beta }}+r\Vert \sum _{(\alpha ,\beta )\in \mathcal{S}}{\overline{w}}_{\mathrm{\alpha \beta }}\Vert \le s,}\\ & & {\overline{\alpha }}_{\mathrm{\alpha \beta }}-r\Vert {\overline{w}}_{\mathrm{\alpha \beta }}\Vert \ge {\gamma }_{\mathrm{lb}},(\alpha ,\beta )\in \mathcal{T},\\ & & {\overline{\alpha }}_{\mathrm{\alpha \beta }}+r\Vert {\overline{w}}_{\mathrm{\alpha \beta }}\Vert \le {\gamma }_{\mathrm{ub}}^c,(\alpha ,\beta )\in \mathcal{C},\\ & & {\overline{\alpha }}_{\mathrm{\alpha \beta }}+r\Vert {\overline{w}}_{\mathrm{\alpha \beta }}\Vert \le {\gamma }_{\mathrm{ub}}^n,(\alpha ,\beta )\in \mathcal{N},\\ & & [({\rho }_0{)}_{\mathrm{\alpha \beta }}+\mathrm{\eta t}+\mathrm{\Gamma t}]z^t\le {\overline{\alpha }}_{\mathrm{\alpha \beta }}^t-r\Vert {\overline{w}}_{\mathrm{\alpha \beta }}^t\Vert ,\\ & & (\alpha ,\beta )\in \mathcal{H},t=0,1,\dots ,t_k-1,\\ & & [({\rho }_{\mathrm{\alpha \beta }}^{\mathrm{observed}}+\eta (t-t_k)+\mathrm{\Gamma }(t-t_k)+\nu ]z^t\\ & & \le {\overline{\alpha }}_{\mathrm{\alpha \beta }}^t-r\Vert {\overline{w}}_{\mathrm{\alpha \beta }}^t\Vert ,(\alpha ,\beta )\in \mathcal{H},t=t_k,\dots ,T,\\ & & {\displaystyle \sum _{t=0}^T{z}^t\ge {\gamma }_{\mathrm{lb}},z^t\ge 0,}\\ & & {\displaystyle {\overline{\alpha }}_{\mathrm{\alpha \beta }}=\sum _{t=0}^T\sum _{p=1}^n{\eta }_p^t{\overline{d}}_{\mathrm{\alpha \beta }}^{p,t},{\overline{\alpha }}_{\mathrm{\alpha \beta }}^t=\sum _{p=1}^n{\eta }_p^t{\overline{d}}_{\mathrm{\alpha \beta }}^{p,t},}\\ & & {\displaystyle {\overline{w}}_{\mathrm{\alpha \beta }}=\sum _{t=0}^T\sum _{p=1}^n{\eta }_p^t{\overline{v}}_{\mathrm{\alpha \beta }},{\overline{w}}_{\mathrm{\alpha \beta }}^t=\sum _{p=1}^n{\eta }_p^t{\overline{v}}_{\mathrm{\alpha \beta }}.}\end{array}\end{array}$

Figure A1. PTV Dose Volume Histogram.

Figure A2. Heart Dose Volume Histogram.