Numerical investigation for a boundary optimal control of reaction-advection-diffusion equation using penalization technique

1.   Introduction

In this paper, we aim to study the optimal control problem for the reaction-advection-diffusion equation with non-sufficiently smooth Dirichlet condition on the boundary on a regular bounded domain of $\mathbb{R}^2$. We start by proving the existence and uniqueness of the solution of the state equation. Since the use of Dirichlet boundary condition poses practical difficulties, we will apply a penalization technique that approaches the state equation by another having a boundary condition of type Robin, which is more easy to implement. This method was introduced by Babu${\breve {\rm{s}}}$ka ^[1] for elliptic equations and considered later in several works, see for example ^[2,3,4] where several authors obtained the expected convergence rate and in ^[5] where an asymptotic expansion of the penalized solution is studied. For the boundary (Dirichlet) control problems of partial differential equations, where the control variable is generated in $L^2$, Robin penalization thus presents a simple and efficient way to solve such problems (see ^{[6,7,8,9,10,11,12,13]}). We therefore, propose in this paper to apply this approach to the problem of boundary control of the reaction-advection-diffusion equation.

2.   Reaction-advection-diffusion equation

Let $\Omega$ be an open regular bounded domain of $\mathbb{R}^2$, with a Lipschitz-continuous boundary $\Gamma$, and we set

and

where $T$ is a positive constant. The optimal control problem consists of finding optimal Dirichlet data $g\in L^2(\Pi)$ which minimizes a cost function depending on the solution $u$ of the state Eq (2.1) given hereafter and on the control $g$. Therefore, it is necessary to define a suitable functional framework which allows us to obtain results of existence, uniqueness, and stability.

Our goal in this work is to mathematically analyze the penalization technique applied on a Dirichlet condition. More precisely, we will replace the Dirichlet condition $u = g$ on $\Pi$ by a Robin-like condition, as follows:

where $\varepsilon$ is a small enough positive constant tending to zero. The advantage of this technique is to define the state variable $u_{\varepsilon}$ in $L^2(0, T, L^2(\Omega))$ for all boundary conditions $g$ in $L^2(\Pi)$. Therefore, we will consider an optimal control problem defined on $L^2(\Pi)$. Furthermore, numerically, taking into account the Dirichlet boundary condition is not made exactly, but it is rather approached by using a penalization term. Since the state variable is used in the expression of the cost function for the optimal control, the convergence of the penalized variable $u_{\varepsilon}$ towards the state variable $u$ should be proved. However, due to some difficulties, this task is in progress and we will give only a numerical investigation regarding the convergence in this paper.

2.1.   Problem position

Let $f \in L^2(\Upsilon)$, $u_0 \in L^2(\Omega)$, and $g \in L^2(\Pi)$, and we aim to find $u$ the solution of

where $\mu$ is the diffusion coefficient, $\sigma$ is the reaction coefficient, and

represents a vector field in $(L^{\infty}(\Upsilon))^2$ describing the advection coefficient. The initial condition is given by $u_0 \in L^2(\Omega)$. $\mu$ and $\sigma$ are assumed to be positive. Equation (2.1) is called a reaction-advection-diffusion equation and has several applications in science and engineering. Transport of heat, momentum and energy, solid mechanics, CO$_2$ sequestration, computational fluid dynamics, and biophysics are a few research fields that need to solve (2.1) numerically as part of the solution strategy. For the sake of simplification, we suppose that:

and

The resolution of this problem relies on the duality technique; but, to explain it we need to recall the essential properties of Eq (2.1) with a homogeneous Dirichlet data. Problem (2.1) admits a weak formulation as follows: Find $u \in L^2(0, T, H_0^1(\Omega))$ such that, for all $v \in H_0^1(\Omega)$, we have:

where

Problem (2.2) admits (according to ^[14]) a unique solution in $L^2(0, T, H_0^1(\Omega))$ with

To write the weak formulation of system (2.1) with $g \in L^2(\Pi)$, we apply the duality technique, which involves looking for $u \in L^2(\Upsilon)$ such that, for all functions $\phi$ in

with $\phi(T) = 0$, we have

In order to prove that problem (2.3) is well-posed, we need to use the Lax-Milgram lemma. It suffices to observe that (2.3) is equivalent to: finding $u \in L^2(\Upsilon)$ such that:

where

is the unique solution of the following problem:

satisfying

Therefore, the operator

is bounded since

Therefore, the problem (2.4) (and then the problem (2.3)) admits a unique solution $u \in L^2(\Upsilon)$ with the stability relation

2.2.   Penalization technique

To overcome the difficulties of the numerical implementation of the Dirichlet condition, we will use the penalization technique using a small enough parameter $\varepsilon$ by introducing the set of problems

to approximate problem (2.1). The penalization technique allows us to work in the Sobolev space $H^1(\Omega)$. The weak formulation associated with (2.2) consists of finding $u_{\varepsilon} \in L^2(0, T, H^1(\Omega))$ such that, for all functions $v \in H^1(\Omega)$, we have

with

Problem (2.5) admits a unique solution in $L^2(0, T, H^1(\Omega))$ and satisfies the following result:

where the constant $C$ is independent of $\varepsilon$.

Concerning the convergence of $u_{\varepsilon}$ to $u$, some tools developed in ^[4] and a fixed point technique used in ^[15], could make it possible to establish a rate of convergence comparable to that demonstrated in the case of the pure diffusion equation. This conjecture is supported by some numerical experiments.

3.   Boundary optimal control problem

The goal is to examine the problem of optimal boundary control of the reaction-advection-diffusion equation. This involves minimizing a functional criterion

which depends on $g$ as well as on $u$ the solution of (2.1). We show the existence and uniqueness of the control problem, which leads us allows to write the first-order optimality conditions. Then, we introduce the sequence of control problems obtained by penalizing the Dirichlet condition and we give an existence result of the optimal solution.

3.1.   Boundary non-penalized control problem

We consider the functional cost given by

where ${u_d \in L^2(\Upsilon)}$ is the predefined state (desired profile). The dependent term of $g$ is the marginal cost.

Consider (2.1) where $g \in L^2(\Pi)$. The boundary optimal control problem is thus expressed by

Let

Then, the function $F$ then verifies the following conditions:

● $F$ is strictly convex;

● $F$ is coercive ($\lim_{\|g\|_{L^2(\Pi)} \rightarrow +\infty} F(g) = +\infty$);

● $F$ is continuous, and this follows from the continuous dependence of $u$ on $g$.

Thus, the minimization problem (3.1) admits a unique solution $\bar g \in L^2(\Pi)$.

We want to act on a parabolic partial differential equation, through the boundary condition of the Dirichlet type, so that the solution is close to a profile fixed beforehand, with the best quality-price ratio. For our purpose, it is a question of minimizing a criterion function which depends on the solution of the problem and on the data at the boundary which represents the control. The difficulty of analyzing this problem strongly depends on the space of possible controls. Indeed, if we consider the space $L^2(0, T, H^{1\over2} (\Gamma))$, we can have a well-posed state equation on $L^ 2(0, T, H^{1}(\Omega))$ that is easy to deal with, while identifying the solution to the optimal control problem presents serious difficulties inherent in the complex form of the norm of $L^2(0, T, H^{1\over2}(\Gamma))$. The alternative that we propose to study therefore amounts to working in $L^2(\Pi)$.

Looking for the necessary conditions of optimality, the triple

formed by the optimal control, optimal state, and associated conjoint state is then the unique solution of the following system:

Conversely, if $(\bar u, \bar p)$ satisfies the previous system, then $(\bar u, \dfrac{1}{\beta} \dfrac{\partial \bar p}{\partial n})$ is the solution of $({\cal{P}})$.

Remark 3.1. Strictly speaking, the advection term in the equation of the adjoint state is of the form $(div({\mathit{\boldsymbol{v}}} p))$, but thanks to the assumptions made on ${\mathit{\boldsymbol{v}}}$ (namely $div\, {\mathit{\boldsymbol{v}}} = 0$) in $\Omega$ and ${\mathit{\boldsymbol{v}}}\cdot n = 0$ on $\Pi$, we obtain the form given above.

3.2.   Penalized control problem

We are interested in this part in the penalization of the control problem $({ \cal P })$ given by (3.3). The sequence of penalized control problems, where the Dirichlet condition is taken into account by the penalization technique, is of the form

which represents an approach to problem $({\cal{P}})$. For the same conditions verified by the functional $J$, we show that the penalized control problem $({\cal{P}} _{\varepsilon})$ has a unique solution

It will then constitute an approximation of the optimal solution $(\bar u, \bar g)$ of $({\cal{P}})$. The triple $(\bar g_{\varepsilon}, \bar u_{\varepsilon}, \bar p_{\varepsilon})$ formed by the optimal control, optimal state, and associated adjoint state in the space

is the unique solution of the system of the following first-order optimality conditions:

Reciprocally, if $(\bar u_{\varepsilon}, \bar p_{\varepsilon})$ satisfies the previous system, then $(\bar u_{\varepsilon}, - \dfrac{1}{\beta} \bar p_{\varepsilon})$ is the solution of $({\cal{P}}_{\varepsilon})$. The convergence of the optimal solution $(\bar u_{\varepsilon}, \bar g_{\varepsilon})$ towards the optimal solution of the Dirichlet problem $(\bar u, \bar g)$ is in progress. We will present in the following section some numerical experiments illustrating the convergence results.

4.   Numerical investigation

In this section, we present the numerical approximation, using the finite element method, for both direct and optimal control problems as well as the penalized problems. We recall that an effective and simple way to deal with the Dirichlet condition consists of using the penalization technique presented above. Based on this procedure, the solution of Eq (2.1) and the associated control problem are respectively approximated by those of the penalized Eq (2.5) and of the penalized control problem (3.3). For the approximation of the state equations, we adopt the finite element method for the discretization in space and the implicit Euler scheme for the discretization in time. For control problems, we use the conjugate gradient method. Finally, relying on numerical tests, we confirm the desired-state convergence results.

4.1.   State equation discretization

Before numerically solving the penalized control problem (3.3), we describe in this section the numerical resolution of the reaction-advection-diffusion Eq (2.1) using the penalization technique in order to take into account the Dirichlet condition. For the discretization in space, we consider for all $h > 0$ a triangulation $\tau_h$ of $\Omega$. This triangulation of step $h$ is a partition into a finite number of triangles $\kappa$ called elements, such that two neighboring triangles have in common a vertex or an entire side. The size $h_\kappa$ of each $\kappa\in\tau_h$ is less than or equal to $h$ (providing that size is the length of the longest side of $\kappa$). The vertices of the triangle $\kappa$ are denoted ${(s_i)}_{1 \leq i \leq 3}$. The triangulation $\tau_h$ is assumed to be regular, i.e., there exists a real $\alpha > 0$ verifying $\sigma_\kappa\geq \alpha h_\kappa$, $\forall \kappa \in \tau_h$, where $\sigma_\kappa$ is the diameter of the circle inscribed in $\kappa$, and $h_\kappa$ is the size of $\kappa$. The finite element discrete space intended to approximate $H^1(\Omega)$ is made up of piecewise affine functions:

where $\mathbb {P}_1$ describes the set of polynomial of first degree. We introduce the space

where ${\tau_h}^{\Gamma}$ describes the triangulation of $\Gamma$ and $S$ is a segment of $\Gamma$. $V_h$ has a finite dimension $( = m)$. It admits a basis ${(\phi_i)}_{i = 1, m}$. The one commonly used is composed of shape functions, which are defined as follows:

The basis of $V_h^{\Gamma}$ is formed by the set of functions

which represent the image by the operator trace of functions ${{\phi_i}}$ (basic functions of $V_h$).

4.2.   Discretization of the reaction-advection-diffusion equation

By injecting into the variational formulation (2.2) the decomposition of the solution $u_{\varepsilon, h}$ in the space $V_h$, given by:

with

and choosing as test functions the functions

we get the semi-discrete system

with $U(t)$ the unknown vector of components $(u_i(t))_{1 \leq i \leq m}$, $G_{{\mathit{\boldsymbol{v}}}}(t)$ the vector of components

$G(t)$ the vector of components

if $s_i$ is a vertex on $\Gamma$ and $g_i(t) = 0$ else, where $g_h$ denotes the space approximation of $g$ in $V_h^\Gamma$, and $F(t)$ the vector of components

et

Matrices $M$, $A$, and $M^\Gamma$ are defined by:

For the temporal discretization, we consider the subdivision of the time interval

and we use the implicit Euler scheme. By denoting $U^{(n)}$ as an approximation of $U(t_n)$, we obtain the total discrete system associated with the reaction-advection-diffusion equation:

where

et

The term $G_{{\mathit{\boldsymbol{v}}}}^{(n)}$, coming from the discretization of the advection term, is the vector of components

with

and

Note that, for reasons of computational simplification, the "advection" term is treated explicitly (see ^[16,17]), which allows us to have a symmetric matrix of the system. One can choose, in principle, to also discretize the advection term implicitly, but this will complicate and make more expensive the resolution of the problem, even if this scheme has better stability.

Solving the algebraic system (4.2) requires the inversion of the matrix $\left(\dfrac{M}{dt }+ A\right)$ which is symmetric definite positive. Among the different algorithms that can be adopted, we have chosen the one that uses the Cholesky factorization. In addition, for the calculation of the integrals appearing in the expressions of the matrices intervening in (4.2), one uses the formulas of quadrature. More precisely, the integrals on the boundary are approximated by Simpson's formula, and the surface integrals are calculated either by the degree 1 formula using the vertices of the triangles or by the degree 2 formula using the midpoints of the edges.

4.3.   Discretization of the control problem

In this subsection, we study the discretization of the control problem associated with the reaction-advection-diffusion equation. The optimal solution $(\bar u, \bar g)$ of the control problem (3.1) is approximated by the optimal solution $(\bar u_\varepsilon, \bar g_\varepsilon)$ of the penalized problem (3.3). Let $(u_{d, h}, g_{d, h})$ to be an approximation of the predefined state (desired profile) on the triangulation $\tau_h$. The discrete equivalent of the control problem (3.3) is given by:

where the discrete version $J_h(., .)$ of the objective function $J(., .)$ is given by:

Using the right rectangle method for time integration, the discrete control problem can be written as:

where the expression of $J_{h, dt }$ is given by:

with

and

The calculation of the optimality conditions associated with the discrete control problem (4.3) involves the discrete adjoint state $P$ which leads to a determination of the optimal discrete control $G$. The triple $(G, U, P)$ formed by the discrete optimal control, the discrete optimal state, and the discrete adjoint state is the unique solution of the system:

where

In order to use the conjugate gradient method for solving the discrete control problem (4.3), we need to calculate the gradient of the functional $J_{h, dt }$, which is given, thanks to the adjoint state technique, by:

4.4.   Optimization algorithm

A problem like (4.3) is usually solved by iterative methods that construct a sequence ${\{u_i\}}^k_{i = 0} \; (k \in \mathbb {N})$ to gradually move towards a better approximation of the "solution", i.e, the point satisfying the optimality condition. Now that we have given a discretization or an approximation of the cost function and the deputy state, we use to evaluate the minimum of the problem (4.3) the gradient method which belongs to a family of iterative methods called the descent method since the construction of the sequence ${\{u_i\}}^k_{i = 0}$ is such that $F(u_{i+1}) < F(u_i)$ $\forall i \in \mathbb {N}$, where $F$ denotes the function to be minimized, which, in our case, is given by $J_{h, dt }$. The following outlines the gradient method:

(1) Available data:

● Initialization $u_0$, a tolerance $\gamma > 0$.

● The function $F$ to be minimized.

(2) Construction:

● A sequence ${({u_i})}_{1\leq i \leq n }$ such that $F(u_{i+1}) < F(u_i)$ $\forall i$.

● If ${u_i}$ satisfies $\|\nabla F (u_i) \| \leq \gamma$ then ${u_i}$ is an approximation of the minimum sought.

Algorithm 1. Optimization algorithm.

1: Initialization $u_0$.

2: Solving the state equation corresponding to $u_0$ and calculation of $F$.

3: Calculation of the adjoint state and the $\nabla F$.

4: If the stop condition is satisfied, then go to step 9. Otherwise, continue with the following steps:

5: Calculation of the direction of descent $d_n$.

6: Exact search for the linear step ${\alpha}_n$ (the criterion is quadratic).

7: Construction of $u_{n+1} = u_n+ {\alpha}_n d_n$.

8: Repeat steps 2 and 3 for $u_{n+1}$ then go to step 4.

9: $u_{n+1}$ is an approximation of the desired minimum.

4.5.   Exact linear search

The descent step ${\alpha}_n$ plays a very important role in various optimization methods. The search for this coefficient then requires the resolution of a one-dimensional minimization problem. We consider the merit function $\alpha \longmapsto q(\alpha)$ representing $F(g_n+{\alpha} d_n)$ which is strictly convex. So, the optimal step ${\alpha}_n$ is the unique minimum of the merit function $q$, it satisfies the following problem:

In the case of the criterion functional used in this work, we show that $\alpha_n$ is given by:

where $z_{d_n}$ is the solution of the following problem:

4.6.   Results

We present here a series of numerical tests to illustrate the results of convergence of the penalized solutions (equations of state and problems of associated control) to solutions of non-penalized problems when the penalty parameter $\varepsilon$ tends to zero. For this, a computer code has been implemented which consists, solving the penalized control problem $({\cal P} _\varepsilon)$. We recall that, for the numerical resolution of the partial differential equations whose state and the adjoint state are solutions, the matrix of the system is obtained by the Cholesky process, and its inversion is calculated by the descent-ascent method. The boundary condition considered in the problems treated here represents a penalization of the Dirichlet condition. Since the functional considered here is quadratic, the used optimization method is the conjugate gradient combined with an exact linear search. Let us recall that for the numerical resolution, the discretization in space is given by the finite element method (FEM) and the discretization in time is given by the implicit Euler scheme. In the following we present some numerical tests illustrating the convergence of the penalized solution towards the non-penalized solution when the penalization parameter $\varepsilon$ tends to zero, first for the equation of state and then for the associated control problem.

4.6.1.   State equation

The parameter values are used to illustrate the theoretical results, but they have no physical significance. The simulations are done using MATLAB software.

In the first example, we take the domain as the square

(see Figure 1) and consider the function

Let

where

and

and ${\mathit{\boldsymbol{v}}}$ verify the assumptions

and

One can easily verify that

is a solution (Figure 2) of:

in $\Upsilon,$ $u = \phi$ on $\Pi$, $u(0) = \phi(0)$ in $\Omega$, where $\mu = 1$, and

The associated penalized problem is given by:

In Figure 3, and in order to evaluate the rate of convergence, we represent the error $(u_\varepsilon -u)$ in norm $L^2(\Upsilon)$ as a function of $\varepsilon$ for time steps $dt \in \{0.01, 0.001\}$. We can observe that the error converges linearly to zero when the penalty parameter $\varepsilon$ tends to zero. This is confirmed by the calculation of the slope, in the region far from the plateau (Figure 4). The slope is worth $0.92$ when $dt = 0.01$ and $0.90$ when $dt = 0.001$. In Figure 3, we display also the curve of the error $(u_\varepsilon -u)$ in norm $L^2(0, T, H^1(\Omega)$ for the two time steps $dt \in \{ 0.01, 0.001 \}$.

In the second example, we get a non-regular exact solution by considering the non-convex domain

as shown in Figure 5. The exact solution given here is

where $(r, \theta)$ are the polar coordinates of a point $M(x_1, x_2)$, $r$ being the distance between $M(x_1, x_2)$ and $O({1 \over 2}, {1 \over 2})$, and $\theta$ the angle between $(OA)$ and $(OM)$ where $A({1 \over 2}, 0)$.

We use the same numerical experiments as in Example 1. In order to better take into account the regularity of the solution and to obtain a better precision calculation, we consider here a mesh twice as fine as that used in the regular case ($h = 0.00625$). For the temporal discretization, we choose $dt \in \{ 0.01, 0.001 \}$. For different values of $\varepsilon$, we calculated the error $(u_\varepsilon - u)$ in norm $L^2(\Upsilon)$ and $L^2(0, T, H^1(\Omega))$, where $u$ and $u_\varepsilon$ are the solutions of Eqs (2.1) and (2.5), respectively (Figure 6).

Note that the theoretical results given in ^[4] for a pure diffusion predict a linear behavior in $\varepsilon$ for the norm error $L^2(\Upsilon)$ and an $\varepsilon^{\frac{2}{3}}$ behavior for the norm error $L^2(0, T, H^1(\Omega))$. In Figure 7, and for the same data as above, we represent the error $(u_\varepsilon - u)$ in norm $L^2(\Upsilon)$ and $L^2(0, T, H^1(\Omega)$ as a function of $\varepsilon$, where this time $u$ and $u_\varepsilon$ are the solutions of the reaction-advection-diffusion Eqs (2.1) and (2.5), respectively, where it is necessary to replace the expression of $f$ by

We observe the same error behavior as in the regular case, namely a convergence toward zero (Figure 6). We also observe a saturation of the error from $\varepsilon = 10^{-3}$. The calculated slopes are given by $0.92$ and $0.77$ for the norms $L^2$ and $L^2(H^1)$, respectively. A more advanced numerical study using higher-order time schemes (e.g., the Crank-Nicholson scheme) and quadrature formulas more precise for the integration of the advection term could give better precision on the calculation of the slopes.

4.6.2.   Control problem

In this subsection, we present some numerical results on the convergence of the optimal solution of the control problem associated with the reaction-advection-diffusion equation. For this, we consider two examples, the first concerns a very regular exact optimal solution and the second an exact solution containing a singular part. We then compare the calculated optimal solution, using the penalty technique, to the exact solution in order to evaluate the convergence rates as a function of the penalty parameter $\varepsilon$ intervening in the Robin condition penalizing the Dirichlet condition.

We start first with the case where the exact optimal solution to be approximated is given by:

on the domain

Then, let

and consider the optimal control problem

where $u$ is the solution of:

with

in $\Omega,$ $\mu = 1$,

and ${\mathit{\boldsymbol{v}}}$ the function used previously. One can easily verify that the minimum of $J$ is reached in $(u_d, g_d)$ and is equal to zero. The discretization of the control problem is the one described in Section 4.3 with a mesh step $h = 0.0125$ and a time step $dt = 0.01$. The minimum of the cost function is calculated by the gradient method combined with a tolerance of $10^{-5}$ and an exact linear search. At each iteration, the discrete algebraic systems associated with the Robin condition are inverted by the Cholesky method, knowing that the factorization is performed only once in pre-processing. For several values of the penalty parameter, $\varepsilon$, the discrete optimal state $\overline u_{\varepsilon, h}^{dt }$ and the discrete optimal control $\overline g_{\varepsilon, h}^{dt }$ were calculated. In Figure 8, we represent, on a logarithmic scale, the relative errors:

in norms $L^2(\Upsilon)$ and $L^2(0, T, H^1(\Omega))$ for the optimal state and in $L^2(\Pi)$ norm for the optimal control.

One can observe on Figure 8 a decrease of the calculated errors towards zero, which confirms the results of convergences of the optimal penalized solution (Robin) to the exact optimal solution (Dirichlet) when the penalty parameter $\varepsilon$ tends to zero.

In addition, we calculated the slopes (convergence rate) corresponding to the different curves, and for values of $\varepsilon$ far from the saturation zone, which, as explained in the examples preceding numeric concerning the equation of state, is due to the fact that the error in $h$ and $dt$ outweighs that coming from the penalty parameter, $\varepsilon$. The calculated slopes are given by $0.92$, $0.73$, and $0.94$ for the error on the optimal state in $L^2(\Upsilon)$ and $L^2(0, T, H^1(\Omega))$ norms and for the error in $L^2(\Pi)$ norm on the optimal control, respectively.

The last example concerns an exact optimal solution containing a singular part. For this, we take the non-convex domain

used below (see Figure 5). We give the functions:

and

and the same objective function $J(u, g)$ considered in the previous example, where $u$ is the solution of (4.6), and the second member $f$ is given by (4.5). The optimal solution is given by $(\overline u = u_d, \overline g = g_d)$. By carrying out the same numerical experiment as that of the preceding example, we show the various errors associated with the control of the reaction-advection-diffusion equation in the Figure 9.

The evaluation of the slope for the different errors leads to the following results: the slopes for the norms $L^2(L^2)$ and $L^2(H^1)$ on the optimal state $\overline u$ and for the norm $L^2$ on the optimal control $\overline g$ are $0.92$, $0.76$, and $0.91$, respectively. The results obtained here allow us to conclude that the penalized optimal solution (Robin) converges to the optimal solution associated with the Dirichlet condition when the penalization parameter $\varepsilon$ tends to $0$.

5.   Conclusions

We discussed the boundary optimal control problem for the reaction-advection-diffusion equation for not very regular Dirichlet boundary condition. We checked the well-posedness of the state equation. Then, we introduced the penalization technique of the boundary condition and we considered, then, the corresponding boundary optimal control problem and gave the optimality conditions. Finally, we conducted a numerical analysis of a time scheme/finite element discretization of the penalized control system. We investigated the convergence by some numerical experiments. This work should be confirmed by some theoretical results on the convergence of the optimal penalized solution towards the optimal solution of the Dirichlet problem, which is in progress.

Acknowledgments

The authors extend their appreciation to the Deanship of Scientific Research at Northern Border University, Arar, KSA for funding this research work through the project number "NBU-SAFIR-2024". The authors are also grateful to the unknown referees for the many constructive suggestions, which helped to improve the presentation of the paper.

Conflict of interest

All authors declare no conflicts of interest in this paper.