# adjoint method for reaction-diffusion problem

I'm trying to code a parameter estimation for a reaction-diffusion problem. Namely, knowing the distribution of tumor density $u$ at time $0$ and $T_f$ ($u^0$ and $u^f$), what are the best distributions of the $k$ and $D$ parameters explaining the observed tumor evolution, assuming

$$\partial{u}/\partial{t} = D(x)u + k(x)u(1-u). \;\; \mbox{(1)}$$

I have already coded the forward simulation code, giving the final distribution, knowing $D$, $k$ and $u^0$. I first tried to just set evenly distributed control points and simply used the matlab fmincon function to try to solve the problem by interpolating the $D$ and $k$ maps between the control points (thus limiting the number of free parameters to ~200, and implicitly enforcing smoothness of the $D$ and $k$ fields) but unsuccessfully (I still don't understand why but that's another question)

I'm slowly understanding that using the adjoint method would be a (the?) correct way to solve my problem, but while I'm reading some tutorial papers about this method I'm stucked on how to apply it to my problem.

Let say that my goal is to find the set of distribution $D$ and $k$ which minimize the function

$$G(u,D,k)= \sum_i{(\tilde{u}_i-u^f_i)^2} + \sum_i{D_i}^2 + \sum_i{k_i}^2 \;\; \mbox{(2)},$$

subject to an evolution function $$F(u^f,D,k) = 0.$$

The goal is to find the set of $D$ and $k$ minimizing $G$, which requires to evaluate $\partial{G}/\partial{D}=\partial_D{G}$ (and $\partial_k{G}$) multiple times per iteration if simply using iterative minimization algorithm on the function defined in (2).

If I understand well the concept of the adjoint method (as explained here), the idea is to replace the evaluation of $\partial_D{G}$ by observing that the optimal parameters are those for which $dG/dD$ and $dG/dk$ are zeros. First,

$$dG/dD = \partial_D{G} + \partial_{u^f}{G}.\partial_D{u^f}.$$

Since one has $F(u^f,D,k)=0$ everywhere, one also has the relation $$dF/dD = \partial_D{F} + \partial_{u^f}{F}.\partial_D{u^f} = 0$$

$$\Rightarrow \partial_D{u^f} = - (\partial_{u^f}{F})^{-1} . \partial_D{F}$$

$$\Rightarrow dG/dD = \partial_D{G} - \partial_{u^f}{G}. (\partial_{u^f}{F})^{-1} . \partial_D{F}$$

So by defining $\lambda$ by $$^T(\partial_{u^f}{F})\lambda = \partial_{u^f}{G}$$

one gets $$dG/dD = \partial_D{G} - \lambda . \partial_D{F}$$

the trick being that $\lambda$ does not depend on $D$ and has to be evaluated only once, and that $\partial_D{G}$ and $\partial_D{F}$ usually have an analytical expression, and do not require costly minimization. This should allow for fast evaluation of the derivative of $G$ and consequently faster minimization of $G$.

Now I'm trying to apply it to my specific problem of reaction diffusion, and the only expression I can find to fit with the above definition of $F$ would be

$$F(u^f,D,k) = u^0 + \int_0^{Tf}{((D\nabla u + k.u(1-u)) dt)}-u^f = 0$$

but then I don't understand how to determine $\partial_D{F}$ and $\partial_k{F}$ (I mean, not numerically but analytically).

Clearly, I missed some points... but I'm lost at this stage... If someone can give me a hand on this topic, I'd be very grateful!

Some remarks:

Notation: $[t_1,t_2]$: the time interval; $\Omega$: the spatial domain; $\bar{u}^i(x)$: the known tumour profile at $t_i$; $\left\| \cdot \right\|_\Omega$: a suitably chosen norm on $\Omega$. I assume that the design variables $D$ and $k$ do not depend on $t$.

1. The cost functional should be a scalar (real-valued), something like:

$${\cal G}(u; D, k) := \frac{1}{2} \sum_{i=1,2} \left\| u^i(x) - \bar{u}^i(x) \right\|^2_\Omega + \frac{1}{2} \left\| D(x) \right\|^2_{\Omega} \ + \frac{1}{2} \left\| k(x) \right\|^2_{\Omega}.$$

2. Your reaction-diffusion PDE problem (the "forward model", including appropriate boundary and initial conditions) acts as a constraint for the optimization problem:

$${\cal F}(u; D,k) = 0.$$

3. The Lagrangian functional for the constrained problem reads something like:

$${\cal L} (u; \lambda; D,k) := {\cal G}(u;D,k) + \left< \lambda , {\cal F}(u; D,k) \right>_{[t_1,t_2] \times \Omega}.$$

Here $\lambda$ is the adjoint variable defined on $[t_1,t_2] \times \Omega$ and $<\cdot>$ an inner product. It's up to you to define the latter.

4. The first-order necessary conditions for a local optimum (the so-called first-order KKT conditions, see chapter 12 of the book by Nocedal and Wright) give you the optimality system:

• the adjoint equation: ${\cal L}_u (u; \lambda; D,k)(\delta u) = 0 \;\; \forall \; \delta u$;
• the primal equation (the forward model): ${\cal L}_\lambda (u; \lambda; D,k)(\delta \lambda) = 0 \;\; \forall \; \delta \lambda$;
• the gradient (or optimal) equations: ${\cal L}_D (u; \lambda; D,k)(\delta D) = 0 \;\; \forall \; \delta D$ and ${\cal L}_k (u; \lambda; D,k)(\delta k) = 0 \;\; \forall \; \delta k$.
5. The choice of cost functional and the primal (forward) PDE have to be compatible such that the dual (adjoint) problem is well posed. Compatibility is discussed further here, here, and here.

6. The proper way to derive the analytic form of these equations (using integration by parts) is shown in chapters 2 and 3 of Max Gunzburger's optimal control book.

• Many thanks for your comments, I really appreciate ! It's clear that I have to work a little bit more on the fundamentals of the methods in order to apply it effectively, so I'll buy the book you advised to me. May be will return here later to ask precise questions on technical points if so required (I'm actually physicist and not mathematician... so I'm sometime get lost when it comes to real mathematical explanation, but I have to try :-) ) – david guez Dec 11 '17 at 15:25
• Gunzurger's book is great -- you might also get a lot out of Weinstock's Calculus of Variations. – Daniel Shapero Dec 11 '17 at 17:09