Boundary Control for a Certain Class of Reaction-Advection-Diffusion System

: There are physical phenomena, involving diffusion and structural vibrations, modeled by partial differential equations (PDEs) whose solution reﬂects their spatial distribution. Systems whose dynamics evolve on an inﬁnite-dimensional Hilbert space, i.e., inﬁnite-dimensional systems, are modeled by PDEs. The aim when designing a controller for inﬁnite-dimensional systems is similar to that for ﬁnite-dimensional systems, i.e., the control system must be stable. Another common goal is to design the controller in such a way that the response of the system does not be affected by external disturbances. The controller design for ﬁnite-dimensional systems is not an easy task, so, the controller design for inﬁnite-dimensional systems is even more challenging. The backstepping control approach is a dominant methodology for boundary feedback design. In this work, we try with the backstepping design for the boundary control of a reaction-advection-diffusion (R-A-D) equation, namely, a type parabolic PDE, but with constant coefﬁcients and Neumann boundary conditions, with actuation in one of these latter. The heat equation with Neumann boundary conditions is considered as the target system. Dynamics of the open- and closed-loop solution of the PDE system are validated via numerical simulation. The MATLAB-based numerical algorithm related with the implementation of the control scheme is here included.


Introduction
There are a large number of dynamic systems that can be modeled by partial differential equations (PDEs). Many of these systems can be found in the study of mechanical, electrical, or civil engineering as well as in the study of physics and chemistry sciences since all these disciplines deal with the study of fluid flow, heat transfer, electromagnetic field, beam deflection, acoustics, and many other physical phenomena. In the study of health and social sciences, PDEs can also be used to model spread of diseases and viruses as well as population growth. One of the main features of the phenomena mentioned above is that these are listed as distributed parameters systems (DPSs), which are characterized by having two or more independent variables where the combination of one spatial variable with one temporal variable is the most common. Another feature from the DPSs is that these are infinite-dimensional systems. PDEs can be categorized as of parabolic type, elliptic type and hyperbolic type. Each one of these categories describes the behavior of different phenomena, e.g., parabolic PDEs are mainly used to describe heat transfer systems and diffusion systems, elliptic PDEs are used to describe steady state behavior for systems such as electrostatic fields, and hyperbolic PDEs are associated with vibration or wave phenomena. To get the general solution (dynamics) from a PDE system, it is necessary to set not only initial conditions but boundary conditions which can be from three types, namely, the so-called Dirichlet, Neumann, and Robin conditions [1].
From control theory, the backstepping methodology can be used to force a finite-dimensional nonlinear dynamical system, mainly modeled by ordinary differential equations (ODEs), to behave like a linear system in a new set of coordinates [2,3]. The backstepping methodology has been extended from control of finite-dimensional ODEs systems to control of infinite-dimensional PDEs systems. Control of PDEs systems can be performed in two ways, depending on the location of sensors and actuators, either in the domain dimension or at the boundary.
One of the first attempts about backstepping control design for linear parabolic PDEs was reported in [4]. The backstepping methodology applied to the boundary control of PDEs, from considering Dirichlet and Neumann boundary value problems on the control of an unstable heat equation, was introduced in [5]. In [6], from the structure for a class of one-dimensional (1-D) linear parabolic partial integro-differential equations (P(I)DEs), with one Neumann boundary condition and control signal applied through either Dirichlet and Neumann boundary actuation, the backstepping methodology was used to establish an explicit expression for the gain kernel, namely the gain kernel P(I)DE system, of the boundary control law. Furthermore, the backstepping approach was extended to the design of inverse optimal controllers. From the gain kernel, explicit solutions for the boundary state feedback laws to the control of an unstable heat equation, a heat equation with destabilizing boundary condition, a reaction-diffusion (R-D) equation with spatially varying coefficient for the reaction term, and for a solid propellant rocket model, were obtained. In [7], the backstepping methodology was translated to the design of state observers with boundary sensing. Adaptive boundary controllers for parabolic PDEs were derived in [8][9][10][11]. The adaptive approaches were extended to reaction-advection-diffusion (R-A-D) systems in higher dimensions (2-D and 3-D) and to systems with spatially varying unknown parameters [12,13]. Most of the work cited above was compiled in [12,14]. Adaptive control of linear hyperbolic PDEs following the backstepping (Volterra integral transformation) strategy can be found in [15]. Furthermore, through the use of Volterra and Fredholm integral transformations, the input-to-state stability (ISS) approach has been extended to the control of 1-D PDEs in [16].
Boundary control design for a R-A-D equation via backstepping approach seems like a challenging task. Because of the inclusion of partial derivatives in its structure, the searching for the solution of the gain kernel PDE represents a tedious work when trying with the solution for an integral equation equivalent to the kernel PDE. In our work, we focus our attention in the design of a Neumann stabilizing controller for a R-A-D equation, namely, a linear parabolic PDE, but with constant coefficients and Neumann boundary conditions, with actuation in one of these latter. A heat equation with Neumann boundary conditions is considered as target system. This work is motivated from [14], where the Neumann stabilizing controller design to this class of PDE system is left as an exercise to the reader.
The manuscript is summarized as follows. The description of functions spaces is given in Section 2. In Section 3, the R-A-D PDE with constant parameters and boundary conditions of the Neumann type is shown. The R-A-D equation is reduced through a proper change of variables to a R-D system in order to implement the backstepping methodology for boundary control. In Section 4, once that the original system has been reduced to the form of a R-D equation, a target system is proposed in order to, invoking the Lyapunov approach, guarantee exponential stability. In Section 5, the main idea of the backstepping approach, a Volterra integral transformation in order to map the original PDE system into a reduced one, namely, the target system, is applied. As a result, a new PDE system is obtained in function of the so-called kernel gain. Then, by a change of variables, the solution for the integral equation is obtained. By using the successive approximations method for integral equations, as a result, the value of the kernel gain is obtained at Section 6. Once the kernel gain value is obtained, in Section 7 we back to the original variables getting so the stabilizing controller to be implemented on the boundary actuated from the original system. In Section 8, the numerical algorithm used here, in order to implement and replicate our work, is given. Simulation results are discussed at Section 9. Conclusions are drawn at the end of the paper.

Function Spaces
Let us denote by Ω an open set of R n with boundary Γ. Let us assume that either Ω is Lipschitz or Ω is of class C r for r ≥ 1. The generic point of R n is denoted by x = {x 1 , . . . , x n }. The Lebesgue measure on R n is denoted by dx = dx 1 . . . dx n .
Let us denote by C(Ω) (resp. C k (Ω), k ∈ N or k = ∞) the space of real continuous functions on Ω (resp. the space of k times continuously differentiable functions on Ω), C(Ω) (resp. C k (Ω)) represents the space of real continuous functions on Ω (resp. the space of k times continuously differentiable functions on Ω). The spaces of real C ∞ functions on Ω with a compact support in Ω are denoted by C ∞ 0 . For p ∈ [1, ∞), L p (Ω) is the space of real functions on Ω which are L p for the Lebesgue measure. It is a Banach space for the norm For p = ∞, L ∞ (Ω) is the space of measurable and bounded real functions on Ω. It is also a Banach space for the norm u L ∞ (Ω) = ess sup x∈Ω |u(x)|.
For p = 2, L 2 (Ω) is a Hilbert space for the scalar product with the corresponding norm denoted by Let us denote by W s,p (Ω), s ∈ N, p ∈ [1, ∞], the Sobolev space of functions u in L p (Ω) whose distribution derivatives of order less than or equal to s are in L p (Ω). It is a Banach space for the norm where D i u = ∂u/∂x i , i ∈ [1, n], is the notation for the partial differential derivatives of a function u, D α u = D α 1 1 . . . D α n n u = (∂ α 1 +···+α n u)/(∂x α 1 1 . . . ∂x α n n ), α = {α 1 , . . . , α n } ∈ N n , and [α] = α 1 + · · · + α n . For p = 2, W s,2 (Ω) = H s (Ω) and this is a Hilbert space for the scalar product Let us consider the Sobolev spaces H s (Ω) which contain C ∞ (Ω) and C s (Ω). The closure of C ∞ 0 (Ω) in H s (Ω) (resp. W s,p (Ω)) is denoted by H s 0 (Ω) (resp. W s,p 0 (Ω)). In particular, which are Hilbert spaces for the scalar product For a bounded Ω, the Poincaré inequality implies that H 1 0 (Ω) is a Hilbert space for the scalar product and that the corresponding norm is equivalent to the norm induced by H 1 (Ω) [17].

Reaction-Advection-Diffusion (R-A-D) System
In order to show how to apply the backstepping methodology to the boundary control of a certain class of R-A-D system, we will illustrate the procedure that must be systematically followed to reach the goal.
Let us consider a R-A-D system of the form for a function u(x, t) defined in the spatial variable x ∈ [0, 1] and time t ∈ [0, ∞), b and λ are arbitrary constant coefficients, U (t) is the actuation signal (Neumann actuation) and Neumann boundary conditions are imposed. u xx (x, t) is the diffusion term, bu x (x, t) is the advection term, and λu(x, t) is the reaction term. First of all, let us consider the following change of variables v(x, t) = u(x, t)e (b/2)x (16) in order to eliminate the advection term u x (x, t). To do this, we calculate the temporal derivative from (16) which results Now, calculating the spatial derivatives from (16) these are given by Then, the R-A-D Equation (13) becomes with λ 0 = λ − b 2 4 . Next, the boundary conditions are calculated by deriving (16), so, which evaluated at the boundary From (14), Now, evaluating (21) In this way, we have obtained (20) as a reduced model from (13) given in the form of a R-D equation with boundary conditions (22) and (23).

Target System
At this point, we must define through a target system the desirable dynamics to be described by (20) in order to that it will behave like a stable system. A R-D system with negative magnitude to the reaction term can be a good choice as an exponentially stable one [5,12]. However, let us consider the heat equation as target system, as it is usual in boundary control of PDEs [14] due to its physical motivation [4], but with Neumann boundary conditions for a function w(x, t) defined in the spatial variable x ∈ [0, 1] and time t ∈ [0, ∞). In the following, we must prove that the target system (24)-(26) is exponentially stable in the sense of the L 2 -norm of the state w(x, t) with respect to the spatial variable limited to the unit domain. The L 2 -norm is defined by Let us consider the Lyapunov function whose time derivative is given bẏ From (24), the time derivativeV of V can be rewritten aṡ By using integration by parts we obtaiṅ from which taking into account the boundary conditions (25) and (26) resultṡ In order to continue with our analysis, the following lemma which establishes the relationship between the L 2 -norms of w(x, t) and w x (x, t) will be of utility.
Multiplying both sides of (31) by a constant we get From Lemma 1, this last identity can be written aṡ which, considering (28), can be rewritten aṡ Let us try to solve this last relation by integration, where (35) can be equivalently written as So, separating variables and considering the integral symbol we get Then, solving each integral it yields Since the natural logarithm ln(x) is defined as the inverse function to the exponential function e x , this last inequality can be expressed by i.e., from which, setting t = 0 results Additionally, we need to prove that V(x, t) ≤ V(x, 0) in order to guarantee that the Lyapunov function is bounded. So, from the definition for which is desirable thatV(x, t) be negative semidefinite, i.e., integrating both sides from such inequality, i.e., and evaluating limits, it results Consequently, or, by using (28) in terms of the L 2 -norm (27), Then, taking square root at both sides from (40) it yields where w 0 (x) stands for w(x, 0). From the analysis above, we have demonstrated that the PDE system (24)-(26) is exponentially stable in the sense of the L 2 -norm of the state w(x, t).

Backstepping Transformation
From the backstepping methodology to the bounded control of PDEs [14], the coordinate transformation, also called the Volterra integral transformation, is used to transform the system (20), (22) and (23) into the target system (24)-(26). First of all, the requirement about to find a kernel gain k(x, y), that makes that the plant behaves like the target system, must be fulfilled. It is well known that (42) is invertible, so, the smoothness of the kernels of the direct and inverse transformation in x and y establishes the equivalence of norms of v and w in both L 2 and H 1 spaces. Thus, from the properties of the heat Equations (24)-(26) can be concluded that the closed-loop system is exponentially stable in L 2 and H 1 .
Let k(x, t) be a continuous function such that its partial derivative is continuous in x ∈ [0, 1] for t ∈ [0, ∞). From the Leibniz differentiation rule [18] d dx the temporal derivative from (42) is given as which, from (20), can be written as From this last equation, applying integration by parts twice for the first term into the integral yields Then, the spatial derivatives from (42), considering again (43), result From (24), subtracting (48) and (49) from (51) yields By setting, from the integral term in (52), it follows that k xx (x, y) − k yy (x, y) = λ 0 k(x, y).
Furthermore, from (52), setting equal to zero the term we have k y (x, 0) = 0.
Moreover, from the first term in (52), by setting from (46), it can be written as Thus, from this last identity, subtracting λ 0 and integrating with the differential dx it yields Hence, in order to zeroing the right-hand side from (52), for all v(x, t), the following identities must be fulfilled. By inspection of (60) and (61) we have a hyperbolic PDE system whose solution will be the kernel gain k(x, y). We can find the kernel by converting (60) and (61) into an integral equation. To this end, let us define the following change of variables and let us denoting Next, by calculating the derivative with respect to the x and y variables from (64) we have Substituting (65)-(68) into (60) it yields By applying the change of variables (63) to the condition (61) we get Then, by applying (70) into (62) it results From all of the above, we arrive to the PDE (69) with conditions (70) and (71).
By integrating (69) respect to η, with limits from 0 to η, results Evaluating limits from the integral on the left side of (72) it yields From (70), we calculate the derivative respect to ξ to then replace it in (73) yielding Then, by integrating (74) with respect to ξ, for limits from η to ξ, we have Evaluating limits of both simple integrals in (75) yields Now, we need to express G(η, η) from this last equation in terms of an integral function. By the identity (71), we can express By applying the identity (46) into (77) results By integrating (79) respect to η, with limits from 0 to η, we get Evaluating the limits of the integral on the left side from this last equation we write and, from (70), if ξ = 0 then G(0, 0) = 0. By the identity (71), if ξ = η we can say that both map into the same domain. So, we can express (81) as given by From (74), the integral term is given as So, substituting this last equation in (82) and then expanding it we get and evaluating the limits of the first integral on the right side it results By substituting (85) into (76) it yields and adding similar terms in (86), this last expression can be written as Finally, we have arrived to the integral equation (87) which is equivalent to the PDE (60) with boundary conditions (61) and (62).

Integral Equation Solution
Our next goal is to find the solution for (87). By using the method of successive approximations, let us begin with the initial guess Here, we have to set a recursive formula for (87) that allow us to approximate the step ahead solution G n+1 (ξ, η) by using the previous solution with the initial guess as the first solution for G(ξ, η). This formula is set up as Let us denote the difference between two consecutive terms as so, Under the assumption that (89) converges to a limit, the solution G(ξ, η) can be written as or, by using (90), it can be alternatively written as Setting n = 0 and taking into account (88)-(91) we calculate We have already calculated ∆G 1 (ξ, η), so using (91) now we get At this point, it must be clear that we are obtaining the term (91) for every new value of n. So, for the case when n = 2 we only need to calculate ∆G 3 (ξ, η), this because ∆G 2 (ξ, η) must be calculated for n = 1. Hence, Now, when setting n = 3 it yields Thus, we could still evaluating n for any time because from (94)-(98) the pattern to follow is given as Then, the solution (93) can be expressed as In order to simplify (100) for software implementation, we use the first-order modified Bessel function of the first kind, i.e., By setting m = 1, from (101) we have To express (102) in the form (100), rearranging terms in (100) it can be written as Moreover, separating terms in (102) it can be written as x 2 Matching the second term from (103) and (104), i.e., it is easy to see that By the knowledge of (105), the Bessel function can be rewritten as Now, all the terms appearing in (106) must appear in (100). So, multiplying (103) by then it is possible to express (102) in the form (100) as is shown next, i.e., Taking into account (63) we get so, substituting (110) and (111) into (109) yields So, in the manner described above, the kernel gain (112), key piece in the design of the controller, has been determined. Thus, we have found the solution for k(x, y) in order to achieve the boundary control of the PDE system (20), (22) and (23).

Neumann Stabilizing Controller
The control action is applied through the boundary condition (15), i.e., u x (1, t) is actuated, so by using the first spatial derivative (50) from the Volterra integral transformation (42) and setting x = 1 we have Taking into account (26) and (42) and separating v x (1, t) from the last equation we get Considering (62) and from the knowledge of k(1, 1), simplifying (114) it can be written as It is easy to get k(1, y) from the kernel gain (112) by setting x = 1. In order to get k x (1, y), then is necessary to know the derivative properties for the Bessel function. The Bessel function must be dependant of a single variable but in (112) it is depending of two variables. Let us consider the next change of variables in order to obtain a Bessel function dependant of a single variable. So, we can rewrite (112) as k(q(x, y)) = −λ 0 xQ, with Q = q −1 I 1 (q). The derivative of a Bessel function is given by Thus, the derivative of (117) results By applying the chain rule to the second term into parentheses it yields where and Substituting (120) and (121) into (119) results and returning to the original variables x, y we get Setting x = 1 in (112) and (122), and replacing them in (115) results At this point, we need to express the controller in terms of the original variable (16). To this, we set x = 1 in (16) and (18) respectively. Substituting (123) and (124) into (125) we obtain which is the Neumann stabilizing controller, i.e., the controller to implement through (15), to the control of the R-A-D system (13)-(15).

MATLAB Code
In order to validate the dynamic response for the closed-loop system, in what follows it is described how to implement the whole control system through MATLAB [19]. To this end, we use the function pdepe. This function is able to solve initial-boundary value problems for parabolic and elliptic PDEs in the one space variable x and time t. This function solves PDEs of the form c x, t, u, ∂u ∂x By comparing (13) with (127), it can be seen that From MATLAB, the syntax for the pdepe function is as follows sol=pdepe(m,pdefun,icfun,bcfun,xlinspace,tlinspace); where m is defined in (131), pdefun is a function that defines the components of the PDE to be solved, icfun is a function that defines the initial condition, bcfun is a function that defines the boundary conditions, xlinspace is a vector with elements [x 0 , x 1 , . . . , x n ] concerning with n specific points (defined by the user) for which a solution is required for every value of tlinspace, also a vector with elements [t 0 , t 1 , . . . , t f ] for f specific points (also defined by the user). For all t 0 and x, the solution components are satisfied for all initial conditions of the form Furthermore, for all t and either x = xlinspace(1)=a or x = xlinspace(end)=b, where a represents the left boundary condition and b represents the right boundary condition, the solution components must satisfy a boundary condition of the form For the PDE (13), based on (128)-(130), pdefun function is declared as it is shown below.
function u0=icfun(x) global initial global i u0 = initial(i); i=i+1; For the boundary conditions (14)- (15), bcfun function is declared following the form given by (133). cntrllr is for the control law (actuation signal). The spatial vector x is defined below.
x0=0; xf=1; xn=40; x=linspace(x0,xf,xn+1); To define the temporal vector, it must be aware of the next explanation; the solution of a PDE has a mesh where for every spatial value and every temporal value there is a solution value, so we can imagine that for each temporal value there is a spatial vector that we can refer to it as a slice of the solution. We will use two meshes for software implementation: one for a slice of the solution which we call t and the another one for the complete solution of the PDE which we will referring it as T. We need that the slice of the solution has at least three data, so we need to define the step that all slices will have in order to get the solution from the PDE. To this end, the code used is given next.
global initial; initial=5*(1-2*sin(3*pi*x/2)); u=initial; global i; i=1; pdepe function returns the solution in a 3-D array called sol, where sol(i,j,k) approximates the k-th component of the solution u k evaluated at t(i) and x(j). The size of sol is ttn-by-xn-by-initial and the command to extract the solution of the PDE for the defined size is as follows. u = sol(:,:,1); ND=5; tg=zeros(1,ttn/ND+1); ug=zeros(ttn/ND+1,xn+1); ug(1,:)=initial; A for cycle is implemented to compute the whole solution of the PDE system. The code shown at the bottom is used to plot the solution.

Simulation Results
In this section, we exhibit simulation results for both open-and closed-loop dynamics of the PDE system (13)- (15). The aim is to show that the Neumann stabilizing controller (126) applied at the boundary condition (15) achieves boundary control of the R-A-D system. The open-loop dynamic behavior is that for which U (t) = 0. The parameters for the PDE (13), and for which the code was executed, are λ = 30 and b = 10 with initial condition u 0 (x) = 5 1 − 2 sin( 3πx 2 ) . Figure 1 shows the open-loop dynamic behavior of the R-A-D system from which it can be seen that the system is unstable. From Figure 2, it can be seen that the closed-loop response from the boundary control of the PDE system the controller makes of this PDE system a stable one. Figure 3 shows the open-loop response behavior at the right boundary for u x (1, t) = 0. Figure 4 shows the control effort when applying the Neumann stabilizing controller (126) at the right boundary of the R-A-D system in order to implement the closed-loop system. From Figure 4, it can be seen that the control effort is bounded.

Conclusions
In this work, a Neumann stabilizing controller to the boundary control of a certain class of reaction-advection-diffusion (R-A-D) system, with constant parameters to the PDE and Neumann boundary conditions, via backstepping approach is designed. The performance of the stabilizing controller is validated via numerical simulation. Simulation results show that the control goal has been achieved. Furthermore, our numerical algorithm based on MATLAB is provided in order to replicate the simulation results. It should be noticed that our numerical algorithm also has been validated for others problems about boundary control of PDEs but with the restriction that the PDEs must be of parabolic or elliptic type. As future work, from the R-A-D system used here, the approaches proposed in [20][21][22][23] could be explored.