Robust Optimal Power Distribution for Hyperthermia Cancer Treatment

Abstract

We consider an optimization problem for spatial power distribution generated by an array of transmitting elements. Using ultrasound hyperthermia cancer treatment as a motivating example, the signal design problem consists of optimizing the power distribution across the tumor and healthy tissue regions, respectively. The models used in the optimization problem are, however, invariably subject to errors. To combat such unknown model errors, we formulate a robust signal design framework that can take the uncertainty into account using a worst-case approach. This leads to a semi-infinite programming (SIP) robust design problem, which we reformulate as a tractable convex problem that potentially has a wider range of applications.

1. Introduction

Local hyperthermia is a noninvasive technique for cancer treatment, in which targeted body tissue is exposed to high temperatures to damage cancer cells, leaving surrounding tissue unharmed. This technique is used both to kill-off cancer cells in tumors and as a means to enhance other treatments such as radiotherapy and chemotherapy. Hyperthermia has the potential to treat many types of cancer, including sarcoma, melanoma, and cancers of the head and neck, brain, lung, esophagus, breast, bladder, rectum, liver, appendix, cervix, etc. [1, 2, 3].

Hyperthermia treatment planning involves modeling patient-specific tissue, using medical imaging techniques such as microwave, ultrasound, magnetic resonance or computed tomography, and calculating the spatial distribution of power deposited in the tissue to heat it [4]. There exist two major techniques to concentrate the power in a well-defined tumor region: electromagnetic and ultrasound, each with its own limitations. The drawback of electromagnetic microwaves is its poor penetration in biological tissue; while for ultrasound, the short acoustic wavelength renders the focal spot very small. Using signal design methods, however, one can improve the spatial power deposition generated by an array of acoustic transducers. Specifically, standard phased array techniques do not make use of combining a diversity of signals transmitted at each transducer. When this diversity is exploited, it is possible to dramatically improve the power distribution in the tumor tissue, thus improving the effectiveness of the method and reducing treatment time [5, 6]. Given a set of spatial coordinates that describe the tumor region and the healthy tissue, respectively, the transmitted waveforms can be designed to optimize the spatial power distribution while subject to certain design constraints.

One critical limitation, however, is the assumption of an ideal wave propagation model from the transducers to a given point in the tissue. Specifically, model mismatches may arise from hardware imperfections, tissue inhomogeneities, inaccurately specified propagation velocities, etc. Thus, the actual power distribution may differ substantially from the ideal one designed by an assumed model. This results in suboptimal clinical outcome due to loss of power in the tumor region and safety issues due to the possible damage of healthy tissue. These considerations motivate developing robust design schemes that take such unknown errors into account.

In this chapter, we derive a robust optimization method that only assumes the unknown model errors to be bounded. The power is then optimized with respect to “worst-case” model errors. By using a worst-case model, we provide an optimal signal design scheme that takes into account all possible, bounded model errors. Such a conservative approach is warranted in signal design for medical applications due to safety and health considerations. Our method further generalizes the approach in [5] by obviating the need to specify a fictitious tumor center point. The framework developed here has potential use in wider signal design applications where the resulting transmit power distributions are subject to model inaccuracies. More specifically, the design problem formulated in this chapter and the proposed robust scheme can be exploited to robustify the spatial power distribution for applications that an array equipped with multiple elements is used to emit waveforms in order to deliver power to an area of interest in a controlled manner.

The core of this study is built upon exploiting waveform diversity which has been introduced in multiple-input multiple-output (MIMO) radar literature [7], and later has been applied for local hyperthermia cancer treatment improvement in [5]. In the MIMO radar field, robustness studies have been carried out in different applications under varying design parameter uncertainties, cf., [8, 9]. Recently, in [10], we have studied the robustification of the waveform diversity methodology for MIMO radar applications. It should be highlighted that, in this chapter, a more generic problem formulation has been studied with respect to those of [10], where a new application area is considered to illustrate the performance of our proposed robust design. In the array processing literature, beamforming under array model errors has also spawned extensive work, cf., [11, 12, 13, 14].

For hyperthermia therapy, the need for robust solutions when optimizing for phase and amplitude of conventional phased array has been investigated in [15], considering perfusion uncertainties, and in [16] considering dielectric uncertainties. The authors emphasize on the role of uncertainty in such designs (hyperthermia planning) since it influences the calculation of power distribution, and correspondingly temperature distribution.

The chapter is organized as follows: in Section 2, we describe the system model and the relevant variables. In Section 3, the signal design problem is presented. First, we consider the state-of-the-art method based on “waveform diversity” [5, 7, 17], then we generalize the design problem by introducing a deterministic and bounded set of possible model errors which results in an infinite number of constraints. Importantly, we show that this seemingly intractable problem can be equivalently formulated as a tractable convex optimization problem. In Section 4, we evaluate the design scheme. We evaluate the performance of our proposed robust power distribution scheme specifically for local hyperthermia breast cancer treatment. This example application is motivated by the alarming statistics pointing to breast cancer as one of the leading causes of death among women worldwide [18, 19, 20].¹ The case of no model mismatch is investigated first, and then the robust design scheme is applied, where its power distribution in the worst-case model is evaluated and compared to the nonrobust formulation.

Notation: Boldface (lower case) is used for column vectors, $\mathbf{x}$, and (upper case) for matrices, $\mathbf{X}$. $\parallel \mathbf{a}{\parallel }_{\mathbf{W}}\triangleq \sqrt{{\mathbf{a}}^H\mathbf{Wa}}$, where $\mathbf{W}\succ \mathbf{0}$. ${\mathbf{x}}^T$ and ${\mathbf{x}}^H$ denote transpose and Hermitian transpose, respectively. $\mathbf{R}\succcurlyeq \mathbf{0}$ signifies a positive semi-definite matrix and ${\mathbf{R}}^{1/2}$ signifies a matrix square-root, e.g., Hermitian. The set of complex numbers is denoted by $\mathcal{C}$.

Abbreviations: semi-infinite programming (SIP); multiple-input multiple-output (MIMO); semidefinite program (SDP); linear matrix inequality (LMI).

2. System model

We consider an array of $M$ acoustic transducers to heat target points. These transducers are located at known positions ${\mathit{\theta }}_m$, for $m=1,2,\dots ,M$, around the tissue at risk, cf., [5, 10]. We parameterize an arbitrary point in 3D space using Cartesian coordinates $\mathbf{r}={\left[xyz\right]}^T$.

Let $x_m\left(n\right)$ denote the baseband representation of narrowband discrete-time signal transmitted at the $m$th transducer, at sample $n=1,\dots ,N$. Then, the baseband signal received at a generic location $\mathbf{r}$ equals the superposition of signals from all $M$ transducers, i.e.,

$$\begin{array}{c}y\left(\mathbf{r},n\right)=\sum _{m=1}^M{a}_m\left(\mathbf{r}\right)x_m\left(n\right),n=1,\dots ,N\\ ={\mathbf{a}}^H\left(\mathbf{r}\right)\mathbf{x}\left(n\right),n=1,\dots ,N,\end{array}$$E1

where the $m$th signal is attenuated by a factor $a_m\left(\mathbf{r}\right)$ which depends on the properties of the transducers, the carrier wave, and the tissue. This factor is modeled as

$$a_m\left(\mathbf{r}\right)=\frac{e^{-j2\pi f_c{\tau }_m\left(\mathbf{r}\right)}}{\parallel {\mathit{\theta }}_m-\mathbf{r}{\parallel }^{\frac{1}{2}}},$$E2

where $f_c$ is the carrier frequency, and ${\tau }_m\left(\mathbf{r}\right)=\frac{\parallel {\mathit{\theta }}_m-\mathbf{r}\parallel }{c}$ is the required time for any signal to arrive at location $\mathbf{r}$ where $c$ is the sound speed inside the tissue. Note that the root-squared term in the denominator of (2) represents the distance-dependent propagation attenuation of the acoustic waveforms. In (1), the narrowband signals are represented in vector form $\mathbf{x}\left(n\right)={\left[x_1\left(n\right)\dots x_m\left(n\right)\dots x_M\left(n\right)\right]}^T\in {\mathcal{C}}^{M\times 1}$ and $\mathbf{a}\left(\mathbf{r}\right)\triangleq {\left[a_1\left(\mathbf{r}\right)\dots a_m\left(\mathbf{r}\right)\dots a_M\left(\mathbf{r}\right)\right]}^T\in {\mathcal{C}}^{M\times 1}$ is the array steering vector as a function of $\mathbf{r}$.

At a generic location $\mathbf{r}$ in the tissue, the power of the transmitted signal, i.e., the transmit beampattern, is given by

$$p\left(\mathbf{r}\right)=\mathbb{E}\left\{{\left|y(\mathbf{r},n)\right|}^2\right\}={\mathbf{a}}^H\left(\mathbf{r}\right)\mathbf{Ra}\left(\mathbf{r}\right),$$E3

where

$$\mathbf{R}\triangleq \mathbb{E}\left\{\mathbf{x}\left(n\right){\mathbf{x}}^H\left(n\right)\right\}$$E9000

is the $M\times M$ covariance matrix of the signal $\mathbf{x}\left(n\right)$. As Eq. (3) suggests, the transmit beampattern is dependent on the waveform covariance matrix $\mathbf{R}$ and the array steering vector $\mathbf{a}\left(\mathbf{r}\right)$. In the following, we analyze how one can form and control the beampattern by optimizing the covariance matrix $\mathbf{R}$, so as to heat up the tumor region of the tissue while keeping the power deposition in the healthy tissue minimal. In this work, we consider schemes which allow for the lowest possible power leakage to the healthy area.

Once an optimal covariance matrix $\mathbf{R}$ has been determined, the waveform signal $\mathbf{x}\left(n\right)$ can be synthesized accordingly. One simple approach is $\mathbf{x}\left(n\right)={\mathbf{R}}^{1/2}\mathbf{w}\left(n\right)$, where $\mathbf{w}\left(n\right)$ is a sequence of independent random vectors with mean zero and covariance matrix $\mathbf{I}$. For detailed discussion see [21, 22, 23, Ch. 14].

A significant challenge to this approach, however, is that the true steering vector $\mathbf{a}\left(\mathbf{r}\right)$ in (3) does not exactly match the model in (2) for a host of reasons: array calibration imperfections, variations in transducing elements, tissue inhomogeneities, inaccurately specified propagation velocity, etc. We will therefore consider the aforementioned design problem subject to model uncertainties in the array steering vector at any given point $\mathbf{r}$. We refer to this approach as robust waveform diversity.

3. Problem formulation

The waveform-diversity-based technique [5, 7, 10, 22, 24] have been used for designing beampatterns (3) subject to practical constraints. In general, we aim to control and shape the spatial power distribution at a set of target points while simultaneously minimizing power leakage in the remaining area. By exploiting a combination of different waveforms in (1), the degrees of freedom increase for optimizing the beampattern under constraints.

After reviewing the standard waveform diversity approach, we focus on the practical scenario where the assumed array steering vector model is subject to perturbations. In the subsequent section, the proposed robust technique is evaluated by numerical simulations, comparing the performance with and without robustified solution under perturbed steering vectors.

3.1. Waveform-Diversity-based Ultrasound System

In the MIMO radar literature, sidelobe minimization is a beampattern design problem that has been addressed by using the waveform diversity methodology, cf., [7, 21, 22, 24]. This design problem can be thought of as an optimization problem, where the probing waveforms covariance matrix $\mathbf{R}$ is the optimization variable to be chosen under positive semi-definiteness assumption and with a constraint on the total power. The waveform-diversity-based scheme for ultrasound system has been introduced and explained in detail in [5] based on the transmit beampattern design technique for MIMO radar systems [7, 24].

In the following, we consider the practical power constraint, where all array elements have the same power. Therefore, the covariance matrix $\mathbf{R}$ belongs to the following set $\mathcal{R}$:

$$\mathcal{R}\triangleq \left\{\mathbf{R}|\mathbf{R}\succcurlyeq \mathbf{0},R_{\mathit{mm}}=\frac{\gamma }{M},m=1,2,\dots ,M\right\},$$E4

where $\gamma$ is the total transmitted power and $R_{\mathit{mm}}$ is the $m$th diagonal element of $\mathbf{R}$ corresponding to the power emitted by $m$th transducer. The healthy tissue and the tumor regions are represented by two sets of discrete control points $\mathbf{r}$:

$${\mathrm{\Omega }}_S=\left\{{\mathbf{r}}_1,{\mathbf{r}}_2,\dots ,{\mathbf{r}}_{N_S}\right\}$$E5

$${\mathrm{\Omega }}_T=\left\{{\mathbf{r}}_1,{\mathbf{r}}_2,\dots ,{\mathbf{r}}_{N_T}\right\},$$E6

where $N_S$ and $N_T$ denote the number of points in the healthy tissue region and the tumor regions, respectively. Without loss of generality, let ${\mathbf{r}}_0$ be a representative point which is taken to be the center of the tumor region ${\mathrm{\Omega }}_T$. The objectives for this optimization problem can be summarized as follows: design the waveform covariance matrix $\mathbf{R}$ so as to

• maximize the gap between the power at the tumor center ${\mathbf{r}}_0$ and the power at the control points $\mathbf{r}$ in the healthy tissue region ${\mathrm{\Omega }}_S$;

• while guaranteeing a certain power level for control points $\mathbf{r}$ in the tumor region ${\mathrm{\Omega }}_T$.

Mathematically, this problem is formulated as (see [5])

$$\begin{array}{lll}\underset{\mathbf{R},t}{\text{max}}& t& \\ \mathrm{s}.\mathrm{t}.& {\mathbf{a}}^H\left({\mathbf{r}}_0\right)\mathbf{Ra}\left({\mathbf{r}}_0\right)-{\mathbf{a}}^H\left(\mathbf{r}\right)\mathbf{Ra}\left(\mathbf{r}\right)\ge t,\forall \mathbf{r}\in {\mathrm{\Omega }}_S& \\ & {\mathbf{a}}^H\left(\mathbf{r}\right)\mathbf{Ra}\left(\mathbf{r}\right)\ge \left(1-\delta \right){\mathbf{a}}^H\left({\mathbf{r}}_0\right)\mathbf{Ra}\left({\mathbf{r}}_0\right),\forall \mathbf{r}\in {\mathrm{\Omega }}_T& \\ & {\mathbf{a}}^H\left(\mathbf{r}\right)\mathbf{Ra}\left(\mathbf{r}\right)\le \left(1+\delta \right){\mathbf{a}}^H\left({\mathbf{r}}_0\right)\mathbf{Ra}\left({\mathbf{r}}_0\right),\forall \mathbf{r}\in {\mathrm{\Omega }}_T& \\ & \mathbf{R}\in \mathcal{R}& \end{array}$$E7

where $t$ denotes the gap between the power at ${\mathbf{r}}_0$ and the power at the control points $\mathbf{r}$ in the healthy region ${\mathrm{\Omega }}_S$. The parameter $\delta$ is introduced here to control the required certain power level at the control points in the tumor region. For instance, if we set $\delta =0.1$, then we aim for having power at the tumor region ${\mathrm{\Omega }}_T$ to be within 10% of $p\left({\mathbf{r}}_0\right)$, i.e., the power at the tumor center. This is an SDP problem which can be solved efficiently in polynomial time using any SDP solver, e.g., CVX [25, 26].

3.2. Robust waveform-diversity-based ultrasound system

The convex optimization problem (7) and consequently its optimal solution, i.e., the optimal covariance matrix $\mathbf{R}$, are functions of the steering vectors $\mathbf{a}\left(\mathbf{r}\right)$. In practice, however, the assumed steering vector model used to optimize $\mathbf{R}$ is inaccurate. Hence, using nominal steering vectors $\widehat{\mathbf{a}}\left(\mathbf{r}\right)$ based on an ideal model, in lieu of the unknown true steering vectors $\mathbf{a}\left(\mathbf{r}\right)$ in (7), may result in undesired beampatterns with low power at the tumor region and damaging power deposition in the healthy tissue region. Such health considerations in medical applications motivate an approach that is robust with respect to the worst-case model uncertainties.

In order to formulate the robust design problem mathematically, we parameterize the steering vector uncertainties as follows. Let the true steering vector for the transducer array be $\mathbf{a}\left(\mathbf{r}\right)=\widehat{\mathbf{a}}\left(\mathbf{r}\right)+\tilde{\mathbf{a}}\left(\mathbf{r}\right)$, where $\tilde{\mathbf{a}}\left(\mathbf{r}\right)$ is an unknown perturbation from the nominal steering vector. The deterministic perturbation at any generic point $\mathbf{r}$ belongs to the uncertainty set ${\mathcal{E}}_{\mathbf{r}}$ that is bounded

$${\mathcal{E}}_{\mathbf{r}}\triangleq \left\{\tilde{\mathbf{a}}\left(\mathbf{r}\right)|\parallel \tilde{\mathbf{a}}\left(\mathbf{r}\right){\parallel }_{\mathbf{W}}^2\le {\epsilon }_{\mathbf{r}}\right\},$$E8

where $\mathbf{W}$ is a $M\times M$ diagonal weight matrix with positive elements. The weight matrix $\mathbf{W}$ can be derived based on the type of uncertainty. Using $\mathbf{W}$, the set ${\mathcal{E}}_{\mathbf{r}}$ indicates an ellipsoidal region. The bound ${\mathrm{ϵ}}_{\mathbf{r}}$ for the set can be a constant or a function of $\mathbf{r}$, i.e., ${\mathrm{ϵ}}_{\mathbf{r}}=f\left(\mathbf{r}\right)$. This set enables parameterization of element-wise uncertainties in the nominal steering vector $\widehat{\mathbf{a}}\left(\mathbf{r}\right)$ at each $\mathbf{r}$.

Besides this consideration, we generalize the problem formulation (7) further by setting a uniform bound (power level) $P$ across the tumor region ${\mathrm{\Omega }}_T$ as an optimization variable to which the power of all the control points in the healthy region ${\mathrm{\Omega }}_S$ are compared. This is in contrast to (7) and the robust formulation in [10], where the power levels of all the healthy grid points ${\mathrm{\Omega }}_S$ are compared with the power of only a single reference point at fictitious tumor center ${\mathbf{r}}_0$. There is no need to limit our problem to a single point as a reference power level. Rather, the desired tightness of the power level across ${\mathrm{\Omega }}_T$ is specified by the parameter $0\le \delta <1$. This generalization also improves the efficiency when it comes to solving the robust design problem.

With these considerations, the robust beampattern design problem can be formulated as

$$\begin{array}{l}\underset{\mathbf{R},t,P}{\text{max}}t\text{subject to}\hfill \\ P-{\left(\widehat{\mathbf{a}}\left(\mathbf{r}\right)+\tilde{\mathbf{a}}\left(\mathbf{r}\right)\right)}^H\mathbf{R}\left(\widehat{\mathbf{a}}\left(\mathbf{r}\right)+\tilde{\mathbf{a}}\left(\mathbf{r}\right)\right)\ge t,\forall \tilde{\mathbf{a}}\left(\mathbf{r}\right)\in {\mathcal{E}}_{\mathbf{r}},\mathbf{r}\in {\mathrm{\Omega }}_S\hfill \\ {\left(\widehat{\mathbf{a}}\left(\mathbf{r}\right)+\tilde{\mathbf{a}}\left(\mathbf{r}\right)\right)}^H\mathbf{R}\left(\widehat{\mathbf{a}}\left(\mathbf{r}\right)+\tilde{\mathbf{a}}\left(\mathbf{r}\right)\right)\ge \left(1-\delta \right)P,\forall \tilde{\mathbf{a}}\left(\mathbf{r}\right)\in {\mathcal{E}}_{\mathbf{r}},\mathbf{r}\in {\mathrm{\Omega }}_T\hfill \\ {\left(\widehat{\mathbf{a}}\left(\mathbf{r}\right)+\tilde{\mathbf{a}}\left(\mathbf{r}\right)\right)}^H\mathbf{R}\left(\widehat{\mathbf{a}}\left(\mathbf{r}\right)+\tilde{\mathbf{a}}\left(\mathbf{r}\right)\right)\le \left(1+\delta \right)P,\forall \tilde{\mathbf{a}}\left(\mathbf{r}\right)\in {\mathcal{E}}_{\mathbf{r}},\mathbf{r}\in {\mathrm{\Omega }}_T\hfill \\ \mathbf{R}\in \mathcal{R},\hfill \end{array}$$E9

where $t$ is the gap between the desired power levels set across ${\mathrm{\Omega }}_T$ and power deposition in the healthy tissue ${\mathrm{\Omega }}_S$, similar to (7). Note that we take into account every possible perturbation $\tilde{\mathbf{a}}\left(\mathbf{r}\right)\in {\mathcal{E}}_{\mathbf{r}}$.

In contrast to the optimization problem (7), which is a tractable convex problem, the robust problem (9) is an SIP problem. For a given $\mathbf{R}$ in (9), there are infinite number of constraints in terms of $\tilde{\mathbf{a}}\left(\mathbf{r}\right)$ to satisfy which makes the problem non-trivial. However, in the following theorem, extending the approach in [10], we reformulate the robust power deposition problem as a convex SDP problem whose solution is the optimally robust covariance matrix.

Theorem 1. The robust power deposition for an M-element transducer array with the probing signal covariance matrix $\mathbf{R}\in \mathcal{R}$ and the perturbation vector $\tilde{\mathbf{a}}\left(\mathbf{r}\right)\in {\mathcal{E}}_{\mathbf{r}}$, i.e., the solution of (9), is given as a solution to the following SDP problem

$$\begin{array}{lll}& \underset{\mathbf{R},t,P,{\beta }_i,{\beta }_{j,1},{\beta }_{j,2}}{\text{max}}t\text{subject to}& \\ & {\mathrm{\Omega }}_S:\left[\begin{array}{ll}{\beta }_i\mathbf{W}-\mathbf{R}& -\mathbf{R}\widehat{\mathbf{a}}\left({\mathbf{r}}_i\right)\\ -\widehat{\mathbf{a}}{\left({\mathbf{r}}_i\right)}^H\mathbf{R}& P-t-\widehat{\mathbf{a}}{\left({\mathbf{r}}_i\right)}^H\mathbf{R}\widehat{\mathbf{a}}\left({\mathbf{r}}_i\right)-{\beta }_i{\epsilon }_{{\mathbf{r}}_i}\end{array}\right]\succcurlyeq \mathbf{0},& \\ & {\mathrm{\Omega }}_T:\left[\begin{array}{ll}{\beta }_{j,1}\mathbf{W}+\mathbf{R}& \mathbf{R}\widehat{\mathbf{a}}\left({\mathbf{r}}_j\right)\\ \widehat{\mathbf{a}}{\left({\mathbf{r}}_j\right)}^H\mathbf{R}& \widehat{\mathbf{a}}{\left({\mathbf{r}}_j\right)}^H\mathbf{R}\widehat{\mathbf{a}}\left({\mathbf{r}}_j\right)-\left(1-\delta \right)P-{\beta }_{j,1}{\epsilon }_{{\mathbf{r}}_j}\end{array}\right]\succcurlyeq \mathbf{0},& \\ & {\mathrm{\Omega }}_T:\left[\begin{array}{ll}{\beta }_{j,2}\mathbf{W}-\mathbf{R}& -\mathbf{R}\widehat{\mathbf{a}}\left({\mathbf{r}}_j\right)\\ -\widehat{\mathbf{a}}{\left({\mathbf{r}}_j\right)}^H\mathbf{R}& \left(1+\delta \right)P-\widehat{\mathbf{a}}{\left({\mathbf{r}}_j\right)}^H\mathbf{R}\widehat{\mathbf{a}}\left({\mathbf{r}}_j\right)-{\beta }_{j,2}{\epsilon }_{{\mathbf{r}}_j}\end{array}\right]\succcurlyeq \mathbf{0},& \\ & \mathbf{R}\in \mathcal{R},{\beta }_i,{\beta }_{j,1},{\beta }_{j,2}\ge 0,i=1,\dots ,N_S,j=1,\dots ,N_T.& \end{array}$$E10

Proof: See Appendix A.

Observe that the notations ${\mathrm{\Omega }}_S$ and ${\mathrm{\Omega }}_T$ indicate that the corresponding linear matrix inequalities (LMIs) should be satisfied for the points ${\mathbf{r}}_i\in {\mathrm{\Omega }}_S$ and ${\mathbf{r}}_j\in {\mathrm{\Omega }}_T$, respectively. Note that the robust SDP problem in this chapter, which is stated in Theorem 1, can be solved more efficiently than the SDP problem in [10] since the matrices $\mathbf{R}$ and $\mathbf{W}$ in the current formulation have half of the size of the matrices involved in the latter problem. This occurs due to the generalization of the robust problem by using the uniform power level as a benchmark.

Note that other robust problems with similar objectives can also be addressed using the above approach which are outlined in the following subsection.

3.3. Alternative robust formulations

Similar robust problems to that of (9) can be formulated in many different ways. For example, by restricting the power level outside the tumor in a weighted fashion.

$$\begin{array}{l}\underset{t,R}{\text{min}}t\text{subject}\text{to}\hfill \\ {\left(\widehat{\mathbf{a}}\left(\mathbf{r}\right)+\tilde{\mathbf{a}}\left(\mathbf{r}\right)\right)}^{H}R\left(\widehat{\mathbf{a}}\left(\mathbf{r}\right)+\tilde{\mathbf{a}}\left(\mathbf{r}\right)\right)\le \mathit{tw}\left(\mathbf{r}\right),\forall \tilde{\mathbf{a}}\left(\mathbf{r}\right)\in {\mathcal{E}}_{\mathbf{r}},\\ \mathbf{r}\in {\mathrm{\Omega }}_{\mathrm{S}}\hfill \\ {\left(\widehat{\mathbf{a}}\left(\mathbf{r}\right)+\tilde{\mathbf{a}}\left(\mathbf{r}\right)\right)}^{H}R\left(\widehat{\mathbf{a}}\left(\mathbf{r}\right)+\tilde{\mathbf{a}}\left(\mathbf{r}\right)\right)\ge \left(1-\delta \right)P,\forall \tilde{\mathbf{a}}\left(\mathbf{r}\right)\in {\mathcal{E}}_{\mathbf{r}},\\ \mathbf{r}\in {\mathrm{\Omega }}_{\mathrm{T}}\hfill \\ {\left(\widehat{\mathbf{a}}\left(\mathbf{r}\right)+\tilde{\mathbf{a}}\left(\mathbf{r}\right)\right)}^{H}R\left(\widehat{\mathbf{a}}\left(\mathbf{r}\right)+\tilde{\mathbf{a}}\left(\mathbf{r}\right)\right)\le \left(1+\delta \right)P,\forall \tilde{\mathbf{a}}\left(\mathbf{r}\right)\in {\mathcal{E}}_{\mathbf{r}},\mathbf{r}\in {\mathrm{\Omega }}_{\mathrm{T}}\\ R\in \mathcal{R}\end{array}$$E11

where $P,\delta$ are fixed and $w\left(\mathbf{r}\right)$ is a weighting function constructed, e.g., so that the energy bound close to the tumor is less restrictive.

One could also construct problems that minimize the sum of the energy in the non-tumor area, where $t\left(\mathbf{r}\right)$ denotes the energy at $\mathbf{r}$:

$$\begin{array}{l}\underset{t\left(\mathbf{r}\right),R}{\text{min}}\sum _{\mathbf{r}\in {\mathrm{\Omega }}_{\mathrm{S}}}t\left(\mathbf{r}\right)\text{subject}\text{to}\hfill \\ {\left(\widehat{\mathbf{a}}\left(\mathbf{r}\right)+\tilde{\mathbf{a}}\left(\mathbf{r}\right)\right)}^{H}R\left(\widehat{\mathbf{a}}\left(\mathbf{r}\right)+\tilde{\mathbf{a}}\left(\mathbf{r}\right)\right)\le t\left(\mathbf{r}\right),\forall \tilde{\mathbf{a}}\left(\mathbf{r}\right)\in {\mathcal{E}}_{\mathbf{r}},\\ \mathbf{r}\in {\mathrm{\Omega }}_{\mathrm{S}}\hfill \\ {\left(\widehat{\mathbf{a}}\left(\mathbf{r}\right)+\tilde{\mathbf{a}}\left(\mathbf{r}\right)\right)}^{H}R\left(\widehat{\mathbf{a}}\left(\mathbf{r}\right)+\tilde{\mathbf{a}}\left(\mathbf{r}\right)\right)\ge \left(1-\delta \right)P,\forall \tilde{\mathbf{a}}\left(\mathbf{r}\right)\in {\mathcal{E}}_{\mathbf{r}},\\ \mathbf{r}\in {\mathrm{\Omega }}_{\mathrm{T}}\hfill \\ {\left(\widehat{\mathbf{a}}\left(\mathbf{r}\right)+\tilde{\mathbf{a}}\left(\mathbf{r}\right)\right)}^{H}R\left(\widehat{\mathbf{a}}\left(\mathbf{r}\right)+\tilde{\mathbf{a}}\left(\mathbf{r}\right)\right)\le \left(1+\delta \right)P,\forall \tilde{\mathbf{a}}\left(\mathbf{r}\right)\in {\mathcal{E}}_{\mathbf{r}},\mathbf{r}\in {\mathrm{\Omega }}_{\mathrm{T}}\\ R\in \mathcal{R}.\end{array}$$E12

Both of the alternative formulations described above can be addressed following the steps derived in Appendix A by using $\mathcal{S}$-lemma, since we are still dealing with quadratic constraints.

In the next section, we illustrate the reference performance of a nominal scenario where the steering vectors are perfectly known. Then, we observe how much power can leak to the healthy tissue and cause damages when subject to uncertain steering vectors. Finally, we evaluate the proposed robust scheme in terms of improving the power deposition along our stated design goals.

4. Numerical results

To illustrate the performance of the proposed robust scheme, we consider a 2D model of the organ at risk. Here, similar to [5], we focus on the ultrasonic hyperthermia treatment for breast cancer where a 10-cm-diameter semicircle is assumed to model breast tissues with a 16-mm-diameter tumor embedded inside. The tumor center is located at ${\mathbf{r}}_0={\left[034\right]}^T$ mm. Figure 1 shows this schematic model. We consider a curvilinear array with $M=51$ acoustic transducers and half wavelength element spacing. Acoustic waveforms used to excite the array have the carrier frequency of $500$ kHz. The acoustic wave speed for the breast tissue is considered as $1500$ m/s.

Figure 1.

A schematic 2-D breast model with a $16$-mm embedded tumor at $\left(0,34\right)$ as a reference geometry. A curvilinear ultrasonic array with $51$ transducers is located near to the organ at risk. The ultrasonic array is used for hyperthermia treatment.

To characterize (discretize) the healthy tissue region ${\mathrm{\Omega }}_S$ and the tumor region ${\mathrm{\Omega }}_T$, two grid sets with the spacing $4$ mm are considered. For optimization, a rectangular surface of the dimension $64\times 42$ in mm is assumed symmetric around the tumor to model the healthy region ${\mathrm{\Omega }}_S$, while the grid points belonging to the circular tumor region are excluded from this surface and they model ${\mathrm{\Omega }}_T$. Overall, $174$ and $13$ number of control points are considered to characterize ${\mathrm{\Omega }}_S$ and ${\mathrm{\Omega }}_T$ in order to optimize the array beampattern.

The total transmitted power is constrained to $\gamma =1$. For simplicity, the uncertainty set ${\mathcal{E}}_{\mathbf{r}}$ is modeled with $\mathbf{W}={\mathbf{I}}_M$ and with ${\mathrm{ϵ}}_{\mathbf{r}}\equiv \mathrm{ϵ}$ for all $\mathbf{r}$, where $\mathrm{ϵ}=0.25$. Furthermore, the tightness of the desired power level in the across tumor region, $\delta$, is set to $0.7$. Note that for the small values of $\delta$ and/or large values of $\mathrm{ϵ}$, the problem may turn infeasible. In general, the feasibility of the problem depends on the value of the tightness bound $\delta$ relative to the size of the existing uncertainty in the system, i.e., the volume of the uncertainty set $\mathrm{ϵ}$, and the number of grid points $N_S$ and $N_T$ used to control the beampattern at the area of interest. When $\delta$ is too small, the desired power level across ${\mathrm{\Omega }}_T$ is close to uniform and there may not exist enough degrees of freedom for the design problem to have a solution.

For reference, the optimal covariance matrix when no uncertainty is taken into account, ${\mathbf{R}}_{\mathit{nr}}$, is obtained by solving problem (9) using only nominal steering vectors $\widehat{\mathbf{a}}\left(\mathbf{r}\right)$, i.e., $\tilde{\mathbf{a}}\left(\mathbf{r}\right)\equiv \mathbf{0}$. The optimal robust covariance matrix, denoted ${\mathbf{R}}^{\star }$, is obtained by solving (10), where $\tilde{\mathbf{a}}\left(\mathbf{r}\right)\in {\mathcal{E}}_{\mathbf{r}}$. For performance evaluation, we consider the power deposition in the tissue under the worst-case perturbations of the steering vectors. This scenario provides a lower bound to the achievable performances of all steering vector perturbations $\tilde{\mathbf{a}}\left(\mathbf{r}\right)$, which belong to the deterministic uncertainty set ${\mathcal{E}}_{\mathbf{r}}$. In other words, for the points $\mathbf{r}$ in the healthy region ${\mathrm{\Omega }}_S$, the worst-case performance is rendered by the steering vectors which provide the highest power; whereas, for the points $\mathbf{r}$ in the tumor region ${\mathrm{\Omega }}_T$, those steering vectors which attain the lowest power are the ones which contribute in the worst-case performance. They are collectively referred to as the worst steering vectors. Therefore, for a given $\mathbf{R}$, either ${\mathbf{R}}_{\mathit{nr}}$ or ${\mathbf{R}}^{\star }$, the worst steering vectors for the control points $\mathbf{r}$ in ${\mathrm{\Omega }}_S$ and ${\mathrm{\Omega }}_T$ are obtained by maximizing and minimizing the transmit beampattern (3), respectively. Observe that finding the worst steering vectors for the points in the tumor region ${\mathrm{\Omega }}_T$ equals solving the following convex minimization problem at each $\mathbf{r}\in {\mathrm{\Omega }}_T$, i.e.,

$$\underset{\parallel \tilde{\mathbf{a}}\left(\mathbf{r}\right){\parallel }^2\le \epsilon }{\text{min}}{\left(\widehat{\mathbf{a}}\left(\mathbf{r}\right)+\tilde{\mathbf{a}}\left(\mathbf{r}\right)\right)}^H\mathbf{R}\left(\widehat{\mathbf{a}}\left(\mathbf{r}\right)+\tilde{\mathbf{a}}\left(\mathbf{r}\right)\right)$$E13

using CVX [25, 26]. Whereas, for finding the worst steering vectors for the points in the healthy region ${\mathrm{\Omega }}_S$, we obtain a local optimum for the following nonconvex maximization problem at each $\mathbf{r}\in {\mathrm{\Omega }}_S$, i.e.,

$$\underset{\parallel \tilde{\mathbf{a}}\left(\mathbf{r}\right){\parallel }^2\le \epsilon }{\text{max}}{\left(\widehat{\mathbf{a}}\left(\mathbf{r}\right)+\tilde{\mathbf{a}}\left(\mathbf{r}\right)\right)}^H\mathbf{R}\left(\widehat{\mathbf{a}}\left(\mathbf{r}\right)+\tilde{\mathbf{a}}\left(\mathbf{r}\right)\right),$$E14

using semidefinite relaxation techniques from [27].

We evaluate the designed beampatterns (3) plotting the spatial power distribution in decibel scale, i.e., $20{\text{log}}_{10}\left(p\left(\mathbf{r}\right)\right)$. Two different scenarios are considered, namely nominal and perturbed, to evaluate the proposed robust power distribution scheme for the ultrasonic array. In the first scenario, nominal, we assume that the array steering vectors are precisely modeled, i.e., $\tilde{\mathbf{a}}\left(\mathbf{r}\right)=\mathbf{0}$. In Figure 2, the beampattern generated by the array is plotted for the nominal scenario. This figure represents how power is spatially distributed over the organ at risk in an idealistic situation. Here, the covariance matrix of the waveforms is optimized under the assumption that the steering vectors are accurately modeled by (2), and the performance is evaluated using exactly the same steering vectors without any perturbations. The power is noticeably concentrated in the tumor region and importantly the power in the healthy tissue is several decibels lower.

Figure 2.

Power distribution (transmit beampattern in dB) for the nominal scenario, i.e., using ${\mathbf{R}}_{\mathit{nr}}$ and $\tilde{\mathbf{a}}\left(\mathbf{r}\right)\equiv 0$.

In the second scenario, perturbed, the idealistic assumptions are relaxed and model uncertainties and imperfections are taken into account. The second scenario represents the case, where the true steering vectors are perturbed versions of the nominal steering vectors $\widehat{\mathbf{a}}\left(\mathbf{r}\right)$, i.e., the true steering vector equals $\widehat{\mathbf{a}}\left(\mathbf{r}\right)+\tilde{\mathbf{a}}\left(\mathbf{r}\right)$ where $\tilde{\mathbf{a}}\left(\mathbf{r}\right)\in {\mathcal{E}}_{\mathbf{r}}$. The perturbation vectors $\tilde{\mathbf{a}}\left(\mathbf{r}\right)$ are unknown but deterministically bounded. In the following, we illustrate the worst-case performance, i.e., using the worst steering vectors to calculate the power distribution at each point. We start by illustrating the beampattern for the nonrobust covariance matrix ${\mathbf{R}}_{\mathit{nr}}$ under the worst steering vectors. Figure 3 shows how steering vector errors can degrade the array performance. Notice that in the worst-case, there is a substantial power leakage that occurs in the healthy tissue surrounding the tumor compared to Figure 2. While, in Figure 4, the robust optimal covariance matrix ${\mathbf{R}}^{\star }$, i.e., the solution to (10), is used to calculate the power for the worst steering vectors. Comparing Figures 3 and 4, we see that by taking model uncertainties into account it is possible to obtain a noticeable increase in power in the tumor region for the worst-case, and importantly, dramatic reductions of power deposited in the healthy tissue.

Figure 3.

Power distribution (transmit beampattern in dB) for the perturbed scenario, i.e., using ${\mathbf{R}}_{\mathit{nr}}$ and $\tilde{\mathbf{a}}\left(\mathbf{r}\right)\in {\mathcal{E}}_{\mathbf{r}}$.

Figure 4.

Power distribution (transmit beampattern in dB) for the perturbed scenario, i.e., using ${\mathbf{R}}^{\star }$ and $\tilde{\mathbf{a}}\left(\mathbf{r}\right)\in {\mathcal{E}}_{\mathbf{r}}$.

To finalize the numerical analysis, we provide a quantitative description for the performance of our proposed scheme summarized in Table 1. It shows the average power calculated in dB received at the tumor region ${\mathrm{\Omega }}_T$ and at the healthy region ${\mathrm{\Omega }}_S$.

Scenarios	${\mathrm{\Omega }}_T$	${\mathrm{\Omega }}_S$
Nominal, ${\mathbf{R}}_0$	$-16.54$	$-29.78$
Perturbed, ${\mathbf{R}}_0$	$-36.40$	$-11.69$
Perturbed, ${\mathbf{R}}^{\star }$	$-27.17$	$-17.43$

Table 1.

Average power for different regions.

5. Conclusion

The robust transmit signal design for optimizing spatial power distribution of a multi-antenna array is investigated. A robustness analysis is carried out to combat against inevitable uncertainty in model parameters which results in performance degradation. Such degradation occurs in practice quite often due to relying on imperfect prior and designs based upon them. Particularly, in this chapter, the transmit signal design is based on exploiting the waveform diversity property, but where errors in the array steering vector are taken into account. These errors are modeled as belonging to a deterministic set defined by a weighted norm. Then, the resulting robust signal covariance optimization problem with infinite number of constraints is translated to a convex problem which can be solved efficiently by using the $\mathcal{S}$-procedure.

Designs that are robust with respect to the worst-case are particularly vital in biomedical applications due to health risks and possible damage. Herein, we have focused on local hyperthermia therapy as one of the cancer treatments to be used either individually or along with other treatments such as radio/chemotherapy. Specifically, we consider hyperthermia treatment of breast cancer motivated by the fact that breast cancer is a major global health concern. The proposed robust signal design scheme aims to reduce unwanted power leakage into the healthy tissue surrounding the tumor while guaranteeing certain power level in the tumor region itself.

We should emphasize on the fact that the robust design problem formulation and the analysis carried out herein yielding to the robust waveforms are general enough to be exploited whenever spatial power distribution is a concern to be addressed in real world scenarios dealing with uncertainties, e.g., for radar applications.

Numerical examples representing different scenarios are given to illustrate the performance of the proposed scheme for hyperthermia therapy. We have observed significant power leakage into the healthy tissue that can occur if the design is based on uncertain model parameters. Importantly, we have shown how such damaging power deposition can be avoided using the proposed robust design for optimal spatial power distribution.

Acknowledgments

The authors would like to acknowledge Prof. Jian Li for providing an implementation of examples from [4].

A.1. Proof of Theorem 1

We start the proof by first stating the $\mathcal{S}$-Procedure lemma, which helps us to turn the optimization problem (6) with infinitely many quadratic constraints into a convex problem with finite number of LMIs.

Lemma 1. ($\mathcal{S}$-Procedure ([28], Lemma 4.1): Let $f_k\left(\mathbf{x}\right):{\mathbb{C}}^n\to \mathbb{R}$, $k=0,1$, be defined as $f_k\left(\mathbf{x}\right)={\mathbf{x}}^H{\mathbf{A}}_k\mathbf{x}+2\mathit{Re}\left\{{\mathbf{b}}_k^H\mathbf{x}\right\}+c_k$, where ${\mathbf{A}}_k={\mathbf{A}}_k^H\in {\mathbb{C}}^{n\times n},{\mathbf{b}}_k\in {\mathbb{C}}^n$, and $c_k\in \mathbb{R}$. Then, the statement (implication) $f_0\left(\mathbf{x}\right)\ge 0$ for all $\mathbf{x}\in {\mathbb{C}}^n$ such that $f_1\left(\mathbf{x}\right)\ge 0$ holds if and only if there exists $\beta \ge 0$ such that²

$$\left[\begin{array}{ll}{\mathbf{A}}_0& {\mathbf{b}}_0\\ {\mathbf{b}}_0^H& c_0\end{array}\right]-\beta \left[\begin{array}{ll}{\mathbf{A}}_1& {\mathbf{b}}_1\\ {\mathbf{b}}_1^H& c_1\end{array}\right]\succcurlyeq 0,$$E15

if there exists a point $\widehat{\mathbf{x}}$ with $f_1\left(\widehat{\mathbf{x}}\right)>0$.

The constraints in the optimization problem (9) can be rewritten as the following functions of $\tilde{\mathbf{a}}\left(\mathbf{r}\right)$ for $\mathbf{r}\in {\mathrm{\Omega }}_S$ and $\mathbf{r}\in {\mathrm{\Omega }}_T$. For notation simplicity, we only specify the set from which the control points are drawn, and we also drop $\mathbf{r}$.

$${\mathrm{\Omega }}_S:\left\{\begin{array}{l}{f}_0=-{\tilde{\mathbf{a}}}^H\mathbf{R}\tilde{\mathbf{a}}-2\mathit{Re}\left({\widehat{\mathbf{a}}}^H\mathbf{R}\tilde{\mathbf{a}}\right)-{\widehat{\mathbf{a}}}^H\mathbf{R}\widehat{\mathbf{a}}-t+P\ge 0\\ f_1=-{\tilde{\mathbf{a}}}^H\mathbf{W}\tilde{\mathbf{a}}+{\epsilon }_{\mathbf{r}}\ge 0\end{array}\right.$$

$${\mathrm{\Omega }}_T:\left\{\begin{array}{l}{f}_0={\tilde{\mathbf{a}}}^H\mathbf{R}\tilde{\mathbf{a}}+2\mathit{Re}\left({\widehat{\mathbf{a}}}^H\mathbf{R}\tilde{\mathbf{a}}\right)+{\widehat{\mathbf{a}}}^H\mathbf{R}\widehat{\mathbf{a}}-\left(1-\delta \right)P\ge 0\\ f_1=-{\tilde{\mathbf{a}}}^H\mathbf{W}\tilde{\mathbf{a}}+{\epsilon }_{\mathbf{r}}\ge 0\end{array}\right.$$

$${\mathrm{\Omega }}_T:\left\{\begin{array}{l}{f}_0=-{\tilde{\mathbf{a}}}^H\mathbf{R}\tilde{\mathbf{a}}-2\mathit{Re}\left({\widehat{\mathbf{a}}}^H\mathbf{R}\tilde{\mathbf{a}}\right)-{\widehat{\mathbf{a}}}^H\mathbf{R}\widehat{\mathbf{a}}+\left(1+\delta \right)P\ge 0\\ f_1=-{\tilde{\mathbf{a}}}^H\mathbf{W}\tilde{\mathbf{a}}+{\epsilon }_{\mathbf{r}}\ge 0\end{array}\right.$$E16

Now, according to the $\mathcal{S}$-Procedure lemma, each pair of the quadratic constraints above is replaced with an LMI for each grid points in the pre-defined sets. In other words, all these quadratic constraints are satisfied simultaneously if we find ${\beta }_i$ for $i=1,\dots ,N_S$, ${\beta }_{j,1}$ and ${\beta }_{j,2}$ for $j=1,\dots ,N_T$, for which the mentioned LMIs in Theorem 1 holds. Thus, the problem boils down to the SDP problem (10) with $2N_T+N_S$ LMIs of the size $\left(M+1\right)\times \left(M+1\right)$ as the constraints. $\square$