Optimal Control for Chemotaxis Systems and Adjoint-Based Optimization with Multiple-Relaxation-Time Lattice Boltzmann Models

TL;DR

This paper develops adjoint-based optimization methods for controlling complex chemotaxis systems using coupled lattice Boltzmann models, enabling efficient and accurate solutions for biochemical and biological applications.

Contribution

It introduces a novel adjoint MRT lattice Boltzmann model for optimal control of nonlinear chemotaxis systems, with detailed derivation and validation.

Findings

The MRT model accurately recovers the chemotaxis system via Chapman-Enskog analysis.

The adjoint AMRT model is simple, efficient, and highly parallelizable.

The combined approach effectively solves complex control problems in biochemical systems.

Abstract

This paper is devoted to continuous and discrete adjoint-based optimization approaches for optimal control problems governed by an important class of Nonlinear Coupled Anisotropic Convection-Diffusion Chemotaxis-type System (NCACDCS). This study is motivated by the fact that the considered complex systems (with complex geometries) appear in diverse biochemical, biological and biosocial criminology problems. To solve numerically the corresponding nonlinear optimization problems, the primal problem NCACDCS is discretised by a coupled Lattice Boltzmann Method with a general Multiple-Relaxation-Time collision operators (MRT) while for the adjoint problem, an Adjoint Multiple-Relaxation-Time lattice Boltzmann model (AMRT) is proposed and investigated. First, the optimal control problems are formulated and first-order necessary optimality conditions are established by using sensitivity and adjoint calculus. The resulting problems are discretised by the coupled MRT and AMRT models and solved via gradient descent methods. First of all, an efficient and stable modified MRT model for NCACDCS is developed, and through the Chapman-Enskog analysis we show that NCACDCS can be correctly recovered from the proposed MRT model. For the adjoint problem, the discretisation strategy is based on AMRT model, which is found to be as simple as MRT model with also highly-efficient parallel nature. The derivation of AMRT model and the discrete cost functional gradient are derived mathematically in detail using the developed MRT model. The obtained method is reliable, efficient, practical to implement and can be easily incorporated into any existing MRT code.

Equations

$$\begin{array}[]{lcr}\displaystyle\frac{\partial u}{\partial t}-div(D_{u}(\mathbf{x},t;{\bf u})\nabla u)+div(\mathfrak{T}_{u}(\mathbf{x},t;u))=\Phi_{u}(\mathbf{x},t,{\bf f}_{1};{\bf u})\\ \displaystyle\hskip 56.9055pt-div(D_{u,1}(\mathbf{x},t;{\bf u})\nabla v)-div(D_{u,2}(\mathbf{x},t;{\bf u})\nabla w)~{}~{}~{}on~{}{\cal Q},\\ \displaystyle\frac{\partial v}{\partial t}-div(D_{v}(\mathbf{x},t;{\bf u})\nabla v)+div(\mathfrak{T}_{v}(\mathbf{x},t;v))=\Phi_{v}(\mathbf{x},t,{\bf f}_{2};{\bf u})\\ \displaystyle\hskip 56.9055pt-div(D_{v,1}(\mathbf{x},t;{\bf u})\nabla u)-div(D_{v,2}(\mathbf{x},t;{\bf u})\nabla w)~{}~{}~{}on~{}{\cal Q},\\ \displaystyle\frac{\partial w}{\partial t}-div(D_{w}(\mathbf{x},t;{\bf u})\nabla w)+div(\mathfrak{T}_{w}(\mathbf{x},t;w))=\Phi_{w}(\mathbf{x},t,{\bf f}_{3};{\bf u})\\ \displaystyle\hskip 56.9055pt-div(D_{w,1}(\mathbf{x},t;{\bf u})\nabla u)-div(D_{w,2}(\mathbf{x},t;{\bf u})\nabla v)~{}~{}~{}on~{}{\cal Q},\\ \text{with the initial condition}\\ (u,v,w)(.,0)=(u_{0},v_{0},w_{0})~{}~{}~{}on~{}\Omega,\end{array}$$

$$\begin{array}[]{lcr}\displaystyle(D_{u}\nabla u).{\bf n}=(D_{u,1}\nabla v).{\bf n}+(D_{u,2}\nabla w).{\bf n}~{}~{}on~{}\Sigma,\\ \displaystyle(D_{v}\nabla v).{\bf n}=(D_{v,1}\nabla u).{\bf n}+(D_{v,2}\nabla w).{\bf n}~{}~{}on~{}\Sigma,\\ \displaystyle(D_{w}\nabla w).{\bf n}=(D_{w,1}\nabla u).{\bf n}+(D_{w,2}\nabla v).{\bf n}~{}~{}on~{}\Sigma,\end{array}$$

$$\mathcal{J}(\mathbf{f})={\int }_{\Omega }{J}_0(\mathbf{u}(T))dx+{\int }_0^T{\int }_{\Omega }\mathcal{Z}(\mathbf{u},\mathbf{f})dxdt,$$

$$\begin{array}[]{lcr}\text{Find ${\bf f}^{*}\in{\cal K}_{ad}$ such that}\\ \displaystyle{\cal J}({\bf f}^{*})=\inf_{{\bf f}\in{\cal K}_{ad}}{\cal J}({\bf f}).\end{array}$$

$$\begin{array}[]{lcr}\displaystyle{\cal J}({\bf f})=\frac{a_{1}}{2}\int_{\Omega}\!\!\mid u(T)-u_{f}\mid^{2}dx+\frac{a_{2}}{2}\int_{\Omega}\!\!\mid v(T)-v_{f}\mid^{2}dx+\frac{a_{3}}{2}\int_{\Omega}\!\!\mid w(T)-w_{f}\mid^{2}dx\\ \displaystyle+\frac{b_{1}}{2}\int_{0}^{T}\!\!\!\int_{\Omega}\!\!\mid u-u_{d}\mid^{2}dxdt+\frac{b_{2}}{2}\int_{0}^{T}\!\!\!\int_{\Omega}\!\!\mid v-v_{d}\mid^{2}dxdt+\frac{b_{3}}{2}\int_{0}^{T}\!\!\!\int_{\Omega}\!\!\mid w-w_{d}\mid^{2}dxdt\\ \displaystyle+\sum_{i=1}^{3}\sum_{j=1}^{d_{i}}\alpha_{i,j}\int_{0}^{T}\!\!\!\int_{\Omega}\!\!\mid f^{(j)}_{i}-f^{(j)}_{r,i}\mid^{p_{i,j}}dxdt\end{array}$$

$$\begin{array}[]{lcr}\displaystyle\frac{\partial u}{\partial t}-div(\sigma\nabla u)=-u+f^{(1)}_{1}uv+f^{(2)}_{1}~{}~{}~{}on~{}{\cal Q},\\ \displaystyle\frac{\partial v}{\partial t}-div(\nabla v)=-div(2\frac{v}{u}\nabla u)+f_{2}-uv~{}~{}~{}on~{}{\cal Q}.\\ \end{array}$$

$$\begin{array}[]{lcr}\displaystyle\frac{\partial u}{\partial t}-div(\sigma_{1}(u)\nabla u)+div(\vec{\omega}u)=-div(\chi_{1}(.;u,v)\nabla v)-div(\chi_{2}(.;u,w)\nabla w)+F(u)~{}~{}~{}on~{}{\cal Q},\\ \displaystyle\frac{\partial v}{\partial t}-div(\sigma_{2}\nabla v)+div(\vec{\omega}v)=f^{(1)}_{2}u-f^{(2)}_{2}v-f^{(3)}_{2}uv~{}~{}~{}on~{}{\cal Q},\\ \displaystyle\frac{\partial w}{\partial t}-div(\sigma_{3}\nabla w)+div(\vec{\omega}w)=f^{(1)}_{3}u-f^{(2)}_{3}w+f^{(3)}_{3}~{}~{}~{}on~{}{\cal Q},\end{array}$$

$$\begin{array}[]{lcr}\displaystyle\frac{\partial u}{\partial t}-div(\sigma_{1}(u)\nabla u)+div(\vec{\omega}u)=-div(\chi_{1}(.;u,w)\nabla w)+F_{1}(u,v)~{}~{}~{}on~{}{\cal Q},\\ \displaystyle\frac{\partial v}{\partial t}-div(\sigma_{2}(v)\nabla v)+div(\vec{\omega}v)=-div(\chi_{2}(.;v,w)\nabla w)+F_{2}(u,v)~{}~{}~{}on~{}{\cal Q},\\ \displaystyle\frac{\partial w}{\partial t}-div(\sigma_{3}\nabla w)+div(\vec{\omega}w)=F_{3}(u,v,w;{\bf f}_{3})~{}~{}~{}on~{}{\cal Q},\end{array}$$

$$\begin{array}[]{lcr}\displaystyle\frac{\partial u}{\partial t}-div(\sigma_{1}(u)\nabla u)+div(\vec{\omega}u)=\mu_{1}\frac{uw}{w+\theta}-\mu_{2}u+\mu_{3}u(1-m-u)\\ \displaystyle\hskip 85.35826pt-div(u\chi_{1}(.;m)\nabla m)-div(u\chi_{2}(.;v)\nabla v)~{}~{}~{}on~{}{\cal Q},\\ \displaystyle\frac{\partial v}{\partial t}-div(\sigma_{2}\nabla v)+div(\vec{\omega}v)=r_{1}u(1-u)-r_{2}v+f~{}~{}~{}on~{}{\cal Q},\\ \displaystyle\frac{\partial w}{\partial t}-div(\sigma_{3}\nabla w)+div(\vec{\omega}w)=\gamma_{1}m-\gamma_{2}w-\gamma_{3}\frac{uw}{w+\theta}+g~{}~{}~{}on~{}{\cal Q},\\ \displaystyle\displaystyle\frac{\partial m}{\partial t}=-\nu_{1}wm+\nu_{2}u(1-m-u)~{}~{}~{}on~{}{\cal Q},\end{array}$$

$${\mathcal{J}}^{\prime }({\mathbf{f}}^{\ast }).(\mathbf{f}-{\mathbf{f}}^{\ast })\ge 0,\forall \mathbf{f}\in {\mathcal{K}}_{ad},$$

$$\begin{array}[]{lcr}\displaystyle\frac{\partial p}{\partial t}-div(D_{u}(.;{\bf u})\nabla p)-div((\tilde{H}_{u;u}(.;{\bf u},\nabla{\bf u})-\frac{\partial\mathfrak{T}_{u}}{\partial u}(.;u))p)-\frac{\partial\Phi_{u}}{\partial u}(.;{\bf f}_{1},{\bf u})p\\ \displaystyle\hskip 42.67912pt=div(\tilde{H}_{u;v}(.;{\bf u},\nabla{\bf u})q)+div(\tilde{H}_{u;w}(.;{\bf u},\nabla{\bf u})\rho)+\frac{\partial\Phi_{u}}{\partial v}(.;{\bf f}_{1},{\bf u})q+\frac{\partial\Phi_{u}}{\partial w}(.;{\bf f}_{1},{\bf u})\rho\\ \displaystyle\hskip 56.9055pt+\frac{\partial\Phi_{u}}{\partial{\bf f}_{1}}(.;{\bf f}_{1},{\bf u}){\bf h}_{1}-div(D_{u,1}(.;{\bf u})\nabla q)-div(D_{u,2}(.;{\bf u})\nabla\rho),\\ \displaystyle\frac{\partial q}{\partial t}-div(D_{v}(.;{\bf u})\nabla q)-div((\tilde{H}_{v;v}(.;{\bf u},\nabla{\bf u})-\frac{\partial\mathfrak{T}_{v}}{\partial v}(.;v))q)-\frac{\partial\Phi_{v}}{\partial v}(.;{\bf f}_{2},{\bf u})q\\ \displaystyle\hskip 42.67912pt=div(\tilde{H}_{v;u}(.;{\bf u},\nabla{\bf u})p)+div(\tilde{H}_{v;w}(.;{\bf u},\nabla{\bf u})\rho)+\frac{\partial\Phi_{v}}{\partial u}(.;{\bf f}_{2},{\bf u})p+\frac{\partial\Phi_{v}}{\partial w}(.;{\bf f}_{2},{\bf u})\rho\\ \displaystyle\hskip 56.9055pt+\frac{\partial\Phi_{v}}{\partial{\bf f}_{2}}(.;{\bf f}_{2},{\bf u}){\bf h}_{2}-div(D_{v,1}(.;{\bf u})\nabla p)-div(D_{v,2}(.;{\bf u})\nabla\rho),\\ \displaystyle\frac{\partial\rho}{\partial t}-div(D_{w}(.;{\bf u})\nabla\rho))-div((\tilde{H}_{w;w}(.;{\bf u},\nabla{\bf u})-\frac{\partial\mathfrak{T}_{w}}{\partial w}(.;w))\rho)-\frac{\partial\Phi_{w}}{\partial w}(.;{\bf f}_{3},{\bf u})\rho\\ \displaystyle\hskip 42.67912pt=div(\tilde{H}_{w;u}(.;{\bf u},\nabla{\bf u})p)+div(\tilde{H}_{w;v}(.;{\bf u},\nabla{\bf u})q)+\frac{\partial\Phi_{w}}{\partial u}(.;{\bf f}_{3},{\bf u})p+\frac{\partial\Phi_{w}}{\partial v}(.;{\bf f}_{3},{\bf u})q\\ \displaystyle\hskip 56.9055pt+\frac{\partial\Phi_{w}}{\partial{\bf f}_{3}}(.;{\bf f}_{3},{\bf u}){\bf h}_{3}-div(D_{w,1}(.;{\bf u})\nabla p)-div(D_{w,2}(.;{\bf u})\nabla q),\end{array}$$

$$\begin{array}[]{lcr}(u,v,w)(.,0)=0~{}~{}~{}on~{}\Omega,\\ \displaystyle(D_{u}\nabla p).{\bf n}+\left(p\tilde{H}_{u;u}(.;{\bf u},\nabla{\bf u})+q\tilde{H}_{u;v}(.;{\bf u},\nabla{\bf u})+\rho\tilde{H}_{u;w}(.;{\bf u},\nabla{\bf u})\right).{\bf n}\\ \displaystyle\hskip 42.67912pt=(D_{u,1}\nabla q).{\bf n}+(D_{u,2}\nabla\rho).{\bf n},~{}~{}~{}on~{}\Sigma\\ \displaystyle(D_{v}\nabla q).{\bf n}+\left(p\tilde{H}_{v;u}(.;{\bf u},\nabla{\bf u})+q\tilde{H}_{v;v}(.;{\bf u},\nabla{\bf u})+\rho\tilde{H}_{v;w}(.;{\bf u},\nabla{\bf u})\right).{\bf n}\\ \displaystyle\hskip 42.67912pt=(D_{v,1}\nabla p).{\bf n}+(D_{v,2}\nabla\rho).{\bf n},~{}~{}~{}on~{}\Sigma\\ \displaystyle(D_{w}\nabla\rho).{\bf n}+\left(p\tilde{H}_{w;u}(.;{\bf u},\nabla{\bf u})+q\tilde{H}_{w;v}(.;{\bf u},\nabla{\bf u})+\rho\tilde{H}_{w;w}(.;{\bf u},\nabla{\bf u})\right).{\bf n}\\ \displaystyle\hskip 42.67912pt=(D_{w,1}\nabla p).{\bf n}+(D_{w,2}\nabla q).{\bf n},~{}~{}~{}on~{}\Sigma\end{array}$$

$$\begin{array}[]{lcr}\displaystyle\tilde{D}_{u;z}(.;{\bf u},\nabla u)=\frac{\partial D_{u}}{\partial z}(.;{\bf u})\nabla u,~{}~{}\displaystyle\tilde{D}_{v;z}(.;{\bf u},\nabla v)=\frac{\partial D_{v}}{\partial z}(.;{\bf u})\nabla v,\\ \displaystyle\tilde{D}_{w;z}(.;{\bf u},\nabla w)=\frac{\partial D_{w}}{\partial z}(.;{\bf u})\nabla w,\\ \displaystyle\tilde{T}_{u;z}(.;{\bf u},\nabla v,\nabla w)=\frac{\partial D_{u,1}}{\partial z}(.;{\bf u})\nabla v+\frac{\partial D_{u,2}}{\partial z}(.;{\bf u})\nabla w,\\ \displaystyle\tilde{T}_{v;z}(.;{\bf u},\nabla u,\nabla w)=\frac{\partial D_{v,1}}{\partial z}(.;{\bf u})\nabla u+\frac{\partial D_{v,2}}{\partial z}(.;{\bf u})\nabla w,\\ \displaystyle\tilde{T}_{w;z}(.;{\bf u},\nabla u,\nabla v)=\frac{\partial D_{w,1}}{\partial z}(.;{\bf u})\nabla u+\frac{\partial D_{w,2}}{\partial z}(.;{\bf u})\nabla v,\\ \displaystyle\tilde{H}_{u;z}(.;{\bf u},\nabla{\bf u})=\tilde{D}_{u;z}-\tilde{T}_{u;z},~{}~{}\tilde{H}_{v;z}(.;{\bf u},\nabla{\bf u})=\tilde{D}_{v;z}-\tilde{T}_{v;z},\\ \displaystyle\tilde{H}_{w;z}(.;{\bf u},\nabla{\bf u})=\tilde{D}_{w;z}-\tilde{T}_{w;z}.\end{array}$$

$$\begin{array}[]{lcr}\displaystyle-\frac{\partial\tilde{u}^{*}}{\partial t}-div(D_{u}(.;{\bf u}^{*})\nabla\tilde{u}^{*})-\frac{\partial\Phi_{u}}{\partial u}(.;{\bf f}^{*}_{1},{\bf u}^{*})\tilde{u}^{*}+(\tilde{H}_{u;u}(.;{\bf u}^{*},\nabla{\bf u}^{*})-\frac{\partial\mathfrak{T}_{u}}{\partial u}(.;u^{*}))\nabla\tilde{u}^{*}\\ \displaystyle\hskip 28.45274pt+\tilde{H}_{v;u}(.;{\bf u}^{*},\nabla{\bf u}^{*})\nabla\tilde{v}^{*}-\frac{\partial\Phi_{v}}{\partial u}(.;{\bf f}^{*}_{2},{\bf u}^{*})\tilde{v}^{*}\displaystyle+\tilde{H}_{w;u}(.;{\bf u}^{*},\nabla{\bf u}^{*})\nabla\tilde{w}^{*}-\frac{\partial\Phi_{w}}{\partial u}(.;{\bf f}^{*}_{3},{\bf u}^{*})\tilde{w}^{*}\\ \displaystyle\hskip 48.36958pt=-div(D_{v,1}(.;{\bf u}^{*})\nabla\tilde{v}^{*})-div(D_{w,1}(.;{\bf u}^{*})\nabla\tilde{w}^{*})+\frac{\partial{\cal Z}}{\partial u}({\bf u}^{*},{\bf f}^{*}),\\ \displaystyle-\frac{\partial\tilde{v}^{*}}{\partial t}-div(D_{v}(.;{\bf u}^{*})\nabla\tilde{v}^{*})-\frac{\partial\Phi_{v}}{\partial v}(.;{\bf f}^{*}_{2},{\bf u}^{*})\tilde{v}^{*}+(\tilde{H}_{v;v}(.;{\bf u}^{*},\nabla{\bf u}^{*})-\frac{\partial\mathfrak{T}_{v}}{\partial v}(.;v^{*}))\nabla\tilde{v}^{*}\\ \displaystyle\hskip 28.45274pt+\tilde{H}_{u;v}(.;{\bf u}^{*},\nabla{\bf u}^{*})\nabla\tilde{u}^{*}-\frac{\partial\Phi_{u}}{\partial v}(.;{\bf f}^{*}_{1},{\bf u}^{*})\tilde{u}^{*}+\tilde{H}_{w;v}(.;{\bf u}^{*},\nabla{\bf u}^{*})\nabla\tilde{w}^{*}-\frac{\partial\Phi_{w}}{\partial v}(.;{\bf f}^{*}_{3},{\bf u}^{*})\tilde{w}^{*}\\ \displaystyle\hskip 48.36958pt=-div(D_{u,1}(.;{\bf u}^{*})\nabla\tilde{u}^{*})-div(D_{w,2}(.;{\bf u}^{*})\nabla\tilde{w}^{*})+\frac{\partial{\cal Z}}{\partial v}({\bf u}^{*},{\bf f}^{*}),\\ \displaystyle-\frac{\partial\tilde{w}^{*}}{\partial t}-div(D_{w}(.;{\bf u}^{*})\nabla\tilde{w}^{*})-\frac{\partial\Phi_{w}}{\partial w}(.;{\bf f}^{*}_{3},{\bf u}^{*})\tilde{w}^{*}\displaystyle+(\tilde{H}_{w;w}(.;{\bf u}^{*},\nabla{\bf u}^{*})-\frac{\partial\mathfrak{T}_{w}}{\partial w}(.;w^{*}))\nabla\tilde{w}^{*}\\ \displaystyle\hskip 28.45274pt+\tilde{H}_{u;w}(.;{\bf u}^{*},\nabla{\bf u}^{*})\nabla\tilde{u}^{*}-\frac{\partial\Phi_{u}}{\partial w}(.;{\bf f}^{*}_{1},{\bf u}^{*})\tilde{u}^{*}\displaystyle+\tilde{H}_{v;w}(.;{\bf u}^{*},\nabla{\bf u}^{*})\nabla\tilde{v}^{*}-\frac{\partial\Phi_{v}}{\partial w}(.;{\bf f}^{*}_{2},{\bf u}^{*})\tilde{v}^{*}\\ \displaystyle\hskip 48.36958pt=-div(D_{u,2}(.;{\bf u})\nabla\tilde{u}^{*})-div(D_{v,2}(.;{\bf u}^{*})\nabla\tilde{v}^{*})+\frac{\partial{\cal Z}}{\partial w}({\bf u}^{*},{\bf f}^{*}),\\ \text{under the boundary conditions }\\ \displaystyle(\tilde{u}^{*}\frac{\partial\mathfrak{T}_{u}}{\partial u}(.;u^{*})-\tilde{u}^{*}\tilde{H}_{u;u}(.;{\bf u}^{*},\nabla{\bf u}^{*})+D_{u}(.;{\bf u}^{*})\nabla\tilde{u}^{*}).{\bf n}\\ \displaystyle\hskip 42.67912pt=(D_{v,1}(.;{\bf u}^{*})\nabla\tilde{v}^{*}+D_{w,1}(.;{\bf u}^{*})\nabla\tilde{w}^{*}).{\bf n}=0,\\ \displaystyle(\tilde{v}^{*}\frac{\partial\mathfrak{T}_{v}}{\partial v}(.;v^{*})-\tilde{v}^{*}\tilde{H}_{v;v}(.;{\bf u}^{*},\nabla{\bf u}^{*})+D_{v}(.;{\bf u}^{*})\nabla\tilde{v}^{*}).{\bf n}\\ \displaystyle\hskip 42.67912pt=(D_{u,1}(.;{\bf u}^{*})\nabla\tilde{u}^{*}+D_{w,2}(.;{\bf u}^{*})\nabla\tilde{w}^{*}).{\bf n}=0,\\ \displaystyle(\tilde{w}^{*}\frac{\partial\mathfrak{T}_{w}}{\partial v}(.;w^{*})-\tilde{w}^{*}\tilde{H}_{w;w}(.;{\bf u}^{*},\nabla{\bf u}^{*})+D_{w}(.;{\bf u}^{*})\nabla\tilde{w}^{*}).{\bf n}\\ \displaystyle\hskip 42.67912pt=(D_{u,2}(.;{\bf u}^{*})\nabla\tilde{u}^{*}-D_{v,2}(.;{\bf u}^{*})\nabla\tilde{v}^{*}).{\bf n}=0\\ \text{and the final condition}\\ \displaystyle\tilde{\bf u}^{*}(T)=\frac{\partial J_{0}}{\partial{\bf u}}({\bf u}^{*}(T))\end{array}$$

$$\begin{array}[]{lcr}\displaystyle\int_{0}^{T}\!\!\!\int_{\Omega}(\frac{\partial{\cal Z}}{\partial{\bf f}_{1}}({\bf u}^{*},{\bf f}^{*})+\tilde{u}^{*}\frac{\partial\Phi_{u}}{\partial{\bf f}_{1}}(.;{\bf f}^{*}_{1},{\bf u}^{*})).({\bf f}_{1}-{\bf f}^{*}_{1})dxdt\geq 0,\\ \displaystyle\int_{0}^{T}\!\!\!\int_{\Omega}(\frac{\partial{\cal Z}}{\partial{\bf f}_{2}}({\bf u}^{*},{\bf f}^{*})+\tilde{v}^{*}\frac{\partial\Phi_{v}}{\partial{\bf f}_{2}}(.;{\bf f}^{*}_{2},{\bf u}^{*})).({\bf f}_{2}-{\bf f}^{*}_{2})dxdt\geq 0,\\ \displaystyle\int_{0}^{T}\!\!\!\int_{\Omega}(\frac{\partial{\cal Z}}{\partial{\bf f}_{3}}({\bf u}^{*},{\bf f}^{*})+\tilde{w}^{*}\frac{\partial\Phi_{w}}{\partial{\bf f}_{3}}(.;{\bf f}^{*}_{3},{\bf u}^{*})).({\bf f}_{3}-{\bf f}^{*}_{3})dxdt\geq 0.\end{array}$$

$$\begin{array}[]{lcr}\displaystyle{\cal J}^{\prime}({\bf f}).{\bf h}=\lim_{\theta\longrightarrow 0^{+}}\frac{({\cal J}({\bf f}+\theta{\bf h})-{\cal J}({\bf f}))}{𝜖\theta}\\ \hskip 56.9055pt\displaystyle=\int_{\Omega}\frac{\partial J_{0}}{\partial{\bf u}}({\bf u}(T)).\wp(T)dx+\int_{0}^{T}\!\!\int_{\Omega}\frac{\partial{\cal Z}}{\partial{\bf u}}({\bf u},{\bf f}).\wp~{}dxdt\\ \displaystyle\hskip 56.9055pt+\int_{0}^{T}\!\!\int_{\Omega}\frac{\partial{\cal Z}}{\partial{\bf f}}({\bf u},{\bf f}).{\bf h}~{}dxdt.\end{array}$$

$$\begin{array}[]{lcr}\displaystyle\int_{\Omega}p(T)\tilde{u}(T)dx+\int_{0}^{T}\!\!\!\int_{\Gamma}p(\frac{\partial\mathfrak{T}_{u}}{\partial u}(.;u).{\bf n})\tilde{u}d\Gamma dt+\int_{0}^{T}\!\!\!\int_{\Gamma}p(D_{u}(.;{\bf u})\nabla\tilde{u}).{\bf n}d\Gamma dt\\ \displaystyle-\int_{0}^{T}\!\!\!\int_{\Gamma}p(\tilde{H}_{u;u}.{\bf n})\tilde{u}d\Gamma dt\displaystyle-\int_{0}^{T}\!\!\!\int_{\Gamma}q(D_{u,1}(.;{\bf u})\nabla\tilde{u}).{\bf n}d\Gamma dt-\int_{0}^{T}\!\!\!\int_{\Gamma}\rho(D_{u,2}(.;{\bf u})\nabla\tilde{u}).{\bf n}d\Gamma dt\\ \displaystyle+\int_{0}^{T}\!\!\!\int_{\Omega}p\left[-\frac{\partial\tilde{u}}{\partial t}-div(D_{u}(.;{\bf u})\nabla\tilde{u})-\frac{\partial\Phi_{u}}{\partial u}(.;{\bf f}_{1},{\bf u})\tilde{u}+(\tilde{H}_{u;u}-\frac{\partial\mathfrak{T}_{u}}{\partial u}(.;u))\nabla\tilde{u}\right]dxdt\\ \displaystyle=\int_{0}^{T}\!\!\!\int_{\Omega}q\left[-\tilde{H}_{u;v}\nabla\tilde{u}+\frac{\partial\Phi_{u}}{\partial v}(.;{\bf f}_{1},{\bf u})\tilde{u}\right]dxdt+\int_{0}^{T}\!\!\!\int_{\Omega}\rho\left[-\tilde{H}_{u;w}\nabla\tilde{u}+\frac{\partial\Phi_{u}}{\partial w}(.;{\bf f}_{1},{\bf u})\tilde{u}\right]dxdt\\ \displaystyle+\int_{0}^{T}\!\!\!\int_{\Omega}\left[\frac{\partial\Phi_{u}}{\partial{\bf f}_{1}}(.;{\bf f}_{1},{\bf u}){\bf h}_{1}\right]\tilde{u}dxdt\\ \displaystyle-\int_{0}^{T}\!\!\!\int_{\Omega}div(D_{u,1}(.;{\bf u})\nabla\tilde{u})qdxdt-\int_{0}^{T}\!\!\!\int_{\Omega}div(D_{u,2}(.;{\bf u})\nabla\tilde{u})\rho dxdt,\end{array}$$

$$\begin{array}[]{lcr}\displaystyle\int_{\Omega}q(T)\tilde{v}(T)dx+\int_{0}^{T}\!\!\!\int_{\Gamma}q(\frac{\partial\mathfrak{T}_{v}}{\partial v}(.;v).{\bf n})\tilde{v}d\Gamma dt+\int_{0}^{T}\!\!\!\int_{\Gamma}q(D_{v}(.;{\bf u})\nabla\tilde{v}).{\bf n}d\Gamma dt\\ \displaystyle-\int_{0}^{T}\!\!\!\int_{\Gamma}q(\tilde{H}_{v;v}.{\bf n})\tilde{v}d\Gamma dt\displaystyle-\int_{0}^{T}\!\!\!\int_{\Gamma}p(D_{v,1}(.;{\bf u})\nabla\tilde{v}).{\bf n}d\Gamma dt-\int_{0}^{T}\!\!\!\int_{\Gamma}\rho(D_{v,2}(.;{\bf u})\nabla\tilde{v}).{\bf n}d\Gamma dt\\ \displaystyle+\int_{0}^{T}\!\!\!\int_{\Omega}q\left[-\frac{\partial\tilde{v}}{\partial t}-div(D_{v}(.;{\bf u})\nabla\tilde{v})-\frac{\partial\Phi_{v}}{\partial v}(.;{\bf f}_{2},{\bf u})\tilde{v}+(\tilde{H}_{v;v}-\frac{\partial\mathfrak{T}_{v}}{\partial v}(.;v))\nabla\tilde{v}\right]dxdt\\ \displaystyle=\int_{0}^{T}\!\!\!\int_{\Omega}p\left[-\tilde{H}_{v;u}\nabla\tilde{v}+\frac{\partial\Phi_{v}}{\partial u}(.;{\bf f}_{2},{\bf u})\tilde{v}\right]dxdt+\int_{0}^{T}\!\!\!\int_{\Omega}\rho\left[-\tilde{H}_{v;w}\nabla\tilde{v}+\frac{\partial\Phi_{v}}{\partial w}(.;{\bf f}_{2},{\bf u})\tilde{v}\right]dxdt\\ \displaystyle+\int_{0}^{T}\!\!\!\int_{\Omega}\left[\frac{\partial\Phi_{v}}{\partial{\bf f}_{2}}(.;{\bf f}_{2},{\bf u}){\bf h}_{2}\right]\tilde{v}dxdt\\ \displaystyle-\int_{0}^{T}\!\!\!\int_{\Omega}div(D_{v,1}(.;{\bf u})\nabla\tilde{v})pdxdt-\int_{0}^{T}\!\!\!\int_{\Omega}div(D_{v,2}(.;{\bf u})\nabla\tilde{v})\rho dxdt\end{array}$$

$$\begin{array}[]{lcr}\displaystyle\int_{\Omega}\rho(T)\tilde{w}(T)dx+\int_{0}^{T}\!\!\!\int_{\Gamma}\rho(\frac{\partial\mathfrak{T}_{w}}{\partial w}(.;w).{\bf n})\tilde{w}d\Gamma dt+\int_{0}^{T}\!\!\!\int_{\Gamma}\rho(D_{w}(.;{\bf u})\nabla\tilde{w}).{\bf n}d\Gamma dt\\ \displaystyle-\int_{0}^{T}\!\!\!\int_{\Gamma}\rho(\tilde{H}_{w;w}.{\bf n})\tilde{w}d\Gamma dt\displaystyle-\int_{0}^{T}\!\!\!\int_{\Gamma}p(D_{w,1}(.;{\bf u})\nabla\tilde{w}).{\bf n}d\Gamma dt-\int_{0}^{T}\!\!\!\int_{\Gamma}q(D_{w,2}(.;{\bf u})\nabla\tilde{w}).{\bf n}d\Gamma dt\\ \displaystyle+\int_{0}^{T}\!\!\!\int_{\Omega}\rho\left[-\frac{\partial\tilde{v}}{\partial t}-div(D_{w}(.;{\bf u})\nabla\tilde{w})-\frac{\partial\Phi_{w}}{\partial w}(.;{\bf f}_{3},{\bf u})\tilde{w}+(\tilde{H}_{w;w}-\frac{\partial\mathfrak{T}_{w}}{\partial w}(.;w))\nabla\tilde{w}\right]dxdt\\ \displaystyle=\int_{0}^{T}\!\!\!\int_{\Omega}p\left[-\tilde{H}_{w;u}\nabla\tilde{w}+\frac{\partial\Phi_{w}}{\partial u}(.;{\bf f}_{3},{\bf u})\tilde{w}\right]dxdt+\int_{0}^{T}\!\!\!\int_{\Omega}q\left[-\tilde{H}_{w;v}\nabla\tilde{w}+\frac{\partial\Phi_{w}}{\partial v}(.;{\bf f}_{3},{\bf u})\tilde{w}\right]dxdt\\ \displaystyle+\int_{0}^{T}\!\!\!\int_{\Omega}\left[\frac{\partial\Phi_{w}}{\partial{\bf f}_{3}}({\bf f}_{3},{\bf u}){\bf h}_{3}\right]\tilde{w}dxdt\\ \displaystyle-\int_{0}^{T}\!\!\!\int_{\Omega}div(D_{w,1}(.;{\bf u})\nabla\tilde{w})pdxdt-\int_{0}^{T}\!\!\!\int_{\Omega}div(D_{w,2}(.;{\bf u})\nabla\tilde{w})qdxdt.\end{array}$$

$$\begin{array}[]{lcr}\displaystyle(\tilde{u}\frac{\partial\mathfrak{T}_{u}}{\partial u}(.;u)-\tilde{u}\tilde{H}_{u;u}+D_{u}(.;{\bf u})\nabla\tilde{u}-D_{v,1}(.;{\bf u})\nabla\tilde{v}-D_{w,1}(.;{\bf u})\nabla\tilde{w}).{\bf n}=0,\\ \displaystyle(\tilde{v}\frac{\partial\mathfrak{T}_{v}}{\partial v}(.;v)-\tilde{v}\tilde{H}_{v;v}+D_{v}(.;{\bf u})\nabla\tilde{v}-D_{u,1}(.;{\bf u})\nabla\tilde{u}-D_{w,2}(.;{\bf u})\nabla\tilde{w}).{\bf n}=0,\\ \displaystyle(\tilde{w}\frac{\partial\mathfrak{T}_{w}}{\partial w}(.;w)-\tilde{w}\tilde{H}_{w;w}+D_{w}(.;{\bf u})\nabla\tilde{w}-D_{u,2}(.;{\bf u})\nabla\tilde{u}-D_{v,2}(.;{\bf u})\nabla\tilde{v}).{\bf n}=0\end{array}$$

$$\begin{array}[]{lcr}\displaystyle\tilde{\bf u}(T)=\frac{\partial J_{0}}{\partial{\bf u}}({\bf u}(T)).\end{array}$$

$$\begin{array}[]{lcr}\displaystyle\int_{\Omega}\frac{\partial J_{0}}{\partial{\bf u}}({\bf u}(T)).\wp(T)dx\displaystyle-\int_{0}^{T}\!\!\!\int_{\Omega}\tilde{u}\frac{\partial\Phi_{u}}{\partial{\bf f}_{1}}(.;{\bf f}_{1},{\bf u}){\bf h}_{1}dxdt\\ \displaystyle-\int_{0}^{T}\!\!\!\int_{\Omega}\tilde{v}\frac{\partial\Phi_{v}}{\partial{\bf f}_{2}}(.;{\bf f}_{2},{\bf u}){\bf h}_{2}dxdt\displaystyle-\int_{0}^{T}\!\!\!\int_{\Omega}\tilde{w}\frac{\partial\Phi_{w}}{\partial{\bf f}_{3}}(.;{\bf f}_{3},{\bf u}){\bf h}_{3}dxdt\\ \displaystyle+\int_{0}^{T}\!\!\!\int_{\Omega}p\left[-\frac{\partial\tilde{u}}{\partial t}-div(D_{u}(.;{\bf u})\nabla\tilde{u})-\frac{\partial\Phi_{u}}{\partial u}(.;{\bf f}_{1},{\bf u})\tilde{u}+(\tilde{H}_{u;u}-\frac{\partial\mathfrak{T}_{u}}{\partial u}(.;u))\nabla\tilde{u}\right]dxdt\\ \displaystyle-\int_{0}^{T}\!\!\!\int_{\Omega}p\left[-\tilde{H}_{v;u}\nabla\tilde{v}+\frac{\partial\Phi_{v}}{\partial u}(.;{\bf f}_{2},{\bf u})\tilde{v}\right]dxdt\displaystyle-\int_{0}^{T}\!\!\!\int_{\Omega}p\left[-\tilde{H}_{w;u}\nabla\tilde{w}+\frac{\partial\Phi_{w}}{\partial u}(.;{\bf f}_{3},{\bf u})\tilde{w}\right]dxdt\\ \displaystyle+\int_{0}^{T}\!\!\!\int_{\Omega}div(D_{v,1}(.;{\bf u})\nabla\tilde{v})pdxdt+\int_{0}^{T}\!\!\!\int_{\Omega}div(D_{w,1}(.;{\bf u})\nabla\tilde{w})pdxdt\\ \displaystyle+\int_{0}^{T}\!\!\!\int_{\Omega}q\left[-\frac{\partial\tilde{v}}{\partial t}-div(D_{v}(.;{\bf u})\nabla\tilde{v})-\frac{\partial\Phi_{v}}{\partial v}(.;{\bf f}_{2},{\bf u})\tilde{v}+(\tilde{H}_{v;v}-\frac{\partial\mathfrak{T}_{v}}{\partial v}(.;v))\nabla\tilde{v}\right]dxdt\\ \displaystyle-\int_{0}^{T}\!\!\!\int_{\Omega}q\left[-\tilde{H}_{u;v}\nabla\tilde{u}+\frac{\partial\Phi_{u}}{\partial v}(.;{\bf f}_{1},{\bf u})\tilde{u}\right]dxdt-\int_{0}^{T}\!\!\!\int_{\Omega}q\left[-\tilde{H}_{w;v}\nabla\tilde{w}+\frac{\partial\Phi_{w}}{\partial v}(.;{\bf f}_{3},{\bf u})\tilde{w}\right]dxdt\\ \displaystyle+\int_{0}^{T}\!\!\!\int_{\Omega}div(D_{u,1}(.;{\bf u})\nabla\tilde{u})qdxdt+\int_{0}^{T}\!\!\!\int_{\Omega}div(D_{w,2}(.;{\bf u})\nabla\tilde{w})qdxdt\\ \displaystyle+\int_{0}^{T}\!\!\!\int_{\Omega}\rho\left[-\frac{\partial\tilde{w}}{\partial t}-div(D_{w}(.;{\bf u})\nabla\tilde{w})-\frac{\partial\Phi_{w}}{\partial w}(.;{\bf f}_{3},{\bf u})\tilde{w}\displaystyle+(\tilde{H}_{w;w}-\frac{\partial\mathfrak{T}_{w}}{\partial w}(.;w))\nabla\tilde{w}\right]dxdt\\ \displaystyle-\int_{0}^{T}\!\!\!\int_{\Omega}\rho\left[-\tilde{H}_{u;w}\nabla\tilde{u}+\frac{\partial\Phi_{u}}{\partial w}(.;{\bf f}_{1},{\bf u})\tilde{u}\right]dxdt\displaystyle-\int_{0}^{T}\!\!\!\int_{\Omega}\rho\left[-\tilde{H}_{v;w}\nabla\tilde{v}+\frac{\partial\Phi_{v}}{\partial w}(.;{\bf f}_{2},{\bf u})\tilde{v}\right]dxdt\\ \displaystyle+\int_{0}^{T}\!\!\!\int_{\Omega}div(D_{u,2}(.;{\bf u})\nabla\tilde{u})\rho dxdt+\int_{0}^{T}\!\!\!\int_{\Omega}div(D_{v,2}(.;{\bf u})\nabla\tilde{v})\rho dxdt=0.\end{array}$$

$$\begin{array}[]{lcr}\displaystyle-\frac{\partial\tilde{u}}{\partial t}-div(D_{u}(.;{\bf u})\nabla\tilde{u})-\frac{\partial\Phi_{u}}{\partial u}(.;{\bf f}_{1},{\bf u})\tilde{u}+(\tilde{H}_{u;u}-\frac{\partial\mathfrak{T}_{u}}{\partial u}(.;u))\nabla\tilde{u}\\ \displaystyle+\tilde{H}_{v;u}\nabla\tilde{v}-\frac{\partial\Phi_{v}}{\partial u}(.;{\bf f}_{2},{\bf u})\tilde{v}\displaystyle+\tilde{H}_{w;u}\nabla\tilde{w}-\frac{\partial\Phi_{w}}{\partial u}(.;{\bf f}_{3},{\bf u})\tilde{w}\\ \displaystyle=-div(D_{v,1}(.;{\bf u})\nabla\tilde{v})-div(D_{w,1}(.;{\bf u})\nabla\tilde{w})+\frac{\partial{\cal Z}}{\partial u}({\bf u},{\bf f}),\\[5.69046pt] \displaystyle-\frac{\partial\tilde{v}}{\partial t}-div(D_{v}(.;{\bf u})\nabla\tilde{v})-\frac{\partial\Phi_{v}}{\partial v}(.;{\bf f}_{2},{\bf u})\tilde{v}+(\tilde{H}_{v;v}-\frac{\partial\mathfrak{T}_{v}}{\partial v}(.;v))\nabla\tilde{v}\\ \displaystyle+\tilde{H}_{u;v}\nabla\tilde{u}-\frac{\partial\Phi_{u}}{\partial v}(.;{\bf f}_{1},{\bf u})\tilde{u}+\tilde{H}_{w;v}\nabla\tilde{w}-\frac{\partial\Phi_{w}}{\partial v}(.;{\bf f}_{3},{\bf u})\tilde{w}\\ \displaystyle=-div(D_{u,1}(.;{\bf u})\nabla\tilde{u})-div(D_{w,2}(.;{\bf u})\nabla\tilde{w})+\frac{\partial{\cal Z}}{\partial v}({\bf u},{\bf f}),\\[5.69046pt] \displaystyle-\frac{\partial\tilde{w}}{\partial t}-div(D_{w}(.;{\bf u})\nabla\tilde{w})-\frac{\partial\Phi_{w}}{\partial w}(.;{\bf f}_{3},{\bf u})\tilde{w}\displaystyle+(\tilde{H}_{w;w}-\frac{\partial\mathfrak{T}_{w}}{\partial w}(.;w))\nabla\tilde{w}\\ \displaystyle+\tilde{H}_{u;w}\nabla\tilde{u}-\frac{\partial\Phi_{u}}{\partial w}(.;{\bf f}_{1},{\bf u})\tilde{u}\displaystyle+\tilde{H}_{v;w}\nabla\tilde{v}-\frac{\partial\Phi_{v}}{\partial w}(.;{\bf f}_{2},{\bf u})\tilde{v}\\ \displaystyle=-div(D_{u,2}(.;{\bf u})\nabla\tilde{u})-div(D_{v,2}(.;{\bf u})\nabla\tilde{v})+\frac{\partial{\cal Z}}{\partial w}({\bf u},{\bf f}),\\[5.69046pt] \text{under the final and boundary conditions \eqref{ACF}-\eqref{ACB}}.\end{array}$$

$$\begin{array}[]{r}\displaystyle\int_{\Omega}\frac{\partial J_{0}}{\partial{\bf u}}({\bf u}(T)).\wp(T)dx+\int_{0}^{T}\!\!\int_{\Omega}\frac{\partial{\cal Z}}{\partial{\bf u}}({\bf u},{\bf f}).\wp dxdt=\int_{0}^{T}\!\!\!\int_{\Omega}\tilde{u}\frac{\partial\Phi_{u}}{\partial{\bf f}_{1}}(.;{\bf f}_{1},{\bf u}){\bf h}_{1}dxdt\\ \displaystyle\hskip 34.14322pt+\int_{0}^{T}\!\!\!\int_{\Omega}\tilde{u}\frac{\partial\Phi_{v}}{\partial{\bf f}_{2}}(.;{\bf f}_{2},{\bf u}){\bf h}_{2}dxdt+\int_{0}^{T}\!\!\!\int_{\Omega}\tilde{u}\frac{\partial\Phi_{w}}{\partial{\bf f}_{3}}(.;{\bf f}_{3},{\bf u}){\bf h}_{3}dxdt\end{array}$$

$$\begin{array}[]{lcr}\displaystyle\frac{\partial{\cal J}}{\partial{\bf f}_{1}}({\bf f}).{\bf h}_{1}=\int_{0}^{T}\!\!\!\int_{\Omega}(\frac{\partial{\cal Z}}{\partial{\bf f}_{1}}({\bf u},{\bf f})+\tilde{u}\frac{\partial\Phi_{u}}{\partial{\bf f}_{1}}(.;{\bf f}_{1},{\bf u})){\bf h}_{1}dxdt,\\ \displaystyle\frac{\partial{\cal J}}{\partial{\bf f}_{2}}({\bf f}).{\bf h}_{2}=\int_{0}^{T}\!\!\!\int_{\Omega}(\frac{\partial{\cal Z}}{\partial{\bf f}_{2}}({\bf u},{\bf f})+\tilde{v}\frac{\partial\Phi_{v}}{\partial{\bf f}_{2}}(.;{\bf f}_{2},{\bf u})){\bf h}_{2}dxdt,\\ \displaystyle\frac{\partial{\cal J}}{\partial{\bf f}_{3}}({\bf f}).{\bf h}_{3}=\int_{0}^{T}\!\!\!\int_{\Omega}(\frac{\partial{\cal Z}}{\partial{\bf f}_{3}}({\bf u},{\bf f})+\tilde{w}\frac{\partial\Phi_{w}}{\partial{\bf f}_{3}}(.;{\bf f}_{3},{\bf u})){\bf h}_{3}dxdt.\end{array}$$

$$\begin{array}[]{lcr}\displaystyle\frac{\partial{\cal J}}{\partial{\bf f}_{1}}({\bf f})=\frac{\partial{\cal Z}}{\partial{\bf f}_{1}}({\bf u},{\bf f})+\tilde{u}\frac{\partial\Phi_{u}}{\partial{\bf f}_{1}}(.;{\bf f}_{1},{\bf u}),\\ \displaystyle\frac{\partial{\cal J}}{\partial{\bf f}_{2}}({\bf f})=\frac{\partial{\cal Z}}{\partial{\bf f}_{2}}({\bf u},{\bf f})+\tilde{v}\frac{\partial\Phi_{v}}{\partial{\bf f}_{2}}(.;{\bf f}_{2},{\bf u}),\\ \displaystyle\frac{\partial{\cal J}}{\partial{\bf f}_{3}}({\bf f})=\frac{\partial{\cal Z}}{\partial{\bf f}_{3}}({\bf u},{\bf f})+\tilde{w}\frac{\partial\Phi_{w}}{\partial{\bf f}_{3}}(.;{\bf f}_{3},{\bf u}).\end{array}$$

$$\begin{array}[]{lcr}\displaystyle\int_{0}^{T}\!\!\!\int_{\Omega}(\frac{\partial{\cal Z}}{\partial{\bf f}_{1}}({\bf u}^{*},{\bf f}^{*})+\tilde{u}^{*}\frac{\partial\Phi_{u}}{\partial{\bf f}_{1}}(.;{\bf f}^{*}_{1},{\bf u}^{*})).({\bf f}_{1}-{\bf f}^{*}_{1})dxdt\geq 0,\\ \displaystyle\int_{0}^{T}\!\!\!\int_{\Omega}(\frac{\partial{\cal Z}}{\partial{\bf f}_{2}}({\bf u}^{*},{\bf f}^{*})+\tilde{v}^{*}\frac{\partial\Phi_{v}}{\partial{\bf f}_{2}}(.;{\bf f}^{*}_{2},{\bf u}^{*})).({\bf f}_{2}-{\bf f}^{*}_{2})dxdt\geq 0,\\ \displaystyle\int_{0}^{T}\!\!\!\int_{\Omega}(\frac{\partial{\cal Z}}{\partial{\bf f}_{3}}({\bf u}^{*},{\bf f}^{*})+\tilde{w}^{*}\frac{\partial\Phi_{w}}{\partial{\bf f}_{3}}(.;{\bf f}^{*}_{3},{\bf u}^{*})).({\bf f}_{3}-{\bf f}^{*}_{3})dxdt\geq 0,\end{array}$$

$$\begin{array}[]{lcr}\displaystyle 0=\frac{\partial{\cal Z}}{\partial{\bf f}_{1}}({\bf u}^{*},{\bf f}^{*})+\tilde{u}^{*}\frac{\partial\Phi_{u}}{\partial{\bf f}_{1}}(.;{\bf f}^{*}_{1},{\bf u}^{*}),\\ \displaystyle 0=\frac{\partial{\cal Z}}{\partial{\bf f}_{2}}({\bf u}^{*},{\bf f}^{*})+\tilde{v}^{*}\frac{\partial\Phi_{v}}{\partial{\bf f}_{2}}(.;{\bf f}^{*}_{2},{\bf u}^{*}),\\ \displaystyle 0=\frac{\partial{\cal Z}}{\partial{\bf f}_{3}}({\bf u}^{*},{\bf f}^{*})+\tilde{w}^{*}\frac{\partial\Phi_{w}}{\partial{\bf f}_{3}}(.;{\bf f}^{*}_{3},{\bf u}^{*}).\end{array}$$

$$\begin{array}[]{lcr}\displaystyle\frac{\partial{\cal J}}{\partial{\bf f}_{1}}({\bf f}^{(k)})=\frac{\partial{\cal Z}}{\partial{\bf f}_{1}}({\bf u}^{(k)},{\bf f}^{(k)})+\tilde{u}^{(k)}\frac{\partial\Phi_{u}}{\partial{\bf f}_{1}}({\bf f}_{1}^{(k)},{\bf u}^{(k)}),\\ \displaystyle\frac{\partial{\cal J}}{\partial{\bf f}_{2}}({\bf f}^{(k)})=\frac{\partial{\cal Z}}{\partial{\bf f}_{2}}({\bf u}^{(k)},{\bf f}^{(k)})+\tilde{v}^{(k)}\frac{\partial\Phi_{v}}{\partial{\bf f}_{2}}({\bf f}_{2}^{(k)},{\bf u}^{(k)}),\\ \displaystyle\frac{\partial{\cal J}}{\partial{\bf f}_{3}}({\bf f}^{(k)})=\frac{\partial{\cal Z}}{\partial{\bf f}_{3}}({\bf u}^{(k)},{\bf f}^{(k)})+\tilde{w}^{(k)}\frac{\partial\Phi_{w}}{\partial{\bf f}_{3}}({\bf f}_{3}^{(k)},{\bf u}^{(k)}),\\ \end{array}$$

$$\begin{array}[]{lcr}\displaystyle\zeta_{k}=\min_{\lambda>0}{\cal J}({\bf f}^{(k)}-\lambda E_{k}).\end{array}$$

Full text

First, the optimal control problems are formulated and first-order necessary optimality conditions are established by using sensitivity and adjoint calculus. The resulting problems are discretised by the coupled MRT and AMRT models and solved via gradient descent methods. First of all, an efficient and stable modified MRT model for NCACDCS is developed, and through the Chapman-Enskog analysis we show that NCACDCS can be correctly recovered from the proposed MRT model. For the adjoint problem, the discretisation strategy is based on AMRT model, which is found to be as simple as MRT model with also highly-efficient parallel nature. The derivation of AMRT model and the discrete cost functional gradient are derived mathematically in detail using the developed MRT model. The obtained method is reliable, efficient, practical to implement and can be easily incorporated into any existing MRT code.

keywords Optimal control, nonlinear coupled anisotropic convection-diffusion chemotaxis-type system, self and cross diffusion tensors, adjoint system, adjoint-based multiple-relaxation-time lattice Boltzmann method, multiscale Chapman-Enskog expansion, optimization.

MSC 2010: 92B05, 49J20, 35Q92, 49M25, 92C50, 65K15, 35K57, 65M99

1 Introduction and mathematical setting of the problem

The mathematical problem considered in this paper derives from the modeling of directed motion of biological organisms in response to changes of a number of chemical and physical cues in their environment. During the last few years, the effect of chemical gradients on organism migration has been widely studied, and complex systems have been developed to analyze organism behavior under simultaneous chemical, electrical, and mechanical stimuli.

In contrast to random diffusion without orientation, chemotaxis describes directed movements and migrations of bacteria, biological cells and organisms under the effect of a chemical stimulus (produced by a substance therein inhomogeneously distributed), in which organisms move either toward or away from the stimulus (i.e., attractive or repulsive movements); it is represented by the gradients of concentration of chemical substances (which have been characterized at molecular level). It is also by far the best understood mode of directed organism movements. The chemical signals are either secreted by the organisms themselves (autogenous) or come from foreign sources (xenogenous), and lead to aggregation of organisms and to formation of complex patterns.

Chemotaxis plays a crucial fundamental role, in physiology and pathological conditions, in many biological processes (see e.g., [31]) such as tissue formation, regeneration or immune functions, and pathological processes like cancer metastasis, as well as in many social activities of micro-organisms, e.g. virulence, social motility, fruiting body formation and biofilm formation.

Recently, numerous and various problems concerning chemotaxis have been intensively studied either from a theoretical or from a numerical point of view (the literature is extensive in this field). Many mathematical models have been developed to describe the chemotactic organism. The most generalized and widely-used model of chemotaxis is in connection with Keller-Segel models [45, 46] which are a coupled nonlinear/degenerate parabolic equations. Well-known examples are the dynamics and aggregation of Escherichia coli colonies and of the slime mould amoebae Dictyostelium discoideum. Also, endothelial cells react to vascular endothelial growth factor to form the primary blood vessels through aggregation (phenomenon of vasculogenesis) [2, 68]. Keller-Segel model has been successfully used by many scientists for the analysis of various biological phenomena. Others, after evaluations of this classical (minimal) chemotaxis model in specific situations, have concluded that many of the hypotheses (foundational to the model) are not sufficient to describe fundamental phenomena (in chemotaxis) in a satisfactory way. Then these latter ones have modified and generalized the model to adequate systems by incorporating additional informations (such as describing aggregation and blow up of attractive chemotaxis and collapse of repulsive chemotaxis, tumor angiogenesis and invasion, vasculogenesis, organism pattern formation, morphogenesis, attraction and repulsion between mutually interacting species in ecology and epidemiology, criminal behavior, etc) and obtain more complicated and more realistic models, see e.g., [2, 3, 5, 6, 14, 15, 16, 17, 28, 31, 36, 38, 39, 40, 41, 43, 52, 53, 54, 57, 58, 60, 61, 62, 64, 65, 68, 69, 75, 77, 78, 80], and the references therein. Furthermore, well-posedness results for these types of models, which are highly nonlinear systems, include, for example, existence, uniqueness and regularity of solutions (local, global, weak or radial), qualitative properties of positive solutions, uniform boundedness/blow-up for finite/infinite times and convergence of solutions to an equilibrium, and asymptotic behavior, see e.g., [8, 19, 21, 29, 30, 35, 44, 56, 59, 60, 63, 66, 72, 73, 76, 89], and the references therein. Recently, several numerical algorithms based on finite-volume and finite-element methods, upwind difference methods, operator-splitting methods, discontinuous Galerkin methods or simple-relaxation time lattice Boltzmann methods are used to solve these types of models (see e.g., [4, 25, 26, 32, 33, 34, 67, 70, 71, 81, 86], and the references therein).

Construction of quantitative models and predict optimal behaviors are among the main objectives of quantitative biology, which aim to understand complex cellular and organismal phenomena, at the system level lie, in a systematic, quantitative, and optimal manner. Optimal behavior is defined as an action that optimizes and evaluates the costs and benefits that influence the outcome of a decision, and contributes to an understanding of the system. The biologists control, in general, their experimental devices by using a certain number of functions or parameters of control which enable them to optimize and to stabilize the system, in accordance with the desired performance and targets.

To predict the dynamic response of interconnected biological systems from given parameters, data and source terms, requires a mathematical model of the behavior of process under investigation and physical and biological theory linking the state variables of the model to data and parameters. This prediction of the observation constitutes the so-called direct problem (primal problem, or also forward problem).

If any of the conditions necessary to define a direct problem are unknown or rather badly known, a control problem (or inverse problem) results, typically when modeling physico-chemical biology situations where the model parameters (intervening either in the boundary conditions, in initial conditions or in equations model itself) or material properties are unknown, difficult to measure biologically or partially known. Certain parameters or data can influence considerably the material behavior or modify phenomena in biological or medical matter; then their knowledge is an invaluable help to, for example, medical, biological and biosocial criminology applications. The resolution of these inverse problems thus provides them essential informations which are necessary to comprehension of various processes which can intervene in these models. In addition to mathematical model (primal problem), the formulation of control problem in form of a known cost (evaluation or objective) function ${\cal J}$ is a pre-requirement for achieving desired optimal conditions. This resolution, in order to ensure sufficiently smoothness of cost function with respect to the variation of design parameters and data, needs some regularity and additional conditions, and partial informations of some unknown parameters and fields (observations) given, for example, by experiment measurements. For optimal control in chemotaxis-type systems, we can mention e.g., [9, 13, 50, 55, 74], and the references therein.

The problem studied in this paper derives from nonlinear optimal control of a class of dynamic systems which are characterized by chemotaxis. In order to target a wide range of applications, we have written the continuous chemotaxis-type models in a generalized form by considering a class of coupled nonlinear/degenerate anisotropic convection-diffusion chemotaxis-type systems with Neumann-type boundary conditions (for simplicity here we use null flux conditions). We have also opted to place the control variables in all equations of the coupled systems. To facilitate the presentation, we consider a system with only three equations (the developed analysis will remains valid if we consider a system with $N$ equations ($N\geq 4$) or a coupled systems of reaction-diffusion and ordinary differential equations, as for example the system modeling atherogenesis or habitat-destruction tumor model). More precisely, we consider the following class of system

[TABLE]

subject to the following homogeneous Neumann boundary conditions

[TABLE]

where $(\mathbf{x},t)$ are the space-time variables, ${\bf n}$ being the outward normal to $\Gamma$, ${\bf f}=({\bf f}_{1},{\bf f}_{2},{\bf f}_{3})$, where ${\bf f}_{i}:(\mathbf{x},t)\longrightarrow{\bf f}_{i}(\mathbf{x},t)\in\mbox{IR}^{d_{i}}$ with $d_{i}\geq 1$, represents the input of the system (source terms, parameters and others) and ${\bf u}=(u,v,w)$ (where $u$,$v$ and $w$ are a scalar function) represents the state or output of the system (e.g., attractant concentrations, attractiveness values, density of the bacterial/cellular/criminal species, nutrient concentrations, etc), $(\Phi_{u},\Phi_{v},\Phi_{w})$ represents the kinetics/source term. The operators $D_{\Theta}$, for $\Theta=u,v$ or $w$ are the self diffusion tensors and $D_{\Theta,k}$, for $\Theta=u,v$ or $w$ and $k=1,3$ are cross-diffusion tensors. These last tensors are assumed to be symmetric and sufficiently regular. The domain $\Omega$ is an open bounded subset of $\mbox{IR}^{m}$, $m\leq 3$, with the Lipschitz-continuous boundary $\Gamma=\partial\Omega$ and the time $T>0$ is a fixed constant (a given final time). Furthermore we define ${\cal Q}=\Omega\times(0,T)$ and $\Sigma=\Gamma\times(0,T)$.

The nonlinear operators, $D$ (tensor-valued), $\mathfrak{T}$ (vector-valued) and $\Phi$ (scalar-valued) in (1) are assumed to be Carathéodory functions.

Remark 1**.**

If the tensor operators $D_{\Theta}$ and $D_{\Theta,k}$ (for $\Theta=u,v$ or $w$ and $k=1,2$) are of the form $D_{\Theta}=\sigma_{\Theta}{\bf I}_{d}$ and $D_{\Theta,k}=\sigma_{\Theta,k}{\bf I}_{d}$, with ${\bf I}_{d}$ identity matrix and $\sigma_{\Theta}$, $\sigma_{\Theta,k}$ scalar operators, to close the system (1) we can impose the following boundary conditions: $\displaystyle\nabla u.{\bf n}=0,~{}~{}\nabla v.{\bf n}=0\text{~{}~{}~{}and~{}~{}~{}}\nabla w.{\bf n}=0~{}~{}on~{}\Sigma.$ $\Box$

Remark 2**.**

To act in accordance with a desired realistic situation, one can change the boundary conditions (2) by considering, for example, mixed boundary conditions (Dirichlet, Neumann, linear/nonlinear Robin-type, etc) on different portions $\Gamma_{i}$ (i=1,p) of $\displaystyle\Gamma=\partial\Omega=\cup_{i=1,p}\Gamma_{i}$ (such that $\displaystyle\cap_{i=1,p}\Gamma_{i}=\emptyset$). $\Box$

In order to solve numerically the considered optimal control problems, this work proposes and develops an adjoint-based coupling multiple-relaxation-time lattice Boltzmann model, reliable, efficient, stable and easy to implement in the context of general anisotropic reaction-diffusion systems, with Neumann (Robin) type boundary conditions in complex geometry boundaries, coupled to gradient-based algorithms. While a huge number of papers has been devoted to adjoint optimisation approach for PDE-constrained optimal control problems, there are few papers dedicated to combination with adjoint problem and lattice Boltzmann model, and most of them relates to the fluid dynamics with the classical and simple-relaxation-time lattice Boltzmann method (which is limited to the description of isotropic diffusion problems), see e.g., [37, 47, 49, 83] and the references therein. Recently, in [7] we have developed a new mathematical framework which models the competition between tumor and normal cells under chemotherapy constraints in treatment of brain tumors. The goal is to predict the distribution and necessary quantity of drugs delivered in drug-therapy by using a control problem. The control estimates simultaneously blood perfusion rate, reabsorption rate of drug and drug dosage administered, which affect the effects of brain-tumor chemotherapy. In order to solve numerically this problem, we have proposed and investigated an adjoint multiple-relaxation-time lattice Boltzmann method for a class of nonlinear coupled anisotropic convection-diffusion system (which includes the developed model for brain tumor targeted drug delivery system). One of the purposes of this paper is to extend the numerical method developed in [7] to systems derive from the modeling of directed motion of biological organisms in response to changes of a number of chemical and physical cues in their environment. The main focus of this method is to express the gradient of cost functional from the adjoint of coupled multiple-relaxation-time lattice Boltzmann system.

Lattice Boltzmann model (LBM) was originated from Boltzmann’s kinetic theory of gases (70s), and attracts more and more attentions for simulating complex fluid flows since 90s. More recently, LBM has been extended successfully to simulate different types of parabolic reaction-diffusion equation in biological model as monodomain and bidomain model in cardiac electrophysiology, see [25, 23] and Keller-Segel chemotaxis model see [81] and the references therein. In order to take into account anisotropic diffusion problems, multiple-relaxation-time LBM has been developed for nonlinear anisotropic convection-diffusion equations, see e.g. [7, 42, 88] and the references therein. Some LBM models have been developed to be kinetic models (LBKM) so that the non-equilibrium flow systems can be better accessed (the kinetic moment relations used in developing this method are physical in comparison with thus used for LBM), see e.g., [51, 79] and the references therein.

The concept of LBM, which is efficient in terms of parallelization, is based on Boltzmann equations which describe the evolution of particles in kinetic theory. Indeed, traditional numerical methods as finite difference method or finite element method directly solve governing equations for deriving macroscopic variable, whereas LBM is based on the particle (the discrete) distribution function and numerical solving of continuous Boltzmann transport equation. Then the macroscopic variables of the considered system can be recovered from the discrete equations through multi-scaling Chapman-Enskog expansion procedure. LBM has two main phases: the first one, local, models the collision between particles and the second one, along each direction of interpolation, models the transport phase.

The rest of paper is organized as follows. In the next section, we start by formulating the control problem in terms of a general control configuration. Afterwards, we establish the (primal-dual) first order optimality conditions for the optimal solution by using sensitivity and adjoint calculus. The optimality system requires calculation of gradients of the given cost functional which are also necessary to develop numerical optimization algorithms. For ease of understanding, some examples are presented. Some generalities about a possible strategy for numerical realization of optimal control problems based on gradient optimization algorithms (by using adjoint variables) are also given. In Section 3, first, we develop our coupled multiple-relaxation-time LBM (MRT) for numerical simulations of the primal problem and investigate its asymptotic behavior. Next, we describe the adjoint-based strategy, derive a multiple-relaxation-time adjoint system associated to MRT and then gradients of the cost functional. Finally, comments are presented in Section 4.

2 Study of control framework

2.1 Formulation of the problem

We suppose that there exist ${\cal V}_{0}(\Omega)$, ${\cal U}={\cal U}_{1}\times{\cal U}_{2}\times{\cal U}_{3}$ (where ${\cal U}_{i}=\prod_{j=1}^{d_{i}}{\cal U}^{j}_{i}$, for $i=1,3$) and ${\cal W}({\cal Q})$ three Banach spaces such that for initial variables and input variables, respectively, in ${\cal V}_{0}(\Omega)$ and ${\cal U}$, the problem (1) with boundary conditions (2) (under some hypotheses for the data and some regularity of the nonlinear operators) is well-posed in the sense of Hadamard and the unique solution ${\bf u}$ is in ${\cal W}({\cal Q})$. Moreover, we assume that the nonlinear operator are Fréchet differentiable in order to ensure the Fréchet differentiability of the control-to-state operator in suitable Banach spaces.

Assume now that the initial condition ${\bf u}_{0}$ is a given function in ${\cal V}_{0}(\Omega)$ (fixed) and introduce the mapping ${\cal F}$ which maps the input ${\bf f}\in{\cal U}$ into the corresponding solution ${\bf u}={\cal F}({\bf f})\in{\cal W}({\cal Q})$ of (1)-(2). We denote by ${\cal K}_{ad}={\cal K}^{ad}_{1}\times{\cal K}^{ad}_{2}\times{\cal K}^{ad}_{3}\subset{\cal U}$ (where ${\cal K}^{ad}_{i}=\prod_{j=1}^{d_{i}}{\cal K}^{ad,j}_{i}$, for $i=1,3$) the set of admissible controls which is a given non-empty, closed and convex subset of ${\cal U}$.

The nonlinear optimal problem that we are going to study may be present in a very large class of problems, which can be formally described as follows : find ${\bf f}^{*}\in{\cal K}_{ad}$ such that the optimal criteria which is expressed by the following cost functional ${\cal J}$ (in the reduced form)

[TABLE]

is minimized with respect to ${\bf f}$ subject to system (1)-(2). In an other way, we will study

[TABLE]

The Fréchet differential operator $J_{0}$ and ${\cal Z}$ are assumed to be sufficiently regular to ensure sufficiently smoothness of cost function ${\cal J}$ with respect to the variation of design ${\bf f}$.

Remark 3**.**

The cost functional ${\cal J}$ describing the control problem depends on function ${\bf f}$ and state function ${\bf u}$ in domain $\Omega$ over the time interval under consideration $[0,T]$. It also depends on a given measurement data $\mathfrak{m}_{data}({\bf x},t)\in\mbox{IR}^{s}$ and reference data $\mathfrak{r}_{data}({\bf x},t)\in\mbox{IR}^{r}$. In order to simplify the presentation, we have used in expression (3), the reduced form of the functional i.e., ${\cal J}({\bf f})$ in the place of classical form ${\cal J}({\bf f},{\cal F}({\bf f}))$. More classical cost functional $\cal J$ is given by

[TABLE]

where $p_{i,j}\geq 2$, the functions $f^{(j)}_{r,i}$, for $i=1,3$ and $j=1,d_{i}$, are given reference data, the data $u_{f},v_{f},w_{f},u_{d},v_{d}$ and $w_{d}$ are the target profile for state functions (observations) and the coefficients $a_{k}\geq 0$, $b_{k}\geq 0$ ($k=1,3$) and $\alpha_{i,j}\geq 0$ ($i=1,3$* and $j=1,d_{i}$) such that $\sum_{k=1}^{3}(a_{k}+b_{k})>0$ and $\sum_{i=1}^{3}\sum_{j=1}^{d_{i}}\alpha_{i,j}>0$ are chosen as constants to establish the relative weight of the terms that appear in the definition of the functional.* 2. 2.

In accordance with practical applications and available clinical or experimental observations, it is clear that, we can consider other observations and/or controls, which can appear in boundary conditions, in initial conditions or in cross-diffusion tensors (see e.g., **[10, 11, 12]**). The study carried out in this paper, by changing the boundary conditions accordingly, will remain valid without major modifications. $\Box$

Remark 4**.**

If ${\cal K}^{ad,j}_{i}$ is for example of the form $\{g\in L^{2}({\cal Q})~{}:~{}m_{1}\leq g\leq m_{2}~{}~{}a.e.~{}in~{}{\cal Q}\}$ (where $m_{1}$ and $m_{2}$ are two constants), although it is a subset of $L^{\infty}({\cal Q})$, we prefer to use the standard norms of the space $L^{2}({\cal Q})$. The reason is that we would like to take advantages of the differentiability of the latter norm away from the origin to perform our variational analysis. $\Box$

2.2 Some examples

To illustrate the considered general problem, in this paragraph, we give some application examples.

2.2.1 Crime problem

In this application we consider a reaction-diffusion system involving mobile criminal offenders as

[TABLE]

The state $u$ is the attractiveness field, the state $v$ is the expectation value of the criminal density, and the number of burglaries being committed at time $t$ and location ${\bf x}$ is given by $u({\bf x},t)v({\bf x},t)$. The parameter $f^{(1)}_{1}$ is an attractive force pulling offenders back to locations where they have successfully committed crimes. The criminal agents are being created at a rate $f_{2}$ and are removed from the model when a burglary is committed. The function $f^{(2)}_{1}$ represents for example a static, spatially varying, component of the attractiveness. The function ${\bf f}=(f^{(1)}_{1},f^{(2)}_{1},f_{2})$ influences considerably the dynamics of criminals and of attractiveness field, and then their control will be very beneficial and great help to predict accurately several regimes of aggregation, including hotspots of high criminal activity.

2.2.2 Attraction-repulsion chemotaxis type model with logistic source

Our second example concerns an attraction-repulsion chemotaxis-type model with fast diffusion and logistic source

[TABLE]

where $F$ can be given by $F(u)=a_{1}u-a_{2}u^{\kappa}$ (with $\kappa\geq 1$ and the positive parameters $a_{1}$, $a_{2}$ describing respectively the organism growth rate and the carrying capacity), the states $u$, $v$, $w$ are respectively the cell density, the concentration of an attractive chemical signal, and the concentration of a repulsive chemical signal, and the vector field $\vec{\omega}$ is the flow velocity which is generated by a force due to aggregation of cells. The positive parameters representing the chemotactic behavior are $f^{(1)}_{2}$ and $f^{(1)}_{3}$ which represent the chemical production rate, $f^{(2)}_{2}$ and $f^{(2)}_{3}$ which represent the chemical degradation rate, $f^{(3)}_{2}$ which represents the chemical consumed rate, $f^{(3)}_{3}$ which corresponds source function, and the matrix-valued operators $\chi_{k}$, $k=1,2$ which represent chemotactic signal of attraction and repulsion, respectively. The operator $\sigma_{1}$ can be described for example by $\sigma_{1}(u)=\frac{1}{(1-u)^{\alpha}}$ with $\alpha\geq 2$, by $\sigma_{1}(u)=\frac{1+u}{(1-u+ln(u))}$ or by $\sigma_{1}(u)=u^{p-1}$ (with $p>1$). The operators $\sigma_{2}$ and $\sigma_{3}$ are assumed to be constant. Finally, the operator $\chi_{k}$ can be described for example by $\chi_{k}({\bf x},t;u,\zeta)=\mathbb{X}({\bf x},t)u{\bf I}_{d}$, or by $\chi_{k}({\bf x},t;u,\zeta)=\mathbb{X}({\bf x},t)\frac{u}{1+\zeta^{2}}{\bf I}_{d}$, where $\mathbb{X}$ is a sufficiently regular function.

We can control some of parameters/source $f^{(j)}_{i}$ from par example observations given by measurements.

2.2.3 Two species chemotaxis blow up

The third example deals with two-species chemotaxis system with competitive kinetics

[TABLE]

where $F_{i}$ (for $i=1,3$) can be given by $F_{1}(u,v)=\mu_{1}u(1-u-a_{1}v)$, $F_{2}(u,v)=\mu_{2}v(1-v-a_{2}u)$ and $F_{3}(u,v,w,;{\bf f}_{3})=-w+f^{(1)}_{3}u+f^{(2)}_{3}v+f^{(3)}_{3}$or $-w(f^{(1)}_{3}u+f^{(2)}_{3}v)+f^{(3)}_{3}$ with ${\bf f}_{3}=(f^{(1)}_{3},f^{(2)}_{3},f^{(3)}_{3})$, the state $u$ and $v$ denote the cell densities of the first and second species, respectively, and $w$ presents the concentration of the chemical signal. The parameters present in the model are positive and the operators are of the same type as those of the previous application.

2.2.4 Tumor invasion into surrounding environment and habitat-destruction type model

In our last example, we propose a modified version of the haptotaxis model introduced in [22]. The model, which takes into account the competition for space between cancer cells and extracellular matrix macromolecules (ECM), consists of an ODE reflecting the degradation of ECM, coupled with a reaction-diffusion chemotaxis-haptotaxis system describing the evolution of cancer cell density (which uses ECM for movement), the matrix degrading enzyme concentration and the oxygen concentration (which is produced by the matrix molecules). The system is of the form

[TABLE]

where the states $u,v,w$ and $m$ are cancer cell density, matrix degradation enzyme concentration, oxygen concentration and extracellular matrix density, respectively. The parameters present in the model are positive.

We can identify (estimate) the best optimal prognostic values of some parameters appearing in the model or control the functions $f$ and $g$ from observations (desired target).

2.3 Continuous adjoint-based optimality conditions

Assume that the nonlinear control problem (4) admits an optimal solution; the necessary conditions for this optimum is given by the following theorem (see [13]).

Theorem 1**.**

If ${\cal J}$ attains a (local) minimum at a point ${\bf f}^{*}\in{\cal K}_{ad}$, then the following first optimality conditions hold:

[TABLE]

where ${\cal J}^{\prime}$ is the directional derivative of ${\cal J}$.

In order to solve numerically the optimal control problem, it is necessary to derive the gradient of the cost functional ${\cal J}$𝐽 with respect to the control ${\bf f}=({\bf f}_{1},{\bf f}_{2},{\bf f}_{3})$. For this, we suppose that the operator solution ${\cal F}$ is continuously differentiable on ${\cal K}_{ad}$ and its derivative,𝑤 at ${\bf f}$ in direction ${\bf h}=({\bf h}_{1},{\bf h}_{2},{\bf h}_{3})$, $\displaystyle{\wp}=(p,q,\rho)={\cal F}^{\prime}({\bf f}).{\bf h}=\lim_{\theta\longrightarrow 0^{+}}\frac{({\cal F}({\bf f}+\theta{\bf h})-{\cal F}({\bf f}))}{𝜖\theta}$ is the unique solution of the following system (which is called the tangent linear model (TLM) or sensitivity problem)

[TABLE]

under the initial and Robin-type boundary conditions:

[TABLE]

where

[TABLE]

We can now show the first-order necessary conditions (optimality conditions) and calculate the gradient of ${\cal J}$ by using TLM and by introducing an intermediate adjoint (dual ou costate) model.

Theorem 2**.**

Suppose that $({\bf u}^{*},{\bf f}^{*})\in{\cal W}({\cal Q})\times{\cal K}_{ad}$, is an optimal solution such that ${\bf f}^{*}=({\bf f}^{*}_{1},{\bf f}^{*}_{2},{\bf f}^{*}_{3})$ is defined by (4) and ${\bf u}^{*}={\cal F}({\bf f}^{*})$ is the solution of (1)-(2). Then there exists an adjoint variable $\tilde{\bf u}^{*}=(\tilde{u}^{*},\tilde{v}^{*},\tilde{w}^{*})$ satisfying the following adjoint problem

[TABLE]

and the following variational inequality (for all ${\bf f}=({\bf f}_{1},{\bf f}_{2},{\bf f}_{3})\in{\cal K}_{ad}$)

[TABLE]

Proof. By using the same technique as in [13], we start by calculating the variation of ${\cal J}$.

According to the regularity of operator solution ${\cal F}$ and the nature of cost function ${\cal J}$ (which is the composition of Fréchet differentiable mappings), we have that ${\cal J}$ is differentiable and the directional derivative of ${\cal J}$ at point ${\bf f}$ along the direction ${\bf h}=({\bf h}_{1},{\bf h}_{2},{\bf h}_{3})$ can be given by

[TABLE]

Now, we simplify the directional derivative of ${\cal J}$. For this we multiply the first part of (11) by some regular function $\tilde{{\bf u}}=(\tilde{u},\tilde{v},\tilde{w})$, integrating over ${\cal Q}$, using Green’s formula and integrating by parts in times, we obtain (according to the boundary and initial conditions for ${\wp}$)

[TABLE]

[TABLE]

and

[TABLE]

Assume that the function $\tilde{\bf u}$ satisfies the following boundary conditions on $\Sigma$ :

[TABLE]

and the final condition

[TABLE]

Then, by summing the three relations of previous system, we obtain (according to (16) and (17))

[TABLE]

In order to simplify (15), according to (18), we suppose that $\tilde{\bf u}=(\tilde{u},\tilde{v},\tilde{w})$ is the solution of the following adjoint system

[TABLE]

Then the previous relation becomes

[TABLE]

and then

[TABLE]

Otherwise (in weak sense)

[TABLE]

Since ${\bf f}^{*}$ is an optimal solution we have ($\forall{\bf f}=({\bf f}_{1},{\bf f}_{2},{\bf f}_{3})\in{\cal K}_{ad}$)

[TABLE]

with ${\bf u}^{*}={\cal F}({\bf f}^{*})$ and $\tilde{\bf u}^{*}$ the solution of (19) (corresponding to $({\bf f}^{*},{\bf u}^{*})$). This completes the proof. $\Box$

Remark 5**.**

*1. The solution $\tilde{\bf u}$ of (19) is said to be the costate, adjoint or dual solution corresponding to the primal solution ${\bf u}={\cal F}({\bf f})$ and will be denoted in the sequel by $\tilde{\bf u}=\tilde{\cal F}({\bf f},{\bf u})$.

2. If ${\cal K}_{ad}={\cal U}$ (i.e., optimal control without constraints) the optimality condition (14) become*

[TABLE]

Remark 6**.**

In order to solve the nonlinear control problem numerically, by using the adjoint variables, we can combine, for example, the obtained optimal necessary conditions and the gradient-iterative algorithm or Newton-iterative algorithm (see e.g., [13]). This resolution requires, at each iteration of optimization algorithm, the numerical resolution of direct problem and its corresponding adjoint problem. The direct, sensitivity and adjoint problems, which are coupled systems of anisotropic convection-diffusion chemotaxis-type equations, can be solved by using a multiple-relaxation-time lattice Boltzmann method. $\Box$

We end this section by a description of a gradient algorithm to solve Problem (4). The method is formulated in terms of continuous variables which are independent of a specific numerical discretization.

2.4 Gradient-iterative algorithm and continuous adjoint models

For a given function $J_{0}$ and ${\cal Z}$ (and then measurement data $\mathfrak{m}_{data}$ and reference data $\mathfrak{r}_{data}$) and initial states $(u_{0},v_{0},w_{0})$, we present a gradient-iterative algorithm where the descent direction is calculated by using adjoint variables, particularly by choosing an admissible step size. The gradient algorithm is given bellows (for $k=0,\ldots$, (iteration index) we denote by ${\bf f}^{(k)}=({\bf f}^{(k)}_{1},{\bf f}^{(k)}_{2},{\bf f}^{(k)}_{3})$ the approximation of control variables at the $kth$ iteration of the algorithm):

Step 1: Initialization: ${\bf f}^{(0)}=({\bf f}^{(0)}_{1},{\bf f}^{(0)}_{2},{\bf f}^{(0)}_{3})$ (given initial guess).

Step 2: Solve problem (1)-(2) with source term ${\bf f}^{(k)}$, gives ${\bf u}^{(k)}=(u^{(k)},v^{(k)},w^{(k)})={\cal F}({\bf f}^{(k)})$.

Step 3: Solve problem (13) (based on ${\bf u}^{(k)}$), gives adjoint state $\tilde{\cal F}({\bf f}^{(k)},{\bf u}^{(k)})=(\tilde{u}^{(k)},\tilde{v}^{(k)},\tilde{w}^{(k)})$.

Step 4: Gradient of ${\cal J}$, $E_{k}={\cal J}^{\prime}({\bf f}^{(k)})=(\frac{\partial{\cal J}}{\partial{\bf f}_{1}}({\bf f}^{(k)}),\frac{\partial{\cal J}}{\partial{\bf f}_{2}}({\bf f}^{(k)}),\frac{\partial{\cal J}}{\partial{\bf f}_{3}}({\bf f}^{(k)}))$ at point ${\bf f}^{(k)}$ is given by

[TABLE]

Step 5: Determine ${\bf f}^{(k+1)}$ : $\displaystyle{\bf f}^{(k+1)}={\bf f}^{(k)}-\zeta_{k}E_{k}$,

where $\zeta_{k}$ with $0<m_{0}\leq\zeta_{k}\leq M_{0}$ are the sequences of step lengths.

Step 6: If gradient of ${\cal J}$ is sufficiently small: end; else set $k:=k+1$ and goto Step 2.

Optimal Solution: ${\bf f}^{*}:={\bf f}^{(k)}$ and $(u^{*},v^{*},w^{*}):=(u^{(k)},v^{(k)},w^{(k)})$.

Remark 7**.**

1. In order to obtain an algorithm which is numerically efficient, the best choice of $\zeta_{k}$ will be the result of a line minimization algorithm. Otherwise, at each iteration step $k$ of the previous algorithm, we solve the one-dimensional optimization problem of parameter $\zeta_{k}$:

[TABLE]

To derive an approximation for $\zeta_{k}$ we can use a purely heuristic approach, e.g., by taking $\zeta_{k}=min(1,\parallel E_{k}\parallel^{-1}_{\infty})$ or by using the linearization of ${\cal F}({\bf f}^{(k)}-\lambda E_{k})$ at ${\bf f}^{(k)}$ by : $\displaystyle{\cal F}({\bf f}^{(k)}-\lambda E_{k})\approx{\cal F}({\bf f}^{(k)})-\lambda{\wp}^{(k)}={\bf u}^{(k)}-\lambda{\wp}^{(k)}$, where ${\wp}^{(k)}={\cal F}^{\prime}({\bf f}^{(k)}).E_{k}=(p^{(k)},q^{(k)},\rho^{(k)})$ is the solution of the sensitivity problem (11). For example in the case of quadratic cost function ${\cal J}$, according to the previous approximation, we can approximate (25) by

[TABLE]

*where $P$ is a polynomial function of degree $2$ (since ${\cal J}$ is quadratic), then problem (26) can be solved exactly. Consequently, we obtain explicitly the value of $\zeta_{k}$.

2. We can write similarly the conjugate gradient algorithm combined with the Wolfe-Powell line search procedure for computing admissible step-sizes along the descent direction. The advantage of this method, compared to the gradient method, is that it performs a soft reset whenever the GC search direction yields no significant progress. In general, the method has the following form:*

[TABLE]

where $G_{k}$ denotes gradient of the functional to be optimized at point $z_{k}$, $\zeta_{k}$ is a step length obtained by a line search, $D_{k}$ is the search direction and $r_{k}$ is a constant. Several varieties of this method differ in the way of selecting $r_{k}$. Some well-known formula for $r_{k}$ are given by Fletcher-Reeves (FR), Polak-Ribière (PR), Hestenes-Stiefel (HS) and Dai-Yuan (DY). $\Box$

3 Adjoint-based numerical optimization approaches

For given functions $J_{0}$ and ${\cal Z}$, and initial states ${\bf u}_{0}=(u_{0},v_{0},w_{0})$, in numerical treatments of nonlinear optimal control problems, the direct system, sensitivity system, adjoint system and objective functional must be discretized (reduction of the infinite-dimensional dynamics to finite-dimensional problems). The discretized formulation for direct, sensitivity and adjoint problems can be performed by combining the classical Galerkin and finite element methods, to the variational formulations, for the space discretization and Euler method for the time discretization, or by using lattice Boltzmann methods. It is clear that the amount of computation is crucial in solution of control problems, because the computation of discrete gradient of objective functional by discrete adjoint methods requires one forward solve of discrete state system and one backward solve of adjoint systems in which the state trajectory is an input. In order to reduce the implementation cost to evaluate discrete gradients, we propose and develop an adjoint multiple-relaxation-time coupled lattice Boltzmann method which is efficient, stable and easy to implement and to parallelize. In objective functional, the integrals with respect to time can be approximated by the composition trapezoidal rules.

Based on work presented in [7, 81], in next section we develop the multiple-relaxation-time lattice Boltzmann method for (1)-(2) and show that the method can recover (1) through multi-scaling Chapman-Enskog procedure.

3.1 Multiple-relaxation-time LBM for general coupled anisotropic convection-diffusion chemotaxis-type equations

In order to develop multiple-relaxation-time LBM for (1) and show that this method can recover (1), let us introduce the evolution equation of the LBM (see [20]) for solving a reaction-diffusion equation, with the macroscopic variable $\Theta$,

[TABLE]

with

[TABLE]

where $\displaystyle{\cal J}(\mathbf{x},t;\mathbf{e})=Q^{col}(\varphi(\mathbf{x},t;\mathbf{e})-\varphi^{eq}(\mathbf{x},t;\mathbf{e}))+P(\mathbf{x},t;\mathbf{e})$, $\varphi(\mathbf{x},t;\mathbf{e})$ is the distribution function of the single particle moving with velocity $\mathbf{e}$ at position $\mathbf{x}$ and time $t$, $P$ is the distribution type function of particle of macroscopic external force $\Phi$ moving with velocity $\mathbf{e}$, $Q^{col}$ is a Bhatnagar-Gross-Krook (BGK) type collision operator (see [18]) and $\varphi^{eq}$ is the Maxwell-Boltzmann equilibrium distribution function.

LBM leads us to approximate (27) to recover reaction-diffusion equation (1) with multi-scaling Chapman-Enskog expansion (see e.g., [24]). In method of Chapman-Enskog, the distribution function $\varphi$ is expanded as : $\varphi=\varphi^{(0)}+\varepsilon\varphi^{(1)}+\varepsilon^{2}\varphi^{(2)}+....=\sum_{k=0}^{\infty}\varepsilon^{k}\varphi^{(k)}$, which can be regarded as a power series in a small variable $\varepsilon$ or an expression that keeps track of the relative orders of magnitude of the different terms through scaling parameter $\varepsilon$ (the so-called Knudsen number). The distribution functions $\varphi^{(0)}$, $\varphi^{(1)}$ and $\varphi^{(2)}$ represent the zero, first and second approximation to the distribution function $\varphi$, and so on.

The numerical solution of Boltzmann equation (27) requires to discretize the configuration spaces, velocities and time. For that, we discretize $\overline{\Omega}\times(0,T)$ in time and space. For a given space discretization parameters $h$ and lattice time step size $\tau={T}/{N}$, $N\in\mbox{IN}^{*}$, the domain $\overline{\Omega}=\Omega\cup\partial\Omega$ is approximated by a uniform, regular grid with spacing $h$ (which is commonly called the lattice), and the interval $[0,T]$ is split by the equidistant partitioning using the following points $0=t_{0}\leq t_{1}\leq\cdots\leq t_{N}=T$ with $t_{n}=n\tau$, for $n=0,\ldots,N$. The steps $h$ and $\tau$ are chosen sufficiently small to guarantee both the time accuracy and convergence of solution.

Finally, we define the discrete time interval $\daleth_{\tau}=\{t\in[0,T]:t=n\tau,~{}~{}n\in\mbox{IN}\}$, streaming lattice speed $c=h/\tau$ and lattice sound speed $C_{s}$ by $3C_{s}^{2}=c^{2}$.

We note $(\mathbf{e}_{i})_{i=0,q-1}$ the discrete set of $q$ admissible particle velocities (the magnitude of each $\mathbf{e}_{i}=(e_{i}^{(k)})_{k=1,m}\in\mbox{IR}^{m}$ depends on speed $c$) and we assume that for each node x of lattice ${\cal L}_{h}$, and each velocity $\mathbf{e}_{i}$, the point $\mathbf{x}+\mathbf{e}_{i}\tau$ is also a node of the lattice ${\cal L}_{h}$. For the latter, we use the so-called $DmQq$ LBM scheme (i.e., m-dimensional and $q$ velocity vectors). First, we consider the finite discrete-velocity system of Boltzmann equation with finite discrete velocity $\mathbf{e}_{i}$, for $i=0,q-1$ (by discretizing the velocity space by discrete set of microscopic velocities), for system (1) (for $\Theta=u,v$ or $w$)

[TABLE]

where, for each particle on lattice, we associate discrete functions, for $i=0,q-1$ and $\Theta=u,v$ or $w$, $\varphi_{\Theta,i}(\mathbf{x},t)=\varphi_{\Theta}(\mathbf{x},t;\mathbf{e}_{i})$, $\varphi^{eq}_{\Theta,i}(\mathbf{x},t)=\varphi_{\Theta}^{eq}(\mathbf{x},t;\mathbf{e}_{i})$ and $\displaystyle{\cal J}_{\Theta,i}=\sum_{j}\hat{\Lambda}^{(\Theta)}_{ij}(\varphi_{\Theta,j}-\varphi^{eq}_{\Theta,j})+\Phi_{\Theta,i}+{\cal S}_{\Theta,i}$, where $\Phi_{\Theta,i}$ (resp. ${\cal S}_{\Theta,i}$) is the discrete operator of $\Phi_{\Theta}$ (resp. the discrete operator corresponding to cross-diffusion operators), which describes the probability of streaming in one particular direction and $\hat{\bf\Lambda}_{\Theta}$, defined by $\hat{\bf\Lambda}_{\Theta}(\mathbf{x},t)=(\hat{\Lambda}^{(\Theta)}_{ij}(\mathbf{x},t))_{i,j=0,q-1}$, is the collision matrix (see e.g., [27]).

The aim of the so-called multiple-relaxation-time LBM scheme is to precise the distribution function of particle $(\varphi_{\Theta,i}(\mathbf{x},t))_{i=0,q-1}$ for $\mathbf{x}\in{\cal L}_{h}$ and $t\in\daleth_{\tau}$ by solving the discretization of Boltzmann equation (29) in two steps (for each time): collision and streaming processes. Then, by using a second order time integration scheme to approximate (29) (by limiting physical space to a lattice and velocity space to discrete set of microscopic velocities), we can obtain

[TABLE]

where the discrete cross-diffusion operators are defined by

[TABLE]

and the discrete source terms are defined by ($\Theta=u,v$ or $w$)

[TABLE]

where the differential tensors $A_{\Theta,k}$ and ${\vec{V}_{\Theta}}$ (for $\Theta=u,v$ or $w$, $k=1,2$) will be determined below and $\mathfrak{h}_{\Theta,k}$ is the identity operator, $\kappa_{x}$ and $\kappa_{t}$ are two parameters to be determined later on and $\Pi_{\Theta}=\tau{\cal M}\hat{\bf\Lambda}_{\Theta}{\cal M}^{-1}$ is the relaxation matrix in velocity space. ${\cal M}$ is a transformation matrix (constant) which is constructed from discrete velocities via the Gram-Schmidt orthogonalization procedure. Then the collision process is carried out in moment space and ${\cal M}$ can be used to project $\varphi_{i}$ and $\varphi^{eq}_{\Theta,i}$ in discrete velocity space onto macroscopic variables in moment space as

[TABLE]

where ${\bf F}_{\Theta}=(\varphi_{\Theta,i})_{i=0,q-1}$ and ${\bf F}^{eq}_{\Theta}=(\varphi^{eq}_{\Theta,i})_{i=0,q-1}$. Therefore, the collision process in moment space can be expressed as

[TABLE]

Introduce also the following family of matrix as (for $k=1,m$)

[TABLE]

Then the equilibrium distribution functions $\varphi^{eq}_{\Theta,i}(\mathbf{x},t)$, for $\Theta=u,v$ or $w$, of the LBM for system (1), associated with discrete velocity $\mathbf{e}_{i}$ at position $\mathbf{x}$ and time $t$, are defined as (for $i=0,q-1$)

[TABLE]

where ${\bf I}_{d}$ is the unit matrix, $d_{\Theta}>0$ is a parameter related to the diffusion tensor $D_{\Theta}$, the matrix operator $(\vec{\nu}\otimes\vec{\mu})$, for $\vec{\nu}=(\nu_{i})_{i=1,m}$ and $\vec{\mu}=(\mu_{i})_{i=1,m}$, is defined by $(\vec{\nu}\otimes\vec{\mu})_{ij}=\nu_{i}\mu_{j}$, for $i,j=1,m$, the contraction product between two matrices $A=(a_{ij})_{i,j=1,m}$ and $B=(b_{ij})_{i,j=1,m}$ is defined by $A:B=\sum_{i,j}a_{ij}b_{ij}$ and $\omega_{i}$ is a weight coefficient in direction $\mathbf{e}_{i}$. Finally the second-order differential tensor $\mathfrak{C}_{\Theta}$ is such that (for all $({\bf x},t)$ and $\Theta$)

[TABLE]

and the matrix $\mho_{i}$ is given by

[TABLE]

Remark 8**.**

For different lattice models $DmQq$, the velocities $\mathbf{e}_{i}$ and weight coefficients $\omega_{i}$, for $i=0,\dots,q$ can be defined as (for example)

For $D2Q9$:

$\displaystyle\displaystyle\mathbf{e}_{0}=c\begin{pmatrix}0\\ 0\end{pmatrix},~{}\mathbf{e}_{1}=c\begin{pmatrix}1\\ 0\end{pmatrix},~{}\mathbf{e}_{2}=c\begin{pmatrix}0\\ 1\end{pmatrix},~{}\mathbf{e}_{3}=-\mathbf{e}_{1}~{}\mathbf{e}_{4}=-\mathbf{e}_{2},\displaystyle\displaystyle\mathbf{e}_{5}=c\begin{pmatrix}1\\ 1\end{pmatrix},~{}\mathbf{e}_{6}=c\begin{pmatrix}-1\\ 1\end{pmatrix},$

$~{}\mathbf{e}_{7}=-\mathbf{e}_{5},~{}\mathbf{e}_{8}=-\mathbf{e}_{6}$* and *

[TABLE]

For $D3Q19$:

$\displaystyle\displaystyle\mathbf{e}_{0}=c\begin{pmatrix}0\\ 0\\ 0\end{pmatrix},~{}\mathbf{e}_{1}=c\begin{pmatrix}1\\ 0\\ 0\end{pmatrix},~{}\mathbf{e}_{2}=-\mathbf{e}_{1},~{}\mathbf{e}_{3}=c\begin{pmatrix}0\\ 1\\ 0\end{pmatrix},~{}\mathbf{e}_{4}=-\mathbf{e}_{3},~{}\displaystyle\mathbf{e}_{5}=c\begin{pmatrix}0\\ 0\\ 1\end{pmatrix},~{}\mathbf{e}_{6}=-\mathbf{e}_{5},$

$\mathbf{e}_{7}=c\begin{pmatrix}1\\ 1\\ 0\end{pmatrix},~{}\mathbf{e}_{8}=-\mathbf{e}_{7},~{}\mathbf{e}_{9}=c\begin{pmatrix}1\\ -1\\ 0\end{pmatrix},~{}\mathbf{e}_{10}=-\mathbf{e}_{9},~{}\mathbf{e}_{11}=c\begin{pmatrix}1\\ 0\\ 1\end{pmatrix},~{}\mathbf{e}_{12}=-\mathbf{e}_{11},$

$~{}\mathbf{e}_{13}=c\begin{pmatrix}-1\\ 0\\ 1\end{pmatrix},~{}\mathbf{e}_{14}=-\mathbf{e}_{13},~{}\mathbf{e}_{15}=c\begin{pmatrix}0\\ 1\\ 1\end{pmatrix},~{}\mathbf{e}_{16}=-\mathbf{e}_{15},~{}\mathbf{e}_{17}=c\begin{pmatrix}0\\ 1\\ -1\end{pmatrix},~{}\mathbf{e}_{18}=-\mathbf{e}_{17}$* and*

[TABLE]

$\Box$**

Finally, in order to derive precisely the macroscopic variable $\Theta$ in terms of microscopic particle distribution function governed by lattice Boltzmann equation (30), the following conditions also need to be satisfied

[TABLE]

Remark 9**.**

During streaming and collision processes, in order to satisfy boundary conditions, boundary nodes need special treatments on distribution functions, which are essential to stability and accuracy of the method. We can use extrapolation methods (see e.g., [87]) and modified bounce-back methods (see e.g.,[85, 82]). $\Box$

In the following proposition, we prove formally that governing macroscopic equation (1) is recovered correctly from the developed multiple-relaxation-time LBM model.

Proposition 1**.**

The LB equations (30) with (31) and (32) approximate solutions to continuous macroscopic system (1) via Chapman Enskog asymptotic expansion with errors of order $O(\varepsilon^{3})$, where $\varepsilon$ is the so-called non-uniformity parameter $\varepsilon$ (which can be interpreted as the Knudsen number that controls the order of the approximation).

Proof: To derive the macroscopic equation from developed multiple-relaxation-time lattice BGK model, the Chapman-Enskog expansion is applied under the assumption of small Knudsen number $\varepsilon$ (which can be a ratio between a characteristic length $L$ and a particular mean free path $l$) to determine equilibrium functions, source terms and link functions. Each function is decomposed around its equilibrium state with different scale of perturbations. Then, the difference between associated LBE and Taylor expansion of the distribution function lead us to recover the governing macroscopic equation. For every distribution functions $\varphi_{\Theta,i}$, the Chapman-Enskog expansion is applied with the small parameter $\varepsilon$ as:

[TABLE]

where $\varphi^{(k)}_{\Theta,i}$ are assumed to be sufficiently smooth functions in $({\bf x},t)$ and independent of $\varepsilon$, $t_{i}$ is the time scale (for $i=1,2$), ${\bf x}_{1}$ is the space scale, $\nabla_{1}=\nabla_{{\bf x}_{1}}$ is the gradient with respect to ${\bf x}_{1}$ and $\varphi_{\Theta,i}^{(0)}$ is associated to unperturbed state. Correspondingly, macroscopic quantity $\Theta$ is expanded as

[TABLE]

We assume also that every term $\Phi_{\Theta,i}$ and ${\cal S}_{\Theta,i}$ take the form

[TABLE]

The previous expansions are inserted into LBE (30) and Taylor-series expansion222Taylor-series expansion with truncation to the $p$-order derivative terms and remainder term, $p\geq 2$ (depending on the regularity of $\varphi_{\Theta,i}$) can be applied. is applied to $\varphi_{\Theta,i}(\mathbf{x}+\mathbf{e}_{i}\tau,t+\tau)$ ($\varphi_{i}$ moving with velocity $\mathbf{e}_{i}$) at point $({\bf x},t)$ as

[TABLE]

where

[TABLE]

with $\displaystyle(\mathbf{e}_{i}\cdot\nabla)^{n}=(\mathbf{e}_{i}\cdot\nabla)\left(\mathbf{e}_{i}\cdot\nabla\right)^{n-1}.$

So, by using the expressions (42)-(45), the expansion (46) and LBE (30) become ($\forall i=0,\dots,q-1$)

[TABLE]

and

[TABLE]

By calculating the difference between relations (47) and (48), we can obtain

[TABLE]

which can be rewritten as follows (according to different scales of $\varepsilon$)

[TABLE]

where probability density functions $\Psi_{\Theta,i}^{(n)}$ regroup all terms of (49) in order of $\varepsilon^{n}$, $n=0,...,$ respectively (which are partial differential equations for $n\neq 0$)

[TABLE]

To treat second order terms, we rewrite $\displaystyle\frac{\tau^{2}}{2}\left(\frac{\partial}{\partial t_{1}}+\mathbf{e}_{i}\cdot\nabla_{1}\right)^{2}\varphi_{\Theta,i}^{(0)}$ with using (52) as:

[TABLE]

and then $\Psi_{i}^{(2)}$ becomes

[TABLE]

If we rewrite the previous equations (51), (52) and (54) in a vector form and we multiply on both sides of them, we can deduce by retaining terms up to $O(\varepsilon^{3})$

[TABLE]

where $({\cal M}^{-1}{\cal H}^{(n)}_{\Theta})_{i}=\Phi_{\Theta,i}^{(n)}$ and $({\cal M}^{-1}{\cal L}^{(n)}_{\Theta})_{i}={\cal S}_{\Theta,i}^{(n)}$, for $n=1,...$, and ${\cal E}\nabla_{1}=\sum_{k=1,m}{\cal E}_{k}\nabla_{1,k}$.

In order to facilitate the presentation we consider, in the rest of this section, the popular two-dimensional $D2Q9$ lattice model (see Fig.1), which involves 9 velocity vectors as an example (the arguments developed below can be extended to a 3-dimensional $D3Q19$ lattice without any substantial difficulties, because this last has essentially the same procedure as that for the $D2Q9$ one).

According to the expression of microscopic velocities $\mathbf{e}_{i}$, we prove easily that

[TABLE]

The transportation matrix ${\cal M}$ can be written as ${\cal M}=E_{0}{\cal M}_{0}$ with ${\cal M}_{0}$ given by (see e.g. [48])

[TABLE]

and $E_{0}=\mathrm{diag}(1,c^{2},c^{4},c,c^{3},c,c^{3},c^{2},c^{2})$ a diagonal matrix. We can deduce that ${\bf m}^{(eq)}_{\Theta}=(m^{(eq)}_{\Theta,i})_{i=0,8}$ is given by

[TABLE]

where

[TABLE]

In order to handle anisotropic diffusion problems, the relaxation matrix $\Pi_{\Theta}$ in the moment space is not a diagonal matrix, as in the isotropic case, but can be expressed as a block-diagonal invertible matrices : $\Pi_{\Theta}=\mathrm{diag}(\Lambda_{\Theta,1},\Lambda_{\Theta,2},\Lambda_{\Theta,3})$, with

[TABLE]

where the diagonal elements $s_{\Theta,j}$ ($j=0,2$), $k_{\Theta,11}$, $s_{\Theta,3}$ and $k_{\Theta,22}$, $s_{\Theta,j}$ ($j=4,6$) are the relaxation parameters corresponding to jth moment ${m}_{j}$, $j=0,8$, respectively, and the terms $k_{\Theta,12}$ and $k_{\Theta,21}$ correspond to the rotation of the principal axis of anisotropic diffusion (see e.g.[84]).

Remark 10**.**

We prove easily that $E_{0}\Pi_{\Theta}=\Pi_{\Theta}E_{0}$ and then ${\cal M}^{-1}\Pi_{\Theta}{\cal M}={\cal M}_{0}^{-1}\Pi_{\Theta}{\cal M}_{0}$. Then from (30)

[TABLE]

Based on the second relation of (55) we can deduce that333we consider only the first, fourth and sixth components which are useful in order to cover the state equation. (according to (41))

[TABLE]

where the invertible matrix ${\bf K}_{\Theta}$ is given by

[TABLE]

and on the third relation of (55) (i.e. the second order relation in $\varepsilon$) we can deduce (by considering only the first which is useful in order to cover the state equation)

[TABLE]

Then, from the second equation of (59), the equation (60) becomes

[TABLE]

and then

[TABLE]

Calculate now the value $\frac{\partial\mathfrak{T}_{\Theta}}{\partial t_{1}}+div_{1}({\bf B}_{\Theta})$. From (56), we can deduce

[TABLE]

According to (37), we can deduce

[TABLE]

Since $div_{1}(\mathfrak{T}_{\Theta}(\mathbf{x},t;\Theta(\mathbf{x},t))=div_{1}(\mathfrak{T}_{\Theta})(\mathbf{x},t;\Theta)+(\frac{\partial\mathfrak{T}_{\Theta}}{\partial\Theta}(\mathbf{x},t;\Theta))\nabla_{1}\Theta,$ we can deduce that

[TABLE]

Denotes now by $DV$ the following operator

[TABLE]

Then (from the first equation of (59))

[TABLE]

Consequently, equation (61) becomes

[TABLE]

By taking $\kappa_{t}=1$, $\vec{V}_{\Theta}=({\bf K}_{\Theta}^{-1}-\frac{\kappa_{x}}{2}{\bf I}_{d})^{-1}({\bf K}_{\Theta}^{-1}-\frac{1}{2}{\bf I}_{d})\frac{\partial\mathfrak{T}_{\Theta}}{\partial\Theta}$, $A_{\Theta,k}=\frac{1}{\tau}{\bf K}_{\Theta}D_{\Theta,k}$ (for $k=1,2$), and the matrix ${\bf K}_{\Theta}$ such the diffusion tensors $D_{\Theta}=d_{\Theta}C_{s}^{2}\tau({\bf K}_{\Theta}^{-1}-\frac{1}{2}{\bf I}_{d})$, we obtain

[TABLE]

In equation (65) there exists the deviation term $\displaystyle div_{1}\left(({\bf K}_{\Theta}^{-1}-\frac{1}{2}{\bf I}_{d})DV(\mathbf{x},t;\Theta)\right).$ If this term is null (e.g., $\mathfrak{T}_{\Theta}$ is independent on $(\mathbf{x},t)$ or $\mathfrak{T}_{\Theta}$ is independent on times with $div(\mathfrak{C}_{\Theta})(\mathbf{x},t;\Theta)=0$ and $div(\mathfrak{T}_{\Theta})(\mathbf{x},t;\Theta)=0$) or is negligible compared to the other terms, then (65) becomes

[TABLE]

By combining the first equation of (59) (multiplying by $\varepsilon$) and (66) (multiplying by $\varepsilon^{2}$), we can recover equation (1). This completes the proof of the recovery of (1) from the developed LBM method. $\Box$

3.2 Discrete adjoint-based optimization

Since the state ${\bf u}$ depends on the microscopic distribution function ${\cal G}=(G_{i})_{i=0,q-1}$ with $G_{i}=(\varphi_{u,i},\varphi_{v,i},\varphi_{w,i})$, the continuous lattice Boltzmann type equations (LBE) with force (for $i=0,q-1$), which is associated to the three equations of system (1), can be reformulated as follows (since the gradient of state functions can be computed locally as function of distributions functions)

[TABLE]

where ${\mathfrak{F}}_{\Theta,i}$ (resp. $\mathfrak{S}_{\Theta,i}$), for $\Theta=u,v$ or $w$ are corresponding to $\Phi_{\Theta,i}$ (resp. to ${\cal S}_{\Theta,i}$).

The corresponding discrete form of (67) is (at $(\mathbf{x},t)\in{\cal L}_{h}\times\daleth_{\tau}$)

[TABLE]

where, for $\Theta=u,v$ or $w$ and for $i=0,q-1$, the equilibrium distribution function $\varphi^{eq}_{\Theta,i}(\mathbf{x},t)$, associated with the discrete velocity $\mathbf{e}_{i}$ at position $\mathbf{x}$ and time $t$, is defined by (36) and the discrete source terms ${\mathfrak{F}}_{\Theta,i}$ and ${\mathfrak{S}}_{\Theta,i}$ are defined, respectively, by (32) and (31). After each time step, we obtain new particular repartitions which lead us to recover the macroscopic state ${\bf u}$ with the following summation of microscopic particle distribution functions $G_{i}$

[TABLE]

Remark 11**.**

In order to calculate $G_{i}$ by (67), for a given ${\bf f}$, it is necessary to approach the third terms of the right hand side of system (67) which depend on the parameter $\kappa_{x}$. For simplicity, in the case when $\kappa_{x}=0$ or 1, we can use an explicit characteristic method to approximate these terms i.e.,

[TABLE]

We can also use an implicit characteristic method (see e.g., [7, 88]). $\Box$

Remark 12**.**

Since the developed method is explicit time-stepping scheme, then it is necessary to impose some conditions on discretization parameters, in the spirit of the CFL condition (as e.g. of type $\tau\leq C_{CFL}~{}h^{2}$), which guarantee the convergence of the method. $\Box$

According to the previous development and using similar approach as in [7], we can now formulate and develop the adjoint-based multiple-relaxation-time lattice Boltzmann. The method consists in using the adjoint of the multiple-relaxation-time LBE (67) and its time discretization form (68).

3.2.1 Adjoint multiple-relaxation-time lattice Boltzmann method and algorithm

Since the cost functional ${\cal J}$ depends on ${\bf f}$ and ${\bf u}=(u,v,w)$, and from LBM approximation of solution $\displaystyle{\bf u}=\sum_{i=0}^{q-1}G_{i}$, with $G_{i}=(\varphi_{u,i},\varphi_{v,i},\varphi_{w,i})$, then ${\cal J}$ can be written as the form

[TABLE]

where ${\cal G}=(G_{i})_{i=0,q-1}$ and $G_{i}$ is the solution of the continuous LBE with force (67) (for $i=0,q-1$). Since the distribution function ${\cal G}$ depends on control variable ${\bf f}$ implicitly, then a variation $\delta{\bf f}$ of ${\bf f}$ causes a variation $\delta{\cal G}$ of ${\cal G}$. Consequently a variation $\delta{\cal J}$ of cost functional ${\cal J}$ can be given as

[TABLE]

We can now show the gradient of ${\cal J}$ by introducing intermediate sensitivity and adjoint LB equations.

Theorem 3**.**

The components of the gradient of ${\cal J}$ with respect to control variable ${\bf f}$ is

[TABLE]

with $(\psi_{u,i},\psi_{v,i},\psi_{w,i})$ the solution of the following adjoint LB equations

[TABLE]

where adjoint equilibrium distribution functions $\psi^{eq}_{\Theta,i}$ (for $\Theta=u$, $v$ or $w$ and $i=0,q-1$)

[TABLE]

and

[TABLE]

Proof: Since $G_{i}=(\varphi_{u,i},\varphi_{v,i},\varphi_{w,i})$ is a solution of (67), then the variation $\delta G_{i}=(\delta\varphi_{u,i},\delta\varphi_{v,i},\delta\varphi_{w,i})$ is a solution of the following sensitivity lattice Boltzmann equation

[TABLE]

Let $H_{i}=(\psi_{u,i},\psi_{v,i},\psi_{w,i})$ be sufficiently regular such that $\displaystyle H_{i}(T)=(\frac{\partial{\cal I}_{0}}{\partial\varphi_{u,i}},\frac{\partial{\cal I}_{0}}{\partial\varphi_{v,i}},\frac{\partial{\cal I}_{0}}{\partial\varphi_{w,i}})$. Multiplying (76) by $H_{i}$, summing on $i$ and integrating over $\Omega\times(0,T)$ this gives (since $\delta G_{i}(t=0)=(0,0,0)$)

[TABLE]

[TABLE]

and

[TABLE]

By setting for $\Theta=u,v$ or $w$, $({\bf e}_{i}.{\bf n})\psi_{\Theta,i}=0$ on $\partial\Omega$, and noting by

[TABLE]

we can deduce

[TABLE]

Introduce the following adjoint equilibrium distribution functions ($\text{for~{}}\Theta=u,~{}v~{}\text{or~{}}w$ and $i=0,q-1$)

[TABLE]

Otherwise

[TABLE]

where $\Psi_{\Theta}=(\psi_{\Theta,i})_{i}$, $\Psi^{eq}_{\Theta}=(\psi^{eq}_{\Theta,i})_{i}$, $\tilde{\bf m}_{\Theta}=\hat{\bf\Lambda}^{*}_{\Theta}\Psi_{\Theta}$, $\tilde{\bf m}^{eq}_{\Theta}=\hat{\bf\Lambda}^{*}_{\Theta}\Psi^{eq}_{\Theta}$. Then (78) becomes

[TABLE]

Assuming that $(\psi_{u,i},\psi_{v,i},\psi_{w,i})$ is a solution of the following adjoint lattice Boltzmann equation

[TABLE]

where (for $\Theta=u$, $v$ or $w$) : ${\mathfrak{R}}_{\Theta,k}={\mathfrak{F}}_{\Theta,k}+{\mathfrak{S}}_{\Theta,k}.$

By adding the first part, second part and third part of (81) we have (according to (82))

[TABLE]

Then (71) becomes

[TABLE]

Consequently the components of the gradient of ${\cal J}$ with respect to variable ${\bf f}$ is then

[TABLE]

This complete the proof. $\Box$

To close the adjoint lattice Boltzmann equation (82) and sensitivity lattice Boltzmann equation (76), it is necessary to evaluate the adjoint distribution function $\psi^{eq}_{\Theta,i}$ (fo $\Theta=u,v~{}\text{or~{}}w$). For this, it is sufficient (according to (80) and (77)) to derive $\displaystyle\frac{\partial\varphi^{eq}_{\Theta,j}}{\partial\varphi_{\Theta,i}}$. According to (41) and (36), we can compute each $\displaystyle\frac{\partial\varphi^{eq}_{\Theta,j}}{\partial\varphi_{\Theta,i}}$ as follows:

[TABLE]

The expression of $\aleph_{j,\Theta}$ can be derived easily from the expressions of $\frac{\partial\mathfrak{C}_{\Theta}}{\partial\Theta}$ and $\mho_{j}$. Consequently, from (80), the adjoint equilibrium function $\tilde{\bf m}^{eq}_{\Theta}$ satisfies

[TABLE]

where $\displaystyle{\cal T}_{\Theta}=\sum_{j}\omega_{j}\tilde{m}_{j}\aleph_{j,\Theta}$.

We can now present numerical algorithm to solve the optimal control problem. The continuous adjoint LBE equation (73) in time and space with force can be written as follows

[TABLE]

where

[TABLE]

This problem is similar the continuous LBE in time and space with force (67). The main difference of these systems is their streaming steps. In LBE (67), the distributions $\varphi_{\Theta,i}$ for $\Theta=u,v$ and $w$ in a lattice node streams to the adjacent node along $i$ lattice line with ${\bf e}_{i}$ as forward in time. Conversely, in adjoint LBE (88), the adjoint distributions $\psi_{\Theta,i}$ for $\Theta=u,v$ and $w$ in a lattice node streams to the adjacent node along $i$ lattice line with $-{\bf e}_{i}$ as backward in time. Consequently we can derive, in same way as to obtain (68), the discrete adjoint LBE corresponding to (88) (when ${\bf x}\in{\cal L}_{h}$ and $t=t_{n}$, n=N,N-1,…,1 (backward in time)).

In the sequel we denote by $G_{i}^{(h,n)}=(\varphi^{h,n}_{u,i},\varphi^{h,n}_{v,i},\varphi^{h,n}_{w,i})$ (respectively $H_{i}^{(h,n)}=(\psi^{h,n}_{u,i},\psi^{h,n}_{v,i},\psi^{h,n}_{w,i})$) the solution of the discrete LBE (68) (respectively of the discrete adjoint LBE corresponding to (88)), ${\cal G}^{h,n}=(G_{i}^{(h,n)})_{i=0,q-1}$, ${\cal H}^{h,n}=(H_{i}^{(h,n)})_{i=0,q-1}$ and ${\bf u}^{h,n}=(u^{h,n},v^{h,n},w^{h,n})=\sum_{i}G_{i}^{(h,n)}$.

In the objective functional ${\cal J}$, the integrals with respect to time can be approximated by the composite trapezoidal rule, where we adopt the point $t_{n}$ from the time discretization of the state and adjoint systems, as follows (for ${\bf f}^{h}\in{\cal K}^{h}_{ad}$)

[TABLE]

where finite dimensional subspace ${\cal K}^{h}_{ad}$ is the internal approximation of ${\cal K}_{ad}$. Here $\theta_{0}=\theta_{N}={1}/{2}$ and $\theta_{n}=1$ for $n\neq 0$ and $n\neq N$. The gradient of ${\cal J}_{h}$ at point ${\bf f}^{h}=({\bf f}_{1}^{h},{\bf f}_{2}^{h},{\bf f}_{3}^{h})$ in direction of ${\bf g}^{h}=({\bf g}_{1}^{h},{\bf g}_{2}^{h},{\bf g}_{3}^{h})$ is given by (analogous to the one of ${\cal J}$)

[TABLE]

According to the previous discrete formula (90) and (91) we can now present the optimization process using the adjoint LBE in order to solve the lattice Boltzmann discrete optimal control problem (in ${\cal K}^{h}_{ad}$).

For $k=0,...$ (iteration index) we denote by $X_{k}$ the numerical approximation of all functions $X$ (involved in algorithm) at $k$th iteration of algorithm.

Initialization : ${\bf f}^{h}_{k}$ (given) 2. 2.

Compute primal problem (1) with source term ${\bf f}^{h}_{k}$ by using LBE (68) (with initial conditions), gives the distribution function ${\cal G}^{h,n}_{k}$ and primal solution ${\bf u}_{k}^{h,n}$, for $n=0,...,N$ (forward in time). 3. 3.

Compute discrete adjoint LBE corresponding to (88) (with final conditions), based on ${\cal G}^{h}_{k}$, ${\bf u}_{k}^{h}$ and ${\bf f}^{h}_{k}$, gives the adjoint distribution function ${\cal H}_{k}^{h,n}$, for $n=N,...,0$ (backward in time). 4. 4.

Compute ${\cal J}^{\prime}_{h}({\bf f}_{k}^{h})$ from (91) (based on ${\cal H}_{k}^{h}$, ${\cal G}^{h}_{k}$, ${\bf u}_{k}^{h}$ and ${\bf f}^{h}_{k}$). 5. 5.

Perform an appropriate optimization algorithm (Gradient, Conjugate Gradient,…), gives the new control variable ${\bf f}^{h}_{k+1}$. 6. 6.

If gradient of ${\cal J}_{h}$ is sufficiently small (convergence) then end; else set k:=k+1 and goto step 2. Discrete approximation of optimal solution $({\bf u}^{*},{\bf f}^{*})$ is $({\bf u}_{k}^{h},{\bf f}^{h}_{k})$.

3.2.2 Continuous adjoint problem and multiple-relaxation-time LBM

The treated control problem, can be solved also numerically from the discretization of the continuous adjoint problem (19) by using the proposed multiple-relaxation-time lattice Boltzmann method. Indeed, the adjoint problem (19) (with final conditions), which is backward in time, can be transformed into an initial-boundary value problem by reversing the sense of time i.e., $t:=T-t$ and the derived system is similar to (1) (with initial conditions). The algorithm can be formulated as follows.

For $k=0,...$ (the iteration index) we denote by $X_{k}$ the numerical approximation of all functions $X$ (involved in the algorithm) at $k$th iteration of the algorithm.

Initialization : ${\bf f}^{h}_{k}$ (given) 2. 2.

Compute primal problem (1) with the source term ${\bf f}^{h}_{k}$ by using LBE (68) (with initial conditions), gives the primal solution ${\bf u}_{k}^{h,n}$, for $n=0,...,N$. 3. 3.

Compute adjoint (13), based on ${\bf u}_{k}^{h}$ and ${\bf f}^{h}_{k}$, gives adjoint solution $\tilde{\bf u}_{k}^{h,n}$, for $n=N,...,0$ (by reversing the sense of time and using the LBE (68)). 4. 4.

Compute ${\cal J}^{\prime}_{h}({\bf f}_{k}^{h})$ from (24) (based on $\tilde{\bf u}_{k}^{h}$, ${\bf u}_{k}^{h}$ and ${\bf f}^{h}_{k}$). 5. 5.

Perform an appropriate optimization algorithm (Gradient, Conjugate Gradient,…), gives the new control variable ${\bf f}^{h}_{k+1}$. 6. 6.

If gradient of ${\cal J}_{h}$ is sufficiently small (convergence) then end; else set k:=k+1 and goto step 2. Discrete approximation of optimal solution $({\bf u}^{*},{\bf f}^{*})$ is $({\bf u}_{k}^{h},{\bf f}^{h}_{k})$.

4 Conclusions

This paper presents a complete formulation for the continuous and discrete adjoint approach for automatic control of some complex real-life systems, which appear in diverse biochemical, biological and biosocial criminology problems, governed by a class of nonlinear coupled anisotropic convection-diffusion chemotaxis-type system (NCACDCS). The optimal control problem subject to primal system NCACDCS has been formulated within a general framework. First-order necessary optimality conditions are established by using sensitivity and adjoint calculus. These conditions allows easy access to gradients of objective functions (even with a large number of estimate parameters and/or constraints), that is very useful in numerical optimization algorithms. A general adjoint-based sensitivity methodology for these inverse optimization problems using a general multiple-relaxation-time Lattice Boltzmann (MRT) is developed and an adjoint multiple-relaxation-time lattice Boltzmann model (AMRT), which is found to be as simple as MRT model with also highly-efficient parallel nature, is derived.

Despite its apparent complexity, the developed discrete algorithm is quite easy to implement with low implementation cost (even in the case of large time horizons $T$ or fine space-time meshes). Due to the local nature of the adjoint and primal collision operators and the fact that nearest neighbor nodes only interact during the streaming stages (allowing for straightforward and efficient parallelization and use of GPU computing), this proposed method is a reliable and efficient approach to solve this type of practical real-life problems (without forgetting that the ability of LBM to approximate complicated geometries with simple algorithmic modifications). Moreover, to help reduce further the computational complexity, CPU and memory costs, we can coupled the developed method and Proper Orthogonal Decomposition Reduced Order Model (POD).

5 Conflict of interest

Author declared that there are no conflicts of interests