Controlling the Kelvin Force: Basic Strategies and Applications to Magnetic Drug Targeting

TL;DR

This paper presents a control strategy to generate a nearly constant Kelvin magnetic force for targeted drug delivery, using optimization of magnetic sources to steer drug concentration effectively along prescribed paths.

Contribution

It introduces a novel optimization-based approach to control magnetic fields for precise drug targeting, considering realistic dipole sources and PDE-governed drug dynamics.

Findings

Effective steering of drug concentration demonstrated

Magnetic field control achieves limited spreading

Optimization approach adapts to prescribed trajectories

Abstract

Motivated by problems arising in magnetic drug targeting, we propose to generate an almost constant Kelvin (magnetic) force in a target subdomain, moving along a prescribed trajectory. This is carried out by solving a minimization problem with a tracking type cost functional. The magnetic sources are assumed to be dipoles and the control variables are the magnetic field intensity, the source location and the magnetic field direction. The resulting magnetic field is shown to effectively steer the drug concentration, governed by a drift-diffusion PDE, from an initial to a desired location with limited spreading.

Equations

$${\partial }_{t}c+\mathrm{div}\left(-A\nabla c+(\mathbf{u}+\gamma \mathbf{f}(\mathbf{h}))c\right)=0\text{mbox}in\tilde{\Omega }\times (0,T)$$

$$(-A\nabla c+(\mathbf{u}+\gamma \mathbf{f}(\mathbf{h}))c)\cdot \mathbf{n}=0\text{mbox}on\partial \tilde{\Omega }\times (0,T)c(x,0)=c_0\text{mbox}in\tilde{\Omega }$$

$$\mathrm{curl}\mathbf{h}=0\text{mbox}in\tilde{\Omega }\mathrm{div}\left(\mu \mathbf{h}\right)=0\text{mbox}in\tilde{\Omega }$$

$$\mathbf{f}(\mathbf{h})=\nu \nabla |\mathbf{h}{|}^2,$$

$$\mathrm{div}\mathbf{f}(\mathbf{h})=\mathrm{div}(\nabla |\mathbf{h}{|}^2)=2\Delta \mathbf{h}+2|\nabla \mathbf{h}{|}^2=2|\nabla \mathbf{h}{|}^2\ge 0\text{in }\tilde{\Omega }$$

$${\int }_{\partial D}\mathbf{f}(\mathbf{h})\cdot \mathbf{n}ds={\int }_D\mathrm{div}\mathbf{f}(\mathbf{h})dx\ge 0.$$

$$\mathbf{f}(\mathbf{h})$$

$$\mathrm{curl}\mathbf{h}$$

$$\mathbf{h}=\sum _{i=1}^{n_p}{\alpha }_i(t)\left(d\frac{{\mathbf{r}}_i{\mathbf{r}}_i^{\top }}{|{\mathbf{r}}_i{|}^2}-\mathbb{I}\right)\frac{\widehat{{\mathbf{d}}_i}}{|{\mathbf{r}}_i{|}^d},$$

$$\underset{(\mathbf{\alpha },\mathbf{\theta })\in {\mathcal{H}}_{ad}\times {\mathcal{V}}_{ad}}{\min }\mathcal{J}(\mathbf{h}(\mathbf{\alpha },\mathbf{\theta }))$$

$$\mathcal{J}(\mathbf{h}(\mathbf{\alpha },\mathbf{\theta }))=\frac{1}{2}{\int }_0^T\Vert \nabla |\mathbf{h}{|}^2-\mathbf{f}{\Vert }_{{\mathrm{L}}^2(D_t)}^2+\frac{\lambda }{2}{\int }_0^T|{\partial }_t\mathbf{\alpha }{|}^2+\frac{\eta }{2}{\int }_0^T|{\partial }_t\mathbf{\theta }{|}^2$$

$$\mathbf{h}(\mathbf{\alpha },\mathbf{\theta })=\sum _{i=1}^{n_p}{\alpha }_i(t)\left(d\frac{{\mathbf{r}}_i{\mathbf{r}}_i^{\top }}{|{\mathbf{r}}_i{|}^2}-\mathbb{I}\right)\frac{\widehat{{\mathbf{d}}_i}({\theta }_i)}{|{\mathbf{r}}_i{|}^d},$$

$${\mathcal{H}}_{ad}\text{mathchar}58=\left\{\mathbf{\alpha }\in [{\mathrm{H}}^1(0,T){]}^{n_p}\text{mathchar}58\mathbf{\alpha }(0)={\mathbf{\alpha }}_0,\text{mbox}and{\mathbf{\alpha }}_{\ast }\le \mathbf{\alpha }(t)\le {\mathbf{\alpha }}^{\ast },\forall t\in [0,T]\right\}$$

$${\mathcal{V}}_{ad}\text{mathchar}58=\left\{\mathbf{\theta }\in [{\mathrm{H}}^1(0,T){]}^{n_p}\text{mathchar}58\mathbf{\theta }(0)={\mathbf{\theta }}_0,\text{mbox}and{\mathbf{\theta }}_{\ast }\le \mathbf{\theta }(t)\le {\mathbf{\theta }}^{\ast },\forall t\in [0,T]\right\}.$$

$$\mathbb{D}({\bm{\theta}})=\big{(}{\bm{D}}_{1}(\theta_{1})\,{\bm{D}}_{2}(\theta_{2})\,\ldots\,{\bm{D}}_{n_{p}}(\theta_{n_{p}})\big{)}\in[\mathrm{H}^{1}(0,T;C^{\infty}(\overline{\Omega}))]^{d\times{n_{p}}}$$

$${\mathbf{D}}_i({\theta }_i)\text{mathchar}58=\left(d\frac{(\mathbf{x}-{\mathbf{x}}_i)(\mathbf{x}-{\mathbf{x}}_i{)}^{\top }}{|\mathbf{x}-{\mathbf{x}}_i{|}^2}-\mathbb{I}\right)\frac{\widehat{{\mathbf{d}}_i}({\theta }_i)}{|\mathbf{x}-{\mathbf{x}}_i{|}^d}={\mathbb{M}}_i(\mathbf{x})\widehat{{\mathbf{d}}_i}({\theta }_i).$$

$$\nabla |\mathbf{h}{|}^2$$

$$=\left(\begin{array}{c}\mathbf{\alpha }(t{)}^{\top }{\partial }_{x_1}(\mathbb{D}(\mathbf{\theta }{)}^{\top }\mathbb{D}(\mathbf{\theta }))\mathbf{\alpha }(t)\\ ⋮\\ \mathbf{\alpha }(t{)}^{\top }{\partial }_{x_d}(\mathbb{D}(\mathbf{\theta }{)}^{\top }\mathbb{D}(\mathbf{\theta }))\mathbf{\alpha }(t)\end{array}\right)=\left(\begin{array}{c}{\mathbf{\alpha }}^{\top }{\mathbb{B}}_1(\mathbf{\theta })\mathbf{\alpha }\\ ⋮\\ {\mathbf{\alpha }}^{\top }{\mathbb{B}}_d(\mathbf{\theta })\mathbf{\alpha }\end{array}\right)$$

$$\underset{(\mathbf{\alpha },\mathbf{\varphi })\in {\mathcal{H}}_{ad}\times {\mathcal{U}}_{ad}}{\min }\mathcal{F}(\mathbf{h}(\mathbf{\alpha },\mathbf{\varphi }))$$

$$\mathcal{F}(\mathbf{h}(\mathbf{\alpha },\mathbf{\varphi }))=\frac{1}{2}{\int }_0^T\Vert \nabla |\mathbf{h}{|}^2-\mathbf{f}{\Vert }_{{\mathrm{L}}^2(D_t)}^2+\frac{\lambda }{2}{\int }_0^T|{\partial }_t\mathbf{\alpha }{|}^2+\frac{\beta }{2}{\int }_0^T|{\partial }_t\mathbf{\varphi }{|}^2$$

$$\mathbf{h}(\mathbf{\alpha },\mathbf{\varphi })=\sum _{i=1}^{n_p}{\alpha }_i(t)\left(d\frac{{\mathbf{r}}_i({\varphi }_i){\mathbf{r}}_i({\varphi }_i{)}^{\top }}{|{\mathbf{r}}_i({\varphi }_i){|}^2}-\mathbb{I}\right)\frac{\widehat{{\mathbf{d}}_i}}{|{\mathbf{r}}_i({\varphi }_i){|}^d},$$

$${\mathbf{R}}_i({\varphi }_i)\text{mathchar}58=\left(d\frac{{\mathbf{r}}_i({\varphi }_i){\mathbf{r}}_i({\varphi }_i{)}^{\top }}{|{\mathbf{r}}_i({\varphi }_i){|}^2}-\mathbb{I}\right)\frac{\widehat{{\mathbf{d}}_i}}{|{\mathbf{r}}_i({\varphi }_i){|}^d},$$

$$\nabla |\mathbf{h}{|}^2$$

$$=\left(\begin{array}{c}\mathbf{\alpha }(t{)}^{\top }{\partial }_{x_1}(\mathbb{R}(\mathbf{\varphi }{)}^{\top }\mathbb{R}(\mathbf{\varphi }))\mathbf{\alpha }(t)\\ ⋮\\ \mathbf{\alpha }(t{)}^{\top }{\partial }_{x_d}(\mathbb{R}(\mathbf{\varphi }{)}^{\top }\mathbb{R}(\mathbf{\varphi }))\mathbf{\alpha }(t)\end{array}\right)=\left(\begin{array}{c}{\mathbf{\alpha }}^{\top }{\mathbb{G}}_1(\mathbf{\varphi })\mathbf{\alpha }\\ ⋮\\ {\mathbf{\alpha }}^{\top }{\mathbb{G}}_d(\mathbf{\varphi })\mathbf{\alpha }\end{array}\right)$$

$$j=\underset{n\to \infty }{\lim }\mathcal{J}({\mathbf{\alpha }}_n,{\mathbf{\theta }}_n).$$

$${\mathbf{\alpha }}_n\rightharpoonup \bar{\mathbf{\alpha }}\text{mbox}in{\mathrm{H}}^1(0,T),\bar{\mathbf{\alpha }}\in {\mathcal{H}}_{ad}\text{mbox}and{\mathbf{\theta }}_n\rightharpoonup \bar{\mathbf{\theta }}\text{mbox}in{\mathrm{H}}^1(0,T),\bar{\mathbf{\theta }}\in {\mathcal{V}}_{ad}.$$

$${\mathbf{\alpha }}_n\to \bar{\mathbf{\alpha }}\text{mbox}and{\mathbf{\theta }}_n\to \bar{\mathbf{\theta }}\text{mbox}inC([0,T]).$$

$${\int }_0^T\Vert \nabla |\mathbf{h}({\mathbf{\alpha }}_n,{\mathbf{\theta }}_n){|}^2-{\mathrm{f}}_i{\Vert }_{{\mathrm{L}}^2(D_t)}^2\to {\int }_0^T\Vert \nabla |\mathbf{h}(\bar{\mathbf{\alpha }},\bar{\mathbf{\theta }}){|}^2-{\mathrm{f}}_i{\Vert }_{{\mathrm{L}}^2(D_t)}^{2}i=1,\dots ,d,$$

$$\begin{gathered} {\mathbb{B}}_k({\mathbf{\theta }}_n{)}_{i,j}-{\mathbb{B}}_k(\bar{\mathbf{\theta }}{)}_{i,j}=({\hat{\mathbf{d}}}_i^{\top }({\mathbf{\theta }}_n(i))-{\hat{\mathbf{d}}}_i^{\top }(\bar{\mathbf{\theta }}(i)))\left({\partial }_{x_k}({\mathbb{M}}_i(\mathbf{x}{)}^{\top }{\mathbb{M}}_j(\mathbf{x}))\right){\hat{\mathbf{d}}}_j({\mathbf{\theta }}_n(j)) \\ {\hat{\mathbf{d}}}_i^{\top }(\bar{\mathbf{\theta }}(i))\left({\partial }_{x_k}({\mathbb{M}}_i(\mathbf{x}{)}^{\top }{\mathbb{M}}_j(\mathbf{x}))\right)({\hat{\mathbf{d}}}_j({\mathbf{\theta }}_n(j))-{\hat{\mathbf{d}}}_j(\bar{\mathbf{\theta }}(j))). \end{gathered}$$

$$\underset{(\mathbf{\alpha },\mathbf{\theta })\in {\mathcal{H}}_{ad}\times {\mathcal{V}}_{ad}}{\min }\mathcal{J}(\mathbf{\alpha },\mathbf{\theta })=\underset{n\to \infty }{\mathrm{liminf}}\mathcal{J}({\mathbf{\alpha }}_n,{\mathbf{\theta }}_n)\ge \mathcal{J}(\bar{\mathbf{\alpha }},\bar{\mathbf{\theta }}),$$

$$\mathbf{h}(\mathbf{\alpha },\mathbf{\theta })=\sum _{i=1}^{n_p}\mathbf{A}({\alpha }_j){\hat{\mathbf{d}}}_j({\theta }_j)$$

$$\mathbf{A}({\alpha }_i)\text{mathchar}58={\alpha }_i\left(d\frac{(\mathbf{x}-{\mathbf{x}}_i)(\mathbf{x}-{\mathbf{x}}_i{)}^{\top }}{|\mathbf{x}-{\mathbf{x}}_i{|}^2}-\mathbb{I}\right)\frac{1}{|\mathbf{x}-{\mathbf{x}}_i{|}^d},i=1,\dots ,n_p.$$

$$\nabla |\mathbf{h}{|}^2$$

$$=\sum _{m,n=1}^{n_p}\left(\begin{array}{c}{\hat{\mathbf{d}}}_m^{\top }({\theta }_m){\mathbb{Q}}_1({\alpha }_m,{\alpha }_n){\hat{\mathbf{d}}}_n({\theta }_n)\\ ⋮\\ {\hat{\mathbf{d}}}_m^{\top }({\theta }_m){\mathbb{Q}}_d({\alpha }_m,{\alpha }_n){\hat{\mathbf{d}}}_n({\theta }_n)\end{array}\right)$$

$${\mathcal{J}}^{\prime }(\bar{\mathbf{\alpha }},\bar{\mathbf{\theta }})⟨\mathbf{\alpha }-\bar{\mathbf{\alpha }},\mathbf{\theta }-\bar{\mathbf{\theta }}⟩\ge 0\forall (\mathbf{\alpha },\mathbf{\theta })\in {\mathcal{H}}_{ad}\times {\mathcal{V}}_{ad}$$

Full text

**AMS subject classification **49J20, 35Q35, 35R35, 65N30

1 Introduction

Magnetic drug targeting (MDT) is an important application of ferrofluids where drugs, with ferromagnetic particles in suspension, are injected into the body. The external magnetic field then can control the drug and subsequently the drug can target the relevant areas, for example, solid tumors (see, for instance, [22]). Mathematically this process can be modeled as follows: let the concentration of magnetic nanoparticles confined in a domain $\widetilde{\Omega}\subset\mathbb{R}^{d},d=2,3$. Let $c$ be the drug concentration and ${\bm{h}}$ the magnetic field, then the evolution of $c$ by the applied magnetic field is given by the following convection-diffusion model [14, 20]:

[TABLE]

where $\gamma$ is a constitutive constant, $A$ is a diffusion coefficient matrix, ${\bm{u}}$ is a fixed velocity vector, $\bm{n}$ is the outward unit normal, the magnetic permeability $\mu$ is assumed to be constant, and ${\bm{f}}$ is the Kelvin force applied to a magnetic particle:

[TABLE]

where $\nu>0$ is a constant depending on the volume of the particle and its permeability (see [20] for details). The two fundamental units that determine the evolution of concentration $c$ are transport and diffusion (cf. (1.1)).

The aim of magnetic drug delivery is to move a drug, initially concentrated in a small subdomain $D_{0}$ and described by the initial concentration $c_{0}$, to another suddomain $D_{t}$ (desired location) while minimizing the spreading which is characterized by the concentration $c(t)$ at time $t$. This is a challenge because magnetic fields inherently tend to disperse magnetic particles. It is known that it is not possible to concentrate magnetic nanoparticles with a static magnetic source (Earnshaw’s theorem). To understand this essential obstruction, let $D\subset\widetilde{\Omega}$ be a subdomain and observe that

[TABLE]

implies

[TABLE]

This prevents the magnetic force from focussing ${\bm{f}}({\bm{h}})\cdot\bm{n}\leq 0$ on $\partial D$, a dispersion effect we could compensate for. We are thus confronted with the following question:

[TABLE]

A first approach to tackle this problem was proposed in [20]. A feedback control of an arrangement of magnetic sources (with fixed positions) is proposed to drive a spot of distributed ferrofluid from an initial point to a target position. The authors focus on the center of mass of the concentration and its covariance matrix which measures the spreading. This leads to a constrained minimization problem where the dynamics of (1.1)-(1.2) are approximated by two coupled ordinary differential equations.

Our approach to (1.5) aims to design and study an optimization framework for steering a target subdomain $D_{t}$ from an initial to a desired location via a magnetic force ${\bm{f}}({\bm{h}})$ which is almost constant in space within $D_{t}$ and thus minimizes the drug spreading outside $D_{t}$ due to (1.4). We stress that the relative size of the diffusion coefficient matrix $A$ also influences the spreading. This paper continues the program started in [2] upon adding the following features:

1. $\bullet$

Controls: The control variables are now the magnetic field intensity, the magnetic field direction, and the position of the magnetic dipoles, instead of just the former [2]. This provides greater flexibility and controllability of the process. 2. $\bullet$

Final configurations: We present several non-trivial, though extremely simplified, examples for $d=2$ that allow us to assess the question of feasibility of MDT and of our approach for its development and automation. For instance, we study the issue of moving the concentration around obstacles and of magnetic injection while minimizing the spreading. 3. $\bullet$

Advection-diffusion PDE: For the magnetic injection example, we solve (1.1) with a realistic Neumann boundary condition (1.2), instead of the Dirichlet boundary conditions as in [2]. We also consider a smaller domain than the one considered in [20, 2] which may lead to numerical oscillations as the concentration tends to go to the boundary (see [25]). In order to deal with the advection-dominant case we have adopted the edge-averaged finite element method of [30] but other alternatives such as the implicit scheme with SUPG stabilization or the explicit schemes are equally valid. Even though we do not observe oscillations in the solution, we notice a large amount of diffusion due to the upwinding nature of the scheme. This behavior makes it difficult to control the concentration in complicated geometries, for instance, the one needed to control the flow around an obstacle.

For the second example we intend to move the concentration around an obstacle. Given that a small diffusion parameter is needed to achieve this goal, homogeneous Dirichlet boundary condition (instead of Neumann condition) is meaningful if the concentration stays away from the boundary $\partial\widetilde{\Omega}$ (see Section 6.2 for details). Moreover, in order to deal with the numerical diffusion, we consider finite element method for space discretization and the explicit Euler scheme with mass lumping for time discretization, including a correction to tackle the dispersive effects of mass lumping [15]. The resulting scheme works better than the upwind scheme.

It is worth mentioning that this approach is not only applicable for magnetic drug targeting but for many other applications where magnetic force plays a role such as: gene therapy [10], magnetized stem-cells [28], magnetic tweezers [18], lab-on-a-chip systems that include magnetic particles or fluids [13] or separation of particles [21], just to name a few.

The paper is organized as follows: in section 2 we describe the dipole approximation of the magnetic sources to be used throughout the remaining paper. Our main work begins in section 3 where we establish the optimization problems. In section 4 we discuss the numerical approximation and convergence of our discrete scheme. We provide numerical illustrations for the optimization problem in section 5.1 and conclude with a numerical scheme and several illustrative examples for (1.1)-(1.3) in section 6.

2 Dipole approximation

Let $\Omega\subset\mathbb{R}^{d}$, $d=2,3$ be open bounded fictitious domain where we intent to control the magnetic force and $T>0$ be the final time. We assume that the magnetic sources lies outside $\overline{\Omega}$, then, from the Maxwell equations it follows that the magnetic field ${\bm{h}}$ satisfies:

[TABLE]

where the last equation follows by assuming a linear relation between the magnetic induction and the magnetic field. The magnetic field generated by a current distribution and a permanent magnet can be modeled by the Biot-Savart law, which is a magnetostatic approximation. However, for simplicity, in our case we consider a dipole approximation to the magnetic source, which provides a concise and easily tractable representation of the magnetic field (see [27] for a quantification of the error associated with the dipole approximation). This approximation is commonly used for localization of objects in applications ranging from medical imaging to military. It is also extensively used in real-time control of magnetic devices in medical sciences [12, 23, 26].

From now on we will assume that the magnetic field is modeled by the superposition of a fixed number $n_{p}$ of dipoles, namely

[TABLE]

where $\bm{r}=\bm{x}-\bm{x}_{i}\in\mathbb{R}^{d}$ and $\mathbb{I}\in\mathbb{R}^{d\times d}$ is the identity matrix. In addition, $\widehat{{\bm{d}}_{i}}\in\mathbb{R}^{d}$ and $\bm{x}_{i}\in\mathbb{R}^{d}\setminus\overline{\Omega}$, $i=1,\ldots,n_{p}$, denote unit vectors and dipole positions, respectively. It is straightforward to show that the magnetic field given by (2.2) satisfies (2.1).

3 Model optimization problem

To fix ideas, we assume the time dependent domain $D_{t}$ moves along a pre-specified curve smooth curve $\mathcal{C}$ and is strictly contained in $\Omega$ for every $t\in[0,T]$. Our goal here is then to approximate a vector field $\mathbf{f}\in[\mathrm{L}^{2}(0,T;\mathrm{L}^{2}(\Omega))]^{d}$ in $D_{t}$ by the so–called Kelvin force generated by a configuration of magnetic sources. In particular, for MDT, $\mathbf{f}$ could be considered uniform on $D_{t}$. This vector field will then replace $\nabla|{\bm{h}}|^{2}$ in (1.1) and will drive the concentration to a desired location.

In [2] we propose an optimal control problem in order to study the feasibility of the dipole approximation, namely, whether we can create a Kelvin force field that allows us to control, for instance, the concentration of nanoparticles in a bounded domain $\Omega$. With this aim, a desired vector field $\mathbf{f}$ and a time dependent domain $D_{t}\subset\Omega$, $t\in[0,T]$ are assumed to be known (see Figure 1 (left)). The field $\mathbf{f}$ is approximated in $D_{t}\subset\Omega$ by the Kelvin force generated by a configuration of $n_{p}$ dipoles with fixed directions $\widehat{{\bm{d}}}_{i}$ and positions ${\bm{x}}_{i}$. Here the control is the vector magnetic field’s magnitude ${\bm{\alpha}}(t)\mathrel{\mathop{\mathchar 58\relax}}=(\alpha_{1}(t),\ldots,\alpha_{n_{p}}(t))^{T}\in\mathbb{R}^{n_{p}}$.

In this work we will approximate $\mathbf{f}\mathrel{\mathop{\mathchar 58\relax}}=(\mathrm{f}_{1},\ldots,\mathrm{f}_{d})^{\top}$ with two different dipole configurations.

3.1 Controlling intensities and directions

First we consider fixed dipole positions $\left\{\bm{x}_{i}\right\}_{i=1}^{n_{p}}$ and the controls are the magnitudes $\left\{\alpha_{i}\right\}_{i=1}^{n_{p}}$ and directions \big{\{}\widehat{{\bm{d}}}_{i}\big{\}}_{i=1}^{n_{p}} parametrized by one degree of freedom per dipole (see Figure 1 (center)), namely, if ${\bm{\alpha}}(t)\mathrel{\mathop{\mathchar 58\relax}}=(\alpha_{1}(t),\ldots,\alpha_{n_{p}}(t))^{\top}\in\mathbb{R}^{n_{p}}$ and ${\bm{\theta}}(t)\mathrel{\mathop{\mathchar 58\relax}}=(\theta_{1}(t),\ldots,\theta_{n_{p}}(t))^{\top}\in\mathbb{R}^{n_{p}}$ denote the vectors of magnetic field’s magnitude and “angles”, such that $\widehat{{\bm{d}}}_{i}=\widehat{{\bm{d}}}_{i}(\theta_{i})$ for a regular function $\widehat{{\bm{d}}}_{i}$ (for instance, $\widehat{{\bm{d}}}_{i}(\theta_{i})=\left(\cos(\theta_{i}),\sin(\theta_{i})\right)$ when $d=2$), then we want to study the problem

[TABLE]

and where $({\bm{\alpha}},{\bm{\theta}})\in\mathcal{H}_{ad}\times\mathcal{V}_{ad}$ and $\lambda$ and $\eta$ are the costs of the control. For given constant vectors ${\bm{\alpha}}_{0},{\bm{\theta}}_{0},{\bm{\alpha}}_{*},{\bm{\theta}}_{*},{\bm{\alpha}}^{*},{\bm{\theta}}^{*}\in\mathbb{R}^{n_{p}}$, we seek $({\bm{\alpha}},{\bm{\theta}})$ in the following admissible convex sets:

[TABLE]

and

[TABLE]

Notice that, in view of (3.1c), we can rewrite ${\bm{h}}$ as ${\bm{h}}({\bm{\alpha}},{\bm{\theta}})=\mathbb{D}({\bm{\theta}}){\bm{\alpha}}$ where

[TABLE]

is such that

[TABLE]

Thus, the Kelvin force in terms of ${\bm{\alpha}}$ and ${\bm{\theta}}$ is given by

[TABLE]

with $\mathbb{B}_{k}({\bm{\theta}})\mathrel{\mathop{\mathchar 58\relax}}=\partial_{x_{k}}(\mathbb{D}^{\top}({\bm{\theta}})\mathbb{D}({\bm{\theta}}))\in[\mathrm{H}^{1}(0,T;C^{\infty}(\overline{\Omega}))]^{{n_{p}}\times{n_{p}}}$.

3.2 Controlling intensities and positions

For the second configuration we consider fixed dipole directions \big{\{}\widehat{{\bm{d}}}_{i}\big{\}}_{i=1}^{n_{p}} and the controls are the magnitudes $\left\{\alpha_{i}\right\}_{i=1}^{n_{p}}$ and positions $\left\{\bm{x}_{i}\right\}_{i=1}^{n_{p}}$. In this case, we assume that each dipole $i$ moves along a prescribed trajectory $\mathcal{C}_{i}$ parameterized by a $C^{1}([0,T])$ curve $\bm{\rho}_{i}\mathrel{\mathop{\mathchar 58\relax}}[{\bm{\phi}}_{*}(i),{\bm{\phi}}^{*}(i)]\to\mathcal{C}_{i}$, $i=1,\ldots,n_{p}$ where ${\bm{\phi}}_{*},{\bm{\phi}}^{*}\in\mathbb{R}^{n_{p}}$ (see Figure 1 (right)). Thus, we study the problem

[TABLE]

where $\bm{r}_{i}(\phi_{i})=(\bm{x}-\bm{\rho}_{i}(\phi_{i}))$, ${\bm{\phi}}=(\phi_{1},\ldots,\phi_{n_{p}})$, $({\bm{\alpha}},{\bm{\phi}})\in\mathcal{H}_{ad}\times\mathcal{U}_{ad}$ and, for a given constant vector ${\bm{\phi}}_{0}\in\mathbb{R}^{n_{p}}$, we consider the admissible set \mathcal{U}_{ad}\mathrel{\mathop{\mathchar 58\relax}}=\big{\{}{\bm{\phi}}\in[\mathrm{H}^{1}(0,T)]^{n_{p}}\mathrel{\mathop{\mathchar 58\relax}}{\bm{\phi}}(0)={\bm{\phi}}_{0},\,\,\mbox{and}\,\,{\bm{\phi}}_{*}\leq{\bm{\phi}}(t)\leq{\bm{\phi}}^{*},\,\,\forall t\in[0,T]\big{\}}. We first notice that, in view of (3.3c), we can rewrite ${\bm{h}}$ as ${\bm{h}}({\bm{\alpha}},{\bm{\phi}})=\mathbb{R}({\bm{\phi}}){\bm{\alpha}}$ where $\mathbb{R}({\bm{\phi}})=\left({\bm{R}}_{1}(\phi_{1})\,{\bm{R}}_{2}(\phi_{2})\,\ldots\,{\bm{R}}_{n_{p}}(\phi_{n_{p}})\right)$ belongs to $[\mathrm{H}^{1}(0,T;C^{\infty}(\overline{\Omega}))]^{d\times{n_{p}}}$ and is given by

[TABLE]

therefore, we can write the Kelvin force in terms of ${\bm{\alpha}}$ and ${\bm{\phi}}$ as follows

[TABLE]

with $\mathbb{G}_{k}({\bm{\phi}})\mathrel{\mathop{\mathchar 58\relax}}=\partial_{x_{k}}(\mathbb{R}^{\top}({\bm{\phi}})\mathbb{R}({\bm{\phi}}))\in[\mathrm{H}^{1}(0,T;C^{\infty}(\overline{\Omega}))]^{{n_{p}}\times{n_{p}}}$, $k=1,\ldots,d$. We will show existence of solution to problems (3.1) and (3.3) via a minimizing sequence argument (see, for instance, [29]). The feasibility of the dipole approximation can be studied with the first term of $\mathcal{J}({\bm{h}}({\bm{\alpha}},{\bm{\theta}}))$ and $\mathcal{F}({\bm{h}}({\bm{\alpha}},{\bm{\phi}}))$, because the computed Kelvin force can be used to control a desired concentration. On the other hand, the last two terms in (3.1b) (and (3.3b)) will enforce a smooth evolution of the intensities and angles (and positions).

3.3 Existence and optimality conditions

Next, we focus on the existence of a solution to the minimization problem (3.1), which is nonconvex. For notational simplicity, from now on we denote by $\mathrm{V}$ both a Banach space $\mathrm{V}$ and the Banach tensor product $\mathrm{V}^{n_{p}}$.

Theorem 1** (existence of minimizers).**

There exist at least one solution $(\bar{{\bm{\alpha}}},\bar{{\bm{\theta}}})$ to the minimization problem (3.1).

Proof.

We apply the direct method of the calculus of variations. Given that $\mathcal{J}$ is bounded below by zero, we deduce that $j=\inf_{({\bm{\alpha}},{\bm{\theta}})\in\mathcal{H}_{ad}\times\mathcal{V}_{ad}}\mathcal{J}({\bm{\alpha}},{\bm{\theta}})$ is finite. We can thus construct a minimizing sequence $\{({\bm{\alpha}}_{n},{\bm{\theta}}_{n})\}_{n\in\mathbb{N}}$ such that

[TABLE]

As the sequence $\{({\bm{\alpha}}_{n},{\bm{\theta}}_{n})\}_{n\in\mathbb{N}}$ is uniformly bounded in $\mathcal{H}_{ad}\times\mathcal{V}_{ad}\subset[\mathrm{H}^{1}(0,T)]^{2}$, we can extract a (not relabeled) weakly convergent subsequence $\{{\bm{\alpha}}_{n}\}_{n\in\mathbb{N}}$ such that

[TABLE]

Moreover, according to [5, Theorem 9.16] we have

[TABLE]

To show the optimality of $\bar{{\bm{\alpha}}}$, we first consider (3.6) to get that

[TABLE]

which follows from (3.1) and the fact that

[TABLE]

and $\widehat{{\bm{d}}_{i}}({\bm{\theta}}_{n}(i))\to\widehat{{\bm{d}}_{i}}(\bar{{\bm{\theta}}}(i))$ in $[C([0,T])]^{d}$, $i=1,\ldots,n_{p}$. Notice that here we have used, for instance, ${\bm{\theta}}(i)$ to denote the $i$-th component of vector ${\bm{\theta}}$. This and the fact that the last terms in $\mathcal{J}$ are weakly lower semicontinuous (see [29, Theorem 2.12]) yields

[TABLE]

which concludes the proof. ∎

By applying the same technique as above and using that $\bm{\rho}$ belongs to $C^{1}([0,T])$, we obtain the following existence result for problem (3.3).

Theorem 2** (existence of minimizers).**

There exist at least one solution $(\bar{{\bm{\alpha}}},\bar{{\bm{\phi}}})$ to the minimization problem (3.3).

3.4 First order optimality conditions

To state first order optimality conditions we follow standard arguments because of the fact that $\mathcal{J}\mathrel{\mathop{\mathchar 58\relax}}[\mathrm{H}^{1}(0,T)]^{2}\rightarrow\mathbb{R}$ is Fréchet differentiable. We first notice that, in view of (2.2), we can also rewrite ${\bm{h}}$ as

[TABLE]

where $\mathbb{{\bm{A}}}(\alpha_{i})\in[\mathrm{H}^{1}(0,T;C^{\infty}(\overline{\Omega}))]^{d\times d}$ and is given by

[TABLE]

Thus, the Kelvin force in terms of ${\bm{\alpha}}$ and ${\bm{\theta}}$ can be written as follows

[TABLE]

with $\mathbb{Q}_{k}(\alpha_{m},\alpha_{n})\mathrel{\mathop{\mathchar 58\relax}}=\partial_{x_{k}}(\mathbb{{\bm{A}}}^{\top}(\alpha_{m})\mathbb{{\bm{A}}}(\alpha_{n}))\in[\mathrm{H}^{1}(0,T;C^{\infty}(\overline{\Omega}))]^{d\times d}$ for $k=1,\ldots,d$.

Theorem 3** (first order optimality condition for $\mathcal{J}$).**

If $(\bar{{\bm{\alpha}}},\bar{{\bm{\theta}}})\in\mathcal{H}_{ad}\times\mathcal{V}_{ad}$ denotes an optimal control, given by Theorem 1, then the first order necessary optimality condition satisfied by $(\bar{{\bm{\alpha}}},\bar{{\bm{\theta}}})$ is

[TABLE]

where, for $(\delta{\bm{\alpha}},\delta{\bm{\theta}})=({\bm{\alpha}}-\bar{{\bm{\alpha}}},{\bm{\theta}}-\bar{{\bm{\theta}}})$ we have

[TABLE]

Similarly, we obtain the first order optmimality conditions for $\mathcal{F}$.

Theorem 4** (first order optimality condition for $\mathcal{F}$).**

If $(\bar{{\bm{\alpha}}},\bar{{\bm{\phi}}})\in\mathcal{H}_{ad}\times\mathcal{U}_{ad}$ denotes an optimal control, given by Theorem 2, then the first order necessary optimality condition satisfied by $(\bar{{\bm{\alpha}}},\bar{{\bm{\phi}}})$ is

[TABLE]

where, for $(\delta{\bm{\alpha}},\delta{\bm{\phi}})=({\bm{\alpha}}-\bar{{\bm{\alpha}}},{\bm{\phi}}-\bar{{\bm{\phi}}})$ we have

[TABLE]

Remark 5** (unknown final time).**

In some applications an important quantity to account for is the time it takes for $D_{t}$ to arrive at its final destination. In this case, the final time is an unknown. This problem can be handled, following [20], by reformulating the minimization in terms of the arc length which has a fixed final value. This approach can be applied to problems (3.1) and (3.3), for details see [2].

4 Discretization

In this section we propose and study a numerical approximation of minimization problems (3.1) and (3.3). We introduce a parametrization of the moving domain $D_{t}$ in terms of a fixed domain $\widehat{D}\subset\mathbb{R}^{d}$ (see also [2]). With this in mind, we define a map $\mathbf{X}\mathrel{\mathop{\mathchar 58\relax}}[0,T]\times\overline{\widehat{D}}\to\overline{\Omega}$, such that for all $t\in[0,T]$

[TABLE]

is a one-to-one correspondence which satisfies $\mathbf{X}(t,\widehat{D})=D_{t}$. For simplicity, we assume

[TABLE]

where

[TABLE]

are functions such that $\bm{\varphi}(0)=\bm{0}$ and $\psi(0)=1$, namely $\widehat{D}=D_{0}$. In case $\psi(t)=\widehat{\psi}\in\mathbb{R}^{+}$ for all $t\in[0,T]$, $\bm{\varphi}$ can be viewed as a parameterization of a desired path that the scaled domain $\widehat{\psi}\widehat{D}$ traverses from an initial position to a final position.

4.1 Problem 1: Controlling intensities and directions

We consider the previous parametrization and rewrite the minimization problem (3.1) in terms of $\widehat{D}$. Therefore, for $i=1,\dots,d$

[TABLE]

because $\det(\nabla_{\widehat{\bm{x}}}\mathbf{X}(t,\widehat{\bm{x}}))=\psi(t)^{d}$. Then, we rewrite $\mathcal{J}$ as

[TABLE]

with $\widehat{\mathbb{B}}_{i}({\bm{\theta}},\widehat{\bm{x}})\mathrel{\mathop{\mathchar 58\relax}}=\mathbb{B}_{i}({\bm{\theta}},\mathbf{X}(t,\widehat{\bm{x}}))\psi(t)^{d/2}$ and $\widehat{\mathrm{f}}_{i}(t,\widehat{\bm{x}})=\mathrm{f}_{i}(t,\mathbf{X}(t,\widehat{\bm{x}}))\psi(t)^{d/2}$, $i=1,\dots,d$.

Next, we introduce a time discretization of problem (3.1) upon using $\mathcal{J}$ as defined in (4.1). Let us fix $N\in\mathbb{N}$ and let $\tau\mathrel{\mathop{\mathchar 58\relax}}=T/N$ be the time step. Now, for $n=1,\ldots,N$, we define $t^{n}\mathrel{\mathop{\mathchar 58\relax}}=n\tau$, $\widehat{\mathbb{B}}^{n}_{i}(\cdot)=\widehat{\mathbb{B}}_{i}(\cdot,\mathbf{X}(t^{n},\widehat{\bm{x}}))$ and $\widehat{\mathrm{f}}^{n}_{i}$ to be

[TABLE]

which in turn allows us to incorporate a general $\mathbf{f}$. Then we consider a time discrete version of (3.1): given the initial condition $({\bm{\alpha}}_{0},{\bm{\theta}}_{0})=\mathrel{\mathop{\mathchar 58\relax}}({\bm{\alpha}}_{\tau}(0),{\bm{\theta}}_{\tau}(0))$, find $(\bar{{\bm{\alpha}}}_{\tau},\bar{{\bm{\theta}}}_{\tau})\in\mathcal{H}^{\tau}_{ad}\times\mathcal{V}^{\tau}_{ad}$ solving

[TABLE]

where

[TABLE]

and

[TABLE]

and $\mathcal{H}^{\tau}_{ad}=\mathbb{P}^{1}_{\tau}\cap\mathcal{H}_{ad}$, $\mathcal{V}^{\tau}_{ad}=\mathbb{P}^{1}_{\tau}\cap\mathcal{V}_{ad}$ with $\mathbb{P}^{1}_{\tau}\mathrel{\mathop{\mathchar 58\relax}}=\{{\bm{v}}\in C[0,T]\mathrel{\mathop{\mathchar 58\relax}}{\bm{v}}|_{[t^{n-1},t^{n}]}\in\mathbb{P}^{1},\ n=1,\dots,N\}$. Hereafter $\mathbb{P}^{1}$ is the space of polynomials of degree at most 1. Moreover, by applying the same arguments of Theorem 1, it follows that there exists $(\bar{{\bm{\alpha}}}_{\tau},\bar{{\bm{\theta}}}_{\tau})\in\mathcal{H}^{\tau}_{ad}\times\mathcal{V}^{\tau}_{ad}$ a solution to problem (4.4).

Next we prove the convergence of the discrete problem to a minimizer of the continuous problem. Such a proof is motivated by $\Gamma$-convergence theory [9, 4].

Theorem 6** (Problem 1: convergence to a minimizer).**

The family of minimizers $\left\{(\bar{{\bm{\alpha}}}_{\tau},\bar{{\bm{\theta}}}_{\tau})\right\}_{\tau>0}$ to (4.4) is uniformly bounded in $[\mathrm{H}^{1}(0,T)]^{2}$ and it contains a subsequence that converges weakly to $(\bar{{\bm{\alpha}}},\bar{{\bm{\theta}}})$ in $[\mathrm{H}^{1}(0,T)]^{2}$, a solution to the minimization problem (3.1), and $\lim_{\tau\to 0}\mathcal{J}_{\tau}(\bar{{\bm{\alpha}}}_{\tau},\bar{{\bm{\theta}}}_{\tau})=\mathcal{J}(\bar{{\bm{\alpha}}},\bar{{\bm{\theta}}})$.

Proof.

We proceed is several steps.

1. 1.-

Boundedness of $\left\{(\bar{{\bm{\alpha}}}_{\tau},\bar{{\bm{\theta}}}_{\tau})\right\}_{\tau>0}$ in $[\mathrm{H}^{1}(0,T)]^{2}$: This follows immediately from the fact that $(\bar{{\bm{\alpha}}}_{\tau},\bar{{\bm{\theta}}}_{\tau})$ minimizes $\mathcal{J}_{\tau}$ and $\lambda,\eta>0$: given that the constant function $({\bm{\alpha}}_{0}(t),{\bm{\theta}}_{0}(t))=({\bm{\alpha}}_{0},{\bm{\theta}}_{0})$ belongs to $\mathcal{H}^{\tau}_{ad}$, we have

[TABLE]

This implies the existence of a (not relabeled) weakly convergent subsequence such that $(\bar{{\bm{\alpha}}}_{\tau},\bar{{\bm{\theta}}}_{\tau})\rightharpoonup(\bar{{\bm{\alpha}}},\bar{{\bm{\theta}}})$ in $[\mathrm{H}^{1}(0,T)]^{2}$ and $(\bar{{\bm{\alpha}}},\bar{{\bm{\theta}}})\in\mathcal{H}_{ad}\times\mathcal{V}_{ad}$. It remains to prove that $(\bar{{\bm{\alpha}}},\bar{{\bm{\theta}}})$ solves (3.1) and $\lim_{\tau\to 0}\mathcal{J}_{\tau}(\bar{{\bm{\alpha}}}_{\tau},\bar{{\bm{\theta}}}_{\tau})=\mathcal{J}(\bar{{\bm{\alpha}}},\bar{{\bm{\theta}}})$. 2. 2.-

Lower bound inequality: We show that

[TABLE]

for all $\left\{({\bm{\alpha}}_{\tau},{\bm{\theta}}_{\tau})\right\}_{\tau>0}\subset\mathcal{H}^{\tau}_{ad}\times\mathcal{V}^{\tau}_{ad}$ converging to $({\bm{\alpha}},{\bm{\theta}})$ weakly in $[\mathrm{H}^{1}(0,T)]^{2}$. Consequently $({\bm{\alpha}}_{\tau},{\bm{\theta}}_{\tau})\to({\bm{\alpha}},{\bm{\theta}})$ strongly in $[\mathrm{L}^{2}(0,T)]^{2}$ for a subsequence (not relabeled). Let $\overline{\Pi}_{\tau}$ be the piecewise constant interpolation operator at the nodes $\left\{t^{n}\right\}_{n>0}$. Then by straightforward computations it follows that

[TABLE]

In view of the smoothness of $\widehat{{\bm{d}}}_{i}$, $i=1,\ldots,n_{p}$, then $\widehat{\mathbb{B}}_{j}(\overline{\Pi}_{\tau}{\bm{\theta}}_{\tau})\rightarrow\widehat{\mathbb{B}}_{j}({\bm{\theta}})$ in $\mathrm{L}^{2}(0,T;\mathrm{L}^{2}(\widehat{D}))$, $j=1,\dots,d$. Collecting these results and making use of the bounds $\|{\bm{\alpha}}_{\tau}\|_{\mathrm{L}^{\infty}(0,T)}\leq{\bm{\alpha}}^{*}$ and $\|{\bm{\theta}}_{\tau}\|_{\mathrm{L}^{\infty}(0,T)}\leq{\bm{\theta}}^{*}$ we obtain (cf (3.7))

[TABLE]

which in conjunction with the regularity of $\mathbf{f}$ leads to

[TABLE]

On the other hand, given that $d_{t}{\bm{\alpha}}_{\tau}$ converges weakly to $d_{t}{\bm{\alpha}}$ in $\mathrm{L}^{2}(0,T)$, from the weak lower semi-continuity of the semi-norm it follows that

[TABLE]

From (4.6) and (4.7) we conclude

[TABLE] 3. 3.-

Existence of a recovery sequence: Let $({\bm{\alpha}},{\bm{\theta}})\in\mathcal{H}_{ad}\times\mathcal{V}_{ad}$ be given. Then, the piecewise linear Lagrange interpolant $(\Pi_{\tau}{\bm{\alpha}},\Pi_{\tau}{\bm{\theta}})$ of $({\bm{\alpha}},{\bm{\theta}}$) belongs to $(\mathcal{H}_{ad}^{\tau},\mathcal{V}_{ad}^{\tau})$. Since $d_{t}(\Pi_{\tau}{\bm{\alpha}})\rightarrow d_{t}{\bm{\alpha}}$, $d_{t}(\Pi_{\tau}{\bm{\theta}})\rightarrow d_{t}{\bm{\theta}}$ and $\Pi_{\tau}{\bm{\alpha}}^{\top}_{\tau}\widehat{\mathbb{B}}_{i}(\Pi_{\tau}{\bm{\theta}}_{\tau})\Pi_{\tau}{\bm{\alpha}}_{\tau}\to{\bm{\alpha}}^{\top}\widehat{\mathbb{B}}_{i}({\bm{\theta}}){\bm{\alpha}}$ in $\mathrm{L}^{2}(0,T;\mathrm{L}^{2}(\widehat{D}))$ because of the box constrains, we obtain

[TABLE] 4. 4.-

$(\bar{{\bm{\alpha}}},\bar{{\bm{\theta}}})$* is a minimizer for problem (3.1) :* We need to show

[TABLE]

From step 3 there exists $\{({\bm{v}}_{\tau},{\bm{w}}_{\tau})\}_{\tau>0}$ such that ${\bm{v}}_{\tau}\rightharpoonup{\bm{v}}$ and ${\bm{w}}_{\tau}\rightharpoonup{\bm{w}}$ in $\mathrm{H}^{1}(0,T)$ and

[TABLE]

where we have used that $\bar{{\bm{\alpha}}}_{\tau}$ is a minimizer for $\mathcal{J}_{\tau}$ together with (4.5). 5. 5.-

Convergence: Since $(\bar{{\bm{\alpha}}}_{\tau},\bar{{\bm{\theta}}}_{\tau})$ is a minimizer we deduce the inequality $\mathcal{J}_{\tau}(\bar{{\bm{\alpha}}}_{\tau},\bar{{\bm{\theta}}}_{\tau})\leq\mathcal{J}_{\tau}(\Pi_{\tau}\bar{{\bm{\alpha}}},\Pi_{\tau}\bar{{\bm{\theta}}})$, whence applying first step 3 and next step 2 we see that

[TABLE]

This implies $\lim_{\tau\to 0}\mathcal{J}_{\tau}(\bar{{\bm{\alpha}}}_{\tau},\bar{{\bm{\theta}}}_{\tau})=\mathcal{J}(\bar{{\bm{\alpha}}},\bar{{\bm{\theta}}})$. In addition $\left\{\bar{{\bm{\alpha}}}_{\tau},\bar{{\bm{\theta}}}_{\tau}\right\}_{\tau>0}$ converges $\mathrm{L}^{2}$-strongly and $\mathrm{H}^{1}$-weakly to $(\bar{{\bm{\alpha}}},\bar{{\bm{\theta}}})$, a minimizer of (3.1).

This concludes the proof. ∎

4.2 Problem 2: Controlling intensities and positions

Like in the previous section, we consider the fixed domain $\widehat{D}$ and rewrite $\mathcal{F}$ (cf. (3.3b)) as

[TABLE]

with $\widehat{\mathbb{G}}_{i}({\bm{\phi}},\widehat{\bm{x}})\mathrel{\mathop{\mathchar 58\relax}}=\mathbb{G}_{i}({\bm{\phi}},\mathbf{X}(t,\widehat{\bm{x}}))\psi(t)^{d/2}$, $i=1,\dots,d$.

Next, we consider the previous definition $\mathcal{F}$ in order to introduce a time discretization of problem (3.3). Let us fix $N\in\mathbb{N}$ and let $\tau\mathrel{\mathop{\mathchar 58\relax}}=T/N$ be the time step. Then we consider the following time discrete version of (3.3): given the initial condition $({\bm{\alpha}}_{0},{\bm{\phi}}_{0})=\mathrel{\mathop{\mathchar 58\relax}}(\bar{{\bm{\alpha}}}_{\tau}(0),\bar{{\bm{\phi}}}_{\tau}(0))$, find $(\bar{{\bm{\alpha}}}_{\tau},\bar{{\bm{\phi}}}_{\tau})\in\mathcal{H}^{\tau}_{ad}\times\mathcal{U}^{\tau}_{ad}$ solving

[TABLE]

where $\mathcal{U}^{\tau}_{ad}=\mathbb{P}^{1}_{\tau}\cap\mathcal{U}_{ad}$,

[TABLE]

and

[TABLE]

Moreover, it is straightforward to prove the following convergence result.

Theorem 7** (Problem 2: convergence to minimizers).**

The family of minimizers $\left\{(\bar{{\bm{\alpha}}}_{\tau},\bar{{\bm{\phi}}}_{\tau})\right\}_{\tau>0}$ to (4.10) is uniformly bounded in $[\mathrm{H}^{1}(0,T)]^{2}$ and it contains a subsequence that converges weakly to $(\bar{{\bm{\alpha}}},\bar{{\bm{\phi}}})$ in $[\mathrm{H}^{1}(0,T)]^{2}$, a solution to the minimization problem (3.3), and $\lim_{\tau\to 0}\mathcal{F}_{\tau}(\bar{{\bm{\alpha}}}_{\tau},\bar{{\bm{\phi}}}_{\tau})=\mathcal{F}(\bar{{\bm{\alpha}}},\bar{{\bm{\phi}}})$.

5 Numerical example

In this section we illustrate the performance of the proposed discrete schemes (4.4) and (4.10). We consider the approximation of two constant vector fields $\mathbf{f}$ on $D_{t}$. For all examples $\Omega\subset\mathbb{R}^{2}$ is a ball of unit radius centered at $(0,0)$ and the dipoles belong to a ball of radius 1.2 centered at $(0,0)$.

To solve the minimization problem, we use projected BFGS with Armijo line search [19]; alternative strategies such as semi-smooth Newton [16] can be immediately applied as well. The algorithm terminates when the $l^{2}$-norm of the projected gradient is less or equal to $10^{-6}$. In addition to the initial condition $({\bm{\alpha}}_{0},{\bm{\theta}}_{0})\in\mathbb{R}^{2n_{p}}$ for problem (4.4) or $({\bm{\alpha}}_{0},{\bm{\phi}}_{0})\in\mathbb{R}^{2n_{p}}$ for (4.10), the optimization algorithm requires an initial guess $(\bar{{\bm{\alpha}}},\bar{{\bm{\theta}}})$ and $(\bar{{\bm{\alpha}}},\bar{{\bm{\phi}}})\in\mathbb{R}^{2Nn_{p}}$, for (4.4) and (4.10), respectively. Following [2] we propose below Algorithm 1 in order to compute an initial guess $({\bm{\alpha}}_{int},{\bm{\theta}}_{int})\in\mathcal{H}^{\tau}_{ad}\times\mathcal{V}^{\tau}_{ad}$ for problem (4.4). A similar strategy can be used for problem (4.10).

Given ${\bm{a}}_{*},{\bm{a}}^{*},{\bm{b}}_{*},{\bm{b}}^{*}\in\mathbb{R}^{2Nn_{p}}$ , by $\mbox{Proj}_{[{\bm{a}}_{*},{\bm{a}}^{*},{\bm{b}}_{*},{\bm{b}}^{*}]}$ we denote pointwise projection on the interval $[{\bm{a}}_{*},{\bm{a}}^{*}]\times[{\bm{b}}_{*},{\bm{b}}^{*}]$

[TABLE]

where $\min$ and $\max$ are interpreted componentwise. Notice that the above algorithm, related to model predictive control (MPC), is the so-called finite horizon model predictive or instantaneous optimization algorithm [8, 7, 17, 1]: the minimization takes place over one time step only. This cheaper strategy provide us an “accurate” initial guess to solve (4.4).

5.1 Approximation of a vector field

In this section we focus on the approximation of two vectors fields $\mathbf{f}$ defined on moving domains $D_{t}$. With this aim we solve the discrete minimization problems (4.4) and (4.10). We will achieve our goal by creating a magnetic force which allows us to “magnetically inject” nanoparticles inside the domain. We first consider the approximation of a constant vector field $\mathbf{f}_{1}(\bm{x},t)=(\frac{1}{\sqrt{2}},-\frac{1}{\sqrt{2}})^{\top}$ by controlling the directions and intensities of the $4$ dipoles. For the second example we are interested in the approximation of the field given by $\mathbf{f}_{2}(\bm{x},t)=(1,0)^{\top}$ by controlling the positions and intensities of $3$ dipoles. The moving domains $D_{i,t}$ are such that $D_{i,t}=\bm{x}(t,\widehat{D})$, with $\bm{x}(t,\widehat{\bm{x}})=\bm{\varphi}_{i}(t)+\widehat{\bm{x}}$, $i=1,2$. Here $\bm{\varphi}_{1}(t)=\frac{0.6}{T}(t,-t)^{\top}$ and $\bm{\varphi}_{2}(t)=(t,0)^{\top}$. The trajectories of the barycenters of $D_{1,t}$ and $D_{2,t}$ are shown in Figure 2 (left and right). For both examples we consider final time $T=0.75$. In the first case the reference domain $\widehat{D}$ is a ball of radius $0.2$ centered at $\bm{x}_{I}=(-0.6,0.6)$ whereas for the second one, $\widehat{D}$ is a ball of radius $0.2$ centered at $\bm{x}_{I}=(-0.75,0)$. In both cases, $D_{i,t}$ moves from $\bm{x}_{I}$ to $\bm{x}_{F}=(0,0)$.

For each of these configurations we have solved the discrete problems with $N=100$ time intervals and $\lambda=\eta=10^{-5}$. In order to solve (4.4), we consider an admissible set characterized by initial conditions ${\bm{\alpha}}_{0}=(2,0,0,2),{\bm{\theta}}_{0}=(0,\pi/2,3\pi/2,3\pi/2)\in\mathbb{R}^{4}$, and by upper and lower bounds for the intensities and angles given by: ${\bm{\alpha}}^{*}=(2,\ldots,2),{\bm{\alpha}}_{*}=(-2,\ldots,-2),{\bm{\theta}}^{*}=(2\pi,\ldots,2\pi),{\bm{\theta}}_{*}=(0,\ldots,0)\in\mathbb{R}^{4}$, respectively. The initial guesses $({\bm{\alpha}}_{int},{\bm{\theta}}_{int})$ are computed by Algorithm 1 with $\verb"tol"=10^{-3}$ with a total of 525 iterations. Notice that, at each time step $n$, the iterations of the minimization problem in Algorithm 1 depends on $n_{p}$ unknowns. We recall that, due to the non-convexity of the cost functional $\mathcal{J}$ (and $\mathcal{J}_{\tau}$), we may converge to different local minima depending on the choice of the initial guess.

The solution $\bar{{\bm{\alpha}}}_{\tau}(t)=(\bar{\alpha}_{i,\tau}(t))_{i=1}^{4},\bar{{\bm{\theta}}}_{\tau}(t)=(\bar{\theta}_{i,\tau}(t))_{i=1}^{4}$ of problem (4.4) with initial guesses given by $({\bm{\alpha}}_{int},{\bm{\theta}}_{int})$ is depicted in Figure 3.

From Figure 3 (left) we observe that the optimal intensity ${\bm{\alpha}}_{i}$ due to $i^{th}$ dipole, is smaller when $D_{t}$ is close to the $i^{th}$ dipole and $\mathbf{f}$ is pointing in a direction opposite to the dipole position. This is due to the fact that the magnetic forces on the boundary of the domain $\Omega$ are much higher than inside $\Omega$, thus making it difficult to approximate the magnetic force when $D_{t}$ is close to the boundary. In particular, it can be seen that the intensities of dipole 2 and 3 are small when $t$ is close to 0. Such a behavior is expected because $D_{1,t}$ is close to the boundary of $\Omega$, where the magnetic field generated by these dipoles is large, thus it is difficult for dipoles 2 and 3 to “push” in the $\mathbf{f}_{1}$ direction. We also notice that dipoles 1 and 4, which can create an attractive field in the $\mathbf{f}_{1}$ direction, have the largest intensities at initial times and decrease when the time is close to $0.75$. Figure 3 (right) shows the dynamics of the angles $(\bar{\theta}_{i,\tau}(t))_{i=1}^{4}$.

Figure 4 shows the approximate field and dipole configuration for three time instances $t=0.007,0.375,0.75$. Here, the magnitude of magnetic force $|\nabla|{\bm{h}}|^{2}|$ in logarithmic scale (in the background) and the magnetic force represented by arrows are depicted. The bottom figures show the dipoles direction represented by the arrows outside the domain $\Omega$. The top figures illustrate the normalized magnetic force in $\Omega$, but in the bottom figures the force is restricted to $D_{t}$. It can be seen that the magnetic force is close (i.e., almost constant) to $\mathbf{f}_{1}$ in $D_{1,t}$ as expected whereas it is quite far from constant in the entire domain.

Next, we approximate the vector field $\mathbf{f}_{2}$ defined on $D_{2,t}$. Here, given that we want to create a uniform field $\mathbf{f}_{2}(\bm{x},t)=(1,0)^{\top}$ on a domain which moves from left to the center of $\Omega$ (see Figure 2 (right)), we restrict the position of the dipoles to the right of the domain. Namely, we assume that each dipole $i$ moves along a prescribed trajectory $\mathcal{C}_{i}=1.2[\cos(\phi_{i}),\sin(\phi_{i})],i=1,2,3$, where $\phi_{1}(t)\in(-\pi/90,\pi/90)$, $\phi_{2}(t)\in(\pi/90,3\pi/4)$ and $\phi_{3}(t)\in(-3\pi/4,-\pi/90)$. We solve (4.10) for an admissible set characterized by initial conditions ${\bm{\alpha}}_{0}=(-2,0,0),{\bm{\phi}}_{0}=(0,2\pi/3,4\pi/3)\in\mathbb{R}^{3}$, and by upper and lower bounds for the intensities and angles given by: ${\bm{\alpha}}^{*}=(2,2,2),{\bm{\alpha}}_{*}=(-2,-2,-2),{\bm{\phi}}^{*}=(\pi/90,3\pi/4,2\pi),{\bm{\theta}}_{*}=(-\pi/90,\pi/90,5\pi/4)\in\mathbb{R}^{3}$, respectively.

The vector intensities $(\bar{\alpha}_{i,\tau})_{i=1}^{3}$ and the parametrization of the position $(\bar{\phi}_{i,\tau})_{i=1}^{3}$ solutions to the minimization problem (4.10) are shown in Figure 5 left and right, respectively . We notice that the first dipole attains a fixed position and intensity (the maximum intensity), whereas dipoles 2 and 3 are moving but their intensity is almost constant. Given that we want to approximate a uniform and unitary vector field on $D_{1,t}$, the intensity of dipole 1 decreases as the domain approaches dipole 1. This is due to the fact that the force applied to the edge of $D_{1,t}$, close to dipole 1, is larger than the one applied to the farther edge, and this difference increases as $D_{1,t}$ is close to the source.

The magnetic force and dipole positions are shown for three time instances $t=0.007,0.375$ and $0.75$ in Figures 6. Here dipoles 2 and 3 “follows” the trajectory of $D_{2,t}$ which allows to obtain an almost constant field on the moving domain. From the previous figures it can be seen that our proposed minimization procedure produces a vector field that is close to the target field.

6 Numerical solution of the convection-dominated diffusion equation

From the previous section it follows that an almost uniform magnetic force field can be attained in a small region for two different dipole configurations. In this section we move a concentration of ferrofluid modeled by (1.1)-(1.3) where the advection term is computed either by solving the minimization problem (3.1) or (3.3) depending on the dipole configuration. The aim of these examples is to study the feasibility of the dipole configurations to create a magnetic force that allow us to “push” the concentration inside the domain.

For simplicity, we assume $\gamma=\nu=1$, ${\bm{u}}={\bm{0}}$. To further minimize the spreading we choose a small diffusion coefficient, in particular, we set $A=\varepsilon\mathbb{I}$, with $\varepsilon\ll 1$. As a result the advection diffusion equation (1.1) and (1.2) can be written as

[TABLE]

Multiplying (6.1a) by a test function $v\in\mathrm{H}^{1}(\widetilde{\Omega})$, integrating by parts, and using the boundary condition in (6.1b) leads us to the following problem: Given $c_{0}\in\mathrm{L}^{2}(\widetilde{\Omega})$, find $c\in\mathrm{L}^{2}(0,T;\mathrm{H}^{1}(\widetilde{\Omega}))\cap\mathrm{H}^{1}(0,T;\mathrm{H}^{1}(\widetilde{\Omega})^{\prime})$ such that

[TABLE]

here $\langle\cdot,\cdot\rangle$ is the duality product between $\mathrm{H}^{1}(\widetilde{\Omega})$ and $\mathrm{H}^{1}(\widetilde{\Omega})^{\prime}$. The existence and uniqueness of solution to this problems follows, for instance, from [11, Theorem 6.6]. Notice that in case of the Dirichlet boundary condition $c=0$ we must replace $\mathrm{H}^{1}(\widetilde{\Omega})$ by $\mathrm{H}^{1}_{0}(\widetilde{\Omega})$ in the above calculations.

Below we present two numerical examples. In the first example we consider a simplified $\widetilde{\Omega}$ whose size is comparable to $\Omega$ where the goal is to magnetically inject the concentration of nanoparticles. With the current configuration, the concentration may reach the boundary and, as is mentioned in [25] for a similar problem, boundary layers may appear. On the other hand, given that we are interested in the case $\varepsilon\ll 1$, a proper space discretization has to be considered in order to avoid numerical oscillations. We adopt a monotone scheme for space discretization: edge-averaged finite element (EAFE) [30]. Other alternatives such as SUPG stabilization are equally valid. Given that the EAFE scheme is a type of upwinding scheme, we observe a significant diffusion in all cases. Thus making it difficult to control the concentration in case of complicated geometries, for instance the flow around an obstacle.

As a second example we study the feasibility of moving the concentration around an obstacle. In order to ensure the uniformity of the field we consider a small, but complex, domain $\widetilde{\Omega}$ (see Figure 7 (right)). Moreover, it is shown that Dirichlet boundary condition $c=0$ on $\partial\widetilde{\Omega}$ in (6.1b) is meaningful in this example. For the numerical approximation we consider piecewise linear finite element space discretization and explicit Euler scheme for time discretization. In order to avoid having to solve a linear system at each time step we consider mass lumping with a correction technique to account for the dispersive effects (see [15]). In section 6.1 we discuss the EAFE scheme and in section 6.2 we discuss the explicit scheme. Numerical examples are presented in each case.

6.1 Edge-averaged Finite Element (EAFE) Method and Magnetic Injection

It is well–known that the standard finite element method applied to problem (6.2), yields solution oscillations when $\varepsilon\ll\|\nabla|{\bm{h}}|^{2}\|_{L^{\infty}(0,T;\widetilde{\Omega})}$. However, there is no universal approach to treat such problems. In order to choose a numerical scheme to approximate problem (6.2) first notice an interesting feature of the flux $J(c)$ of (6.1a):

[TABLE]

i.e., it is symmetrizable. Here

[TABLE]

Equations with the above property have been studied by many authors (see, for instance, [6, 24]). They also appear in applications such as semiconductors device simulation. In order to exploit the structure of the equation, for the numerical approximation we consider the edge-average finite element (EAFE), a monotone scheme introduced in [30] for stationary advection-diffusion equations. It was later extended to time-dependent advection-diffusion equations in [3] and compared with SUPG, another well-known stabilization technique. With this in mind we adapt the scheme presented in [3] to our particular case.

We begin by recalling the EAFE scheme of [30] and apply it to the symmetric operator in (6.3). We denote by $\mathcal{T}_{h}$ a shape regular triangulation of $\widetilde{\Omega}$ and by $\mathrm{H}_{h}\subset\mathrm{H}^{1}(\Omega)$ the space of piecewise linear finite elements with hat basis $\{\phi_{i}\}_{i=1}^{N}$. Notice that the scheme proposed below is equally applicable to both the Neumann and the Dirichlet problems. We designate by $E=E_{ij}$ the edge connecting the nodes $x_{i}$ and $x_{j}$. The discrete weak form of $-\textrm{div }(a\nabla u)$ reads

[TABLE]

where $\delta_{E}(w_{h})=w_{i}-w_{j}$ and

[TABLE]

because $\nabla\phi_{i}$ is piecewise constant. We now replace the first factor by the so-called harmonic average over the edge $E$, namely

[TABLE]

Upon setting $\omega_{E}^{T}\mathrel{\mathop{\mathchar 58\relax}}=\int_{T}\nabla\phi_{i}\cdot\nabla\phi_{j}$, we realize that (6.5) can be approximated by

[TABLE]

which is the bilinear form corresponding to the EAFE scheme of [30].

For the time interval $[0,T]$ we introduce a uniform partition $0=t_{0}<t_{1}<\ldots<t_{N}=T$, with $N\geq 1$, and denote by $\Delta t=T/N$ the time step size. By applying the implicit Euler scheme for time discretization, we obtain the following fully discrete scheme of problem $\eqref{eq:weak}$: Given $c_{0}\in\mathrm{L}^{2}(\widetilde{\Omega})$, for $n=1,\ldots,N$ find $c^{n}_{h}\in\mathrm{H}_{h}$ such that

[TABLE]

Error estimates of a similar discrete problem can be found in [3].

In the following subsections we will study the feasibility of the manipulation of nanoparticles concentration by controlling the dipole intensity and direction. To fix ideas, let $\Omega$ be a ball of unit radius centered at $(0,0)$. The configuration of the 4 dipoles is same as in the previous example (cf. Section 5.1).

In our first example we focus on magnetically injecting the initial concentration $c_{0}$ of nanoparticles so that it arrives at a desired location. Instead of just the intensity (cf. [2, 20]) the control variables here are both the intensities and the directions of the dipoles. Additionally, the concentration is close to the boundary for certain time instances. A standard Galerkin scheme in space usually leads to oscillations in this case. We assume that the domain $\widetilde{\Omega}\subset\Omega$ is a $\pi/4$ clockwise rotation of the rectangle $[-0.9,0.9]\times[-0.3,0.3]$ (cf. Figure 7, left). We set the final time $T=0.75$ and $c_{0}(x,y)=\exp\left(3\times 10^{-5}((x+0.53)^{2}+(y-0.53)^{2})\right)$ (cf. Figures 7 and 8 (left)).

We use the optimal magnetic force computed in the first example of Section 5.1 (cf.  Figure 4) as an input to (6.7). Here we have set $\varepsilon=10^{-5}$, $\Delta t=7.5\times 10^{-3}$ and $h=0.0065$. The initial concentration is located at $(-0.53,0.53)$ and the final location is $(0,0)$, see Figure 7 (left). Figure 8 shows the snapshots of concentration at three time instances.

Notice that (cf. Figure 8) the concentration reaches the desired location. Since we are using 4 dipoles instead of 8 dipoles used in (cf. [2, 20]) we find it difficult to control the diffusion, also notice that we have Neumann boundary conditions here. Indeed the spikes in Figure 8 (right) can be explained by the behavior of the magnetic force given in Figures 3 (left) and 4 (right).

6.2 Explicit Scheme and Magnetic Injection Around an Obstacle

Our next example explores the possibility of moving the concentration around an obstacle. Indeed a larger $\varepsilon$ leads to significant diffusion and it is not possible to accomplish our goal. On the other hand, a smaller $\varepsilon$ can computationally lead to an ill-posed problem; recall the Neumann boundary condition:

[TABLE]

So when $\varepsilon$ is small then an appropriate boundary condition is Dirichlet (assuming that $\nabla|{\bm{h}}|^{2}\cdot\bm{n}\neq 0$).

Motivated by this, we consider the advection-diffusion equation (6.1a) but with the Dirichlet boundary condition $c=0$ on $\partial\widetilde{\Omega}\times(0,T)$. The latter is also physically meaningful if the concentration stays away from the boundary at all times. Even for small parameter $\varepsilon$ we observe a relatively large diffusion both in case of EAFE scheme and stabilized implicit (in time) scheme [2]. For that reason and, in order to ensure that the concentration is far from the boundary, we have adopted a piecewise linear finite element scheme for the space approximation and the explicit Euler scheme for the time discretization.

The fully discrete problem is: Given $c_{0}\in\mathrm{L}^{2}(\widetilde{\Omega})$, for $n=1,\ldots,N$ find $c^{n}_{h}\in\mathrm{H}_{h,0}\subset\mathrm{H}_{0}^{1}(\widetilde{\Omega})$ such that

[TABLE]

with $a(c^{n-1}_{h},v_{h})=\int_{\Omega}\varepsilon\nabla c_{h}^{n-1}\cdot\nabla v_{h}-c_{h}^{n-1}\nabla|{\bm{h}}|^{2}\cdot\nabla v_{h}\;dx$.

In the matrix vector notation the system in (6.8) is given by

[TABLE]

where ${\bm{M}}=(m_{ij})_{i,j}$ with $m_{ij}=\int_{\Omega}\phi_{i}\phi_{j}$ and ${\bm{A}}$ are the mass and the stiffness matrices, respectively. Moreover, for every $k\in\{1,\dots,N\}$ the vector ${\bm{c}}^{k}$ contains the finite element nodal values. Clearly (6.9) requires solving a linear system at every time. To avoid such a linear solve we employ the mass lumping (with correction) of [15]. Let $\overline{\bm{M}}=(\overline{m}_{ij})_{i,j}$ denotes a diagonal lumped mass matrix with diagonal entries given by $\overline{m}_{ii}=\sum_{j}m_{ij}$. Then we approximate ${\bm{M}}^{-1}$ by replacing the inverse of the lumped mass matrix by $(\mathbb{I}+{\bm{B}}_{r})\overline{\bm{M}}^{-1}$ where ${\bm{B}}_{r}=\overline{\bm{M}}^{-1}(\overline{\bm{M}}-{\bm{M}})$.

Our domain of interest is $\widetilde{\Omega}=[-0.18,0.18]^{2}\setminus[-0.02,0.02]\times[-0.18,0]$ depicted on Figure 7 (right) where the path traversed by the concentration is given by the dotted line. Our aim is to move the initial concentration

[TABLE]

from the initial location $(0.1,-0.1)$ at $t=0$ to the final location $(-0.1,-0.1)$ at $t=T=0.6$. With this in mind, we first solve the optimal control problem (4.4) with a ball of radius $0.2$ centered at $(0.1,-0.1)$ as a reference domain $\widehat{D}$. We approximate a time dependent vector field given by $\mathbf{f}(\bm{x},t)$ such that: $\mathbf{f}(\bm{x},t)=(0,1)^{\top}$ if $0\leq t<0.2$, $\mathbf{f}(\bm{x},t)=(-1,0)^{\top}$ if $0.2\leq t<0.4$ and $\mathbf{f}(\bm{x},t)=(0,-1)^{\top}$ if $0.4\leq t\leq 0.6$. The moving domains $D_{t}$ is such that $D_{t}=\bm{x}_{I}(t,\widehat{D})$, with $\bm{x}_{I}(t,\widehat{\bm{x}})=\bm{\varphi}(t)+\widehat{\bm{x}}$. Here $\bm{\varphi}(t)$, which represents the trajectory of the barycenter of $D_{t}$, is defined as: $\bm{\varphi}(t)=(0.1,-0.1)+t(0,1)^{\top}$ if $0\leq t<0.2$, $\bm{\varphi}(t)=(0.1,0.1)+(t-0.2)(-1,0)^{\top}$ if $0.2\leq t<0.4$ and $\bm{\varphi}(t)=(-0.1,0.1)+(t-0.4)(0,-1)^{\top}$ if $0.4\leq t\leq 0.6$. Figure 9 (top row) shows the resulting magnetic force at four different times.

In order to steer $c_{0}$ around the obstacle shown in Figure 9, we use the above optimal force as input to (6.8). The bottom row in Figure 9 shows the concentration snapshots at four different times. Here we have set $\varepsilon=10^{-8}$, $\Delta t=3\times 10^{-5}$ and mesh size $h=0.0016$.

The top row in Figure 9 shows the optimal magnetic force obtained by solving (4.4). The bottom row shows the evolution of the concentration at four different time instances. Notice that the magnetic forces prevent the concentration from reaching the boundary $\partial\widetilde{\Omega}$ and enables the concentration to reach the final location at $(-0.1,0.1)$ while avoiding the obstacle.