1. Introduction

The epidemic problem man-environment in [3] is stated as optimal control problem. The functional consists of three members and the state equations are governed by two PDE: parabolic linear equation and nonlinear first order equation with feedback operator on the boundary.

Thus optimal control problem (P) is to minimize

$J(u_{1}, u_{2}, w) =\int _{0}^{T}\int _{\Omega }F(u_{2}(t,x)dxdt +\int _{0}^{T}\int _{ \partial \Omega }h(w(t, x))dxdt +\int _{\Omega }l(u_{2}(T, x))dx \label{fu}$

(1)

over all $(u_{1}, u_{2}, w)$ subject to state system

$\frac{ \partial u_{1}}{ \partial t} - \Delta u_{1} +a_{11} u_{1} =0, \text{in }\text{}\text{}\text{}\text{}Q =(0,T) \times \Omega, \text{ } \label{L}$

(2)

$\frac{ \partial u_{2}}{ \partial t} +a_{22}u_{2} -g(u_{1}) =0, \text{in }Q,\text{ } \label{u2}$

(3)

$u_{1}(0, x) =u_{1}^{0}(x), \text{ }u_{2}\text{ }(0, x) =u_{2}^{0}(x)\text{ for }x \in \Omega, \text{ } \label{inc}$

(4)

$\frac{ \partial u_{1}}{ \partial \nu } +\alpha u_{1} =K \ast u_{2} =\int _{\Omega }K(t, x, \sigma )u_{2}(t, x)dx, \text{ }(t, \sigma ) \in \sum _{1} =(0, T) \times \Gamma _{1}, \label{bon}$

(5)

$\frac{ \partial u_{1}}{ \partial \nu } =0\text{ in }\sum _{2} =(0, T) \times \Gamma _{2}. \label{bon2}$

(6)

where $\Omega$ is a bounded and open subset of $R^{2}$ with a sufficiently smooth boundary $\Gamma =\Gamma _{1} \cup \Gamma _{2}$, $a_{11}, a_{22}$ and $\alpha$ are positive constants, and

$K(t, x, \sigma, w) =\sum\limits_{i =1}^{N}w_{i}(t,\sigma )K_{i}(x, \sigma )\text{ for }t \in [0, T], \text{ }x \in \Omega, \text{}\text{ }\sigma \in \Gamma _{1}, \label{K}$

(7)

$K_{i} \in L^{\infty }(\Omega \times \Gamma _{1}),$ $w_{i} \in L^{\infty }((0, T) \times \Gamma _{1})$ for $i =1,..., N.$ We set $w(t,\sigma ) =(w_{1}(t,\sigma ),..., w_{N}(t,\sigma )),$ it is a control function, control on the boundary$.$ We assume $F, g, l$ to be continuous in $R,$ $F, g \geq 0$ and $h :R^{N} \rightarrow \bar{R} =[ -\infty, +\infty ]$ is convex, lower semicontinuous. Moreover we assume that there exists a bounded closed subset $M \subset R^{N}$ such that $h(w) = +\infty$ for $w \notin M$, i.e. we assume that control $w \in M$. The following particular case is considered as important in the model problem:

$h(w) =\sum\limits_{i =1}^{N}h_{i}(w_{i}), \label{h}$

(8)

where

$h_{i}(r) =\left \{\begin{array}{cc}\lambda /r^{2} & \text{if }0<r \leq a, \\ +\infty & \text{otherwise.}\,\end{array}\right. \label{r=1}$

(9)

The following $K$

$K(t, x, \sigma, w) =w(t)K(\sigma ), \text{ }t \in [0, T], \text{}\text{}\text{ }\sigma \in \Gamma _{1}, \label{K1}$

(10)

with $w_{i} =w,$ $i =1, 2,..., N$ and

$K(\sigma ) =\sum\limits_{i =1}^{m}a_{i}\chi _{i}(\sigma ),$

(11)

where $\chi _{i}$ is the characteristic function of the interval $[x_{i} -\eta _{i}, x_{i} +\eta ],$ $x_{i} \in \Gamma _{1},$ $i =1, 2,..., m,$ $\eta >0$ is also investigated in [3]$.$ The points $x_{i}$ are related to the treatment of the sewage output. In [1] for that problem existence and necessary optimality conditions, as well two gradient type algorithms are derived. In [4] analytical results are given in support of the well posedness of the problem.The essential point in the convergence of gradient algorithm (using the necessary optimality conditions -Pontryagin maximum principle) is that it starts from arbitrary control function and stop when the difference between two computed controls in next two steps is smaller than given $\varepsilon$. However, we do not know whether the calculated sequence of controls converges to optimal control or the values of the cost functional for those controls converge to optimal value. Moreover, we do not know when to stop the proces in order to get near optimal value i.e. whether for calculated controls the cost of the functional is near optimal value (we do not know it a priori). We need sufficient optimal conditions to grasp such an information. In the literature there is not any optimal control theory of sufficient optimality conditions which can be applied to the above control problem. The main reason is that we deal with the state equations having controls on the boundary. In the next section we develop new dual dynamic programming theory to derive verification theorem -sufficient optimality conditions for problem (1)-(5). However the main advantage of this paper is that we also develop sufficient conditions for $\varepsilon$-optimality i.e. we formulate conditions which allow us to assert that for calculated control (e.g. numerically) we know how far we are from optimal value. Just this approximate theory is fundamental for our numerical algorithm. The control $w =(w_{1},..., w_{N})$ on the boundary we call admissible boundary control and a solution $(u_{1}, u_{2})$ corresponding to it we call admissible state. The set of admissible controls and states we denote by $A d$.

---

2. Dual dynamic programming approach for (1)-(6)

The dual approach to dynamic programming was first introduced in [8] and then developed in several papers to different optimal control problems governed by: elliptic, parabolic and wave equations (see e.g. [5], [11], [10], [9]). The essential point in this dual approach is that we do not deal directly with a value function but with some auxiliary function, defined in a dual set, satisfying dual dynamic equation. The auxiliary function allow us to derive sufficient optimality condition for primal value function. The dual approach has some advantages: we do not need any properties of the value function such as smoothness or convexity. However it has also some disadvantages: the auxiliary function must satisfy a kind of generalised transversality conditions which is a little restrictive. The approach we present here was inspired by PhD thesis [7] of the first author where the model of distortion compensation (elliptic system of equations) was investigated. A new challenge in the control problem (1)-(6) is that the problem under consideration has the control on the boundary and with fixed initial conditions. Therefore we need really to construct a new dual dynamic programming type approach for problem (1)-(6). Thus let us start first with the definition of a dual set. Let $P \subset R^{3 +3}$ be an open set of the variables $(t, x, p) =(t, x, y^{0}, y)$, $y \in R^{2}$, $y^{0} \leq 0$, $(t, x) \in Q =\{(s, z) ;$ $z \in \Omega, 0<s<T\}$. Denote by

$\boldsymbol{Y} =\{(y^{0}, y) =p ;(t, x, p) \in P\}.$

(12)

Denote by $c l P$ the closure of $P$ and by $P_{1}$, $P_{2}$, $clY$ its subsets:

$P_{1} =\{(t, x, p) ;(t, x) \in \sum\limits_{1}\},$

$P_{2} =\{(t, x, p);(t, x) \in (0, T) \times \Gamma _{2}\},$

$cl\boldsymbol{Y} =\left \{(y^{0}, y) =p ;(t, x, p) \in clP \right \}\text{.}$

Let $u$ be a vector of pairs coordinates $(u_{1}, u_{2})$, $\Delta _{x}u =(\Delta _{x}u_{1}, 0),$ $f (t, x, u) =$ $(f_{1} (t, x, u)\text{,}$ $f_{2} (t, x, u))$ where

$f_{1} = -a_{11}u_{1, } \label{f1}$

(13)

$f_{2} =g(u_{1}) -a_{22}u_{2}, \label{f2}$

(14)

$g (t, x, u, w) =(g_{1} (t, x, u, w), g_{2} (t, x, u, w))$ where

$g_{1} = -\alpha u_{1} +\text{}\int _{\Omega }K(t, x, \sigma, w)u_{2}(t, x)dx, \label{g1}$

(15)

$g_{2} =0.$

(16)

Let us introduce an auxiliary function $V :clP \rightarrow R$ belonging to $H^{2} (P)$ (Sobolev space of functions having second weak derivatives) and satisfying "transversality condition":

$V (t, x, p) =y^{0} V_{y^{0}} (t, x, p) +y V_{y} (t, x, p)\;\text{for}\;(t, x, p) \in cl P\;\text{,}\; \label{tr}$

(17)

where $V_{y^{0}}\text{,}$ $V_{y}$ are partial derivative of $V$. Denote by $\boldsymbol{u} :clP \rightarrow R^{2}$ a function of six variables $(t, x, y^{0}, y)\text{.}$ In the sequel we shall assume that

$\boldsymbol{u} (t, x, p) = -V_{y} (t, x, p)\;\text{,}\;(t, x, p) \in clP. \label{x}$

(18)

We shall consider not all admissible controls and corresponding to them admissible states but only those which relate to $\boldsymbol{u}$. To this effect we introduce for given fixed $\xi ( \cdot ) \in (H^{2} (\Omega ))^{2}$ the following set

$Ad_{\boldsymbol{u}} = \left \{(u ( \cdot ), w ( \cdot )) \in Ad ;\;\text{exist}\;p (t, x) =(y^{0}, y (t, x))\;\text{,}\;(t, x) \in Q\;\text{,}\;\right. \label{ad} \\ \;\text{}\;y ( \cdot ) \in (H^{2} (Q))^{2}\;\text{,}\;y^{0} \leq 0\;\text{,}\;(t, x, p (t, x)) \in clP, \nonumber \\ y(0, x) =\xi (x), \;\text{}\;u(t, x) =\boldsymbol{u}(t, x, p (t, x)), \;\text{}\;(t, x) \in clQ\}. \nonumber$

(19)

In fact our optimal control problem we shall study just on the set $A d_{\boldsymbol{u}}$. We consider condition (17) and function (18) just on the set $c lP$. The function $p :Q \rightarrow R^{3}$ we call dual trajectory while $u :Q \rightarrow R^{2}$ we call primal trajectory. The function $\boldsymbol{u}$ from (18) builds a relation between dual and primal trajectory. Next define a dual optimal value $S_{D}$

$S_{D} =\inf\limits_{(u, w) \in A d_{\boldsymbol{u}}} -y^{0}J(u, w). \label{SD}$

(20)

Notice that in spite of that our problem depends on time we cannot perturb it with respect to initial data and time (they are fixed) as it is usually done in classical optimal control theory. This is why a dual dynamic approach to the above problem seems to be the only one possible. Thus let us introduce a dual Hamilton-Jacobi equation in $P$ for our problem:

$V_{t} (t, x, p) - \Delta _{x}V (t, x, p) +y f (t, x, -V_{y} (t, x, p)) +y^{0}F( -V_{y_{2}}(t, x, p)) =0 \label{De}$

(21)

and dual Hamilton-Jacobi type equation on $P_{1}$

$\inf\limits_{w \in M}\{\frac{ \partial V (t, x, p)}{ \partial \nu } +y g (t, x, -V_{y} (t, r, p), w) -y^{0}h(w)\} =0, \label{Be}$

(22)

$\frac{ \partial V (t, x, p)}{ \partial \nu } =0, \;\text{for}\;(t, x, p) \in P_{2}, \label{Be1}$

(23)

$-y^{0}\;\text{}V_{y^{0}} (T, x, p) = -y^{0}l( -V_{y_{2}}(T, x, p)). \label{Be2}$

(24)

We should stress that the notion of dual Hamilton-Jacobi equation appears also in convex optimization (see [2]). However the above dual Hamilton-Jacobi equation is completely different than that in [2]. Our problem is nonconvex and we do not use any tools from convex analysis.

---

3. Sufficient optimality conditions

The dual approach to dynamic programming described in the former section allow us to formulate and to prove a kind of verification theorem ensuring sufficient optimality conditions for our problem (20). We would like to stress that we are working now in dual space $c l P$ and with auxiliary function $V$ defining, by (18) the set $A d_{\boldsymbol{u}}$. Define the set

$\mathcal{P} =\{p (t, x) =(y^{0}, y (t, x)), \;\;(t, x) \in Q ;\;\;(t, x, p (t, x)) \in c l P, \;\; \\ \begin{array}{c}y( \cdot ) \in (H^{2}\text{}(Q))^{2}, \;y(0, x) =\xi (x), \text{exist}\;(u ( \cdot ), w ( \cdot )) \in A d_{\boldsymbol{u}}, \;\,\; \\ u(t, x) = -V_{y} (t, x, p (t, x))\;,\;(t, x) \in Q \cup \Gamma \}.\;\end{array}\,$

Theorem 3.1. Assume that there exists $V \in H^{2} (P)$ satisfying (21)-(24), (17). Let $\bar{\boldsymbol{u}}(t, x, p) = -V_{y} (t, x, p)\text{,}$ $(\bar{u} ( \cdot ), \bar{w} ( \cdot )) \in A d_{\bar{\boldsymbol{u}}}\text{,}$ $(\bar{y}^{0}, \bar{y} (t, x)) =\bar{p} (t, x) \in \mathcal{P}$, $\bar{u}(t, x)$ $= -V_{y} (t, x, \bar{p} (t, x))$, $(t, x) \in Q \cup \Gamma$, $\Gamma =\Gamma _{1} \cup \Gamma _{2_,}$ and

$V_{t} (t, x, \bar{p} (t, x)) - \Delta _{x}V (t, x, \bar{p} (t, x)) \label{Deq} \\ \text{}\text{}\text{} +\bar{y} (t, x) f (t, x, -V_{y} (t, x, \bar{p} (t, x))) +\bar{y}^{0}F( -V_{y_{2}}(t, x, \text{}\text{}p\bar{}(t, x))),$

(25)

$\frac{ \partial V (t, x, \bar{p} (t, x))}{ \partial \nu } +\bar{y} (t, x) g (t, x, -V_{y} (t, x, \bar{p} (t, x)), \bar{w} (t, x)) -\bar{y}^{0}h(\bar{w} (t, x)) =0, \label{Beq}$

(26)

$-y^{0}V_{y^{0}}(T, x, \bar{p} (T, x)) = -y^{0}l( -V_{y_{2}}(T, x, \bar{p}(T, x))), \label{vbo}$

(27)

$\frac{ \partial V (t, x, \bar{p} (t, x))}{ \partial n} =0, \;\text{for}\;(t, x, \bar{p} (t, x)) \in P_{2}. \label{Beq1}$

(28)

Then

$(\bar{u} ( \cdot ), \bar{w} ( \cdot ))$

is an optimal pair with respect to all $(u ( \cdot ), w ( \cdot )) \in Ad_{\bar{\boldsymbol{u}}}$ i.e.

$-\bar{y}^{0}J(u\bar{},w\bar{}) \leq -y^{0}J(u, w).$

Proof. We follow the standard way in proofs of verification theorems. Thus take any $(u( \cdot ), w ( \cdot )) \in Ad_{\bar{\boldsymbol{u}}}$ and corresponding to it $p ( \cdot ) \in \mathcal{P}$ such that $u(t, x) = -V_{y} (t, x, p (t, x))$, $(t, x) \in Q \cup \Gamma$, $t \in [0, T]$. From transversality condition (17) we infer that, for $(t, x, p) \in P \text{, (remember } \Delta _{x}V_{y} =( \Delta _{x}V_{y_{1}}, 0$)

$V_{t} (t, x, p (t, x)) - \Delta _{x}V (t, x, p (t, x)) =\bar{y}^{0} (\frac{ d}{ dt}V_{y^{0}} (t, x, p (t, x)) - \Delta _{x}V_{y^{0}} (t, x, p (t, x)))\\ +y (t, x) (\frac{ d}{ dt}V_{y} (t, x, p (t, x)) - \Delta _{x}V_{y} (t, x, p (t, x)))\text{.}$

(29)

From (2)-(5) (see also (13)-(14)) we have

$\frac{ d}{ dt}V_{y} (t, x, p (t, x)) - \Delta _{x}V_{y} (t, x, p (t, x)) = -f(t, x, -V_{y} (t, x, p (t, x))). \label{Veq}$

(30)

Putting (30) into (29) and applying (21) we get equality

$y^{0}(\frac{ d}{ dt}V_{y^{0}} (t, x, p (t, x)) - \Delta _{x}V_{y^{0}} (t, x, p (t, x))) +y^{0}F( -V_{y}(t, x, \overset,{p}(t, x))). \label{vin}$

(31)

Following the same way as above but now using equality (25) we come to the equality

$\bar{y}^{0} (\frac{ d}{ dt}V_{y^{0}} (t, x, \bar{p} (t, x)) - \Delta _{x}\text{}V_{y^{0}} (t, x, \bar{p} (t, x))) +\bar{y}^{0}F( -V_{y}(t, x, \bar{p}(t, x))) =0. \label{veq}$

(32)

Now we consider dual Hamilton-Jacobi type equation on $\Gamma _{1}$ i.e. relations (22)-(24). Considering transversality condition at the points belonging to $P_{1}$ we have

$\frac{ \partial V (t, x, p (t, x))}{ \partial \nu } =y^{0}\frac{ \partial V_{y^{0}} (t, x, p (t, x))}{ \partial \nu } +y (t, x) \frac{ \partial V_{y} (t, x, p (t, x))}{ \partial \nu }. \label{vbtr}$

(33)

From (2)-(5) (see also (15)-(16)) we have, for the same $(u ( \cdot ), w ( \cdot ))$ and $p ( \cdot )$ at $P_{1}$

$\frac{ \partial V_{y} (t, x, p (t, x))}{ \partial \nu } = -g (t, x, -V_{y} (t, x, p (t, x)), w(t, x)). \label{Vbeq}$

(34)

Putting (34) into (33) and applying (22) we get inequality at $P_{1}$

$\bar{y}^{0} \frac{ \partial V_{y^{0}} (t, x, p (t, x))}{ \partial \nu } \geq y^{0}h(w(t, x)) \label{vbin}$

(35)

Similarly we get equality at $P_{1}$

$\bar{y}^{0} \frac{ \partial V_{y^{0}} (t, x, \bar{p} (t, x))}{ \partial \nu } = -\bar{y}^{0}h(\bar{w} (t, x)). \label{vbeq}$

(36)

Let us integrate over $Q$ equality (31) and equality (32). Next we apply boundary conditions (28), (36) and (35), (36), then we get

$y^{0}\int _{\Omega }V_{y^{0}}\text{}(T, x, p(T, x))dx -y^{0}\int _{\Omega }V_{y^{0}}\text{}(0, x, p(0, x))dx\\ +y^{0}\int _{0}^{T}\int _{\Omega }F( -V_{y_{2}}\text{}(t, x, p(t, x)))dxdt +y^{0}\int _{0}^{T}\int _{ \partial \Omega }h(w(t, x))dxdt \leq 0,\\ \bar{y}^{0}\int _{\Omega }V_{y^{0}}\text{}(T, x, \bar{p}(T, x))dx -\bar{y}^{0}\int _{\Omega }V_{y^{0}}\text{}(0, x, \bar{p}(0, x))dx\\ +\bar{y}^{0}\int _{0}^{T}\int _{\Omega }F( -V_{y_{2}}\text{}(t, x, \bar{p}(t, x)))dxdt +\bar{y}^{0}\int _{0}^{T}\int _{ \partial \Omega }h(\bar{w}(t, x))dxdt \leq 0.$

(37)

From the above relations, (24), (27) and taking into account that $\bar{y}(0, x) =y(0, x) =\xi (x)$ we infer that

$\begin{array}{c} -\bar{y}^{0}\int _{0}^{T}\int _{\Omega }F( -V_{y_{2}}(t, x, \bar{p}(t, x)))dxdt -\bar{y}^{0} \int _{0}^{T}\int _{ \partial \Omega }^{0}h(\bar{w}(t, x))dxdt \\ \text{}\text{} -\bar{y}^{0}\int _{\Omega }l( -V_{y_{2}}(T, x, p(T, x)))dx \leq y^{0}\int _{0}^{T}\int _{\Omega }F( -V_{y_{2}}(t, x, p(t, x)))dxdt\\ -y^{0}\int _{0}^{T}\int _{ \partial \Omega }h(w(t, x))dxdt -y^{0}\int _{\Omega }l( -V_{y_{2}}(T, x, p(T, x)))dx.\end{array}$

Directly from (37) and (27) we infer

Corollary 1. The dual optimal value can also be defined with the help of $V_{y^{0}}$ i.e. we have

$\begin{align*}\bar{y}^{0} \int _{\Omega }V_{y^{0}} (0, x, \bar{p}\text{}\text{}(0, x))dx =\text{} -\bar{y}^{0}\int _{0}^{T}\int _{\Omega }F( -V_{y_{2}}(t, x, \bar{p}(t, x)))dxdt \\ -\bar{y}^{0} \int _{0}^{T}\int _{ \partial \Omega }^{0}h(\bar{w}(t, x))dxdt -\bar{y}^{0}\int _{\Omega }l( -V_{y_{2}}(T, x, \bar{p}(T, x)))dx.\end{align*}$

---

4. Dual feedback control

In optimal control theory all what we want to find is to calculate optimal control and optimal value. However, in practice, a feedback control is more important than a value function. It turns out that the dual dynamic programming approach allows to define a kind of a feedback control. In fact with the help of the dual feedback control we can formulate and prove the verification theorem. Surprisingly, the dual feedback control have better properties than the classical one in spite of that it appears on the boundary. First we define general feedback control on the boundary and then optimal feedback control.

Definition 4.1. A function $\boldsymbol{w} =\boldsymbol{w} (t, x, p) =(\boldsymbol{w}_{1}(t, x, p),..., \boldsymbol{w}_{N}(t, x, p))$ defined in $P_{1}$ with values in $M$, we call dual feedback controls, if, for each $p \in c l \boldsymbol{Y}$, there exists any solution $\boldsymbol{u} (t, x, p) =(\boldsymbol{u}_{1}(t, x, p), \boldsymbol{u}_{2}(t, x, p))$ of the equation

$\frac{ \partial u}{ \partial t} - \Delta u_, =f(t, x, u), \text{ }(t, x) \in Q$

with the boundary condition

$\frac{ \partial u_{1}}{ \partial \nu } +\alpha u_{1} =\int _{\Omega }\sum\limits_{i =1}^{N}\boldsymbol{w}_{i}(t, \sigma , p\text{})K_{i}(x, \sigma )u_{2}(t, x)dx,(t, \sigma ) \in \sum\limits_{1}.$

Next step is to define optimal dual feedback control.

Definition 4.2. Dual feedback controls $\bar{\boldsymbol{w}}(t, x, p)$ defined in $P_{1}$ we call optimal dual feedback controls if there exist functions $\bar{\boldsymbol{u}}(t, x, p)$ corresponding to $\bar{\boldsymbol{w}}(t, x, p)$ as in the former definition and there exists a function $\bar{p} ( \cdot ) \in \mathcal{P}$ defined in $\bar{Q}$ such that there exists a pair $(\bar{u} ( \cdot ), \bar{w} ( \cdot )) \in A d_{\bar{\boldsymbol{u}}}$ defined by

$\bar{u}(t, x) =\bar{\boldsymbol{u}}(t, x, \bar{p} (t, x)), \;\text{}\;(t, x) \in \bar{Q},$

$\bar{w} (t, x) =\bar{\boldsymbol{w}} (t, x, \bar{p} (t, x)), \;\text{}\;(t, x) \in \sum\limits_{1}$

with optimal value $S_{D}^{\bar{\boldsymbol{u}}, \bar{y}^{0}}$

Following the same way as in the proof of Theorem 3.1 one can prove the theorem on sufficient optimality conditions for our problem (1)-(5) in terms of optimal dual feedback controls.

Theorem 4.3. Let $\bar{\boldsymbol{w}} (t, x, p),$ defined in $P_{1},$ be dual feedback control and let $\bar{\boldsymbol{u}} (t, x, p)$ be the function corresponding to $\bar{\boldsymbol{w}}$ as in the Definition 4.1. Assume that there exists in $P$ a $H^{2} (P)$ solution $V (t, x, p)$ to (21)-(24) such that $\,$

$V_{y} (t, r, p) = -\bar{\boldsymbol{u}}(t, x, p)$

and that condition (17) in $P$ is satisfied. Let $\bar{p} ( \cdot ) \in \mathcal{P}\text{,}$ defined in $\bar{Q}\text{,}$ be such that there exists a pair $(\bar{u} ( \cdot ), \bar{w} ( \cdot )) \in A d_{\bar{\boldsymbol{u}}}$

$S_{D}^{\bar{\boldsymbol{u}}, \bar{y}^{0}} = -\bar{y}^{0} \int _{\Omega }V_{y^{0}}(0, x, \bar{p}\text{}(0, x))dx$

and that $\bar{u} (t, x) =\bar{\boldsymbol{u}} (t, x, \bar{p} (t, x)),\text{}$ $(t, x) \in \bar{Q},$ $\bar{w} (t, x) =\bar{\boldsymbol{w}} (t, x, \bar{p} (t, x)), \,(t, x) \in \sum _{1}.$ Then $\bar{\boldsymbol{w}} (t, x, p)$ defined in $P_{1}$ is optimal dual feedback control.

---

5. Sufficient conditions for $\varepsilon$-optimality

The theory presented in the last two subsections being in terms of dual dynamic programming gives us a possibility to find at least formally the optimal value. However in practice it is difficult (or even impossible) to solve equations stated there in exact form. In fact we solve such a system using different approximate (numerical) methods. Therefore what we can get then is eventually approximate optimality. This is why in this section we present dual dynamic approach to sufficient conditions for approximate ($\varepsilon$-optimality) optimality. Just dual $\varepsilon$-optimality conditions are base to construct computational method for approximate optimality. Let us recall that for fixed $\bar{y}^{0}$ and $\bar{\boldsymbol{u}}$ the dual optimal value is defined as

$S_{D}^{\bar{\boldsymbol{u}}, \bar{y}^{0}} =\inf\limits_{(w, u) \in Ad_{\bar{\boldsymbol{u}}}\text{}} -y^{0}\int _{0}^{T}F( -V_{y_{2}}(t, x, p(t, x)))dx$

$-y^{0}\int _{0}^{T}\int _{ \partial \Omega }h(w(t, x)dxdt -y^{0}\text{}\int _{\Omega }l( -V_{y_{2}}\text{}(T, x, p(T, x)))dx.$

Dual $\varepsilon$-optimal value for problem (1)-(5) we call each value $S_{\varepsilon D}^{\boldsymbol{u}, \bar{y}_{\varepsilon }^{0}}$ satisfying inequality

$S_{D}^{\bar{\boldsymbol{u}}, \bar{y}^{0}} \leq S_{\varepsilon D}^{\boldsymbol{u}, \bar{y}_{\varepsilon }^{0}}\text{} \leq S_{D}^{\bar{\boldsymbol{u}}, \bar{y}^{0}} -4\varepsilon \bar{y}_{\varepsilon }^{0}.\text{}$

(38)

Let us fix $m >0$. As for $\varepsilon$-optimal value we use in general inequality instead of equality, it suggests that expressions allowing to derive Theorem 3.1 should satisfy also suitable inequalities. Thus we shall use the following system of inequalities for auxiliary function $\tilde{V} :$ dual Hamilton-Jacobi inequality

$\varepsilon \bar{y}_{\varepsilon }^{0} \leq \tilde{V}_{t} (t, x, p) - \Delta _{x}\tilde{V} (t, x, p) +y f (t, x, -\tilde{V}_{y} (t, x, p)) +y^{0}F( -\tilde{V}_{y_{2}}(t, x, p)) \leq 0 \label{ine}$

(39)

and dual Hamilton-Jacobi type inequality on $P_{1}:$

$\varepsilon \bar{y}_{\varepsilon }^{0} \leq \inf\limits_{w \in M}\{\frac{ \partial \tilde{V} (t, x, p)}{ \partial \nu } +y g (t, x, -\tilde{V}_{y} (t, x, p), w) -y^{0}h(w)\} \leq 0, \label{inb}$

(40)

$\frac{ \partial \tilde{V} (t, x, p)}{ \partial \nu } =0, \;\text{for}\;(t, r, p) \in P_{2}, \label{p1}$

(41)

$\; -y^{0}\;\text{}\tilde{V}_{y^{0}} (T, x, p) = -y^{0}l( -\widetilde{V}_{y_{2}}(T, x, p)). \label{v0in}$

(42)

$\tilde{V}_{y}$ satisfies instead of boundary conditions of type (5)-(6) the following inequality in $P_{1}$(each coordinate):

$0 \geq \frac{ \partial \tilde{V}_{y} (t, x, p)}{ \partial \nu } +g (t, x, -\tilde{V}_{y} (t, x, p), w) \geq \frac{\varepsilon }{m} \bar{y}_{\varepsilon }^{0}. \label{vyin}$

(43)

We want to apply our theory to numerical solutions of (2)-(6), therefore instead of system of equations we shall deal with systems of inequalities:

$0 \leq \frac{ \partial u}{ \partial t} -( \Delta u_{1}, 0) -f(t, x, u) \leq -\frac{\varepsilon }{m}\bar{y}_{\varepsilon }^{0} \label{xin}$

(44)

satisfying the boundary condition

$0 \geq -\frac{ \partial u}{ \partial \nu } +g (t, x, u, w (t, x)) \geq \frac{\varepsilon }{m} \bar{y}_{\varepsilon }^{0}\;\text{,} \label{xbin}$

(45)

Thus in this section by the set of admissible controls and states i.e. satisfying (44)-(45) we denote $A d_{\varepsilon }$.

---

6. $\varepsilon$-optimality

Now we are ready to describe the concept of $\varepsilon$-optimal pair, to formulate and to prove sufficient $\varepsilon$-optimality for problem (1)-(6) i.e. $\varepsilon$-version of verification theorem. Assume that there exists $\tilde{V}$ satisfying (17) and (39)-(43). Then we define similarly as in section 2

$\boldsymbol{u}_{\varepsilon } (t, x, p) = -\tilde{V}_{y} (t, x, p)\;\text{,}\;(t, x, p) \in c l P. \label{xe}$

(46)

For $\bar{y}_{\varepsilon }^{0}$ and $\boldsymbol{u}_{\varepsilon }$ we define similarly as in section 2 $A d_{\boldsymbol{u}}$

$A d_{\boldsymbol{u}_{\varepsilon }} = \left \{(u ( \cdot ), w ( \cdot )) \in A d_{\varepsilon } ;\;\text{exist}\;p (t, x) =(y^{0}, y (t, x))\;\text{,}\;(t, x) \in Q\;\text{,}\;\right. \label{ade} \\ \;\text{}\;y ( \cdot ) \in (H^{2} (Q))^{2}\;\text{,}\;y^{0} \leq 0\;\text{,}\;(t, x, p (t, x)) \in c l P, \\ y(0, x) =\xi (x), \;\text{}\;u(t, x) =\boldsymbol{u}_{\varepsilon }(t, x, p (t, x)), \;\text{}\;(t, x) \in clQ\}$

and $\mathcal{P}_{\varepsilon }$

$\begin{array}{c}\mathcal{P}_{\varepsilon } =\left \{p (t, x) =(\bar{y}_{\varepsilon }^{0}, y (t, x)), \;\text{}\;(t, x) \in Q ;\;\text{}\;(t, x, p (t, x)) \in c l P, \;\text{}y ( \cdot ) \in (H^{2}(Q))^{2}\text{,}\right. \\ \left.\sup\limits_{(t, x) \in Q}\left \vert y (t, x)\right \vert _{R^{2}} \leq m\;\text{,}\text{ }y >0\text{, }\text{exist}\;(u ( \cdot ), w ( \cdot )) \in A d_{\boldsymbol{u}_{\varepsilon }}, \;\text{}\;\right. \\ u(t, x) = -\tilde{V}_{y} (t, x, p (t, x))\;\text{,}\;(t, x) \in Q \cup \Gamma \}\text{.}\end{array}$

Now we are ready to define notions of $\varepsilon$-optimal dual feedback control $\bar{\boldsymbol{w}}_{\varepsilon } (t, x, p)$ and of $\varepsilon$-optimal pair $(\bar{u}_{\varepsilon } ( \cdot ), \bar{w}_{\varepsilon } ( \cdot ))$.

Definition 6.1. Dual feedback control $\bar{\boldsymbol{w}}_{\varepsilon } (t, x, p)$ we call $\varepsilon$-optimal dual feedback control if there exist a function $\bar{\boldsymbol{u}}_{\varepsilon } (t, x, p)$ in $P\text{,}$ accordingly to Definition 4.1 and a function $\bar{p}_{\varepsilon } ( \cdot ) \in \mathcal{P}_{\varepsilon }$ defined in $\bar{Q}$, such that the pair defined by

$\bar{u}_{\varepsilon } (t, x) =\bar{\boldsymbol{u}}_{\varepsilon } (t, x, \bar{p}_{\varepsilon } (t, x)), \;\text{}\;(t, x) \in \bar{Q},\\ \bar{w}_{\varepsilon } (t, x) =\bar{\boldsymbol{w}}_{\varepsilon } (t, x, \bar{p}_{\varepsilon } (t, x)), \;\text{}\;(t, x) \in \bar{Q}$

(47)

belongs to $A d_{\bar{\boldsymbol{u}}_{\varepsilon }}$ and that this pair defines $\varepsilon$-optimal value

$S_{\varepsilon D}^{\bar{\boldsymbol{u}}_{\varepsilon }\bar{y}_{\varepsilon }^{0}} = -\bar{y}_{\varepsilon }^{0}\int _{\Omega }\tilde{V}_{y^{0}}(0, x, \bar{p}_{\varepsilon }\text{}\text{}\text{}\text{}(0, x))dx.\,$

(48)

Definition 6.2. For given $\tilde{V} \in H^{2}(P)$ satisfying (17) and (39)-(43) let $\bar{\boldsymbol{u}}_{\varepsilon } (t, x, p)$ in $P$ be defined by (46). Let $\bar{p}_{\varepsilon } ( \cdot ) \in \mathcal{P}_{\varepsilon }$ be defined in $\bar{Q}$ and let $\bar{u}_{\varepsilon }$ be defined by (47). Let $\bar{w}_{\varepsilon } ( \cdot )$ be any admissible control such that $(\bar{u}_{\varepsilon } ( \cdot ), \bar{w}_{\varepsilon } ( \cdot )) \in A d_{\bar{\boldsymbol{u}}_{\varepsilon }}$. The pair $(\bar{u}_{\varepsilon } ( \cdot ), \bar{w}_{\varepsilon } ( \cdot ))$ we call $\varepsilon$-optimal pair with respect to all pairs $(u ( \cdot ), w ( \cdot )) \in A d_{\bar{\boldsymbol{u}}_{\varepsilon }}$ if

$\,\begin{array}{c} -\bar{y}_{\varepsilon }^{0}\int _{0}^{T}\int \text{}_{\Omega }F( -V_{y_{2}}(t, x, \bar{p}_{\varepsilon }(t, x)))dx -\bar{y}_{\varepsilon }^{0}\int _{0}^{T}\int _{ \partial \Omega }h(\bar{w}_{\varepsilon }(t, x))dxdt\\ -\bar{y}_{\varepsilon }^{0}\text{}\int _{\Omega }l( -V_{y_{2}}\text{}(T, x, \bar{p}_{\varepsilon }(T, x)))dx\\ \leq -\bar{y}_{\varepsilon }^{0}\int _{0}^{T}\int _{\Omega }F( -\overset,{V_{y_{2}}}(t, x, p(t, x)))dx -\bar{y}_{\varepsilon }^{0}\int _{0}^{T}\int _{ \partial \Omega }h(w(t, x))dtdx \\ \text{} -\bar{y}_{\varepsilon }^{0}\text{}\int _{\Omega }l( -V_{y_{2}}\text{}(T, x, p(t, x)))dx -4\varepsilon \bar{y}_{\varepsilon }^{0}. \end{array}\begin{array}{c}\end{array}$

Having all the above notions we can formulate the verification theorem for $\varepsilon$-optimality.

Theorem 6.3. Assume that there exists $\tilde{V} \in H^{2} (P)$ satisfying (17) and (39)-(43). Take $\bar{p}_{\varepsilon } ( \cdot ) \in \mathcal{P}_{\varepsilon }$ and $(\bar{u}_{\varepsilon } ( \cdot ), \bar{w}_{\varepsilon } ( \cdot )) \in A d_{\bar{\boldsymbol{u}}_{\varepsilon }}$ such that $\bar{u}_{\varepsilon } (t, x) = -\tilde{V}_{y} (t, x, \bar{p}_{\varepsilon } (t, x))\text{,}$ $(t, x) \in \bar{Q}$. Moreover, assume that the trio $(\bar{u}_{\varepsilon } ( \cdot ), \bar{w}_{\varepsilon } ( \cdot ), \bar{p}_{\varepsilon } ( \cdot ))$ satisfies

$\frac{ d}{ dt}\tilde{V}_{y} (t, x, \bar{p}_{\varepsilon } (t, x)) - \Delta _{x}\tilde{V}_{y} (t, x, \bar{p}_{\varepsilon } (t, x)) \\ +f(t, x, -\tilde{V}_{y} (t, x, \bar{p}_{\varepsilon } (t, x))) \geq \bar{y}_{\varepsilon }^{0}\frac{\varepsilon }{m}, \text{ }(t, x) \in Q,$

(49)

$\,\begin{array}{c} -\frac{\varepsilon }{m} \bar{y}_{\varepsilon }^{0} \geq \tilde{V}_{t} (t, x, \bar{p}_{\varepsilon } (t, x)) - \Delta _{x}\tilde{V} (t, x, \bar{p}_{\varepsilon } (t, x)) \\ \begin{array}{c}\begin{array}{c}\begin{array}{c} +\bar{y}\text{}_{\varepsilon }\text{}(t, x)(\text{}\text{}f (t, x, -\tilde{V}_{y} (t, x, \bar{p}_{\varepsilon } (t, x)) +\bar{y}_{\varepsilon }^{0}F( -\tilde{V}_{y_{2}}(t, x, \bar{p}_{\varepsilon }(t, x, ))), \text{}\text{}\end{array}\end{array}\text{}_,^{}\end{array}\end{array}$

$\frac{ \partial \widetilde{V} (t, x, \bar{p} _{\varepsilon }(t, x))}{ \partial \nu } +\bar{y} _{\varepsilon }(t, x) g (t, x, -\tilde{V}_{y} (t, x, \bar{p}_{\varepsilon } (t, x)), \bar{w}_{\varepsilon }(t, x)) \\ -\bar{y}_{\varepsilon }^{0}h\left (\bar{w}_{\varepsilon }(t, x)\right ) \leq -\text{}\frac{\varepsilon }{m} \bar{y}_{\varepsilon }^{0}, \label{e2}$

(50)

$-\bar{y}_{\varepsilon }^{0}\tilde{V}_{y^{0}}(T, x, \bar{p}_{\varepsilon } (t, x)) = -\bar{y}_{\varepsilon }^{0}l( -\tilde{V}_{y_{2}} (T, x, \bar{p}_{\varepsilon } (T, x))), \\ \frac{ \partial \widetilde{V} (t, x, \bar{p} _{\varepsilon }(t, x))}{ \partial \nu } =0, \text{}\text{for} (t, x,p\bar{}\varepsilon (t,x)) \in P_{2}.$

(51)

Then the pair $(\bar{u}_{\varepsilon } ( \cdot ), \bar{w}_{\varepsilon } ( \cdot ))$ is an $\varepsilon$-optimal with respect to all pairs $(u ( \cdot ), w ( \cdot )) \in A d_{\bar{\boldsymbol{u}}_{\varepsilon }}$

Proof. Take any $(u ( \cdot ), w ( \cdot )) \in A d_{\bar{\boldsymbol{u}}_{\varepsilon }}$ and $p ( \cdot ) \in \mathcal{P}_{\varepsilon }$ such that $u (t, x) = -\tilde{V}_{y} (t, x, p (t, x))$, $(t, x) \in Q$. We follow the same way as in the proof of Theorem 3.1, i.e. from (17) we have, for $(t, x, p) \in P \text{,}$

$\tilde{V}_{t} (t, x, p (t, x)) - \Delta _{x}\tilde{V} (t, x, p (t, x)) =\bar{y}_{\varepsilon }^{0} (\frac{ d}{ dt}\tilde{V}_{y^{0}} (t, x, p (t, x)) - \Delta _{x}\tilde{V}_{y^{0}} (t, x, p (t, x)))\\ +y (t, x) (\frac{ d}{ dt}\tilde{V}_{y} (t, x, p (t, x)) - \Delta _{x}\tilde{V}_{y} (t, x, p (t, x))).$

Similarly, we have by (44)

$-\frac{ d}{ dt}\tilde{V}_{y} (t, x, p (t, x)) + \Delta _{x}\tilde{V}_{y} (t, x, p (t, x)) -f (t, x, -\tilde{V}_{y} (t, x, p (t, x))) \geq 0$

and then applying (39) (having in mind that $y >0$) we get inequality

$\bar{y}_{\varepsilon }^{0} (\frac{ d}{ dt}\tilde{V}_{y^{0}} (t, x, p (t, x)) - \Delta _{x}\tilde{V}_{y^{0}} (t, x, p (t, x)))\\ +\bar{y}_{\varepsilon }^{0}F( -\tilde{V}(t, x, \overset,{p}(t, x))) \geq \varepsilon \bar{y}_{\varepsilon }^{0}$

(52)

and using inequality (25) we come to the inequality

$- \varepsilon \bar{y}_{\varepsilon }^{0} \geq \bar{y_{\varepsilon }}^{0} (\frac{ d}{ dt}\tilde{V}_{y^{0}} (t, x, \bar{p}_{\varepsilon } (t, x)) - \Delta _{x}\tilde{V}_{y^{0}} (t, x, \bar{p}_{\varepsilon } (t, x))) +\bar{y}_{\varepsilon }^{0}F( -\tilde{V}(t, x, \bar{p}_{\varepsilon }(t, x))). \label{eVeq}$

(53)

Considering transversality condition at the points belonging to $\Gamma$ we have

$\frac{ \partial \tilde{V}(t, x, p (t, x))}{ \partial \nu } =\bar{y}^{0} \frac{ \partial \tilde{V}_{y^{0}} (t, x, p (t, x))}{ \partial \nu } +y (t, r) \frac{ \partial \tilde{V}_{y} (t, x, p (t, x))}{ \partial \nu }. \label{Vbtr}$

(54)

From (43) we have, for the same $(u ( \cdot ), w ( \cdot ))$ and $p ( \cdot )\text{,}$ at $P_{1}$

$\frac{\tilde{ \partial V}_{y^{0}}(t, x, p(t, x))}{ \partial \nu } +g(t, x, -\tilde{V_{y}}(t, x, p(t, x), w(t, x))) \leq 0.$

Hence we get inequality at $P_{1}$

$\varepsilon \bar{y}_{\varepsilon }^{0} \leq \bar{y}_{\varepsilon }^{0} \frac{ \partial \overset,{\tilde{V}}_{y^{0}}(t, x, \overset,{p}_{}(t, x))}{ \partial \nu } +\bar{y}_{\varepsilon }^{0}h(w(t, x)). \label{vbineq}$

(55)

Similarly, using (50) we get inequality at $P_{1}$

$\bar{y}_{\varepsilon }^{0} \frac{ \partial \overset,{\tilde{V}}_{y^{0}}(t, x, \overset,{\bar{p}_{\varepsilon }}_{}(t, x))}{ \partial \nu }\text{} +\bar{y}_{\varepsilon }^{0}h(\bar{w}_{\varepsilon }(t, x)) \leq -\varepsilon \text{} \bar{y}_{\varepsilon }^{0}. \label{vbine}$

(56)

Let us integrate over $Q$, inequality (52) and (53). Next we apply boundary conditions (55), (56), respectively, and then we get

$\begin{array}{c}\bar{y}_{\varepsilon }^{0}\int _{\Omega }\tilde{V}_{y^{0}}\text{}(T, x, p(T, x))dx -\bar{y}_{\varepsilon }^{0}\int _{\Omega }\tilde{V}_{y^{0}}\text{}(0, x, p(0, x))dx \\ \bar{y}_{\varepsilon }^{0}\int _{0}^{T}\int _{\Omega }F( -\tilde{V_{y_{2}}}(t, x, \overset,{p}(t, x)))dx +\bar{y}_{\varepsilon }^{0}\int _{0}^{T}\int _{ \partial \Omega }h(\overset,{w}(t, x))dtdx \geq 2\varepsilon \bar{y}_{\varepsilon }^{0},\end{array}$

$\bar{y_{\varepsilon }}^{0} \int _{\Omega }\tilde{V}_{y^{0}} (T, x, \bar{p}_{\varepsilon }(T, x))dx -\bar{y}_{\varepsilon }^{0}\int _{\Omega }\tilde{V}_{y^{0}}(0, x, \bar{p}_{\varepsilon }(0, x))dx\\ +\bar{y}_{\varepsilon }^{0}\int _{0}^{T}\int _{\Omega }F( -\tilde{V_{y_{2}}}(t, x, \bar{p}_{\varepsilon }(t, x)))dx +\bar{y}_{\varepsilon }^{0}\int _{0}^{T}\int _{ \partial \Omega }h(\bar{w}_{\varepsilon }(t, x))dtdx \leq -2\varepsilon \bar{y}_{\varepsilon }^{0}.$

From the above relations we infer that

$\begin{array}{c}\begin{array}{c} -\bar{y}_{\varepsilon }^{0}\int _{0}^{T}\int _{\Omega }F(\bar{u}_{2\varepsilon } (t, x))dx -\bar{y}_{\varepsilon }^{0}\int _{0}^{T}\int _{ \partial \Omega }h(\bar{w}_{\varepsilon }(t, x))dtdx -\bar{y}_{\varepsilon }^{0}\int _{\Omega }l(\bar{u}_{2\varepsilon } (T, x))\text{}\end{array} \\ \text{} \leq -\bar{y}_{\varepsilon }^{0}\int _{0}^{T}\int _{\Omega }F( -\tilde{V_{y_{2}}}(t, x, p(t, x)))dx -\bar{y}_{\varepsilon }^{0}\int _{0}^{T}\int _{ \partial \Omega }h(w(t, x))dtdx\text{} \\ -\bar{y}_{\varepsilon }^{0}\text{}\int _{\Omega }l(u_{2}\text{}(T, x)dx -4\varepsilon \bar{y}_{\varepsilon }^{0}\text{}.\end{array}$

This is just the assertion of the theorem.

---

6.1. Computational algorithm

The sufficient conditions formulated for $\varepsilon$-value function allows us to build numerical approach to calculate suboptimal pair $(\bar{u}_{\varepsilon } ( \cdot ), \bar{w}_{\varepsilon } ( \cdot ))$. The algorithm, we present below, ensures that we find in finite number of steps suboptimal pair.

Algorithm:

1. Fix $m >0\text{,}$ $\varepsilon >0$ and calculate auxiliary function $\tilde{V}$ from (39)-(43).

2. Form $A d_{\boldsymbol{u}_{\varepsilon }}$ as a finite family of $N$ pairs $(u ( \cdot ), w ( \cdot )):$

a) Define controls $w_{n}$ in $\Gamma _{1}$, $n =1,\ldots, N\text{.}$

b) To calculate $u^{n} \text{,} n =1,\ldots, N,$ solve inequalities (44)-(45).

3. Find minimal value of $J(u^{n}, w_{n})\text{,}$ $n =1,\ldots, N$ and corresponding to it pair denote by $(\hat{u} ( \cdot ), \hat{w} ( \cdot ))\text{.}$

4. Assume $\bar{y}_{\varepsilon }^{0} = -1$ and determine $\hat{y} ( \cdot )$ from the relation

$\hat{u}\text{}(t, x) = -\tilde{V}_{y} (t, x, -1, \hat{y} (t, x))\text{.}$

(57)

5. For $\tilde{V}$ and $(\hat{u} ( \cdot ), \hat{w} ( \cdot ), \hat{y} ( \cdot ))$ check the inequalities (53)-(56)

a) If $\tilde{V}$ and $(\hat{u} ( \cdot ), \hat{w} ( \cdot ), \hat{y} ( \cdot ))$ satisfy (53)-(56)

then $(\hat{u} ( \cdot ),\hat{w} ( \cdot ))$ is an $\varepsilon$-optimal pair and $J(\hat{u},\hat{w})$ is an $\varepsilon$-optimal value.

b) If $\tilde{V}$ and $(\hat{u} ( \cdot ), \hat{w} ( \cdot ), \hat{y} ( \cdot ))$ do not satisfy (53)-(56)

then go to 2.

---

6.2. Numerical calculations

Numerical experiments we do for the same data as in [1] i.e. $\Omega =(0, 1) \times (0, 1)$, $T =1$, $\beta =50$, $\gamma =2$, $\Gamma _{1}$ is one of the side of $\Omega$ corresponding to the seashore. It is assumed that control $w$ is only time dependent, one dimensional and then $K(x, \sigma ) =w(t)K_{0}(x, \sigma )$, where $K_{0}(x, \sigma ) =\sum _{i =1}^{m}a_{i}\chi _{i}(\sigma )\tilde{\chi }_{i}(x)$, with $\chi _{i}$ the characteristic function of the interval $[x_{i}, x_{i +1}]$. We denote by $\tilde{\chi }_{i}$ the characteristic function of the rectangular subdomain $\Omega _{i} =[x_{i}, x_{i +1}] \times [0, 1] \subset \Omega$ where the knots $x_{i}$ belong to $\Gamma _{1}$. The constants $a_{i}$ are weights assigned to the subdomains $\Omega _{i}$. For $g$ it is chosen expression $\beta \frac{u_{1}}{1 +\gamma u_{1}}$. The number $m =4$, the knots are equidistant and $a_{i} =0.1$. The initial conditions $u_{1} \neq 0$ on $\Omega _{s} =[0.3, 0.7],$ $u_{1} =0$ on $\Omega \setminus \Omega _{s},$ $u_{2} =0$ on $\Omega$. For the cost functional it is considered $F =0$, $l(u) =u$ and $h(w) =w^{ -2}$. Control $w$ has values in $[1, 2]$. In order to make calculations we apply FreeFem++-cs 14.3 package from the site http://freefem-cs.software.informer.com/14.3/ and we implement in this application the steps from the above algorithm. We divided time interval $[0, 1]$ on $10$ equal intervals and consider differences instead of derivative in time. We choose $N =500$. Next we choose randomly in each step of time controls $w$ from $[1, 2].$ We included to the set of admissible controls as one of the control that which is considered in [1] as optimal control. Next we calculated corresponding $u_{1}, u_{2}$ and the values of $J(u, w)$. Then we found the minimal value among those of $J$. We repeated these procedure $5$ times and chosen minimal value from those five former minimal values and wrote down corresponding $\hat{u}, \hat{w}$. These are at $t =0.0, 0.1, 0.2,..., 0.9, 1$,

$\hat{w}= 1.11941, 1.04171, 1.06656, 1.04117, 1.05112, 1.00715, 1.33937, 1.07568, 1.06864, 1.32913, 1.30369$ and the graph of it is:

$\hat{u}_{1}(t) =\int _{\Omega }\hat{u}_{1}(t, x)dx =0, 0.422815, 0.372519, 0.354227, 0.374098, 0.439481, \\ 0.57045,0.770489, 1.04842, 1.42482, 1.90287$

the graph of it is:

$\hat{u}_{2}(t) =\int _{\Omega 2}(t, x)dx =0,$$1.02899, 1.90447, 2.6737, 3.40012, 4.1412, 4.94705, 5.83374, 6.79449, 7.81148, 8.8582$

the graph of it is:

The $\varepsilon$-value of the functional $J(\hat{u}, \hat{w}) =8.92283$. Next we follow the steps 4. and 5. from the former subsection. It turn out that for $\varepsilon =0.003$ the pair $(\hat{u}, \hat{w})$ is $\varepsilon$-optimal. It differs from [1], the value of our functional is smaller (in [1] it equals $12.38386$) and we found different control $\hat{w\text{.}}$

---

Acknowledgments

The paper was inspired by the lecture given by V. Capasso during the conference Micro and Macro Systems in Life Sciences in Bedlewo 2015.

---