Open AccessArticle

Collocation Method for Optimal Control of a Fractional Distributed System

Wen Cao

^1,*

and

Yufeng Xu

School of Mathematics and Statistics, Hunan University of Finance and Economics, 139 Fenglin 2nd Road, Changsha 410205, China

School of Mathematics and Statistics, Central South University, 932 Lushan South Road, Changsha 410083, China

Author to whom correspondence should be addressed.

Fractal Fract. 2022, 6(10), 594; https://doi.org/10.3390/fractalfract6100594

Submission received: 7 September 2022 / Revised: 8 October 2022 / Accepted: 10 October 2022 / Published: 14 October 2022

(This article belongs to the Special Issue Advances in Fractional Order Derivatives and Their Applications)

Download

Browse Figures

Versions Notes

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.

Keywords:

collocation method; fractional optimal-control problem; fractional variational principle; exponential convergence

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 (x) \in L^{1} ([a, b])

s > 0

be a finite and real number. Then,

(_{a} I_{x}^{s} f) (x) = \frac{1}{Γ (s)} \int_{a}^{x} {(x - t)}^{s - 1} f (t) d t,

(1)

and

(_{x} I_{b}^{s} f) (x) = \frac{1}{Γ (s)} \int_{x}^{b} {(t - x)}^{s - 1} f (t) d t,

(2)

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

Definition 2

([29]). Let

f (x) \in A C ([a, b])

and

_{a} I_{t}^{n - s} f, {_{t} I}_{b}^{n - s} f \in C^{n} ([a, b])

n - 1 < s < n

, n is a natural number. Then,

(_{a} D_{x}^{s} f) (x) = \frac{1}{Γ (n - s)} \frac{d^{n}}{d x^{n}} \int_{a}^{x} {(x - t)}^{n - s - 1} f (t) d t,

(3)

and

(_{x} D_{b}^{s} f) (x) = \frac{{(- 1)}^{n}}{Γ (n - s)} \frac{d^{n}}{d t^{n}} \int_{x}^{b} {(t - x)}^{n - s - 1} f (t) d t,

(4)

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

Definition 3

([29]). Let

f (x) \in A C ([a, b])

n - 1 < s < n

, n is a natural number. Then,

({_{a}^{c} D}_{x}^{s} f) (x) = \frac{1}{Γ (n - s)} \int_{a}^{x} {(x - t)}^{n - s - 1} \{\frac{d^{n}}{d t^{n}} f (t)\} d t,

(5)

and

({_{x}^{c} D}_{b}^{s} f) (x) = \frac{{(- 1)}^{n}}{Γ (n - s)} \int_{x}^{b} {(t - x)}^{n - s - 1} \{\frac{d^{n}}{d t^{n}} f (t)\} d t,

(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:

(_{x} D_{b}^{s} f) (x) = ({_{x}^{c} D}_{b}^{s} f) (x) + f (b) \frac{{(b - x)}^{- s}}{Γ (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 (b) = 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 \leq t \leq 1

0 \leq x \leq 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 (x, t)

, which minimizes the cost functional

J = \int_{0}^{1} \int_{0}^{L} Q z^{2} (x, t) + R u^{2} (x, t) d x d t,

(8)

subject to the system dynamic constraints

_{0}^{c} D_{t}^{s} z (x, t) = \frac{\partial^{2} z (x, t)}{\partial x^{2}} + u (x, t), 0 < s < 1,

(9)

the initial condition

z (x, 0) = z_{0} (x), 0 \leq x \leq L,

(10)

and the boundary condition

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

(11)

where s is the order of Caputo fractional derivative,

z (x, t)

and

u (x, t)

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

λ

\tilde{J} = \int_{0}^{1} \int_{0}^{L} Q z^{2} (x, t) + R u^{2} (x, t) + λ [_{0}^{c} D_{t}^{s} z (x, t) - \frac{\partial^{2} z (x, t)}{\partial x^{2}} - u (x, t)] d x d t .

(12)

Assuming that

z^{*} (x, t)

and

u^{*} (x, t)

are the desired optimum solutions of (8)–(11), define

\begin{matrix} z (x, t) = z^{*} (x, t) + ε_{1} η_{1} (x, t), ε_{1} \in R, \end{matrix}

(13)

\begin{matrix} u (x, t) = u^{*} (x, t) + ε_{2} η_{2} (x, t), ε_{2} \in R, \end{matrix}

(14)

where

ε_{1}

ε_{2}

are perturbation parameters, and

η_{1} (x, t)

and

η_{2} (x, t)

are arbitrary test functions satisfying

η_{1} (x, 0) = η_{1 x} (0, t) = η_{1 x} (L, t) = 0

Substitute Equations (13) and (14) into functional (12), then

\tilde{J}

becomes a scale function of

ε_{1}

and

ε_{2}

. Obviously,

\tilde{J}

has an extremum at

ε_{1} = ε_{2} = 0

. Differentiating

\tilde{J}

with respect to

ε_{1}, ε_{2}

and

λ

, respectively, we have

\begin{matrix} \frac{\partial \tilde{J}}{\partial ε_{1}} & = & \int_{0}^{1} \int_{0}^{L} [2 Q z (x, t) η_{1} (x, t) + λ_{0}^{c} D_{t}^{s} η_{1} (x, t) - λ \frac{\partial^{2} η_{1} (x, t)}{\partial x^{2}}] d x d t, \end{matrix}

(15)

\begin{matrix} \frac{\partial \tilde{J}}{\partial ε_{2}} & = & \int_{0}^{1} \int_{0}^{L} [2 R u (x, t) η_{2} (x, t) - λ η_{2} (x, t)] d x d t, \end{matrix}

(16)

\begin{matrix} \frac{\partial \tilde{J}}{\partial λ} & = & \int_{0}^{1} \int_{0}^{L} [_{0}^{c} D_{t}^{s} z (x, t) - \frac{\partial^{2} z (x, t)}{\partial x^{2}} - u (x, t)] d x d t . \end{matrix}

(17)

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

\begin{matrix} \int_{0}^{1} \int_{0}^{L} λ [{_{0}^{c} D}_{t}^{s} η_{1} (x, t) - \frac{\partial^{2} η_{1} (x, t)}{\partial x^{2}}] d x d t = \int_{0}^{1} \int_{0}^{L} η_{1} (x, t) [{_{t} D}_{1}^{s} λ - \frac{\partial^{2} λ (x, t)}{\partial x^{2}}] d x d t \\ + \int_{0}^{1} [λ_{x} (L, t) η_{1} (L, t) - λ_{x} (0, t) η_{1} (0, t)] d t + \int_{0}^{L} η_{1} {(x, 0)}_{t} I_{1}^{s} λ (x, 0) d x . \end{matrix}

To minimize

\tilde{J}

(the same as to minimize J), it is required that results of (15)–(17) are zeros when

ε_{1} = ε_{2} = 0

. Then, by the arbitrariness of

λ

, one could assume

λ_{x} (0, t) = λ_{x} (L, t) = 0

λ (x, 1) = 0

to have

\begin{matrix} _{t}^{c} D_{1}^{s} λ (x, t) = \frac{\partial^{2} λ (x, t)}{\partial x^{2}} - 2 Q z^{*} (x, t), \end{matrix}

(18)

\begin{matrix} λ (x, t) = 2 R u^{*} (x, t), \end{matrix}

(19)

\begin{matrix} _{0}^{c} D_{t}^{s} z^{*} (x, t) = \frac{\partial^{2} z^{*} (x, t)}{\partial x^{2}} + u^{*} (x, t), \end{matrix}

(20)

since

η_{1} (x, t)

and

η_{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{matrix} _{t}^{c} D_{1}^{s} u (x, t) = \frac{\partial^{2} u (x, t)}{\partial x^{2}} - \frac{Q}{R} z (x, t), \end{matrix}

(21)

\begin{matrix} _{0}^{c} D_{t}^{s} z (x, t) = \frac{\partial^{2} z (x, t)}{\partial x^{2}} + u (x, t), \end{matrix}

(22)

with the initial and boundary conditions

\begin{matrix} z (x, 0) = z_{0} (x), & u (x, 1) = 0, \\ \frac{\partial z (0, t)}{\partial x} = \frac{\partial z (L, t)}{\partial x} = 0, & \frac{\partial u (0, t)}{\partial x} = \frac{\partial u (L, t)}{\partial x} = 0 . \end{matrix}

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{matrix} _{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), \end{matrix}

(23)

\begin{matrix} _{0}^{c} D_{t}^{s} \tilde{z} (x, t) = \frac{\partial^{2} \tilde{z} (x, t)}{\partial x^{2}} + \tilde{u} (x, t) . \end{matrix}

(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} d x + Q \int_{0}^{L} {\tilde{z}}_{0}^{c} D_{t}^{s} \tilde{z} d x = Q \int_{0}^{L} \tilde{z} {\tilde{z}}_{x x} d x + R \int_{0}^{L} \tilde{u} {\tilde{u}}_{x x} d x .

(25)

In view of

\begin{matrix} Q \int_{0}^{L} \tilde{z} {\tilde{z}}_{x x} d x + R \int_{0}^{L} \tilde{u} {\tilde{u}}_{x x} d x \\ = & Q ({\tilde{z} {\tilde{z}}_{x}|}_{0}^{L} - \int_{0}^{L} {\tilde{z}}_{x}^{2} d x) + R ({\tilde{u} {\tilde{u}}_{x}|}_{0}^{L} - \int_{0}^{L} {\tilde{u}}_{x}^{2} d x) \\ \leq & 0, \end{matrix}

we conclude that either one of the following inequalities hold:

R \int_{0}^{L} \tilde{u} \cdot (_{t}^{c} D_{1}^{s} \tilde{u}) d x \leq 0 or Q \int_{0}^{L} \tilde{z} \cdot (_{0}^{c} D_{t}^{s} \tilde{z}) d x \leq 0 .

(26)

If the former case is true, we assume that

\tilde{u}

is approximated by smooth polynomials

φ (x)

and

ϕ (t)

, which is

\tilde{u} = φ (x) ϕ (t)

. Then,

\begin{matrix} \int_{0}^{L} \tilde{u} \cdot (_{t}^{c} D_{1}^{s} \tilde{u}) d x & = & \int_{0}^{L} φ (x) ϕ (t) \cdot φ (x) (_{t}^{c} D_{1}^{s} ϕ (t)) d x \\ = & ϕ (t) \cdot (_{t}^{c} D_{1}^{s} ϕ (t)) \int_{0}^{L} φ^{2} (x) d x . \end{matrix}

(27)

Furthermore, notice that

\begin{matrix} \int_{0}^{1} ϕ (t) \cdot (_{t}^{c} D_{1}^{s} ϕ (t)) \int_{0}^{L} φ^{2} (x) d x d t & = & \int_{0}^{1} ϕ (t) \cdot (_{t}^{c} D_{1}^{s} ϕ (t)) d t \cdot \int_{0}^{L} φ^{2} (x) d x \\ = & \int_{0}^{1} (_{0} D_{t}^{s} ϕ (t)) \cdot ϕ (t) d t \cdot \int_{0}^{L} φ^{2} (x) d x \\ = & \int_{0}^{1} {(_{0} D_{t}^{s / 2} ϕ (t))}^{2} d t \cdot \int_{0}^{L} φ^{2} (x) d x \\ \geq & 0, \end{matrix}

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

\int_{0}^{L} \tilde{u} \cdot (_{t}^{c} D_{1}^{s} \tilde{u}) d x \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) d x = \int_{0}^{L} \frac{\partial^{2} z (x, t)}{\partial x^{2}} d x + \int_{0}^{L} u (x, t) d x .

(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

ζ (t) = \int_{0}^{L} u (x, t) d x < 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) d x = \frac{t^{1 - s}}{Γ (2 - s)} \int_{0}^{L} z_{t} (x, t_{l}^{*}) d x .

(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) d x = \frac{t^{1 - s}}{Γ (2 - s)} \int_{0}^{L} z_{t} (x, t_{l}^{*}) d x .

(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) d x = - \frac{{(1 - t)}^{1 - s}}{Γ (2 - s)} \int_{0}^{L} u_{t} (x, t_{r}^{*}) d x .

(32)

Combining (31) and (32), we obtain

Φ (t, s) : = \int_{0}^{L} z (x, t) d x = C \cdot \frac{{(1 - t)}^{1 - s}}{Γ (1 - s) Γ (2 - s)}, 0 < s, t < 1,

(33)

where

C = - t_{r}^{* - s} \int_{0}^{L} z_{t} (x, t_{l}^{*}) d x \geq 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

α

and

β

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}^{α, β}

with indices

α, β > - 1

of degree n are defined on interval

[- 1, 1]

and can be expressed as

J_{n}^{α, β} (ξ) = \frac{{(- 1)}^{n}}{2^{n} \cdot n!} {(1 - ξ)}^{- α} {(1 + ξ)}^{- β} \frac{d^{n}}{d ξ^{n}} [{(1 - ξ)}^{α + n} {(1 + ξ)}^{β + n}], - 1 \leq ξ \leq 1,

(34)

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

(1 - ξ^{2}) \frac{d^{2} Y (ξ)}{d ξ^{2}} + [β - α - (α + β + 2) ξ] \frac{d Y (ξ)}{d ξ} + n (n + α + β + 1) Y (ξ) = 0 .

In order to use these polynomials on the interval

[0, L]

, we take the change of variable

x = \frac{L}{2} (ξ + 1)

. Therefore,

0 \leq x \leq L

is due to the range of

ξ

. The derivative of

J_{n}^{α, β}

on the interval

[0, L]

can be represented by

\frac{d^{m}}{d t^{m}} J_{n}^{α, β} (\frac{2 x}{L} - 1) = \frac{Γ (n + m + α + β + 1)}{Γ (n + α + β + 1) L^{m}} J_{n - m}^{α + m, β + m} (\frac{2 x}{L} - 1),

(35)

where

m = 1, 2, 3, \dots

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}^{α, β}

can be presented in the form of the Jacobi polynomial. Let

0 < α < 1

β = 0

; we have

_{0} D_{x}^{α} J_{n}^{- α, 0} (\frac{2 x}{L} - 1) = \frac{Γ (n + 1) x^{- α}}{Γ (n - α + 1)} J_{n}^{0, - α} (\frac{2 x}{L} - 1),

(36)

and when

α = 0

0 < β < 1

, it reads that

_{x} D_{1}^{β} J_{n}^{0, - β} (\frac{2 x}{L} - 1) = \frac{Γ (n + 1) {(L - x)}^{- β}}{Γ (n - β + 1)} J_{n}^{- β, 0} (\frac{2 x}{L} - 1) .

(37)

Next, we present the basic idea of the collocation method with Jacobi polynomials. Denoting sets

\{x_{j}^{α, β}\}, 0 \leq j \leq M

and

\{t_{j}^{s}\}, 0 \leq j \leq N

as the nodes of Jacobi–Gauss–Lobatto (JGL) quadratures on the interval

[- 1, 1]

, we set the JGL points

\{x_{L, j}^{α, β}\}, 0 \leq j \leq M

and

\{t_{1, j}^{s}\}, 0 \leq j \leq N

[0, L] \times [0, 1]

\begin{matrix} x_{L, j}^{α, β} & = & \frac{L}{2} (x_{j}^{α, β} + 1), 0 \leq j \leq M, \\ t_{1, j}^{s} & = & \frac{1}{2} (t_{j}^{s} + 1), 0 \leq j \leq N . \end{matrix}

Let

{ψ_{i} (t), ϕ_{j} (x)}

{\hat{ψ_{i}} (t), ϕ_{j} (x)}

be the pairs of Jacobi basis polynomials; the approximations for solutions

z_{N}

and

u_{N}

can be expanded as

\begin{matrix} z_{N} (x, t) = \sum_{i = 0}^{N} \sum_{j = 0}^{M} a_{i j} ϕ_{j} (x) ψ_{i} (t), \end{matrix}

(38)

\begin{matrix} u_{N} (x, t) = \sum_{i = 0}^{N} \sum_{j = 0}^{M} b_{i j} ϕ_{j} (x) {\hat{ψ}}_{i} (t), \end{matrix}

(39)

by using the method of separation of variables. Here,

{ψ_{i} (t), ϕ_{j} (x)}

and

{\hat{ψ_{i}} (t), ϕ_{j} (x)}

are different Jacobi basis polynomials, as follows:

\begin{matrix} ϕ_{j} (x) & = & J_{j}^{α, β} (\frac{2 x}{L} - 1), \\ ψ_{i} (t) & = & J_{i}^{- s, 0} (2 t - 1), \\ {\hat{ψ}}_{i} (t) & = & J_{i}^{0, - s} (2 t - 1) . \end{matrix}

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

\begin{matrix} \sum_{i = 0}^{N} \sum_{j = 0}^{M} b_{i j} ϕ_{j} (x) \cdot {_{t} D}_{1}^{s} {\hat{ψ}}_{i} (t) = \sum_{i = 0}^{N} \sum_{j = 0}^{M} b_{i j} {\hat{ψ}}_{i} (t) \frac{\partial^{2} ϕ_{j} (x)}{\partial x^{2}} - \frac{Q}{R} \sum_{i = 0}^{N} \sum_{j = 0}^{M} a_{i j} ϕ_{j} (x) ψ_{i} (t), \\ \sum_{i = 0}^{N} \sum_{j = 0}^{M} a_{i j} ϕ_{j} (x) \cdot {_{0}^{c} D}_{t}^{s} ψ_{i} (t) = \sum_{i = 0}^{N} \sum_{j = 0}^{M} a_{i j} ψ_{i} (t) \frac{\partial^{2} ϕ_{j} (x)}{\partial x^{2}} + \sum_{i = 0}^{N} \sum_{j = 0}^{M} b_{i j} ϕ_{j} (x) {\hat{ψ}}_{i} (t) . \end{matrix}

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

_{0}^{c} D_{t}^{s} J_{i}^{- s, 0} (2 t - 1) = \frac{Γ (i + 1) t^{- s}}{Γ (i - s + 1)} J_{i}^{0, - s} (2 t - 1) - \frac{J_{i}^{- s, 0} (- 1) t^{- s}}{Γ (1 - s)},

(40)

and

_{t}^{c} D_{1}^{s} J_{i}^{0, - s} (2 t - 1) = \frac{Γ (i + 1) {(1 - t)}^{- s}}{Γ (i - s + 1)} J_{i}^{- s, 0} (2 t - 1) - \frac{J_{i}^{0, - s} (1) {(1 - t)}^{- s}}{Γ (1 - s)} .

(41)

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

\frac{\partial^{2} ϕ_{j} (x)}{\partial x^{2}} = \frac{Γ (j + 3 + α + β)}{Γ (j + 1 + α + β) L^{2}} J_{j - 2}^{α + 2, β + 2} (\frac{2 x}{L} - 1) .

(42)

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

\begin{matrix} \sum_{i = 0}^{N} \sum_{j = 0}^{M} b_{i j} [\frac{Γ (i + 1) {(1 - t)}^{- s}}{Γ (i - s + 1)} J_{i}^{- s, 0} (2 t - 1) - \frac{J_{i}^{0, - s} (1) {(1 - t)}^{- s}}{Γ (1 - α)}] J_{j}^{α, β} (\frac{2 x}{L} - 1) \\ = & \sum_{i = 0}^{N} \sum_{j = 0}^{M} [b_{i j} \frac{Γ (j + 3 + α + β) L^{- 2}}{Γ (j + α + β + 1)} J_{j - 2}^{2 + α, 2 + β} (\frac{2 x}{L} - 1) J_{i}^{0, - s} (2 t - 1) \\ - \frac{Q}{R} a_{i j} J_{j}^{α, β} (\frac{2 x}{L} - 1) J_{i}^{- s, 0} (2 t - 1)], \end{matrix}

(43)

and

\begin{matrix} \sum_{i = 0}^{N} \sum_{j = 0}^{M} a_{i j} [\frac{Γ (i + 1) t^{- s}}{Γ (i - s + 1)} J_{i}^{0, - s} (2 t - 1) - \frac{J_{i}^{- s, 0} (- 1) t^{- s}}{Γ (1 - α)}] J_{j}^{α, β} (\frac{2 x}{L} - 1) \\ = & \sum_{i = 0}^{N} \sum_{j = 0}^{M} [a_{i j} \frac{Γ (j + 3 + α + β) L^{- 2}}{Γ (j + α + β + 1)} J_{j - 2}^{2 + α, 2 + β} (\frac{2 x}{L} - 1) J_{i}^{- s, 0} (2 t - 1) \\ + b_{i j} J_{j}^{α, β} (\frac{2 x}{L} - 1) J_{i}^{0, - s} (2 t - 1)] . \end{matrix}

(44)

Note that JGL points

\{x_{L, j}^{α, β}\}, 0 \leq j \leq M

and

\{t_{1, j}^{s}\}, 0 \leq j \leq N

include

t_{1, 0}^{s} = 0, t_{1, N}^{s} = 1, x_{L, 0}^{α, β} = 0, x_{L, M}^{α, β} = 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

α

β

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{matrix} u (x, t) = (t - 1) cos (x π), \\ z (x, t) = (1 - t) cos (x π) . \end{matrix}

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}{Γ (2 - α)} {(1 - t)}^{1 - s} cos (x π) + cos (x π) [π^{2} (1 - t) + t - 1],

(45)

and

f_{2} (x, t) = \frac{1}{Γ (2 - α)} t^{1 - s} cos (x π) + cos (x π) [π^{2} (t - 1) + t - 1],

(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} (x) = - cos (x π), 0 \leq x \leq 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

α = β = 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 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

$α = β = - 0.5$ (Chebyshev polynomials of the first kind),
$α = β = 0.5$ (Chebyshev polynomial of the second kind),
$α = β = 0$ (Legendre polynomials),
$α = 2$ , $β = 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} (x) = \{\begin{matrix} 10, & 0 \leq x < 0.5, \\ 0, & 0.5 \leq x \leq 1 . \end{matrix}

In computation, we select

N = M = 11

and

α = β = - 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

α = β = 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.

Author Contributions

Conceptualization, W.C. and Y.X.; methodology, W.C.; validation, W.C. and Y.X.; formal analysis, W.C. and Y.X.; writing—original draft preparation, W.C; writing—review and editing, W.C.; and funding acquisition, W.C. All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded by the Hunan Provincial Natural Science Foundation of China (grant number 2019JJ50019) and the Scientific Research Foundation of Hunan Provincial Education Department (grant number 18C0970).

Data Availability Statement

Not applicable.

Acknowledgments

The authors would like to thank the reviewers for their constructive comments to improve the manuscript.

Conflicts of Interest

The authors declare no conflict of interest.

References

Sage, A.P.; White, C.C., III. Optimal Systems Control; Prentice-Hall: Englewood Cliffs, NJ, USA, 1977. [Google Scholar]
Betts, J.T.; Erb, S.O. Optimal low thrust trajectories to the moon. SIAM J. Appl. Dyn. Syst. 2003, 2, 144–170. [Google Scholar] [CrossRef] [Green Version]
Highfill, J.; McAsey, M. An optimal-control problem in economics. Int. J. Mathemetics Math. Sci. 1991, 14, 537–544. [Google Scholar] [CrossRef]
Abu-Rqayiqa, A.; Alayed, H.; Zannon, M. Optimal control of a basic model of oncolytic virotherapy. J. Mathemetics Comput. Sci. 2022, 24, 119–126. [Google Scholar] [CrossRef]
Bittner, L.; Bulirsch, R.; Heier, K.; Schmidt, W. (Eds.) Variational Calculus, Optimal Control and Applications. In Proceedings of the International Conference in Honour of L. Bittner and R. Klötzler, Trassenheide, Germany, 23–27 September 1996.
Hull, D.G. Variational calculus and approximate solution of optimal-control problems. J. Optim. Theory Appl. 2001, 108, 483–497. [Google Scholar] [CrossRef]
Biswas, R.K.; Sen, S. Free final time fractional optimal-control problems. J. Frankl. Inst.-Eng. Appl. Math. 2014, 351, 941–951. [Google Scholar] [CrossRef]
Kamocki, R. On the existence of optimal solutions to fractional optimal-control problems. Appl. Math. Comput. 2014, 235, 94–104. [Google Scholar] [CrossRef]
Almeida, R.; Torres, D.F.M. A discrete method to solve fractional optimal-control problems. Nonlinear Dyn. 2015, 80, 1811–1816. [Google Scholar] [CrossRef] [Green Version]
Xu, Y.; Agrawal, O.P. New fractional operators and application to fractional variational problem. Comput. Math. Appl. 2016; in Press. [Google Scholar] [CrossRef]
Paulen, R.; Kikar, M. Optimal Operation of Batch Membrane Processes; Springer: Berlin/Heidelberg, Germany, 2016. [Google Scholar]
Agrawal, O.P. Formulation of Euler-Lagrange equations for fractional variational problems. J. Math. Anal. Appl. 2002, 272, 368–379. [Google Scholar] [CrossRef] [Green Version]
Agrawal, O.P. A general formulation and solution scheme for fractional optimal-control problem. Nonlinear Dyn. 2004, 38, 323–337. [Google Scholar] [CrossRef]
Agrawal, O.P. Fractional optimal control of a distributed system using eigenfunctions. ASME J. Comput. Nonlinear Dyn. 2008, 3, 021204-1–021204-6. [Google Scholar] [CrossRef]
Baleanu, D.; Defterli, O.; Agrawal, O.P. A central difference numerical scheme for fractional optimal-control problems. J. Vib. Control. 2009, 15, 583–597. [Google Scholar] [CrossRef] [Green Version]
Tricaud, C.; Chen, Y.Q. An approximate method for numerically solving fractional order optimal control problems of general form. Comput. Math. Appl. 2010, 59, 1644–1655. [Google Scholar] [CrossRef] [Green Version]
Lotfi, A.; Dehghan, M.; Yousefi, S.A. A numerical technique for solving fractional optimal-control problems. Comput. Math. Appl. 2011, 62, 1055–1067. [Google Scholar] [CrossRef] [Green Version]
Lotfi, A.; Yousefi, S.A.; Dehghan, M. Numerical solution of a class of fractional optimal-control problems via the Legendre orthonormal basis combined with the operational matrix and the Gauss quadrature rule. J. Comput. Appl. Math. 2013, 250, 143–160. [Google Scholar] [CrossRef]
Özdemir, N.; Karadeniz, D.; İskender, B.B. Fractional optimal-control problem of a distributed system in cylindrical coordinates. Phys. Lett. A 2009, 373, 221–226. [Google Scholar] [CrossRef]
Safaie, E.; Farahi, M.H.; Ardehaie, M.F. An approximate method for numerically solving multi-dimensional delay fractional optimal-control problems by Bernstein polynomials. Comput. Appl. Math. 2015, 34, 831–846. [Google Scholar] [CrossRef]
Özdemir, N.; Agrawal, O.P.; İskender, B.B.; Karadeniz, D. Fractional optimal control of a 2-dimensional distributed system using eigenfunctions. Nonlinear Dyn. 2009, 55, 251–260. [Google Scholar] [CrossRef]
Bhrawy, A.H.; Doha, E.H.; Baleanu, D.; Ezz-Eldien, S.S.; Abdelkawy, M.A. An accurate numerical technique for solving fractional optimal-control problems. Proc. Rom. Acad. Ser. A 2015, 16, 47–54. [Google Scholar]
Darehmiraki, M.; Farahi, M.H.; Effati, S. Solution for the fractional distrubted optimal-control problem by hybrid meshless method. J. Vib. Control 2018, 24, 2149–2164. [Google Scholar] [CrossRef]
Canuto, C.; Hussaini, M.Y.; Quarteroni, A.; Zang, T.A. Spectral Methods: Fundamentals in Single Domains. Springer: Berlin/Heidelberg, Germany, 2006. [Google Scholar]
Yang, Y.; Chen, Y.P.; Huang, Y.Q. Convergence analysis of the Jacobi spectral-collocation method for fractional integro-differential equations. Acta Mathemetica Sci. 2014, 34B, 673–690. [Google Scholar] [CrossRef]
Tohidi, E.; Nik, H.S. A Bessel collocation method for solving fractional optimal-control problems. Appl. Math. Model. 2015, 39, 455–465. [Google Scholar] [CrossRef]
Nikan, O.; Avazzadeh, Z. Numerical simulation of fractional evolution model arising in viscoelastic mechanics. Appl. Numer. Math. 2021, 169, 303–320. [Google Scholar] [CrossRef]
Shojaeizadeh, T.; Mahmoudi, M.; Darehmiraki, M. optimal-control problem of advection-diffusion-reaction equation of kind fractal-fractional applying shifted Jacobi polynomials. Chaos Solitons Fractals 2021, 143, 110568. [Google Scholar] [CrossRef]
Kilbas, A.A.; Srivastava, H.M.; Trujillo, J.J. Theory and Applications of Fractional Differential Equations; Elsevier: Amsterdam, The Netherlands, 2006. [Google Scholar]
Podlubny, I. Fractional Differential Equations; Academic Press: Cambridge, MA, USA, 1999. [Google Scholar]
Hilfer, R. Fractional Calculus in Physics; World Scientific Publishing: Singapore, 2000. [Google Scholar]
Bhrawy, A.H. A Jacobi-Gauss-Lobatto collocation method for solving generalized Fitzhugh-Nagumo equation with time-dependent coefficients. Appl. Math. Comput. 2013, 222, 255–264. [Google Scholar] [CrossRef]
Duan, B.P.; Zheng, Z.S.; Cao, W. Spectral approximation methods and error estimates for Caputo fractional derivative with applications to initial-value problems. J. Comput. Phys. 2016, 3, 108–128. [Google Scholar] [CrossRef]

$Fractalfract 06 00594 g001 550$

Figure 1. Numerical and analytical solution of system (21) and (22) (see (a,b)); the relation between N and error (E) with

s = 0.7

and

α = β = 3

(see (c,d)).

Figure 1. Numerical and analytical solution of system (21) and (22) (see (a,b)); the relation between N and error (E) with

s = 0.7

and

α = β = 3

(see (c,d)).

$Fractalfract 06 00594 g001$

$Fractalfract 06 00594 g002a 550$ $Fractalfract 06 00594 g002b 550$

Figure 2. The state function at different moments with

α = β = - 0.5

(see (a)) and

α = β = 1

(see (b–d)).

Figure 2. The state function at different moments with

α = β = - 0.5

(see (a)) and

α = β = 1

(see (b–d)).

$Fractalfract 06 00594 g002a$ $Fractalfract 06 00594 g002b$

Table 1. Maximum absolute errors (MAE).

M	$α$	$β$	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}$

Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

© 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https://creativecommons.org/licenses/by/4.0/).

Share and Cite

MDPI and ACS Style

Cao, W.; Xu, Y. Collocation Method for Optimal Control of a Fractional Distributed System. Fractal Fract. 2022, 6, 594. https://doi.org/10.3390/fractalfract6100594

AMA Style

Cao W, Xu Y. Collocation Method for Optimal Control of a Fractional Distributed System. Fractal and Fractional. 2022; 6(10):594. https://doi.org/10.3390/fractalfract6100594

Chicago/Turabian Style

Cao, Wen, and Yufeng Xu. 2022. "Collocation Method for Optimal Control of a Fractional Distributed System" Fractal and Fractional 6, no. 10: 594. https://doi.org/10.3390/fractalfract6100594

APA Style

Cao, W., & Xu, Y. (2022). Collocation Method for Optimal Control of a Fractional Distributed System. Fractal and Fractional, 6(10), 594. https://doi.org/10.3390/fractalfract6100594

Article Menu

Collocation Method for Optimal Control of a Fractional Distributed System

Abstract

1. Introduction

2. Mathematical Preparation

3. Analysis of Fractional Optimal-Control Problem

3.1. Formulation of Model

3.2. Fractional Variational Principle and Properties of Model

4. Numerical Algorithm

5. Numerical Experiment

5.1. Numerical Solution for System with Smooth Initial Input Signal

5.2. Numerical Solution for System with Discontinuous Initial Input Signal

6. Conclusions

Author Contributions

Funding

Data Availability Statement

Acknowledgments

Conflicts of Interest

References

Share and Cite

Article Metrics

Article Access Statistics

Further Information

Guidelines

MDPI Initiatives

Follow MDPI