A PSEUDOSPECTRAL METHOD FOR BUDGET-CONSTRAINED INFINITE HORIZON OPTIMAL CONTROL PROBLEMS

Abstract. In this paper a generalization of the indirect pseudo-spectral method, presented in [17], for the numerical solution of budget-constrained infinite horizon optimal control problems is presented. Consideration of the problem statement in the framework of weighted functional spaces allows to arrive at a good approximation for the initial value of the adjoint variable, which is inevitable for obtaining good numerical solutions. The presented method is illustrated by applying it to the budget-constrained linear-quadratic regulator model. The quality of approximate solutions is demonstrated by an example.


Introduction
In the last decades, developing numerical solution methods for infinite horizon optimal control problems has emerged a lot of attention and a plenty of new results were obtained, cf. [5], [7], [9], [17]. Among them, both direct and indirect solution methods were established. In the present paper, we extend the indirect pseudospectral method, introduced in [17], onto the class of budget-constrained infinite horizon optimal control problems. The plausibility of investigating this class of problems can be seen e.g. in [6], [15], where for instance an upper bound on the total amount of available and/or tolerated chemotherapeutic agent or vaccine is postulated and described mathematically by means of an isoperimetric control constraint. This economic and medical evidence is of great importance especially for the low-income countries.
Theoretical results concerning the existence of optimal solutions and necessary optimality conditions in form of a Pontryagin Type Maximum Principle for this class of problems were considered in [11], [13].
A big amount of numerical methods, direct and also indirect ones, consider instead of an infinite time horizon a long but finite horizon (0, T ), i.e. they cut the interval (0, ∞) at t = T . In contrast to "cutting the horizon" methods, the key idea of the presented indirect pseudospectral method consists in expanding the unknown solution in a finite Fourier-Laguerre series and using the fact that the Laguerre polynomials build a complete orthonormal system in the weighted Lebesgue space L 2 (IR + , e −ρt ). The advantages of consistent using the framework of weighted functional spaces has been addressed in [12] and we continue to use this here too.
The paper is structured as follows. Section 1 contains the necessary definitions, whereas Section 2 presents the optimal control problem statement in weighted functional spaces. Section 3 discusses the available necessary optimality conditions for the considered class of problems. Sections 4 and 5 contain the description of the method and, respectively, its application to a budget-constrained linear-quadratic regulator model. Finally, some conclusions finish the paper.

Definitions of weighted functional spaces
Let us write (0, ∞) = IR + . We denote by M n (IR + ), L n p (IR + ) and C 0,n (IR + ) the spaces of all vector functions x : IR + → IR n which are Lebesgue measurable, in the pth power Lebesgue integrable or continuous, respectively, see [3], p. 146 and pp. 285; [4], pp. 228. The Sobolev space W 1,n p (IR + ) is then defined as the space of all functions x : IR + → IR n which belong to L n p (IR + ) and admit distributional derivativeẋ ( [18], p. 49) belonging to L n p (IR + ) as well.

Problem statement
The optimal control problem being considered in the present paper is: where the next assumption has to be fulfilled: Assumption 1. Let U be a compact convex subset of IR m , ν 1 (t) := e −ρt , ρ > 0, i.e. ν 1 is a density function, and let ν 0 be a weight function as defined in Section 2.
With the denotation an equivalent formulation of (P ) B ∞ can be given, wherein only one density function ν 1 (·) appears. We use this equivalent formulation of (P ) B ∞ , namely and introduce the following basic assumptions for (P ) B ∞ and (P ) B ∞ : Assumption 2. Let the matrix functions A : . . , n}) satisfy D ∈ L n 2 (IR + , ν 1 ) and Dν 1 ∈ L n ∞ (IR + ). Further, the functions A, B, C, D are assumed to be continuous. Assumption 3. Let r 0 : IR + × IR n × IR m → IR + , and correspondingly r : IR + × IR n × IR m → IR + , be continuous in the first argument, continuously differentiable in the second and third, and convex in the second and third argument together.
The functions x and u are called the state and the control function respectively. The integral in (2) is understood in Lebesgue sense. The fact that we have to distinguish between different integral types in infinite horizon optimal control problems was discussed in details in [10], [14].

Necessary optimality conditions
We refer some additional assumptions, stated in [13]: Assumption 4. Let ρ > 0 be chosen in such a way that for any T ≥ τ the solution ζ of the initial value where ζ T ∈ IR n is arbitrary, belongs to W 1,n 2 ([τ, ∞), ν 1 ) (dominance condition). Assumption 5. For the optimal solution (x * , u * ) let the inclusions Assumption 6. Let the admissible set be non empty and contain at least one admissible pair besides the optimal one.
Under Assumptions 1 -6, the problem (P ) B ∞ possesses an optimal solution, cf. [11]. Further, the following necessary optimality conditions in form of a Pontryagin Type Maximum Principle hold true, cf. [13]. where This theorem provides the kernel for the numerical method presented in the next section.

Description of the method
We now intend to describe an indirect numerical solution method, i.e. it acts according to the scheme "first optimize, then discretize", and represents an extended version of the indirect pseudospectral method presented in [17], for a class of budget-constrained infinite horizon optimal control problems. For incorporating the method into existing numerical methods compare scheme in Fig. 4. For a density function ν 1 , cf. Assumption 1, considering  as scalar product in the weighted Lebesgue space L 2 (IR + , ν 1 ), a spectral method is developed for solving the control problem (P ) B ∞ numerically. We propose a discrete approximation scheme based on global collocation using the Gauss-Laguerre quadrature formula. For this density with a fixed ρ > 0, we apply the Gram-Schmidt orthogonalization procedure to the system of monomials {1, t, t 2 , ...} with respect to the weighted scalar product (15). As a result we obtain the orthonormal sequence of Laguerre polynomials, which build a complete orthonormal basis in the space L 2 (IR + , ν 1 ), cf. [1]. Remind that the element x a representing the best approximation of element x of a Hilbert space H onto a subspace H N ⊂ H is characterized by the orthogonality condition x a ⊥H ⊥ N . Therefore, similar to Galerkin methods, see e.g. [19], we can write the approximate solution of the control problem as a linear combination of Laguerre polynomials L k (·). Collocation points {t 1 , . . . , t N } are obtained from the Gauss-Laguerre quadrature formula where {t 1 , . . . , t N } are the zeros of the Laguerre polynomial of degree N with corresponding weights {ω 1 , . . . , ω N }, see [1]. Then the integration formula of Gauss-Laguerre type is exact for all polynomials P ∈ P 2N −1 of degree up to 2N − 1 which means that see [2]. Following Lemma 5.5, p. 135 in [1], we obtain the following estimate of the truncation error made when using the Fourier-Laguerre approximation under appropriate assumptions on f : , r ≥ 1.
We now apply the pseudospectral method to the system of necessary optimality conditions from the Pontryagin type Maximum Principle provided by Theorem 1. Let (x(·), u(·), y 0 (·), µ 0 ), ν 1 (t) = e −ρt , ρ > 0, be given, satisfying all the conditions of Theorem 1, cf. Section 3. Further let x(·) ∈ C 1,n (R + ) ∩ W 1,n 2 (IR + , ν 1 ), u(·) ∈ C 0,m (IR + ) ∩ L m 2 (IR + , ν 1 ) and the adjoint y 0 (·) ∈ C 1,n (IR + ) ∩ W 1,n 2 (IR + , ν 1 ) be satisfied. Assuming that the control constraints do not become active 1 we can calculate the control directly from the maximum condition in a unique way as follows whose solution leads to a reduction of the number of optimization variables by elimination of optimization variables corresponding to the unknown control. Now, the state x(·) = (x 1 (·), . . . , x n (·)) and the adjoint y 0 (·) = (y 0 1 (·), . . . , y 0 n (·)) are approximated for each t > 0 component-wisely by the following finite Fourier series expansions: with unknown coefficients a sk , c sk ∈ IR, where L k (·) is the Laguerre polynomial of degree k. Evaluation and differentiation of the series at the collocation points under consideration of L N (t i ) = 0 yield the following expressions: The discretization of the Hamiltonian system dynamics is obtained by evaluating for s ∈ {1, . . . , n} the state equation (12) and the current adjoint equation (K) at the collocation points: which means that the approximate solutions fulfill the state and adjoint equations exactly. Considering the initial time t 0 = 0 as an additional non-collocated point, one arrives at the discretization of the initial condition given by x(0) = x s0 = N k=0 a sk L k (t 0 ). We want to emphasize the following advantage of an indirect method concerning good estimates for the missing initial condition for the adjoint y 0 . Based on the Pontryagin type Maximum Principle, cf. Theorem 1, the transversality condition in form of an inclusion can be used in the framework of the spectral method to conclude the existence of the integral ∞ 0ẏ 0 s (t)ν 1 (t)dt < ∞. Multiplying the differential equation for the current adjoint variable y 0 by the density function ν 1 and integrating this over the interval (0, ∞) one obtains On the other hand, integrating the left hand side of the last equation by parts and having in mind thaṫ Since the right hand sides of the last two equations are equal to each other we arrive at Applying the Gauss-Laguerre quadrature scheme (16) to the right hand side of (26) and taking y 0 as a discretized transversality condition. Thus the budget-constrained infinite horizon optimal control problem is replaced by the following system of 2n(N + 1) + 1 linear equations and 2n(N + 1) + 1 unknowns gathered in the vector a 10 , . . . , a 1N , . . . , a n0 , . . . , a nN ; c 10 , . . . , c 1N , . . . , c n0 , . . . , c nN , µ 0 T : x s0 = N k=0 a sk L k (t 0 ), s = 1, . . . , n;  , a 1N , . . . , a n0 , . . . , a nN ; c 10 , . . . , c 1N , . . . , c n0 , . . . , c nN , µ 0 T ∈ IR 2n(N +1)+1 , if the coefficient matrix of the system is regular.

Application of the method to the budget-constrained linear-quadratic regulator model
Let us consider a linear-quadratic regulator model with an isoperimetrical constraint in order to illustrate the applicability of the described pseudospectral method. Thus the optimal control is to minimize the objective functional with respect to all processes (x, u) ∈ W 1,1 2 (R + , ν 1 ) × L 1 2 (R + , ν 1 ) satisfying the constraintṡ where the density function ν is set as follows: ν 1 (t) := e −ρt , 4 < ρ < 2 + 2 √ 2 and α > − ρ 2 . Here, additional pointwise control constraints in (34) with a compact control set are needed in order to be able to apply the Pontryagin Type Maximum Principle. Therefore, these are chosen in dependence on the initial condition in such manner that they do not become active. We refer to the work [1], p. 193 ff., where this problem is solved in case of active control constraints as well, but without budget-constraints. For the introduced problem, the assumptions of the Pontryagin Type Maximum Principle, cf. Theorem 1, are satisfied for the analytically computed optimal solution and for the given values of ρ and λ 0 = 1. Moreover, the conditions of the existence Theorem, cf. Theorem 1 in [11], are fulfilled as well. Now, with the aim of comparison, we apply the indirect pseudo-spectral method from the previous section to arrive at the numerical solution of the budgetconstrained linear-quadratic optimal control problem. Doing so we arrive Solving the above system of linear equations and substituting the solution vector into the finite Fourier expansions (19) of the optimal solution, we obtain the following numerical solution of the linear-quadratic regulator problem for different numbers of collocation points N and ρ = 5, α = −2, d = 1, cf. Figures 2 -4. The exact optimal solution is given in black, whereas the magenta, green, cyan, blue and red lines stand for the approximate solutions with 3, 5, 8, 10 and 12 collocation points respectively. The numerical values for the objective J ∞ (x N , u N ) as well as the deviation from the optimal value J ∞ (x * (·), u * (·)) = 2( √ 2+3) are given in Table 1 in dependence on the number N of collocation points. Obviously This emphasizes our conjecture about the quality of the approximations caused by the fast convergence of Fourier series, where already several Fourier coefficients are sufficient to deliver a good approximation of a function by a finite Fourier series. Although we do not provide a rigorous convergence result of the suggested pseudospectral method here, the results observed above as well as solutions for other optimal control problems provided by pseudo-spectral method such as inverse pendulum problem or Lotka-Volterra optimal control problem, cf. [1], allow us to assume that the chosen strategy would lead to comparable convergence results also for other budgetconstrained optimal control problems.

Conclusions
We have presented the extension of the indirect pseudospectral method for infinite horizon optimal control problems which arises in a natural way from the method developed in [17] through incorporating an additional isoperimetric, or also called budget-constraint, into the discretization scheme applied to the system of necessary optimality conditions in form of a Pontryagin Type maximum Principle, established in [13]. The quality and the convergence rate of the method presented here is comparable with those of the method in [17] and this fact gives us a reason to presume similar quality of convergence for the whole class of budget-constrained control problems, although we do not present any rigorous convergence proof here. Besides the convergence proof, one of the most important focuses of upcoming research will be the generalization of the method to nonlinear infinite horizon optimal control problems.