Collocation Method for Optimal Control of a Fractional Distributed System

Abstract

In this paper, a collocation method based on the Jacobi polynomial is proposed for a class of optimal-control problems of a fractional distributed system. By using the Lagrange multiplier technique and fractional variational principle, the stated problem is reduced to a system of fractional partial differential equations about control and state functions. The uniqueness of this fractional coupled system is discussed. For spatial second-order derivatives, the proposed method takes advantage of Jacobi polynomials with different parameters to approximate solutions. For a temporal fractional derivative in the Caputo sense, choosing appropriate basis functions allows the collocation method to be implemented easily and efficiently. Exponential convergence is verified numerically under continuous initial conditions. As a particular example, the relation between the state function and the order of the fractional derivative is analyzed with a discontinuous initial condition. Moreover, the numerical results show that the integration of the state function will decay as the order of the fractional derivative decreases.

1. Introduction

Optimal-control problems minimize a function (or a functional) with the states and control inputs of the system over a set of admissible control functions [1], which arises in many scientific and engineering problems, such as aeronautics [2] and economics [3]. In [4], by applying optimal control theory, the dynamics of a basic oncolytic virotherapy model was studied. The constrained dynamics of the optimal-control problem may be divided into two major classes: an ordinary differential equation (ODE) and a partial differential equation (PDE). The fractional differential equation usually came up by replacing the conventional derivatives in ODE/PDE with some sorts of fractional derivatives, whose order of derivative turns out to be some real numbers or even complex numbers. In recent years, with the development of fractional calculus, it has been verified that fractional differential equations model dynamical systems and physical processes more accurately than classic ODEs and PDEs do, and fractional controllers perform better than integer order controllers as well. In this paper, we shall investigate a class of fractional optimal-control problems. The dynamics of the system are described by a partial differential equation with a Caputo fractional derivative on a temporal variable and a second-order derivative on a spatial variable, comprising what is called a fractional optimal-control problem of a distributed system. Although integer-order optimal-control problems and fractional variational problems regarding a single variable have been extensively discussed [5,6,7,8,9,10,11], the optimal-control problems related to the fractional distributed system are not systematically developed yet. Therefore, it is meaningful to further study this topic.

As a matter of course, the above-mentioned optimal-control problems have to be solved to examine their further features. However, only few of them can be handled via analytical techniques. Therefore, tremendous efforts have been made regarding numerical methods. By means of a fractional variational principle and the Lagrange multiplier technique, the underlying fractional optimal-control problem can be converted to a system of partial differential equations with temporal fractional derivatives. As one of the pioneering works, two classes of the fractional variational problem with the Riemann–Liouville derivative are considered in [12], where the necessary condition (i.e., the Euler–Lagrange equation) for both cases are analytical deduced. In [13], a general formulation for a class of fractional optimal-control problems was derived by using techniques of calculus of variations, the Lagrange multiplier, and the formula for fractional integration by parts. Based on the foundations mentioned above, several numerical methods have been applied to solve fractional optimal-control problems. In [14], the dynamic system of fractional optimal-control problems are considered as fractional partial differential equations with a Caputo derivative, and numerical solutions of this problem are obtained by means of eigenfunctions. In [15], the modified Grünwald–Letnikov approach is employed to establish a central difference numerical scheme for both time-invariant and time-variant cases. In [16], the fractional optimal-control problem is converted into a general integer order problem by using Oustaloup’s approximation for the fractional derivative operator. This method is allowed to solve a more complex problem of a fractional-free, final-time optimal-control problem. In [17,18], the considered fractional optimal-control problems are both reduced to a system of algebraic equations by using the normalized Legendre orthogonal basis. Nevertheless, in [17], the fractional integration and multiplication are approximated directly, while in [18], the fractional optimal-control problems are transformed into an equivalent variational problem. In [19], the fractional optimal-control problem of a distributed system is solved in cylindrical coordinates. Applying the method of separation of variables and eigenfunctions, an equivalent problem in terms of generalized state and control variables is defined, and it is computed by the Grünwald–Letnikov approach. Overall, numerical methods based on approximation by smooth polynomials are perfectly suitable for optimal-control problems defined on a bounded domain with regular boundary conditions. For instance, see [20,21,22]. A hybrid meshless method for fractional distributed optimal control was presented in [23].

Motivated by the above results, a collocation method with the Jacobi polynomial is established to solve the fractional optimal-control problem. It is well-known that the collocation method is a global numerical method with high accuracy, which was firstly applied by Slater and Kantorovic in 1934 [24]. The computation of this method at any given point depends on the information of the entire domain. Naturally, it has a great superiority in solving fractional problems since the fractional derivative we used here is defined in the Caputo sense and is non-local. In [25], a spectral Jacobi-collocation approximation for fractional integro-differential equations of the Volterra type is proposed, and the error of this method decaying exponentially is verified theoretically. In [26], the truncated Bessel series is used in a collocation scheme for solving a fractional optimal-control problem with a nonlinear fractional two-point boundary value problem. A local meshless collocation algorithm is applied to approximate the time fractional evolution model in [27]. In [28], by applying shifted Jacobi polynomials, the fractal-fractional derivative in the Atangana–Riemann–Liouville sense and the fractional derivatives in the Caputo and Atangan–Baleau–Caputo concepts have been used to study the optimal-control problem of the advection–diffusion–reaction equation.

In this paper, we will transform the fractional optimal-control problems into a system of fractional partial differential equations through the fractional variational principle, then solve the resulting system by the Jacobi collocation method. The numerical results under both continuous and discontinuous initial conditions will be provided. In the case of the continuous initial condition, the exponential convergence and high accuracy can be arrived with different Jacobi parameters. The influence of the fractional order on the state function will be analyzed in the case of the problem with a discontinuous initial condition. The remainder of this article is organized as follows: In Section 2, we will introduce some necessary preliminaries of fractional calculus to make this paper self-contained. In Section 3, the formulation of the fractional optimal-control problem is given and basic property is discussed. In Section 4, the derivation of the collocation method for solving the fractional optimal-control problem is performed. A numerical example is given in Section 5, and our conclusion is finally drawn in Section 6.

2. Mathematical Preparation

In this section, we introduce two frequently discussed fractional derivatives and their fundamental properties, which are useful in the analysis of the model and computation in what follows. We begin by recalling the Riemann–Liouville fractional integral.

Definition 1

([29]). Let $f\left(x\right)\in L^1\left([a,b]\right)$, $s>0$ be a finite and real number. Then,

$$\left({}_a{I}_x^{s}f\right)\left(x\right)=\frac{1}{\Gamma \left(s\right)}{\int }_a^x{(x-t)}^{s-1}f\left(t\right)dt,$$

(1)

and

$$\left({}_x{I}_b^{s}f\right)\left(x\right)=\frac{1}{\Gamma \left(s\right)}{\int }_x^b{(t-x)}^{s-1}f\left(t\right)dt,$$

(2)

are named as left-sided Riemann–Liouville fractional integral and right-sided Riemann–Liouville fractional integral, respectively.

Definition 2

([29]). Let $f\left(x\right)\in AC\left(\right[a,b\left]\right)$ and ${}_a{I}_t^{n-s}f,{{}_{t}I}_b^{n-s}f\in {\mathcal{C}}^n\left([a,b]\right)$, $n-1<s<n$, n is a natural number. Then,

$$\left({}_a{D}_x^{s}f\right)\left(x\right)=\frac{1}{\Gamma (n-s)}\frac{d^n}{d{x}^n}{\int }_a^x{(x-t)}^{n-s-1}f\left(t\right)dt,$$

(3)

and

$$\left({}_x{D}_b^{s}f\right)\left(x\right)=\frac{{(-1)}^n}{\Gamma (n-s)}\frac{d^n}{d{t}^n}{\int }_x^b{(t-x)}^{n-s-1}f\left(t\right)dt,$$

(4)

are named as the left-sided Riemann–Liouville fractional derivative and right-sided Riemann–Liouville fractional derivative, respectively.

Definition 3

([29]). Let $f\left(x\right)\in AC\left(\right[a,b\left]\right)$, $n-1<s<n$, n is a natural number. Then,

$$\left({{}_a^{c}D}_x^{s}f\right)\left(x\right)=\frac{1}{\Gamma (n-s)}{\int }_a^x{(x-t)}^{n-s-1}\left\{\frac{d^n}{d{t}^n}f\left(t\right)\right\}dt,$$

(5)

and

$$\left({{}_x^{c}D}_b^{s}f\right)\left(x\right)=\frac{{(-1)}^n}{\Gamma (n-s)}{\int }_x^b{(t-x)}^{n-s-1}\left\{\frac{d^n}{d{t}^n}f\left(t\right)\right\}dt,$$

(6)

are named as the left-sided Caputo fractional derivative and right-sided Caputo fractional derivative, respectively.

Theorem 1

([29]). Let $0<s<1$. Then, Riemann–Liouville and Caputo fractional derivatives satisfy the following relation:

$$\left({}_x{D}_b^{s}f\right)\left(x\right)=\left({{}_x^{c}D}_b^{s}f\right)\left(x\right)+f\left(b\right)\frac{{(b-x)}^{-s}}{\Gamma (1-s)}.$$

(7)

Based on Theorem 1, the Caputo fractional derivative is equivalent with the Riemann–Liouville fractional derivative for the continuously differentiable functions satisfying Dirichlet boundary condition $f\left(b\right)=0$. More details on fractional derivatives and their properties can be found in monographes [29,30,31].

3. Analysis of Fractional Optimal-Control Problem

In this part, we focus on the formulation of a class of optimal problems with a distributed system involving the Caputo fractional derivative. The details are described as follows.

3.1. Formulation of Model

Let $0\le t\le 1$, $0\le x\le L$, and $Q,R$ be two positive parameters of the system. Functions z and u depending on t and x (i.e., time and position) simultaneously represent some performance indices. Notice that the time domain may not necessary be $[0,1]$, and it can be generalized to $[0,T]$, $T>0$. Our major task is: finding an optimal control input signal $u\left(x,t\right)$, which minimizes the cost functional

$$J={\int }_0^1{\int }_0^{L}Q{z}^2\left(x,t\right)+R{u}^2\left(x,t\right)dxdt,$$

(8)

subject to the system dynamic constraints

$${}_0^c{D}_t^{s}z\left(x,t\right)=\frac{{\partial }^{2}z\left(x,t\right)}{\partial x^2}+u(x,t),0<s<1,$$

(9)

the initial condition

$$z\left(x,0\right)=z_0\left(x\right),0\le x\le L,$$

(10)

and the boundary condition

$$\frac{\partial z\left(0,t\right)}{\partial x}=\frac{\partial z\left(L,t\right)}{\partial x}=0,0<t<1,$$

(11)

where s is the order of Caputo fractional derivative, $z\left(x,t\right)$ and $u\left(x,t\right)$ are the state and control functions, respectively. In [14], models (8)–(11) was discussed.

3.2. Fractional Variational Principle and Properties of Model

In order to deduce the necessary condition for problem (8)–(11), we impose a modified performance index by using the Lagrange multiplier $\lambda$ as

$$\tilde{J}={\int }_0^1{\int }_0^{L}Q{z}^2(x,t)+R{u}^2(x,t)+\lambda \left[{}_0^c{D}_t^{s}z(x,t)-\frac{{\partial }^{2}z(x,t)}{\partial x^2}-u(x,t)\right]dxdt.$$

(12)

Assuming that $z^{*}(x,t)$ and $u^{*}(x,t)$ are the desired optimum solutions of (8)–(11), define

$$\begin{array}{c}\hfill z(x,t)=z^{*}(x,t)+{\epsilon }_1{\eta }_1(x,t),{\epsilon }_1\in \mathbb{R},\end{array}$$

(13)

$$\begin{array}{c}\hfill u(x,t)=u^{*}(x,t)+{\epsilon }_2{\eta }_2(x,t),{\epsilon }_2\in \mathbb{R},\end{array}$$

(14)

where ${\epsilon }_1$, ${\epsilon }_2$ are perturbation parameters, and ${\eta }_1(x,t)$ and ${\eta }_2(x,t)$ are arbitrary test functions satisfying ${\eta }_1(x,0)={\eta }_{1x}(0,t)={\eta }_{1x}(L,t)=0$.

Substitute Equations (13) and (14) into functional (12), then $\tilde{J}$ becomes a scale function of ${\epsilon }_1$ and ${\epsilon }_2$. Obviously, $\tilde{J}$ has an extremum at ${\epsilon }_1={\epsilon }_2=0$. Differentiating $\tilde{J}$ with respect to ${\epsilon }_1,{\epsilon }_2$ and $\lambda$, respectively, we have

$$\begin{array}{ccc}\hfill \frac{\partial \tilde{J}}{\partial {\epsilon }_1}& =& {\int }_0^1{\int }_0^L\left[2Qz(x,t){\eta }_1(x,t)+\lambda {}_0^c{D}_t^s{\eta }_1(x,t)-\lambda \frac{{\partial }^2{\eta }_1(x,t)}{\partial x^2}\right]dxdt,\hfill \end{array}$$

(15)

$$\begin{array}{ccc}\hfill \frac{\partial \tilde{J}}{\partial {\epsilon }_2}& =& {\int }_0^1{\int }_0^L\left[2Ru(x,t){\eta }_2(x,t)-\lambda {\eta }_2(x,t)\right]dxdt,\hfill \end{array}$$

(16)

$$\begin{array}{ccc}\hfill \frac{\partial \tilde{J}}{\partial \lambda }& =& {\int }_0^1{\int }_0^L\left[{}_0^c{D}_t^{s}z(x,t)-\frac{{\partial }^{2}z(x,t)}{\partial x^2}-u(x,t)\right]dxdt.\hfill \end{array}$$

(17)

Applying fractional integration by parts formulae (see [29] (Page 76, Lemma 2.7)) to (15), we obtain

$$\begin{array}{ccc}\hfill {\int }_0^1{\int }_0^L\lambda \left[{{}_0^{c}D}_t^s{\eta }_1(x,t)-\frac{{\partial }^2{\eta }_1(x,t)}{\partial x^2}\right]dxdt={\int }_0^1{\int }_0^L{\eta }_1(x,t)\left[{{}_{t}D}_1^s\lambda -\frac{{\partial }^2\lambda (x,t)}{\partial x^2}\right]dxdt& & \\ \hfill +{\int }_0^1\left[{\lambda }_x(L,t){\eta }_1(L,t)-{\lambda }_x(0,t){\eta }_1(0,t)\right]dt+{\int }_0^L{\eta }_1{(x,0)}_t{I}_1^s\lambda (x,0)dx.& \end{array}$$

To minimize $\tilde{J}$ (the same as to minimize J), it is required that results of (15)–(17) are zeros when ${\epsilon }_1={\epsilon }_2=0$. Then, by the arbitrariness of $\lambda$, one could assume ${\lambda }_x(0,t)={\lambda }_x(L,t)=0$, $\lambda \left(x,1\right)=0$ to have

$$\begin{array}{ccc}& & {}_t^c{D}_1^s\lambda (x,t)=\frac{{\partial }^2\lambda (x,t)}{\partial x^2}-2Q{z}^{*}(x,t),\hfill \end{array}$$

(18)

$$\begin{array}{ccc}& & \lambda (x,t)=2R{u}^{*}(x,t),\hfill \end{array}$$

(19)

$$\begin{array}{ccc}& & {}_0^c{D}_t^s{z}^{*}(x,t)=\frac{{\partial }^2{z}^{*}(x,t)}{\partial x^2}+u^{*}(x,t),\hfill \end{array}$$

(20)

since ${\eta }_1(x,t)$ and ${\eta }_2(x,t)$ are arbitrary.

Substituting (19) into (18), it is easy to find that the desired optimal solutions $u^{*}(x,t),z^{*}(x,t)$ satisfy the following coupled equations:

$$\begin{array}{ccc}& & {}_t^c{D}_1^{s}u(x,t)=\frac{{\partial }^{2}u(x,t)}{\partial x^2}-\frac{Q}{R}z(x,t),\hfill \end{array}$$

(21)

$$\begin{array}{ccc}& & {}_0^c{D}_t^{s}z(x,t)=\frac{{\partial }^{2}z(x,t)}{\partial x^2}+u(x,t),\hfill \end{array}$$

(22)

with the initial and boundary conditions

$$\begin{array}{ccc}\hfill z\left(x,0\right)=z_0\left(x\right),& & u\left(x,1\right)=0,\hfill \\ \hfill \frac{\partial z\left(0,t\right)}{\partial x}=\frac{\partial z\left(L,t\right)}{\partial x}=0,& & \frac{\partial u\left(0,t\right)}{\partial x}=\frac{\partial u\left(L,t\right)}{\partial x}=0.\hfill \end{array}$$

Now, the problem (8)–(11) has been converted to a fractional coupled system (21) and (22).

Next, we discuss the uniqueness of system (21) and (22) because of the fact that existence result directly follows from the linearity of system itself (see [29] (Chapter 6)). Suppose that $(u_1,z_1)$ and $(u_2,z_2)$ are two pairs of solutions of (21) and (22); then, $\tilde{u}:=u_1-u_2$ and $\tilde{z}:=z_1-z_2$ must be the solutions of (21) and (22) by superposition principle. Therefore,

$$\begin{array}{ccc}& & {}_t^c{D}_1^s\tilde{u}(x,t)=\frac{{\partial }^2\tilde{u}(x,t)}{\partial x^2}-\frac{Q}{R}\tilde{z}(x,t),\hfill \end{array}$$

(23)

$$\begin{array}{ccc}& & {}_0^c{D}_t^s\tilde{z}(x,t)=\frac{{\partial }^2\tilde{z}(x,t)}{\partial x^2}+\tilde{u}(x,t).\hfill \end{array}$$

(24)

Multiplying $\tilde{u}$ and $\tilde{z}$ to (23) and (24), respectively, then subtracting one to another, then integrating from 0 to L with respect to x, we have

$$R{\int }_0^L{\tilde{u}}_t^c{D}_1^s\tilde{u}dx+Q{\int }_0^L{\tilde{z}}_0^c{D}_t^s\tilde{z}dx=Q{\int }_0^L\tilde{z}{\tilde{z}}_{xx}dx+R{\int }_0^L\tilde{u}{\tilde{u}}_{xx}dx.$$

(25)

In view of

$$\begin{array}{ccc}& & Q{\int }_0^L\tilde{z}{\tilde{z}}_{xx}dx+R{\int }_0^L\tilde{u}{\tilde{u}}_{xx}dx\hfill \\ & =& Q\left({\left.\tilde{z}{\tilde{z}}_x\right|}_0^L-{\int }_0^L{\tilde{z}}_x^{2}dx\right)+R\left({\left.\tilde{u}{\tilde{u}}_x\right|}_0^L-{\int }_0^L{\tilde{u}}_x^{2}dx\right)\hfill \\ & \le & 0,\hfill \end{array}$$

we conclude that either one of the following inequalities hold:

$$R{\int }_0^L\tilde{u}·\left({}_t^c{D}_1^s\tilde{u}\right)dx\le 0\mathrm{or}Q{\int }_0^L\tilde{z}·\left({}_0^c{D}_t^s\tilde{z}\right)dx\le 0.$$

(26)

If the former case is true, we assume that $\tilde{u}$ is approximated by smooth polynomials $\phi \left(x\right)$ and $\varphi \left(t\right)$, which is $\tilde{u}=\phi \left(x\right)\varphi \left(t\right)$. Then,

$$\begin{array}{ccc}\hfill {\int }_0^L\tilde{u}·\left({}_t^c{D}_1^s\tilde{u}\right)dx& =& {\int }_0^L\phi \left(x\right)\varphi \left(t\right)·\phi \left(x\right)\left({}_t^c{D}_1^s\varphi \left(t\right)\right)dx\hfill \\ & =& \varphi \left(t\right)·\left({}_t^c{D}_1^s\varphi \left(t\right)\right){\int }_0^L{\phi }^2\left(x\right)dx.\hfill \end{array}$$

(27)

Furthermore, notice that

$$\begin{array}{ccc}\hfill {\int }_0^1\varphi \left(t\right)·\left({}_t^c{D}_1^s\varphi \left(t\right)\right){\int }_0^L{\phi }^2\left(x\right)dxdt& =& {\int }_0^1\varphi \left(t\right)·\left({}_t^c{D}_1^s\varphi \left(t\right)\right)dt·{\int }_0^L{\phi }^2\left(x\right)dx\hfill \\ & =& {\int }_0^1\left({}_0{D}_t^s\varphi \left(t\right)\right)·\varphi \left(t\right)dt·{\int }_0^L{\phi }^2\left(x\right)dx\hfill \\ & =& {\int }_0^1{\left({}_0{D}_t^{s/2}\varphi \left(t\right)\right)}^{2}dt·{\int }_0^L{\phi }^2\left(x\right)dx\hfill \\ & \ge & 0,\hfill \end{array}$$

which implies that the right hand side of (27) is non-negative. Combining (26), one can immediately know that

$${\int }_0^L\tilde{u}·\left({}_t^c{D}_1^s\tilde{u}\right)dx\equiv 0.$$

(28)

That is to say, $\tilde{u}$ is zero. The other argument for $\tilde{z}$ is similar. Thus, the pair of solutions of problem (21) and (22) is unique.

In order to observe the property of state function $z(x,t)$, we integrate both sides of (22) from 0 to L,

$${\int }_0^L{}_0^c{D}_t^{s}z(x,t)dx={\int }_0^L\frac{{\partial }^{2}z(x,t)}{\partial x^2}dx+{\int }_0^{L}u(x,t)dx.$$

(29)

Obviously, the first term on the right side is zero according to boundary conditions. $u(x,t)$ is bounded because the problem is well-posed. That is, there exists a real number M satisfying $\zeta \left(t\right)={\int }_0^{L}u(x,t)dx<M$ for any $t\in [0,1]$. Regarding the left hand side, it is easy to deduce that there exists $t_l^{*}\in [0,1]$, such that

$${\int }_0^L{}_0^c{D}_t^{s}z(x,t)dx=\frac{t^{1-s}}{\Gamma (2-s)}{\int }_0^L{z}_t(x,t_l^{*})dx.$$

(30)

Equations (29) and (30) lead to the relation between the change rate of energy of $z(x,t)$ and time t as

$${\int }_0^{L}u(x,t)dx=\frac{t^{1-s}}{\Gamma \left(2-s\right)}{\int }_0^L{z}_t\left(x,t_l^{*}\right)dx.$$

(31)

On account of the fact that Equations (21) and (22) are coupled, it is easy to arrive at the following identity:

$${\int }_0^{L}z(x,t)dx=-\frac{{\left(1-t\right)}^{1-s}}{\Gamma (2-s)}{\int }_0^L{u}_t(x,t_r^{*})dx.$$

(32)

Combining (31) and (32), we obtain

$$\Phi (t,s):={\int }_0^{L}z(x,t)dx=C·\frac{{(1-t)}^{1-s}}{\Gamma (1-s)\Gamma (2-s)},0<s,t<1,$$

(33)

where $C=-t_r^{*-s}{\int }_0^L{z}_t(x,t_l^{*})dx\ge 0$. This indicates that the volume of $z(x,t)$ in the temporal direction will decay as the order s decreases, due to the non-negativity assumption of $z(x,0)$ and the diffusion property of the considered system. Conversely, if $z_t$ is positive for each fixed x, then solution $z(x,t)$ will blow up eventually, which is physically impossible.

4. Numerical Algorithm

In this section, a collocation method is proposed to the solve fractional optimal control problem (8)–(11). Firstly, we recall the definition and properties of the Jacobi polynomials, occasionally called hypergeometric polynomials, which are used in our collocation method. From now on, we use symbols $\alpha$ and $\beta$ to parameterize the Jacobi polynomials and let s, $0<s<1$ be the order of fractional derivative. Notice that when s approaches one, the corresponding fractional derivatives reduce to conventional integer-order derivatives.

Recently, Jacobi polynomial for solving differential equations has become increasingly popular because of the high accuracy of interpolation for ODE/PDE with certain regularity [32,33]. The well-known Jacobi polynomial $J_n^{\alpha ,\beta }$ with indices $\alpha ,\beta >-1$ of degree n are defined on interval $[-1,1]$ and can be expressed as

$$J_n^{\alpha ,\beta }\left(\xi \right)=\frac{{\left(-1\right)}^n}{2^n·n!}{(1-\xi )}^{-\alpha }{(1+\xi )}^{-\beta }\frac{d^n}{d{\xi }^n}\left[{(1-\xi )}^{\alpha +n}{(1+\xi )}^{\beta +n}\right],-1\le \xi \le 1,$$

(34)

which is a solution of the second order linear homogeneous differential equation

$$(1-{\xi }^2)\frac{d^{2}Y\left(\xi \right)}{d{\xi }^2}+\left[\beta -\alpha -(\alpha +\beta +2)\xi \right]\frac{dY\left(\xi \right)}{d\xi }+n(n+\alpha +\beta +1)Y\left(\xi \right)=0.$$

In order to use these polynomials on the interval $[0,L]$, we take the change of variable $x=\frac{L}{2}(\xi +1)$. Therefore, $0\le x\le L$ is due to the range of $\xi$. The derivative of $J_n^{\alpha ,\beta }$ on the interval $[0,L]$ can be represented by

$$\frac{d^m}{d{t}^m}{J}_n^{\alpha ,\beta }\left(\frac{2x}{L}-1\right)=\frac{\Gamma \left(n+m+\alpha +\beta +1\right)}{\Gamma \left(n+\alpha +\beta +1\right)L^m}{J}_{n-m}^{\alpha +m,\beta +m}\left(\frac{2x}{L}-1\right),$$

(35)

where $m=1,2,3,\cdots$.

One advantage of Jacobi polynomial is that its fractional derivatives still have a recurrence form. For special cases, the Riemann–Liouville fractional derivative of $J_n^{\alpha ,\beta }$ can be presented in the form of the Jacobi polynomial. Let $0<\alpha <1$, $\beta =0$; we have

$${}_0{D}_x^{\alpha }{J}_n^{-\alpha ,0}\left(\frac{2x}{L}-1\right)=\frac{\Gamma (n+1)x^{-\alpha }}{\Gamma (n-\alpha +1)}{J}_n^{0,-\alpha }\left(\frac{2x}{L}-1\right),$$

(36)

and when $\alpha =0$, $0<\beta <1$, it reads that

$${}_x{D}_1^{\beta }{J}_n^{0,-\beta }\left(\frac{2x}{L}-1\right)=\frac{\Gamma (n+1){\left(L-x\right)}^{-\beta }}{\Gamma (n-\beta +1)}{J}_n^{-\beta ,0}\left(\frac{2x}{L}-1\right).$$

(37)

Next, we present the basic idea of the collocation method with Jacobi polynomials. Denoting sets $\left\{x_j^{\alpha ,\beta }\right\},0\le j\le M$ and $\left\{t_j^s\right\},0\le j\le N$ as the nodes of Jacobi–Gauss–Lobatto (JGL) quadratures on the interval $[-1,1]$, we set the JGL points $\left\{x_{L,j}^{\alpha ,\beta }\right\},0\le j\le M$ and $\left\{t_{1,j}^s\right\},0\le j\le N$ on $[0,L]\times [0,1]$ as

$$\begin{array}{ccc}\hfill x_{L,j}^{\alpha ,\beta }& =& \frac{L}{2}\left(x_j^{\alpha ,\beta }+1\right),0\le j\le M,\hfill \\ \hfill t_{1,j}^s& =& \frac{1}{2}\left(t_j^s+1\right),0\le j\le N.\hfill \end{array}$$

Let $\{{\psi }_i\left(t\right),{\varphi }_j\left(x\right)\}$, $\{\widehat{{\psi }_i}\left(t\right),{\varphi }_j\left(x\right)\}$ be the pairs of Jacobi basis polynomials; the approximations for solutions $z_N$ and $u_N$ can be expanded as

$$\begin{array}{c}\hfill z_N(x,t)=\sum _{i=0}^N\sum _{j=0}^M{a}_{ij}{\varphi }_j\left(x\right){\psi }_i\left(t\right),\end{array}$$

(38)

$$\begin{array}{c}\hfill u_N(x,t)=\sum _{i=0}^N\sum _{j=0}^M{b}_{ij}{\varphi }_j\left(x\right){\widehat{\psi }}_i\left(t\right),\end{array}$$

(39)

by using the method of separation of variables. Here, $\{{\psi }_i\left(t\right),{\varphi }_j\left(x\right)\}$ and $\{\widehat{{\psi }_i}\left(t\right),{\varphi }_j\left(x\right)\}$ are different Jacobi basis polynomials, as follows:

$$\begin{array}{ccc}\hfill {\varphi }_j\left(x\right)& =& J_j^{\alpha ,\beta }\left(\frac{2x}{L}-1\right),\hfill \\ \hfill {\psi }_i\left(t\right)& =& J_i^{-s,0}\left(2t-1\right),\hfill \\ \hfill {\widehat{\psi }}_i\left(t\right)& =& J_i^{0,-s}\left(2t-1\right).\hfill \end{array}$$

Inserting (38) and (39) into (21) and (22) leads to

$$\begin{array}{ccc}& & \sum _{i=0}^N\sum _{j=0}^M{b}_{ij}{\varphi }_j\left(x\right)·{{}_{t}D}_1^s{\widehat{\psi }}_i\left(t\right)=\sum _{i=0}^N\sum _{j=0}^M{b}_{ij}{\widehat{\psi }}_i\left(t\right)\frac{{\partial }^2{\varphi }_j\left(x\right)}{\partial x^2}-\frac{Q}{R}\sum _{i=0}^N\sum _{j=0}^M{a}_{ij}{\varphi }_j\left(x\right){\psi }_i\left(t\right),\hfill \\ & & \sum _{i=0}^N\sum _{j=0}^M{a}_{ij}{\varphi }_j\left(x\right)·{{}_0^{c}D}_t^s{\psi }_i\left(t\right)=\sum _{i=0}^N\sum _{j=0}^M{a}_{ij}{\psi }_i\left(t\right)\frac{{\partial }^2{\varphi }_j\left(x\right)}{\partial x^2}+\sum _{i=0}^N\sum _{j=0}^M{b}_{ij}{\varphi }_j\left(x\right){\widehat{\psi }}_i\left(t\right).\hfill \end{array}$$

Combining the relation (7) and the Riemann–Liouville fractional derivative (36) and (37), we obtain

$${}_0^c{D}_t^s{J}_i^{-s,0}\left(2t-1\right)=\frac{\Gamma (i+1)t^{-s}}{\Gamma (i-s+1)}{J}_i^{0,-s}\left(2t-1\right)-\frac{J_i^{-s,0}\left(-1\right)t^{-s}}{\Gamma (1-s)},$$

(40)

and

$${}_t^c{D}_1^s{J}_i^{0,-s}\left(2t-1\right)=\frac{\Gamma (i+1){(1-t)}^{-s}}{\Gamma (i-s+1)}{J}_i^{-s,0}\left(2t-1\right)-\frac{J_i^{0,-s}\left(1\right){\left(1-t\right)}^{-s}}{\Gamma (1-s)}.$$

(41)

The second-order derivative of Jacobi polynomial can be expressed as

$$\frac{{\partial }^2{\varphi }_j\left(x\right)}{\partial x^2}=\frac{\Gamma (j+3+\alpha +\beta )}{\Gamma (j+1+\alpha +\beta )L^2}{J}_{j-2}^{\alpha +2,\beta +2}\left(\frac{2x}{L}-1\right).$$

(42)

Eventually, a numerical scheme for solving (8)–(11) is established as follows:

$$\begin{array}{ccc}& & \sum _{i=0}^N\sum _{j=0}^M{b}_{ij}\left[\frac{\Gamma \left(i+1\right){\left(1-t\right)}^{-s}}{\Gamma \left(i-s+1\right)}{J}_i^{-s,0}\left(2t-1\right)-\frac{J_i^{0,-s}\left(1\right){\left(1-t\right)}^{-s}}{\Gamma \left(1-\alpha \right)}\right]J_j^{\alpha ,\beta }\left(\frac{2x}{L}-1\right)\hfill \\ & =& \sum _{i=0}^N\sum _{j=0}^M\left[b_{ij}\frac{\Gamma \left(j+3+\alpha +\beta \right)L^{-2}}{\Gamma \left(j+\alpha +\beta +1\right)}{J}_{j-2}^{2+\alpha ,2+\beta }\left(\frac{2x}{L}-1\right)J_i^{0,-s}\left(2t-1\right)\right.\hfill \\ & & \left.-\frac{Q}{R}{a}_{ij}{J}_j^{\alpha ,\beta }\left(\frac{2x}{L}-1\right)J_i^{-s,0}\left(2t-1\right)\right],\hfill \end{array}$$

(43)

and

$$\begin{array}{ccc}& & \sum _{i=0}^N\sum _{j=0}^M{a}_{ij}\left[\frac{\Gamma \left(i+1\right)t^{-s}}{\Gamma \left(i-s+1\right)}{J}_i^{0,-s}\left(2t-1\right)-\frac{J_i^{-s,0}\left(-1\right)t^{-s}}{\Gamma \left(1-\alpha \right)}\right]J_j^{\alpha ,\beta }\left(\frac{2x}{L}-1\right)\hfill \\ & =& \sum _{i=0}^N\sum _{j=0}^M\left[a_{ij}\frac{\Gamma \left(j+3+\alpha +\beta \right)L^{-2}}{\Gamma \left(j+\alpha +\beta +1\right)}{J}_{j-2}^{2+\alpha ,2+\beta }\left(\frac{2x}{L}-1\right)J_i^{-s,0}\left(2t-1\right)\right.\hfill \\ & & \left.+b_{ij}{J}_j^{\alpha ,\beta }\left(\frac{2x}{L}-1\right)J_i^{0,-s}\left(2t-1\right)\right].\hfill \end{array}$$

(44)

Note that JGL points $\left\{x_{L,j}^{\alpha ,\beta }\right\},0\le j\le M$ and $\left\{t_{1,j}^s\right\},0\le j\le N$ include $t_{1,0}^s=0,t_{1,N}^s=1,x_{L,0}^{\alpha ,\beta }=0,x_{L,M}^{\alpha ,\beta }=L$, which are necessarily used to satisfy the initial and boundary conditions.

5. Numerical Experiment

In this section, we present two numerical examples to verify the efficiency of our algorithm. The same as the former part s is the order of fractional derivative, and $\alpha$, $\beta$ are Jacobi polynomials parameters. We first discuss an example with an exact solution to illustrate the implementation and efficiency of the proposed numerical method; then, a fractional optimal-control problem with noncontinuous initial input signal is considered. The MATLAB R2015b software is used for all computations in this section.

5.1. Numerical Solution for System with Smooth Initial Input Signal

As the first example, we assume that systems (21) and (22) have an analytical solution given by

$$\begin{array}{c}\hfill u(x,t)=\left(t-1\right)cos\left(x\pi \right),\\ \hfill z(x,t)=\left(1-t\right)cos\left(x\pi \right).\end{array}$$

Without loss of generality, we assume $R=Q=L=1$ to simplify the computing procedures. In order to balance both sides of these equations, let

$$f_1(x,t)=\frac{1}{\Gamma (2-\alpha )}{(1-t)}^{1-s}cos\left(x\pi \right)+cos\left(x\pi \right)\left[{\pi }^2(1-t)+t-1\right],$$

(45)

and

$$f_2(x,t)=\frac{1}{\Gamma (2-\alpha )}{t}^{1-s}cos\left(x\pi \right)+cos\left(x\pi \right)\left[{\pi }^2(t-1)+t-1\right],$$

(46)

which are purposively added to the right hand side parts of Equations (21) and (22), respectively. Therefore, the initial input signal is governed by a cosine waveform

$$z_0\left(x\right)=-cos\left(x\pi \right),0\le x\le 1.$$

(47)

One important advantage of the collocation method is the high accuracy. Firstly, we compare the exact and numerical solutions for $s=0.7$, $\alpha =\beta =3$. Figure 1a,b shows the exact and numerical solutions of state and control functions by severally fixing $x=0.2$ and $t=0.5$ and choosing $N=M=10$ in collocation method. Note that the numerical solutions of both $z(x,t)$ and $u(x,t)$ are highly aligned with the exact solutions.

Table 1. Maximum absolute errors (MAE).

M	$\mathit{\alpha }$	$\mathit{\beta }$	MAE of z	MAE of u
5			1.05 $\times {10}^{-3}$	1.13 $\times {10}^{-3}$
10	−0.5	−0.5	3.36$\times {10}^{-8}$	3.43$\times {10}^{-8}$
15			2.89$\times {10}^{-15}$	3.09$\times {10}^{-14}$
5			2.78$\times {10}^{-3}$	2.96$\times {10}^{-3}$
10	0.5	0.5	1.82$\times {10}^{-7}$	1.86$\times {10}^{-3}$
15			4.51$\times {10}^{-14}$	5.00$\times {10}^{-14}$
5			1.98$\times {10}^{-3}$	2.12$\times {10}^{-3}$
10	0	0	9.45$\times {10}^{-8}$	9.65$\times {10}^{-8}$
15			2.62$\times {10}^{-14}$	2.98$\times {10}^{-14}$
5			5.95$\times {10}^{-3}$	7.90$\times {10}^{-3}$
10	2	3	1.80$\times {10}^{-6}$	2.87$\times {10}^{-6}$
15			2.87$\times {10}^{-13}$	4.13$\times {10}^{-13}$

lists the maximum absolute errors (MAE) of the exact and numerical solution at $x=0.2$ and $t=0.5$ with four kinds of Jacobi polynomials parameters, which are

• $\alpha =\beta =-0.5$ (Chebyshev polynomials of the first kind),
• $\alpha =\beta =0.5$ (Chebyshev polynomial of the second kind),
• $\alpha =\beta =0$ (Legendre polynomials),
• $\alpha =2$, $\beta =3$ (non-specified Jacobi polynomials).

The result shows that the proposed algorithm provides high accuracy in cases of different Jacobi parameters. Figure 1c,d show that the relation between the degree M, N ($M=N$) of the polynomial and the logarithm of error (E) is almost linear, which indicates that the error decays exponentially.

5.2. Numerical Solution for System with Discontinuous Initial Input Signal

In this part, a kind of discontinuous initial condition is considered. Similarly, we take $R=Q=L=1$, and we solve systems (21) and (22) with the initial condition expressed by a step-function signal

$$z_0\left(x\right)=\left\{\begin{array}{cc}10,\hfill & 0\le x<0.5,\hfill \\ 0,\hfill & 0.5\le x\le 1.\hfill \end{array}\right.$$

In computation, we select $N=M=11$ and $\alpha =\beta =-0.5$. The numerical solution for $s=0.9$ is displayed in Figure 2a. One can observe that the state function $z(x,t)$ becomes increasingly smooth after a short time because of the influences of diffusion and Neumann boundary conditions. Meanwhile, the state function will converge to a stable equilibrium that is not identical zero because the control input is not zero.

In Figure 2b–d, we sketched the curves of the control function with three different fractional orders ($s=0.4$, $s=0.6$ and $s=0.9999$) as time goes by. In this case, we set $\alpha =\beta =1$. This shows that when s is relatively small ($s=0.4$), the state function $z(x,t)$ decays much faster than the case of relatively large s ($s=0.6$). This coincides with the memory property of the fractional derivative, which is sometimes called the heavy tail feature. When s approaches 1, a comparatively slow and uniform decay phenomenon can be viewed, which matches our previous analysis.

6. Conclusions

In this paper, a collocation method based on the Jacobi polynomial is presented for solving a class of optimal-control problems with a fractional distributed system. The problem is transformed into a system of fractional diffusion equations by using the fractional variational principle. The major feature of the resulting system is that both left-sided and right-sided fractional derivatives are involved, which makes it more difficult to handle. We approximate the numerical solutions in terms of the Jacobi polynomials in both temporal and spatial directions. The efficiency and exponential convergence are verified numerically. The numerical simulations show that the proposed method is validated for fractional optimal-control problem containing continuous or discontinuous initial conditions. In the case of discontinuous initial conditions, the results shows that the state function decays much faster when the order of fractional derivative s is relatively small, which coincides with the heavy tail feature of the fractional derivative.