Generalized Lagrangian JacobiGaussRadau collocation method for solving a nonlinear 2D optimal control problem with the classical diffusion equation
Abstract
In this paper, a nonlinear 2D Optimal Control Problem (2DOCP) is considered. The quadratic performance index of a nonlinear cost function is endowed with the state and control functions. In this problem, the dynamic constraint of the system is given by a classical diffusion equation. This article is concerned with a generalization of Lagrangian functions. Besides, a Generalized Lagrangian JacobiGaussRadau (GLJGR)collocation method is introduced and applied to solve the aforementioned 2DOCP. Based on initial and boundary conditions, the time and space variables and are considered JacobiGaussRadau points clustered on first or end of interval respectively. Then, to solve the 2DOCP, Lagrange Multipliers are used and the optimal control problem is reduced to a parameter optimization problem. Numerical results demonstrate its accuracy, efficiency, and versatility of the presented method.
Keywords:
Lagrange Multipliers 2D optimal control problemGeneralized Lagrangian functions Generalized Lagrangian Jacobi GaussRadau (GLJGR) collocation method.Msc:
49J20 93C20 34G20∎
1 Introduction
In order to present 2DOCP solved in this manuscript, firstly, we give an introduction to the 2DOCP and provides an explanation of the functions and parameters defined in this problem. A brief review and history of these equations and spectral and Pseudospectral (PS) methods are in the following subsections.
1.1 The governing equations
Optimum control problems rise in the minimization of a functional over a set of admissible control functions subject to dynamic constraints on the state and control functions agrawal2007 (); agrawal2006 (). As the equations of dynamics in the system are reformed by a partial differential equation– with time and space variables– this 2DOCP is known as an optimal control of a distributed system agrawal2007 (). The formulation of this optimal control problem is meme ():
(1) 
subject to
(2) 
with initial and boundary conditions
(3) 
In fact, this is a nonlinear 2D quadratic optimal control problem with the dynamic system of classical diffusion equation.
and are the state and control smooth functions, respectively. and are two arbitrary functions. The upper bound of variable is considered . The parameter is specified in numerical examples as or .
The purpose of solving this problem is to approximate the control and state functions that minimize the .
1.2 The literature of Optimal control problems
Twodimensional (2D) systems and their beneficial applications in many different industrial fields draw the attention of scientists presently. These applications rise in heat transfer, image processing, seismological and geophysical data processing, distributed systems, restoration of noisy images, earthquake signal processing, water stream heating gas absorption, smart materials, and transmission lines meme2015 (); lewis (); mars (); roes (). The miscellaneous chemical, biological, and physical problems are modeled by diffusion processes involving control function mentioned in Eqs. (1)–(3). By the aid of Roesser’s model roes (), Attasi’s model attasi1973 (); attasi1975 () FornasiniMarchesini’s models forn1976 (); forn1978 (), the statespace models of 2D systems are organized meme2015 (). These models are used extensively to analyze controllability, stability, observability of 2D systems.
Remarkable studies have been done in the area of optimal controls, and excellent article are written hereby bryson (); sage (); agrawal1989 (); lotfi2 (); samimi (). Among these studies, numerical techniques have been used to solve optimal control problems agrawal1989 (); gre (). Moreover, Agrawal agrawal2007 () presented a general formulation and a numerical scheme for Fractional Optimal Control for a class of distributed systems. He used eigenfunctions to develop his method. In other works, Manabe manabe (), Bode et al. bode () and Rabiee et al. rabiee () studied fractional order optimal control problems. Additionally, Mamehrashia et al. meme (); meme2015 () and Lotfi et al. lotfi () employed Ritz method to solve optimal control problems. With the Variational method, Yousefi et al. yusefi () found the solution of the optimal control of linear systems, approximately. Li et al. Li () considered a continuous time 2D system and converted it to the discretetime 2D model. In other works, Wei et al. wei () and Zhang et al. zhang () investigated an optimal control problem in continuoustime 2D Roesser’s model. They employed iterative adaptive critic design algorithm and the adaptive dynamic programming method to approximate the solution. Sabeh et al. sabeh () introduced a pseudospectral method for the solution of distributed optimal control problem with the viscous Burgers equation
1.3 The literature of Spectral and PS methods
The main feature of spectral methods is to use different orthogonal polynomials/functions as trial functions. These polynomials/functions are global and infinitely differentiable. These methods are applied to 4 types of problems: periodic problems, nonperiodic problems, whole line problems and half line problems. Trigonometric polynomials for periodic problems; classical Jacobi, ultraspherical, Chebyshev and Legendre polynomials for nonperiodic problems; Hermite polynomials for problems on the whole line; and Laguerre polynomials for problems on the half line doha2012a (). With the truncated series of smooth global functions, spectral methods are giving the solution of a particular problem parand1 (); parand2 (). These methods, with a relatively limited number of degrees of freedom, provide such an accurate approximation for a smooth solution. Spectral coefficients tend to zero faster than any algebraic power of their index bhrawy3 ().
Spectral methods can fall into 3 categories: Tau, Collocation and Galerkin methods bhradoha ().

The Tau spectral method is used to approximate numerical solutions of various differential equations. This method considers the solution as an expansion of orthogonal polynomials/functions. Such coefficients, in this expansion, are set to approximate the solution correctly bhrawy7 ().

Collocation method helps obtain a highly accurate solutions to nonlinear/linear differential equations bhrawy4 (); tal1 (); tal2 (); bhrawy5 (). Two common steps in collocation methods are: First, suitable nodes (Gauss/GaussRadau/GaussLobatto) are selected to restate a finite or discrete form of the differential equations.
Second, a system of algebraic equations from the discretization of the original equation is achieved bhrawy6 (); parand3 (); parand4 (). 
In Galerkin Spectral method, trail and test functions are chosen the same boyd (); This method can result in a highly accurate solution.
It is said that spectral Galerkin methods are similar to Tau methods where in approximation by Tau method, differential equation is enforced bhrawy3 ().
Furthermore, some other numerical methods like Finite difference method (FDM) and Finite element method (FEM) need network construction of data and they perform locally. Although spectral methods are continuous and globally performing, they do not require network construction of data.
As well as spectral methods, PS methods have also attracted researchers recently bhrawy8 (); bhrawy9 (); bhrawy10 (). As mentioned previously, Sabeh et al. sabeh () investigated a PS method to solve optimal control problem. PS methods are also utilized in the solution of other optimal control problems as well moh (); fahroo (); garg (); shamsi (). These methods become popular because of their computational feasibility and efficiency. In fact, in standard PS methods, interpolation operators are used to reducing the cost of computation of the inner product we encounter in some of the spectral methods. For this purpose, a set of distinct interpolation points is considered by which the corresponding Lagrange interpolants are achieved. Besides, when applying collocation points, , the residual function is set to vanish on the same set of points. Notwithstanding, the collocation points do not need to be chosen the same as the interpolation nodes; Indeed, just for having the Kronecker delta property, they are considered to be the same: as a consequence, this property helps reduce computational cost noticeably as well delkhosh (). There are such authors that utilized PS methods for the solution of optimal control as well. William william () introduced a Jacobi PS method for solving an optimal control problem. He reported that significant differences in computation time can be seen for different parameters of the Jacobi polynomial. Garge et al. divya () presented a unified framework for the numerical solution of optimal control problems using collocation at LegendreGauss (LG), LegendreGaussRadau (LGR), and LegendreGaussLobatto (LGL) points and discussed the advantages of each for solving optimal control problems. Chebyshev PS method was utilized by Fahroo et al. fahroo () to provide an optimal solution for optimal control problem.
1.4 The main aim of this paper
To the best of our knowledge, the use of PS methods for solving optimal control problems has been limited in the literature to either Chebyshev or Legendre methods. Noteworthy, the PS method based on Jacobi can encompass a wide range of other PS methods since the Legendre and Chebyshev nodes can be obtained as particular cases of the general Jacobi. This happens when by changing the parameters in the Jacobi polynomial a proper selection of the Jacobi parameters succeed in more accurate realtime solutions to nonlinear optimal control problems. Meanwhile, an arbitrary and not precise selection of nodes may result in a poor interpolation characteristics such as the Runge phenomenon, therefore, the nodes in PS methods are selected as the GaussRadau points william ().
In this paper, we present a general formulation and a suitable numerical method called the GLJGR collocation method to solve 2DOCP for a class of distributed systems. The developed method is exponentially accurate and obtained by generalization of the classical Lagrangian polynomials. Additionally, the equation of the dynamics of optimal control problem is reformed as a partial differential equation.
This paper is arranged as follows: In Section 2, we present some preliminaries and drive some tools for introducing GL function, GLJGR collocation method, and their relevant derivative matrices. In Section 3, we apply the GLJGR collocation method to the solution of the 2DOCP. Section 4 shows numerical examples to demonstrate the effectiveness of the proposed method. Also, a conclusion is given in Section 5.
2 Preliminaries, Conventions and Notations
In this section, we review some necessary definitions and relevant properties of Jacobi polynomials. In the next step, we introduce Generalized Lagrangian (GL) functions. Then, we state and prove the accuracy of GL functions and develop GLJGR collocation method. Finally, in term of GLJGR collocation method, we give a formula that expresses the derivative matrix of the mentioned functions.
2.1 Some properties of Jacobi Polynomials
A basic property of the Jacobi polynomials is that they are the eigenfunctions to a singular Sturm–Liouville problem. Jacobi polynomials are defined on and are of high interest recently bhrawyzaky2015 (); bhrawyAlzaidy (); BhrawyDoha2015 (); BhrawyAHDohaEH2015 (); Doha2015pseudo (). The following recurrence relation generates the Jacobi polynomials dohajacobi ():
(4) 
(5) 
where
The Jacobi polynomials are satisfying the following identities:
(6)  
(7)  
(8) 
(9) 
and its weight function is .
Moreover, the Jacobi polynomials are orthogonal on :
(10) 
where is the Kronecker delta function. We define the weighted space . The inner product and the norm of with respect to the weight function are defined as:
It is noted that the set of Jacobi polynomials forms a complete system.
2.2 Generalized Lagrangian (GL) functions
In this section, generally, the GL functions are introduced and the suitable formulas for the first and secondorder derivative matrices of these functions are presented.
Definition 1
Let , then, the generalized Lagrange (GL) formula is shown as delkhosh ()
(11) 
where , and is a continuous, arbitrary and sufficiently differentiable function, and .
For simplicity and are considered.
Theorem 2.1
Considering the GL functions in Eq. (11), one can exhibit the firstorder derivative matrices of GL functions as
where
Proof
: As the GL functions defined in Eq. (11), the firstorder derivative formula for the case can be achieved as follows:
(12) 
But, when , with Hopital’s rule:
This completes the proof.
Theorem 2.2
Let be the above matrix (first order derivative matrix of GL functions) and define matrix such that , , then, the secondorder derivative matrix of GL functions can be formulated as:
(13) 
Proof
: See Ref. delkhosh ().
For simplicity, from now on is considered .
2.3 Generalized Lagrangian Jacobi GaussRadau (GLJGR) collocation method
It is a wellestablished fact that a proper choice of collocation points is crucial in terms of accuracy and computational stability of the approximation by Lagrangian basis sabeh (). As a good choice of such collocation points, we can refer to the wellknown GaussRadau points in which points lie inside (a,b) and one point is clustered near the endpoints. In this sequel, we use JacobiGaussRadau nodes.
In case of GLJGR collocation method, in Eq. (11) can be considered as two approaches:
and
where is a real constant. As a matter of simplification, we write
(14) 
with the following important properties:
(15) 
(16) 
Lets speak of the first approach. Assume
then, we have:
(17) 
(18) 
(19) 
Recalling that and using formulas in Eqs. (14)–(19), we find the entry of the firstorder derivative matrix of GL functions as:
Similar to this fasion, for the second approach one can write:
(20) 
(21) 
(22) 
Now this is obvious that . Therefore, by the second approach, the entries of defined matrix can be filled as:
More specifically, Legendre, Chebyshev, and ultraspherical polynomials can be obtained as special cases from the proposed method. These cases are summarized in the following corollaries:
Corollary 1: If , we have the all the mentioned formulas of GL functions, , for Gegenbauer (ultraspherical) polynomials (symmetric Jacobi polynomials).
Corollary 2: If , we have the all the mentioned formulas of GL functions, , for Legendre case.
Corollary 3: If , we have all the mentioned formulas of GL functions, , for Chebyshev case (the 1st kind).
Corollary 4: If , we have all the mentioned formulas of GL functions, , for Chebyshev case (the 2nd kind).
Corollary 5: If , , we have all the mentioned formulas of GL functions, , for Chebyshev case (the 3rd kind).
Corollary 6: If , , we have all the mentioned formulas of GL functions, , for Chebyshev case (the 4th kind).
2.4 Operational matrix of GL functions
Defining the onecolumn vectors
for approximation of a function can write
and similarly for the derivative of this function can rewrite it as
(23) 
where is the operational matrix of derivative where
and in other words,
Taking the th row
by collocating the JacobiGaussRadau nodes () in this equation, we obtain:
This means that the operational matrix for these functions are
(24) 
where is defined in Section 2.3. Similarly, for can say
(25) 
where is defined and found in Section 2.3.
3 Numerical Method
The main objective of this section is to develop the GLJGR collocation method to solve 2DOCP. In this section, firstly, a promising function approximation method has been presented. Then, GLJGR collocation is implemented so as to accomplish the introduction of the presented numerical method for 2DOCP of interest.
3.1 Function approximation
If , , where
, is the set of GL functions product and
where is a finite dimensional vector space. For any , one can find the best approximation of in space as such that
Therefore, for any can write
(26) 
in which is a matrix of and are the relevant coefficients. is considered by the first approach mentioned in subsection 2.3 in which , and is based on the second approach and considered as .
3.2 Implementation of GLJGR collocation method for solving the 2DOCP
Now, for approximation of state and control functions
(27) 
(28) 
where
We define residual functions by substituting Eqs. (27) and (28) in Eq. (2)
(29) 
in terms of Eqs. (24) and (25) one can read
(30) 
(31) 
(32) 
, and so the Eq. (29) can be restated as
(33) 
and the initial and boundary conditions of the problem are obtained as
(34) 
and within the assumption of Eq. (28) and Eq. (27)
(35) 
As and are the and , with the aid of the characteristic of Lagrange polynomials, we write
(36) 
where and are the first column and th row of Matrix .
Then, with collocation nodes in space and collocation nodes in space, a set of algebraic equations is constructed by using Eq. (3.2) together with the conditions in Eq. (35) ( simplification of Eq. (36) is also used).
(37) 
in which can be considered as
or
(38) 
The reason why in the second case of (37) is started from 1 is that: in both examples we will consider later, in the first case of (37) we have and this makes a redundancy with the second case of (37) at . To avoid this, is started from 1.
In what follows, we have used and which are considered as the following statements. Firstly,
(39) 
as mentioned , then . Therefore,
(40) 
One can see that this assignment is based on the first approach in Section 2.3.
With the same fashion, for variable we have
(41) 
where and .
(42) 
Again, one can consider this assignment based on the second approach in Section 2.3.
At the next step we approximate the integral existing in the 2DOCP. For this, we exploit Gauss Jacobi quadratures.
For estimating an integral by Gauss Jacobi quadratures, we do as follow: