# Numerical Computation for a Kind of Time Optimal Control Problem for the Tubular Reactor System

Authors: Detang Zeng; Xin Yu; Jingfang Huang; Chunqing Tan
Journal: Mathematical Problems in Engineering
 (2018)

---

## Abstract

This paper is devoted to the study of numerical computation for a kind of time optimal control problem for the tubular reactor system. This kind of time optimal control problem is aimed at delaying the initiation timeτ of the active control as late as possible, such that the state governed by this controlled system can reach the target set at a given ending time T. To compute the time optimal control problem, we firstly approximate the original problem by finite element method and get a new approximation time optimal control problem governed by ordinary differential equations. Then, through the control parameterization method and time-scaling transformation, the approximation problem becomes an optimal parameter selection problem. Finally, we use Sequential Quadratic Program algorithm to solve the optimal parameter selection problem. A numerical simulation is given for illustration.


---

## Body

## 1. Introduction

Tubular reactors have extensive applications in the industrial production. A lot of reaction processes in chemical and biochemical engineering can be described by tubular reactor models, such as chlorine dioxide bleaching model [1]. The dynamic tubular reactors are typically described by nonlinear partial differential equations (PDEs), which include convection, reaction, and diffusion phenomena. For studies of such PDEs system, we would mention some works [2, 3]. Sometimes we regulate the reaction process artificially to improve the reactor rate or reduce the reaction time, by seeking some control strategies such as adding catalyst or placing a heating or cooling jacket. For this purpose, optimal control problems governed by such PDEs system arouse increasing attentions [4–7]. As a kind of important optimal control problem, time optimal control problem plays an important role in many fields of applications.Generally speaking, time optimal control problems can be divided into two types. The first one is to find a minimal timeT and a control belonging to some constraint set which is acted upon from beginning time 0, such that the state governed by the controlled system can arrive at a given target set in the shortest time interval [0,T]. This kind of problem is called the minimal time optimal control problem. The second one is to delay the initiation time τ of the active control as late as possible, such that the state governed by this controlled system can reach the target set at the given ending time T. This is a maximal time optimal control problem. There have been extensive researches on the first kind of time optimal control problem [8–11], but only a few works related to the second kind of time optimal control problem have been studied. For the second kind of time optimal control problem, we would mention the works [12, 13].In this paper, we shall consider the second kind of time optimal control problem with boundary controls governed by a coupled tubular reactor system, which is described by nonlinear partial differential equations and corresponds to the following chemical reaction process [14, 15]: (1)R1+R2⟶P,where R1, R2 are two reactants and P is the product of the chemical reaction. We assume that the kinetics depend on the concentrations of two reactants R1 and R2, which has the form (-r1xy,-r2xy)⊤, where ri,i=1,2, are the kinetic constants and x, y are the concentrations of the reactants R1 and R2, respectively. Moreover, we assume that the length of the reactor is l. A control u is acted upon on the one side of the tubular reactor to control the reactor process, as depicted in Figure 1. More precisely, the above controlled process can be described as follows:(2)xtz,t=β1xzzz,t-vxzz,t-r1xz,tyz,t,0<z<l,0<t≤T,ytz,t=β2yzzz,t-vyzz,t-r2xz,tyz,t,0<z<l,0<t≤T,β1xz0,t=u1t,xzl,t=0,0<t≤T,β2yz0,t=u2t,yzl,t=0,0<t≤T,xz,0=x0z,yz,0=y0z,where x0,y0∈L2(0,l) are the given initial data, T is a given positive number, β1, β2 are the dispersion coefficients, v is the superficial fluid velocity, and u=(u1,u2)⊤∈L2(0,T)×L2(0,T) is the boundary control. We denote by (x(z,t;u,x0),y(z,t;u,y0))⊤ the solution to (2) corresponding to the control u and the initial data (x0,y0)⊤.Figure 1
Tubular reactor.Letai,bi∈R for i=1,2. Define(3)S=x1,x2⊤∣ai≤xi≤bi,i=1,2,U=u,τ∣u∈L20,T×L20,T,τ∈0,T,ut=0fora.e.t∈0,τ,ut∈Sfora.e.t∈τ,T. Let αi,γi∈R for i=1,2. Define (4)K=y1,y2⊤∈L20,l×L20,l∣αi≤yiz≤γifora.e.z∈0,l,i=1,2. The second kind of time optimal control problem governed by (2) can be stated as follows: TPmaxu,τ∈Uτ∣xz,T;u,x0,yz,T;u,y0⊤∈K, where (x(z,t;u,x0),y(z,t;u,y0))⊤ is the solution to (2). We assume that, for given (x0,y0)⊤∈L2(0,l)×L2(0,l), there exists (u,τ)∈U, such that (5)x·,T;u,x0,y·,T;u,y0⊤∈K.u=(u1,u2)⊤∈L2(0,T)×L2(0,T) satisfying the above assumption is called an admissible control for the problem TP. Obviously, the above assumption implies that the admissible control set is not empty. In the problem TP, τ represents the initiation time of the active control. We call (6)τ∗=maxu,τ∈Uτ∣xz,T;u,x0,yz,T;u,y0⊤∈Kthe optimal time for the problem TP, and a control u∗=(u1∗,u2∗)⊤ satisfying (u∗,τ∗)∈U such that (x(z,T;u∗,x0),y(z,T;u∗,y0))⊤∈K is called an optimal control for the problem TP. The goal of this problem is to find a maximal τ∗ and a control u∗ such that (u∗,τ∗)∈U, and the state of (2) driven by this control can reach the given target set K at T with the shortest actively controlled interval [τ∗,T].The time optimal control problems governed by PDEs are infinite dimensional control problems and are difficult in finding out analytical solutions. Thus, the numerical approximations of time optimal control problems attract a lot of attentions. Here, we would like to mention some related works [16–18] on the computations of time optimal control problems. However, the works mentioned above focused on the first kind of time optimal control problem. To the best of our knowledge, it seems that no attempts have been made to develop the computation for the second kind of time optimal control problem governed by nonlinear PDEs systems. In this paper, we shall consider the computation for the second kind of time optimal control problem TP. Firstly, we project the time optimal control problem TP by finite element method into an approximation problem TPh. Although approximation problem TPh is governed by ordinary differential equations, it is difficult to solve directly due to the unknown time variables. Then, we translate problem TPh to an optimal parameter selection problem TPph through the control parameterization method and time-scaling transformation. Finally, we use Sequential Quadratic Program (SQP) algorithm to solve the optimal parameter selection problem.The rest of this paper is organized as follows. In Section2, we give the finite element approximation for the time optimal control problem TP. In Section 3, we provide the control parameter method and time-scaling transformation to reduce the approximation problem to an optimal parameter selection problem. In Section 4, we give the procedure of solving the optimal parameter selection problem. In Section 5, we carry out some numerical experiments to illustrate the effectiveness of our approximation method.


## 2. The Finite Element Approximation of the ProblemTP

In this section, we describe finite element approximation of the problemTP. Firstly, we introduce a standard triangulation Th of [0,l]. We decompose [0,l] into finite number of subintervals [zj-1,zj],j=1,2,…,N, where N is a positive integer and zj,j=0,1,…,N, are grid points which satisfy (7)0=z0<z1<⋯<zN-1<zN=l.h=max⁡{zj-zj-1∣j=1,…,N} is called the mesh size of the triangulation Th. Thus, we can write Th as (8)Th=Ij=zj-1,zj∣j=1,2,…,N. Corresponding to each triangulation Th, we can define a finite dimensional space as follows: (9)Sh=wh∈C0,l∣wh∣Ij∈P1Ijforj=1,2,…,N, where P1(Ij) stands for the space of all linear polynomials defined on the subinterval Ij. Obviously, the basis of Sh can be taken as {ei∣i=0,1,…,N}, where ei∈Sh,i=0,1,…,N, are continuous and piecewise linear polynomials that satisfy (10)eizj=1,i=j,0,i≠j. Thus, it is clear that Sh is a space of N+1 dimensions and Sh⊂H1(0,l).Next we construct the approximation of (2). Let Ph be the L2-projection from L2(0,l) onto Sh, defined by(11)Phv,wh=v,wh,∀v∈L20,l,wh∈Sh.Here and in what follows, (·,·) denotes the usual inner product of L2(0,l). Moreover, if there is no risk to make any confusion, we will use x(t) to denote a function x(·,t)∈L2(0,l) for a.e. t∈[0,T]. Then, the finite element approximation of (2) can be defined as follows: seek functions xh and yh satisfying that xh(t)∈Sh and yh(t)∈Sh for a.e. t∈[0,T], such that(12)∂txht,wh+β1∂zxht,∂zwh+u1twh0+v∂zxht,wh+r1xhtyht,wh=0,∀wh∈Sh,∂tyht,wh+β2∂zyht,∂zwh+u2twh0+v∂zyht,wh+r2xhtyht,wh=0,∀wh∈Sh,xh0=Phx0,yh0=Phy0.Noting that xh(t),yh(t)∈Sh, we can write(13)xhz,t=∑j=0NXjtejz,yhz,t=∑j=0NYjtejz.The initial conditions can be written as(14)xh0=Phx0z=∑j=0NXj0ejz,yh0=Phy0z=∑j=0NYj0ejz.Substituting (13) into (12) and taking wh=ei,i=0,1,…,N, we can obtain that(15)∑j=0Nei,ejX˙jt+β1∑j=0Nei′,ej′Xjt+∑j=0Nvei,ej′Xjt+u1tei0+r1∑j=0N∑k=0NXjtYktejek,ei=0,i=0,1,…,N,∑j=0Nei,ejY˙jt+β2∑j=0Nei′,ej′Yjt+∑j=0Nvei,ej′Yjt+u2tei0+r2∑j=0N∑k=0NXjtYktejek,ei=0,i=0,1,…,N,where ej′ denotes the derivative of ej with respect to the spacial variable z and X˙j and Y˙j denote the derivatives of Xj and Yj with respect to the time variable t, respectively. Moreover, it follows from (11) and (14) that(16)∑j=0Nei,ejXj0=x0,ei,i=0,1,…,N,∑j=0Nei,ejYj0=y0,ei,i=0,1,…,N.To simplify the notation, let (17)Xt=XjtN+1×1∈RN+1,Yt=YjtN+1×1∈RN+1,M=ei,ejN+1×N+1∈RN+1×N+1,K=ei′,ej′N+1×N+1∈RN+1×N+1,W=ei,ej′N+1×N+1∈RN+1×N+1,QXt,Yt=∑j=0N∑k=0NXjtYktejek,eiN+1×1=Xt⊤AiYtN+1×1∈RN+1,where Ai=[(ejek,ei)]N+1×N+1∈R(N+1)×(N+1),i=0,1,…,N.Let(18)B=ei0N+1×1∈RN+1,ζ1=x0,eiN+1×1∈RN+1,ζ2=y0,eiN+1×1∈RN+1. Then, (15) and (16) can be rewritten as(19)MX˙t+β1KXt+vWXt+Bu1t+r1QXt,Yt=0,MY˙t+β2KYt+vWYt+Bu2t+r2QXt,Yt=0,X0=M-1ζ1,Y0=M-1ζ2. Define (20)Zt=XtYt,R=-β1M-1K-vM-1W00-β2M-1K-vM-1W,F=-M-1B00-M-1B,GZt=-r1M-1QXt,Yt-r2M-1QXt,Yt. Thus, (19) is equivalent to the following equation:(21)Z˙t=Ht,Zt,ut,Z0=Z0,where H(t,Z(t),u(t))=RZ(t)+Fu(t)+G(Z(t)), and (22)Z0=M-1ζ1M-1ζ2.Now, we get a new time optimal control problemTPh governed by ordinary differential equations as follows: TPhmaxu,τ∈Uτ∣ZT∈Kh, where (23)Kh=XY∈RN+1×RN+1∣X=X1,X2,…,XN+1⊤,Y=Y1,Y2,…,YN+1⊤,α1≤Xi≤γ1,α2≤Yi≤γ2,i=1,2,…,N+1, and Z(t) is the solution to (21).


## 3. Control Parameterization Method and Time-Scaling Transformation

To solve the problemTPh numerically, which is a time optimal control problem in the finite dimensional space, we apply the classical control parametrization method [19], where the control function can be approximated by the piecewise constant functions. We subdivide the time horizon [0,T] into p subintervals [sk-1,sk],k=1,2,…,p, where p≥2 is a given integer and sk,k=0,1,…,p, satisfy(24)0=s0≤s1=τ≤s2≤⋯≤sp=T.Here each time point sk,k=1,2,…,p-1, is called a switching time. We allow the approximate control to switch at each switching time. To find the optimal switching time points, we assume that these switching points are not prefixed and are the parameters, which need to be chosen optimally in the following set: (25)Γ=ς=s1,…,sp-1⊤∈Rp-1∣0≤s1≤s2≤⋯≤sp-1≤T. Define(26)Up=σ=σ1,σ2,…,σp∈R2×p∣σk=σ1k,σ2k⊤∈R2,σ1=0,ai≤σik≤bi,k=2,…,p,i=1,2.We approximate the control function u(t) by a constant vector on each subinterval: (27)ut≈σk,fort∈sk-1,sk,k=1,…,p. Then the control function can be approximated by the following piecewise constant control of the form (28)ut≈upt=u1pt,u2pt⊤=∑k=1pσkχsk-1,skt,t∈0,T. It is clear that up(t) can switch at times sk,k=1,…,p-1. Substituting (28) into (21), we can obtain that(29)Z˙t=Ht,Zt,0,t∈0,s1,Z˙t=Ht,Zt,σk,t∈sk-1,sk,k=2,…,p,Z0=Z0.Since the switching times are unknown, (29) is not easy to solve numerically. It is also difficult to integrate (29) accurately if some of the subintervals [sk-1,sk],k=1,…,p, are very short. To overcome these difficulties brought about by the unknown switching times, we make use of the time-scaling transformation [20–22]. Let r∈[0,p] be a new time variable, where p is the number of the time subintervals. The time-scaling transformation between r∈[0,p] and t∈[0,T] can be established through the following differential equation:(30)dtrdr=wr,r∈0,p,t0=0,tp=T, where w is a piecewise constant function, which is defined by (31)wr=∑k=1pθkχk-1,kr,r∈0,p, where θk=sk-sk-1,k=1,…,p. Integrating (30) can yield that(32)tr=∑m=1rθm+θr+1r-r,r∈0,p,where ⌊r⌋ denotes the largest integer which is less than or equal to r. It follows from (32) that (33)tk=∑m=1kθm=sk,k=0,1,…,p-1. Noting that t(p)=T in (30), we have t(p)=sp. Thus(34)tk=∑m=1kθm=sk,k=0,1,…,p.Equality (34) shows that the time-scaling transformation maps r=k to the kth switching time t=sk,k=1,…,p-1. Moreover, it also follows from (34) that (35)θ1+θ2+⋯+θp=T. By taking the approximation u¯(r)≈u¯p(r)=up(t(r))=∑k=1pσkχ[k-1,k)(r) and making variable substitution Z¯(r)=Z(t(r)), from (29) and (30), we can derive that(36)Z¯˙r=θ1H¯r,Z¯r,0,r∈0,1,Z¯˙r=θkH¯r,Z¯r,σk,r∈k-1,k,k=2,…,p,Z¯0=Z0, where H¯(r,Z¯(r),σk)=H(t(r),Z(t(r)),σk). Let (37)Θp=θ=θ1,…,θp⊤∈Rp∣θi≥0,i=1,…,p,θ1+θ2+⋯+θp=T.We denote by Z¯(r∣(θ,σ)) the solution to (36) corresponding to (θ,σ)∈Θp×Up, where (38)Z¯r∣θ,σ=X¯r∣θ,σY¯r∣θ,σ.Through the control parameterization method and time-scaling transformation, we can translate the approximation problemTPh into an optimal parameter selection problem TPph as follows: TPphmaxθ,σ∈Θp×Upθ1∣Z¯p∣θ,σ∈Kh, where Z¯(r∣(θ,σ)) is the solution to (36).


## 4. Solving ProblemTPph

In the problemTPph, both θk,k=1,…,p, and σk,k=2,…,p, are decision variables, where θk,k=1,…,p, are the time parameters and σk,k=2,…,p, are the control parameters. The problem TPph is essentially nonlinear mathematical programming problem and hence can be solved numerically with standard mathematical programming algorithms such as Sequential Quadratic Program algorithm. Now we present framework of solving the problem TPph in Algorithm 1.Algorithm 1:Procedure of solving the problem TPph.
( 1 ) Input the number of the subintervals in the space interval [0,l]:N; Input the number of the subintervals in the new time horizon[0,p]:p≥2; Input initial value:x0,y0; Input initial pair:(θ,σ)∈Θp×Up; ( 2 ) Calculate corresponding matrix M,K,W,B, and Ai,i=0,1,…,N; Calculate initial valueZ-(0); ( 3 ) Solve the optimal solution (θ,σ) by making use of SQP algorithm; ( 4 ) Output the value of the optimal parameters (θ∗,σ∗).We make the following finite difference scheme to solveZ¯(r∣(θ,σ)),r∈(0,p], which is necessary in the third step of Algorithm 1. To simplify the notation, we write Z¯(r)=Z¯(r∣(θ,σ)). Firstly, we divide the time horizon [0,p] into Λp subintervals. Let Δr=1/Λ and rj=jΔr,j=0,1,…,Λp. We obtain that(39)X¯j+1Δr=X¯jΔr+ΔrθdM-1×-β1KX¯jΔr-vWX¯jΔr-Bσ1d-r1QX¯jΔr,Y¯jΔr,j=0,1,…,Λp-1,Y¯j+1Δr=Y¯jΔr+ΔrθdM-1×-β1KY¯jΔr-vWY¯jΔr-Bσ2d-r2QX¯jΔr,Y¯jΔr,j=0,1,…,Λp-1, where d=⌊j/Λ⌋+1.  X¯(jΔr) and Y¯(jΔr) at j=0 are known from the initial condition. Therefore, the computation of (39) is straightforward.Following the procedure of Algorithm1, we can solve the parameter selection problem TPph. Let (θ∗,σ∗) be the optimal parameter for the problem TPph. By making use of (32), the approximation of (u∗,τ∗) and (x,y)⊤ for the problem TP can be obtained.


## 5. Numerical Simulation and Discussion

In this section, a numerical simulation for the problemTP is presented. We take (40)β1=1,β2=1,v=1,r1=0.05,r2=0.1,T=1,l=1,x0=sin⁡πz,y0=12sin⁡πz,z∈0,l,in (2). The control bounds are taken as (41)a1=-0.2,b1=0.1,a2=-0.15,b2=0.15,and the target set bounds are (42)α1=0.75,γ1=0.9,α2=0.4,γ2=0.5.In this simulation, we divide the interval[0,1] into N equal subintervals. Thus, we have h=1/N. We choose the basis functions ei(z),i=0,1,…,N, as follows: (43)eiz=z-zi-1h,zi-1≤z<zi,-z-zi+1h,zi≤z<zi+1,0,otherwise,i=1,…,N-1e0z=-z-z1h,z0≤z<z1,0,otherwise,eNz=z-zN-1h,zN-1<z≤zN,0,otherwise.We carry out the numerical simulation experiments within the MATLAB environment (version R2016a). We use personal computer with the following configuration: Intel Core i5-7200 2.50 GHz CPU, 8.00 GB RAM, 64-bit Windows 10 Operating System.The following simulations are conducted by Algorithm1. Using the piecewise constant control parameterization method with p=2,5,10,20, we solve problem TPph for N=10,12,15,20,30,40. The optimal time τ∗ is given in Table 1. From Table 1, we can see that as N increases from 10, 12, 15 to 20, the approximation value of τ∗ varies from 0.3680,0.3682,0.3685 to 0.3686, which means finer finite element triangulation Th can lead to shorter actively controlled interval [τ∗,T]. But when N increases from 20, 30 to 40, the optimal time tends to be stable. Moreover, we see that the change of the value of p has little effect on the approximation results. The optimal controls for N=40 when p=2 and p=20 are presented in Figure 2. It can be found in Figure 2 that the optimal control for the problem TPph is bang-bang control, and the optimal control u∗ only switches one time. The corresponding states without control and with the optimal control acting on the interval (τ∗,T) for N=40 and p=20 are also presented in Figures 3 and 4, respectively.Table 1
Optimal timeτ∗ for different N and different p.
p N 10 12 15 20 30 40 2 0.3680 0.3682 0.3685 0.3686 0.3686 0.3686 5 0.3680 0.3682 0.3685 0.3686 0.3686 0.3686 10 0.3680 0.3682 0.3685 0.3686 0.3686 0.3686 20 0.3680 0.3682 0.3685 0.3686 0.3686 0.3686Figure 2
Optimal controls for differentp.
(a) (b)Figure 3
The state function without control.
(a) (b)Figure 4
The state function with optimal control acting on the interval(τ∗,T).
(a) (b)


## 6. Conclusions

In this paper, we propose a computation method for the second kind of time optimal control problem for the tubular reactors systems, which are widely used in chemical and biochemical engineering. The aim of the second kind of time optimal control problem is to find a control such that the state governed by the tubular reactors system can arrive at a given target set at the ending timeT with the shortest duration of the control action. Currently, the computation of the second kind of time optimal control problem did not attract enough attentions of the researchers. Our approach can be summarized in three steps: (1) A finite dimensional approximation problem can be obtained with finite element method; (2) through control parameterization method and time-scaling transformation, we translate the approximation problem into an optimal parameter selection problem; (3) we solve this optimal parameter selection problem by Sequential Quadratic Program algorithm. Moreover, we give an example for illustration to show the effectiveness of our approximation method.



---

*Source: 9580470-2018-04-04.xml*