Numerical solution of a parabolic optimal control problem in economics

In this paper, we construct and study finite difference approximations of optimal control problems arising in the mathematical modeling of certain management and financial problems. For the definiteness we consider two-dimensional in space variables problems, although all the theoretical results remain valid for any dimension. The numerical calculation results are presented for 1D and 2D problems. To approximate the state equations we use implicit (backward Euler) and fractional steps (operator splitting) methods. The solutions of the constructed mesh state equations are strictly positive and keep an analogue of the mass balance condition. We provide the existence results and first order optimality conditions for the corresponding mesh optimal control problems. We use several iterative methods to implement the constructed nonlinear optimization problems and compare their effectiveness.


Lemma 1.
For given functionsᾱ ∈ V α and m 0 ∈ L 2 (Ω), m 0 0, problem (1) has a unique weak solution u ∈ W (0, T ), which is non-negative: m(x, t) 0 a.e. Q T , and satisfies mass conservation property: Remark 1. We impose the additional constraint ∥ᾱ∥ Vα C α to ensure the boundedness of the set K, because the objective functional J(m,ᾱ) is not coercive.

Approximation
We approximate the optimal control problem using finite difference method on the uniform meshes ω x = {x ij = (ih 1 , jh 2 ), 0 i N 1 , 0 j N 2 } and ω t = {t j = jτ, j = 0, . . . , M, M τ = T }. Let V x be the vector space of mesh functions defined on ω x and (., .) means Euclidian scalar product in V x . The space V 0 Let us define the mesh operators A 01 and A 11 (β) for a β ∈ V 0 x 1 by the following equalities (for all j = 0, 1, . . . N 2 ): Similarly we define the mesh operators A 02 and A 12 (β) for a β ∈ V 0 x 2 . Note that the finite difference expression A 0 m = ( A 01 +A 02 ) m approximates the diffusive part of the state equation, and the finite difference expression ( A 11 (α 1 ) + A 12 (α 2 ) ) m approximates its advective part presented in the form −div(ᾱ + m) + div(ᾱ − m). We consider two variants of the approximation of the state equation: and fractional steps scheme (cf., e.g., [3]) with initial condition m 0 = m 0 ∈ V x for both schemes. We approximate the objective functional J(m,ᾱ) by the discrete objective function where ρ ∈ V x equals to ρ = ρ 1 ρ 2 and the components of vectors has at least one solution. I(m,ᾱ) is not coercive. The properties (8) and (9) of the state function m are crucial for proving this result.

First order optimality conditions and iterative solution methods
Let the function g(x, t, m) be continuously differentiable with respect to m for (x, t) ∈Q T , m > 0.
In the case of implicit approximation (5) of the state equation the adjoint problem is given by the following system of equations: with initial condition λ M +1 = 0.
In the case of fractional steps approximation (6) of the state equation the adjoint problem is given by the following system of equations: for r = 1, 2 and all k = 1, 2, . . . , M .
We use the gradient information provided by the formulas (5), (11), (12) or (6), (13), (14) for implementing the iterative solution methods of the form Aboveᾱ (s) is s-th iteration, H (s) is a preconditioning matrix and θ (s) is the iterative parameter. We consider two variants of the preconditioners: constructed by quasi-Newton method to approximate Hessian of I(m(ᾱ (s) ),ᾱ (s) ) and H (s) ≡ H = diag (ρm) (called descent method below . Iterative parameter θ (s) is searched by line-search algorithm. In both iterative algorithms the multivalued function sign α k r at the points α k r = 0 was taken to be equal 0. As an initial guess we used α (0) = 0, which was found to be the most reasonable after numerical experiments.

Numerical results
We performed numerical experiments for several 1D and 2D problems. Some of them, taken from the articles [1] and [2], are mathematical models of real problems in the field of finance and management. Other optimal control problems were chosen to demonstrate the applicability of the proposed numerical methods for the described class of objective functionals. In 1D problems we took state problem ∂m ∂t − 0.07 (as in [1]) and objective functional J(m, α) = . We approximated the state equation in 1D problems by implicit finite difference scheme on the mesh with h = 1/64, τ = 1/64 and solved by quasi-Newton and descent methods. In the figures fig.1 -fig.3 we represent the rate of decrease of the objective function in relation to the number of iterations. In 2D problems state problem was (1)