Proximal Point Algorithm for an Approximated Stochastic Optimal Control Problem

W. Grecksch^{∗}, F. Heyde^{†} and Chr. Tammer^{‡}

Institut für Optimierung und Stochastik Martin-Luther-Universität Halle-Wittenberg Theodor-Lieser-Str. 5, D-06120 Halle, Germany ∗ E-mail: grecksch@mathematik.uni-halle.de † E-mail: heyde@mathematik.uni-halle.de ‡ E-mail: tammer@mathematik.uni-halle.de

Abstract In this paper a numerical method for solving a stochastic optimal control problem under control restrictions is introduced. For this purpose a special kind of Markov chain approximation is used in order to discretize the problem. For the solution of the discrete Bellman equation a primal dual proximal point algorithm is derived.

1 Introduction

If we consider the deterministic linear-quadratic control problem

F (u) = ∫_{0}^{T} (x(t)^{T} Γx(t) + u(t)^{T} Λu(t))dt + x(T)^{T} Mx(T) → Min!

under constraints x ∈ C^{1}[0, T], u ∈ U ⊆ C[0, T] and state equation

ẋ(t) = Ax(t) + Bu(t), x(0) = x_{0},

then there is a linear operator A : C[0, T] → C^{1}[0, T] and a function a ∈ C^{1}[0, T] such that x = Au − a. So we can write

F (u) = ∥Au − a∥^{2} + ∥u∥^{2}

(if we define the norms in a proper way) and solve this optimization problem with the proximal point algorithm introduced in Benker/Hamel/Tammer [1].

We consider now the stochastic case. Let (Ω, F, P) be a probability space {F_{t}, t ≥ 0} a filtration on this probability space and w a l-dimensional F_{t}-Wiener process. The problem treated in this paper is the following:

(P):

f(x, u) = E [x(τ)^{T} Mx(τ) + ∫_{0}^{τ}(x(t)^{T} Γx(t) + u(t)^{T} Λu(t))dt] → Min! (1)

under constraints:

(i) the control process u(·) is a F_{t}-adapted stochastic process with values in IR^{m} which has to satisfy control restrictions

u(t) ∈ U = {u ∈ IR^{m} : ∥u∥ ≤ C} ∀t ∈ [0, T],  (2)

(ii) the state process xi(·) shall be F_{t}-adapted with values in IR^{d} and a solution of the linear state equation

dx(t) = (Ax(t) + Bu(t))dt + σdw(t) ∀t ∈ [0, T], x(0) = x_{0}  (3)

for given x_{0} (iii) the state process has to satisfy state restrictions x(t) ∈ G ⊂ IR^{d} ∀t ∈ [0, T], (iv) the stopping time τ is defined by

τ := \{ \begin{matrix}T & if x(t) ∈ G ∀t ∈ [0, T] \\ inf\{t ∈ [0, T] : x(t) ∉ G\} & otherwise \end{matrix} .  (4)

We make the following assumptions about the given quantities.

(A1) The matrices M, Γ ∈ IR^{d×d} and Λ ∈ IR^{m×m} are positive definite, A, B, σ are matrices with the dimensions d × d, d × m and d × l respectively, C is a positive constant.

(A2) G is a compact subset of IR^{d} and the closure of its interior.

(A3) τ as a function of the process x is continuous with probability 1 for every admissible control u.

(A4) For each piecewise constant control u(·) with a finite number of values and each initial condition x_{0} there exists a unique solution of (3) in the weak sense.

Definition 1 A F_{t}-adapted control process u which satisfies the control restrictions (2) is called an admissible control. The set of all admissible controls is denoted by U.

With the above described procedure for the deterministic case we can only compute an open loop control which does not take into consideration the previous stochastic influences. For practical purposes it is more desirable to compute feedback controls. But if we consider a feedback control u : IR^{d} → IR^{m} then the operator A which assigns the function x(·) to the feedback control u is not linear, so in this case this procedure does not work. The usual way to get feedback control laws for optimal control problems is to solve Bellman equation for the control problem. Here one has to find an approximation method which leads to a Bellman equation which is in some kind easy to solve.

2 Markov Chain Approximation

In order to solve the problem (P) numerically we have to approximate the given optimal control problem in a suitable way. For this purpose we use the Markov chain approximation method introduced in Kushner/Dupuis [4]. That means we discretize time interval and state space, the discrete time points will be nΔt^{h} (n ∈ IN) and the discrete state space will be the grid IR_{h}^{d} = {nh : n ∈ ℤ}^{d} and we replace the diffusion process x(·) by a Markov chain ξ^{h} = {ξ_{n}^{h}}_{n≥0} with values in IR_{h}^{d} and the control function u(·) by a control sequence u^{h} = {u_{n}^{h}}_{n≥0} with transition probabilities

P (ξ_{n+1}^{h} = y|ξ_{i}^{h}, u_{i}^{h}, i ≤ n, ξ_{n}^{h} = x, u_{n}^{h} = u) =: p(x, y|u).

If the condition:

a_{ii} − ∑_{j:j≠i} |a_{ij}| ≥ 0  (5)

is satisfied for all i where a := σσ^{T} then we choose the transition probabilities of the Markov chain for x ∈ IR_{h}^{d} in the following way:

Q^{h} := ∑_{i} a_{ii} − ∑_{i,j:i≠j} \frac{|a_{ij}|}{2} + 3hC∥B∥_{∞} + h max_{x∈G} |Ax|,  (6)

Δt^{h} := \frac{h^{2}}{Q^{h}} ,  (7)
p(x, x + e_{i}h|u)
:= \frac{1}{Q^{h}} \{\frac{a_{ii}}{2} − ∑_{j:j≠i} \frac{|a_{ij}|}{2} + h[(A_{i}x)^{+} + B_{i}u + C∥B∥_{∞}]\} ,  (8)
p(x, x − e_{i}h|u)
:= \frac{1}{Q^{h}} \{\frac{a_{ii}}{2} − ∑_{j:j≠i} \frac{|a_{ij}|}{2} + h[(A_{i}x)^{−} + C∥B∥_{∞}]\} ,  (9)

p(x, x + e_{i}h + e_{j}h|u) = p(x, x − e_{i}h − e_{j}h|u) := \frac{a_{ij}^{+}}{2Q^{h}} ,  (10)
p(x, x + e_{i}h − e_{j}h|u) = p(x, x − e_{i}h + e_{j}h|u) := \frac{a_{ij}^{−}}{2Q^{h}} .  (11)

All other transition probabilities are zero except

p(x, x|u) := 1 − ∑_{y:y≠x} p(x, y|u)  (12)

In the previous definitions A_{i} and B_{i} denote the i-th row of the matrices A an B, respectively; e_{i} is the i-th unit vector and for a real number z the definitions z^{+} := \frac{1}{2}(|z| + z), z^{−} := \frac{1}{2}(|z| − z) are assumed.

Remark 1 Since the matrix a is nonnegative definite there is an orthogonal matrix D so that â = DaD^{T} is a diagonal matrix satisfying condition (5). If we use coordinate transformation x̂ = Dx then the state equation for x̂ is

dx̂(t) = (DAD^{T} x̂(t) + DBu(t))dt + Dσdw_{t}  (13)

and (Dσ)(Dσ)^{T} = Dσσ^{T} D^{T} = DaD^{T} = â satisfies condition (5).

Theorem 1 Let

V (x̄, t̄) = inf_{u(·)∈U} E [x(τ)^{T} Mx(τ) + ∫ τ ; t̄(x(t)^{T} Γx(t) + u(t)^{T} Λu(t))dt|x(t̄) = x̄] (14)

be the optimal value function for the control problem (P) and let the value function for the discretized problem defined by

V^{h}(x̄, t̄) = inf_{u^{h}} E [ξ_{N^{h}(ξ^{h},t̄)}^{h}^{T} Mξ_{N^{h}(ξ^{h},t̄)}^{h} + ∑_{i=0}^{N^{h}(ξ^{h},t̄)−1}(ξ_{i}^{h}^{T} Γξ_{i}^{h} + u_{i}^{h}^{T} Λu_{i}^{h})Δτ_{i}^{h}|ξ_{0}^{h} = x̄] , (15)

where N^{h}(ξ^{h}, t̄) := inf {i ∈ IN : ξ_{n}^{h} ∉ G_{h}^{0} ∨ t̄+ iΔt^{h} ≥ T } with G_{h}^{0} = IR_{h}^{d} ∩ G^{0} and Δτ_{i}^{h} := τ_{i+1}^{h} − τ_{i}^{h} for random switch times τ_{i}^{h} with the properties τ_{0}^{h} = 0 and

P \{Δτ_{n}^{h} ≤ t|ξ_{i}^{h}, u_{i}^{h}, τ_{i}^{h}, i ≤ n\} = \{ \begin{matrix}1 − exp [\frac{−t}{Δt^{h}} ] & for t ≥ 0 \\ 0 & for t < 0 \end{matrix}

for n > 0 If the assumptions (A1) - (A4) are satisfied then V^{h}(x̄, t̄) → V (x̄, t̄) with h → 0.

The proof of this theorem is based on the results in chapter 10 of Kushner/Dupuis [4]. Kushner/Dupuis work with relaxed controls m(·), which are Borel measures on the σ-algebra B(U × [0, T]) of Borel subsets of U × [0, T] with the property that

m(U × [0, t]) = t ∀t ∈ [0, T]

for each ω ∈ Ω.

A relaxed control is called admissible with respect to some filtration F_{t} if m(A × [0, t]) is F_{t}-adapted for each Borel subset A ∈ U. The set of all admissible relaxed controls is denoted by R.

Now for any t ∈ [0, T] there is a derivative m_{t}(·) so that m_{t}(A) is F̃_{t}-adapted for any Borel subset A ∈ U and with probability one holds

m(B) = ∫_{0}^{T} ∫_{U} I_{B}(α, t)m_{t}(dα)dt ∀B ∈ B(U × [0, T])

(see Kushner [3]).

Given an admissible control u(·) we can define a relaxed control representation m^{u(·)}(·) by

m^{u(·)}(A × [0, t]) = ∫_{0}^{t} I_{A}(u(s))ds

i.e. m_{t}^{u(·)}(A) = I_{A}(u(t)).

The controlled process x(·) : [0, ∞) → IR^{r} considered in that book is a solution of the integral equation

x(t) = x̄ + ∫_{0}^{t} ∫_{U} b(x(s), α)m_{s}(dα)ds + ∫_{0}^{t} σ(x(s))dw(s)
+ ∫_{0}^{t} ∫_{Γ} q(x(s^{−}), ρ)N(dsdρ),  (16)

where w(·) is a F̃_{t}-Wiener process, N(·) is a F̃_{t}-Poisson measure for some filtration F̃_{t} and m(·) is an admissible relaxed control w.r.t this filtration.

The cost function is defined by

W (x̄, m) = Em ; x̄ [∫_{0}^{τ} ∫_{U} e^{−βs}k(x(s), α)m_{s}(dα)ds + e^{−βτ}g(x(τ))];  (17)

where Em ; x̄ means expectation under initial value x̄ and relaxed control m and τ is the first escape time of the process x(·) from the interior G^{0} of a set G. The optimal value function is defined by V(x̄) = min_{m∈R} W (x̄, m)

The functions b, σ, q, k and g shall satisfy the following assumptions:

A1.1 b(·) and σ(·) are bounded and continuous.

A1.2 q(·) is measurable and bounded, q(·, ρ) is continuous for each ρ.

A1.3 k(·) and g(·) are bounded and continuous.

A1.4 Let u(·) be an ordinary admissible control and suppose that u(·) is piecewise constant and takes only a finite number of values. Then, for each initial condition, there exists a solution to (16), where m(·) is a relaxed control representation of u(·) and this solution is unique in the weak sense.

For the set G the following conditions are demanded.

A2.1 The set G is compact and is the closure of its interior G^{0}.

A2.2 The function τ̂(·) : D^{r}[0, ∞) → [0, ∞] defined by

τ̂(ϕ) = \{ \begin{matrix}∞ & if ϕ(t) ∈ G^{0} for all t < ∞ \\ inf\{t : ϕ(t) ∉ G^{0}\} & otherwise \end{matrix}

is continuous with probability one relative to the measure induced by any solution to (16) for the initial condition x of interest.

Let now S_{h} be a discrete subset of IR^{r}, u^{h} an admissible control sequence and ξ^{h} a controlled Markov chain with transition probabilities

P (ξ_{n+1}^{h} = y|ξ_{i}^{h}, u_{i}^{h}, i ≤ n, ξ_{n}^{h} = x, u_{n}^{h} = α) =: p(x, y|α).

Then the Markov chain is called local consistent with a solution of (16) if the following conditions are satisfied:

E_{x,n}^{h,α}Δξ_{n}^{h} = b(x, α)Δt^{h} + o(Δt^{h}),  (18)
E_{x,n}^{h,α}[Δξ_{n}^{h} − E_{x,n}^{h,α}Δξ_{n}^{h}][Δξ_{n}^{h} − E_{x,n}^{h,α}Δξ_{n}^{h}]^{T} = σ(x)σ(x)^{T} Δt^{h} + o(Δt^{h}), (19)
lim_{h→0} sup_{n,ω} |Δξ_{n}^{h}| = 0,  (20)

where Δξ_{n}^{h} = ξ_{n+1}^{h} − ξ_{n}^{h}.

Moreover an Markov chain interpolation is defined by

ψ^{h}(t) = ξ_{n}^{h} for t ∈ [τ_{n}^{h}, τ_{n+1}^{h})

where τ_{n}^{h} are the random switch times defined in Theorem 1.

Let u^{h}(t) = u_{n}^{h} for t ∈ [τ_{n}^{h}, τ_{n+1}^{h}) and τ_{h} = inf{t > 0 : ψ^{h}(t) ∉ G^{0}} the first escape time of ψ^{h} from G^{0} then the cost function for the Markov chain is defined by

W^{h}(x̄, u^{h}) = Eu^{h} ; x̄ [∫_{0}^{τ_{h}} e^{−βs}k(ψ^{h}(s), u^{h}(s))ds + e^{−βτ_{h}}g(ψ^{h}(τ_{h}))];  (21)

and the optimal value function by V^{h}(x̄) = inf_{u^{h}} W^{h}(x̄, u^{h}).

Now the following theorems from chapter 10 of the book of Kushner/Dupuis [4] are of interest.

Theorem 2 (3.1) Assume (A1.1) - (A1.4), (A2.1) and (A2.2’) and use the coast function (17) with β > 0. Fix ϵ_{0} > 0, and let (x(·), m(·), w(·), N(·)) be an ϵ_{0}-optimal solution whose existence is asserted in (A2.2’). Let τ denote the first escape time from G_{0}. Then, for each ϵ > 0, there is a δ > 0 and a probability space on which are defined a pair (w^{ϵ}(·), N^{ϵ}(·)), an ordinary control u^{ϵ}(·) with values in a finite subset U_{ϵ} ∈ U and a solution x^{ϵ}(·) of (16) such that

|W (x̄, m^{u^{ϵ}}) − W (x̄, m)| ≤ ϵ.

If the set of stopping times over all controls is uniformly integrable, then the above conclusions continue to hold for β = 0.

Theorem 3 (5.2) Assume (A1.1) - (A1.4) and (A2.1). Let β > 0. If (A2.2) holds for all controls then V^{h}(x̄) → V(x̄). Let β = 0. Then, the assertion continues to hold if {τ_{h}} is uniformly integrable.

Proof of Theorem 1. Let x(·) be a process satisfying (3) for an admissible control u(·) and x(t̄) = x̄. Then we define a process x(·) with values in IR^{d+1} by

x_{j}(s) := x_{j}(s + t̄) for j = 1...d,  (22)
x_{d+1}(s) := s + t̄,  (23)

and a filtration F̃_{s} := F_{s+t̄}. For x = (x_{1}, ..., x_{d+1})^{T} we introduce the notation x_{(d)} :=

(x_{1}, ..., x_{d})^{T}. If we define x̄ := ( x̄ ; t̄ ), b(x, α) := ( \begin{matrix}Ax_{(d)} + Bα \\ 1 \end{matrix} ), σ(x) := ( \begin{matrix}σ \\ 0 \end{matrix} ),

w(s) := w(s + t̄) and u(s) := u(s + t̄) then w(s) is a F̃_{s}- Wiener process, u(·) is an admissible control w.r.t. F̃_{s} and the process x(·) satisfies the integral equation

x(t) = x̄ + ∫_{0}^{t} b(x(s), u(s))ds + ∫_{0}^{t} σ(x(s))dw(s)  (24)

which is equation (16) with q ≡ 0. Moreover we define k(x, α) := x_{(d)}^{T}Γx_{(d)} + α^{T} Λα, g(x) := x_{(d)}^{T}Mx_{(d)}, β := 0 and G := G × [0, T]. Now

V (x̄, t̄) = inf_{u(·)∈U} W (x̄, m^{u(·)})  (25)

The assumptions (A1.1) - (A1.3) follow directly from the definition of the corresponding functions and the compactness of G. The assumptions (A2), (A3) and (A4) imply the assumptions (A2.1), (A2.2) and (A1.4). Since the set of stopping times over all controls is bounded (due to τ ≤ T) this set is uniformly integrable. Therefore we can apply Theorem 2 which implies

V (x̄, t̄) = inf_{u(·)∈U} W (x̄, m^{u(·)}) = min_{m∈R} W (x̄, m) = V(x̄).  (26)

Moreover we define the discrete state space by S_{h} := IR_{h}^{d} × {nΔt^{h} : n ∈ IN_{0}} Let ξ^{h} = {ξ_{n}^{h}}_{n≥0} be a Markov chain with values in S_{h} and transition probabilities

p(x, y|u) := \{ \begin{matrix}p(x_{(d)}, y_{(d)}|u) & if y_{d+1} = x_{d+1} + Δt^{h} \\ 0 & otherwise \end{matrix}  (27)

Now we proof the local consistency of the Markov chain ξ^{h} with the process x(·)

E_{x,n}^{h,u}Δξ_{n}^{h} = ∑_{y∈S_{h}}(y − x)p(x, y|u)

= ∑_{y_{(d)}∈IR_{h}^{d}} ( \begin{matrix}y_{(d)} − x_{(d)} \\ Δt^{h} \end{matrix} ) p(x_{(d)}, y_{(d)}|u)

= Δt^{h}e_{d+1} + ∑_{i=1}^{d} e_{i}h [p(x_{(d)}, x_{(d)} + e_{i}h|u) − p(x_{(d)}, x_{(d)} − e_{i}h|u)]

+ ∑_{{i,j=1 \atop i<j}}^{d}(e_{i}h + e_{j}h) [p(x_{(d)}, x_{(d)} + e_{i}h + e_{j}h|u) − p(x_{(d)}, x_{(d)} − e_{i}h − e_{j}h|u)]

+ ∑_{{i,j=1 \atop i<j}}^{d}(e_{i}h − e_{j}h) [p(x_{(d)}, x_{(d)} + e_{i}h − e_{j}h|u) − p(x_{(d)}, x_{(d)} − e_{i}h + e_{j}h|u)]

= Δt^{h}e_{d+1} + ∑_{i=1}^{d} e_{i} \frac{h^{2}}{Q^{h}} [A_{i}x_{(d)} + B_{i}u]
= Δt^{h}b(x, u)
cov_{x,n}^{h,u}Δξ_{n}^{h} = E_{x,n}^{h,u}(Δξ_{n}^{h})(Δξ_{n}^{h})^{T} − (E_{x,n}^{h,u}Δξ_{n}^{h})(E_{x,n}^{h,u}Δξ_{n}^{h})^{T}
= ∑_{y∈S_{h}} (y − x)(y − x)^{T} p(x, y|u) − (Δt^{h})^{2}b(x, u)b(x, u)^{T}

= ∑_{y_{(d)}∈IR_{h}^{d}} ( \begin{matrix}y_{(d)} − x_{(d)} \\ Δt^{h} \end{matrix} ) ( \begin{matrix}y_{(d)} − x_{(d)} \\ Δt^{h} \end{matrix} )^{T} p(x_{(d)}, y_{(d)}|u) + o(Δt^{h})

= ∑_{y_{(d)}∈IR_{h}^{d}} ( \begin{matrix}y_{(d)} − x_{(d)} \\ 0 \end{matrix} ) ( \begin{matrix}y_{(d)} − x_{(d)} \\ 0 \end{matrix} )^{T} p(x_{(d)}, y_{(d)}|u)

+ ∑_{y_{(d)}∈IR_{h}^{d}} Δt^{h}e_{d+1} ( \begin{matrix}y_{(d)} − x_{(d)} \\ Δt^{h} \end{matrix} )^{T} p(x_{(d)}, y_{(d)}|u)

+ ∑_{y_{(d)}∈IR_{h}^{d}} ( \begin{matrix}y_{(d)} − x_{(d)} \\ Δt^{h} \end{matrix} ) Δt^{h}e_{d+1}^{T}p(x_{(d)}, y_{(d)}|u)

− ∑_{y_{(d)}∈IR_{h}^{d}}(Δt^{h})^{2}e_{d+1}e_{d+1}^{T}p(x_{(d)}, y_{(d)}|u) + o(Δt^{h})

= ∑_{i=1}^{d} e_{i}e_{i}^{T} h^{2} [p(x_{(d)}, x_{(d)} + e_{i}h|u) + p(x_{(d)}, x_{(d)} − e_{i}h|u)]

+ ∑_{{i,j=1 \atop i<j}}^{d}(e_{i} + e_{j})(e_{i} + e_{j})^{T} h^{2} ·

· [p(x_{(d)}, x_{(d)} + e_{i}h + e_{j}h|u) + p(x_{(d)}, x_{(d)} − e_{i}h − e_{j}h|u)]

+ ∑_{{i,j=1 \atop i<j}}^{d}(e_{i} − e_{j})(e_{i} − e_{j})^{T} h^{2} ·

· [p(x_{(d)}, x_{(d)} + e_{i}h − e_{j}h|u) + p(x_{(d)}, x_{(d)} − e_{i}h + e_{j}h|u)] +(Δt^{h})^{2} [e_{d+1}b(x, u)^{T} + b(x, u)e_{d+1}^{T} − e_{d+1}e_{d+1}^{T}] + o(Δt^{h})

= ∑_{i,j=1}^{d} e_{i}e_{i}^{T} h^{2} [p(x_{(d)}, x_{(d)} + e_{i}h|u) + p(x_{(d)}, x_{(d)} − e_{i}h|u)

+p(x_{(d)}, x_{(d)} + e_{i}h + e_{j}h|u) + p(x_{(d)}, x_{(d)} − e_{i}h − e_{j}h|u) +p(x_{(d)}, x_{(d)} + e_{i}h − e_{j}h|u) + p(x_{(d)}, x_{(d)} − e_{i}h + e_{j}h|u)]

+ ∑_{{i,j=1 \atop i≠j}}^{d} e_{i}e_{j}^{T} h^{2} [p(x_{(d)}, x_{(d)} + e_{i}h + e_{j}h|u) + p(x_{(d)}, x_{(d)} − e_{i}h − e_{j}h|u)

−p(x_{(d)}, x_{(d)} + e_{i}h − e_{j}h|u) − p(x_{(d)}, x_{(d)} − e_{i}h + e_{j}h|u)]

+o(Δt^{h})

= ∑_{i=1}^{d} e_{i}e_{i}^{T} \frac{h^{2}}{Q^{h}} \{a_{ii} + h[|A_{i}x| + B_{i}u + 2C∥B∥_{∞}]\}

+ ∑_{{i,j=1 \atop i≠j}} e_{i}e_{j}^{T} \frac{h^{2}}{Q^{h}} a_{ij} + o(Δt^{h})

= Δt^{h}σ(x)σ(x)^{T} + o(Δt^{h})

Equation (20) follows from |Δξ_{n}^{h}| ≤ \sqrt{2}h.

Moreover from (27) follows

P (ξ_{n+1,(d)}^{h} = y | ξ_{i}^{h}, u_{i}^{h}, i ≤ n, ξ_{n,(d)}^{h} = x, u_{n}^{h} = α) = p(x, y|α).  (28)

That means {ξ_{n,(d)}^{h}}_{n≥0} and {ξ_{n}^{h}}_{n≥0} have the same transition probabilities. Supposing ξ_{0}^{h} = x̄, then

N^{h} := inf{n ∈ IN : ξ_{n}^{h} ∉ G_{h}^{0}}
= inf{n ∈ IN : ξ_{n,(d)}^{h} ∉ G_{h}^{0} ∨ ξ_{n,d+1}^{h} ≥ T }
= inf{n ∈ IN : ξ_{n,(d)}^{h} ∉ G_{h}^{0} ∨ nΔt^{h} + t̄≥ T }

has the same conditional distribution as N^{h}(ξ^{h}, t̄).

Therefore

V^{h}(x̄) = inf_{u^{h}} Eu^{h} ; x̄ [∫_{0}^{τ_{h}} k(ψ^{h}(s), u^{h}(s))ds + g(ψ^{h}(τ_{h}))]

= inf_{u^{h}} Eu^{h} ; x̄ [ ∑_{n=0}^{N^{h}−1} k(ξ_{n}^{h}, u_{n}^{h})Δτ_{n}^{h} + g(ξ_{N^{h}}^{h})]

= inf_{u^{h}} Eu^{h} ; x̄ [ ∑_{n=0}^{N^{h}−1}(ξ_{n,(d)}^{h}^{T} Γξ_{n,(d)}^{h} + u_{n}^{h}^{T} Λu_{n}^{h})Δτ_{n}^{h} + ξ_{N^{h},(d)}^{h}^{T} Mξ_{N^{h},(d)}^{h}]

= inf_{u^{h}} E [ξ_{N^{h}(ξ^{h},t̄)}^{h}^{T} Mξ_{N^{h}(ξ^{h},t̄)}^{h} + ∑_{i=0}^{N^{h}(ξ^{h},t̄)−1}(ξ_{i}^{h}^{T} Γξ_{i}^{h} + u_{i}^{h}^{T} Λu_{i}^{h})Δτ_{i}^{h}|ξ_{0}^{h} = x̄]
= V^{h}(x̄, t̄)

Thus the assertion of Theorem 1 follows from Theorem 3.

3 Bellman Equation for the Approximated Problem

Now we define V_{n}^{h}(x) := V^{h}(x, nΔt^{h}). Then the Bellman equation for the approximated problem yields

V_{n}^{h}(x) = x^{T} Mx  (29)

for x ∉ G_{h}^{0} or n ≥ \frac{T}{Δt^{h}} and

V_{n}^{h}(x) = min_{u∈U} { ∑_{y∈IR_{h}^{d}} p(x, y|u)V_{n+1}^{h}(y) + Δt^{h}(x^{T} Γx + u^{T} Λu)}

= min_{u∈U} \{u^{T} \frac{h^{2}}{Q^{h}} Λu + c(x)u + r(x)\}
= min_{u∈IR^{m}} \{u^{T} \frac{h^{2}}{Q^{h}} Λu + c(x)u + r(x) + χ_{U}(u)\}  (30)

for n < \frac{T}{δt^{h}} and x ∈ G_{h}^{0} with

c(x) = ∑_{i} \frac{h}{Q^{h}} [V_{n+1}^{h}(x + e_{i}h) − V_{n+1}^{h}(x)] B_{i}

r(x) = \frac{1}{2Q^{h}} \{∑_{i} [V_{n+1}^{h}(x + e_{i}h) + V_{n+1}^{h}(x − e_{i}h) − 2V_{n+1}^{h}(x)]

[a_{ii} − ∑_{j:j≠i} |a_{ij}| + h |A_{i}x| + 2hC ∥B∥_{∞}]

+h ∑_{i} [V_{n+1}^{h}(x + e_{i}h) − V_{n+1}^{h}(x − e_{i}h)] A_{i}x

+ ∑_{i,j:i<j} [V_{n+1}^{h}(x + e_{i}h + e_{j}h) + V_{n+1}^{h}(x − e_{i}h − e_{j}h)] (a_{ij})^{+}

+ ∑_{i,j:i<j} [V_{n+1}^{h}(x + e_{i}h − e_{j}h) + V_{n+1}^{h}(x − e_{i}h + e_{j}h)] (a_{ij})^{−}

∑_{i,j:i<j} 2V_{n+1}^{h}(x)|a_{ij}|\} + V_{n+1}^{h}(x) + \frac{h^{2}}{Q^{h}} x^{T} Γx

Thus the optimal feedback control is determined by

u^{∗}(x) ∈ arg min_{u∈IR^{m}} \{u^{T} \frac{h^{2}}{Q^{h}} Λu + c(x)u + r(x) + χ_{U}(u)\} .  (31)

Taking into consideration that the function in braces is convex with respect to u condition (31) is equivalent to

0 ∈ ∂ \{u^{∗T} (x)\frac{h^{2}}{Q^{h}} Λu^{∗}(x) + c(x)u^{∗}(x) + r(x) + χ_{U}(u^{∗})\}
= \frac{2h^{2}}{Q^{h}} Λu^{∗}(x) + c(x) + N_{U}(u^{∗}(x))  (32)

where χ_{U}(u) = \{ \begin{matrix}0 if u ∈ U \\ ∞ otherwise \end{matrix} and N_{U}(u) is the normal cone to U at point u.

4 Proximal Point Algorithm

Now we will construct a proximal point algorithm in order to solve the inclusion (32) using the approach in Benker, Hamel and Tammer [1] which is based on the papers of Michelot and Lefebvre [5] and Idrissi, Lefebvre and Michelot [2].

First we divide the operator T in two parts: Let T_{1}(u) = \frac{2h^{2}}{Q^{h}} Λu + c(x) and T_{2}(u) = N_{U}(u). Then (32) is equivalent to

Find p_{1}^{∗}, p_{2}^{∗}, u^{∗} ∈ IR^{m} with p_{1}^{∗} ∈ T_{1}(u^{∗}), p_{2}^{∗} ∈ T_{2}(u^{∗}) and p_{1}^{∗} + p_{2}^{∗} = 0.  (33)

If we define

A := \{y = ( \begin{matrix}y_{1} \\ y_{2} \end{matrix} ) ∈ IR^{2m} : y_{1} = y_{2}\} ,

B := \{p = ( \begin{matrix}p_{1} \\ p_{2} \end{matrix} ) ∈ IR^{2m} : p_{1} + p_{2} = 0\}

and T : IR^{2m} → 2^{IR^{2m}} by

T(y) := ( \begin{matrix}T_{1}(y_{1}) \\ T_{2}(y_{2}) \end{matrix} ) ,

then (33) is equivalent to

Find y^{∗} ∈ A, p^{∗} ∈ B with p^{∗} ∈ T(y^{∗})  (34)

and A and B are orthogonal complementary subspaces of IR^{2m}. Further the operator T is maximal monotone since its components are subdifferentials of convex functions.

For the numerical algorithm we want to apply we need a condition of the form 0 ∈ T̃(·) with a suitable maximal monotone operator T̃. But if we rewrite condition (34) as 0 ∈ T̃(y^{∗}, p^{∗}) = T(y^{∗}) − p^{∗} the operator T̃ is not maximal monotone. Therefore we use the following idea of Spingarn [7].

Definition 2 Let E be a Hilbert space, T : E → 2^{E} a set valued operator, A and B orthogonal complementary subspaces of E. Then the operator T_{A} : E → 2^{E} with the property

p̂ ∈ T_{A}(ŷ) ⟺ p̂_{A} + ŷ_{B} ∈ T (p̂_{B} + ŷ_{A})  (35)

(v_{A} and v_{B} are the projections of v onto the subspaces A and B, respectively) is called the partial inverse of T with respect to A.

Remark 2 Since A and B are orthogonal complements in IR^{2m}, for every z ∈ IR^{2m} there is exactly one y ∈ A and exactly one p ∈ B with z = y + p. The variables y and p can be computed from z by y = z_{A} and p = z_{B}.

Remark 3 Using the previous remark it is easy to see that (35) is equivalent to

p ∈ T (y) ⟺ p_{A} + y_{B} ∈ T_{A}(p_{B} + y_{A})  (36)

Spingarn [7] proved the following lemma.

Lemma 4 T_{A} is maximal monotone iff T is maximal monotone.

Now (36) implies that (34) is equivalent to

Find y^{∗} ∈ A and p^{∗} ∈ B with 0 ∈ T_{A}(p^{∗} + y^{∗})  (37)

and due to Remark 2 (37) is equivalent to:

Find z^{∗} ∈ IR^{2m} with 0 ∈ T_{A}(z^{∗})
and compute y^{∗} ∈ A and p^{∗} ∈ B by y^{∗} = z_{A}^{∗}, p^{∗} = z_{B}^{∗}.

(38)

Now we can apply the following proximal point algorithm (compare Michelot/Lefebvre [5] and Idrissi/Lefebvre/Michelot [2]) for solving (38):

z^{k+1} = (I + c_{k}T_{A})^{−1}z^{k}, k = 1, 2, ...  (39)

with arbitrary starting point z^{1} and positive real c_{k}. From Rockafellar [6] we know:

If (c_{k}) is bounded away from 0 then either (z^{k}) converges weakly to a zero of T_{A} or ∥z^{k}∥ → ∞ and in this case T_{A} has no zeros.

Now we will apply this algorithm with c_{k} = 1 for all k to our problem (38).

Let z^{k} = y^{k} + p^{k} be the unique decomposition of z^{k} in elements y^{k} ∈ A and p^{k} ∈ B (i.e. y^{k} = z_{A}^{k} and p^{k} = z_{B}^{k}) then one gets

(39) ⇔ z^{k} − z^{k+1} ∈ T_{A}(z^{k+1})
⇔ y^{k} − y^{k+1} + p^{k+1} ∈ T(p^{k} − p^{k+1} + y^{k+1}).

By setting ỹ^{k} := p^{k} − p^{k+1} + y^{k+1} we obtain that (39) is equivalent to

y^{k} + p^{k} − ỹ^{k} ∈ T(ỹ^{k}).

By taking into consideration the concrete form of the operator T and the subspaces A and B we get

y_{1}^{k} + p_{1}^{k} − ỹ_{1}^{k} = \frac{2h^{2}}{Q^{h}} Λỹ_{1}^{k} + c(x)  (40)
y_{2}^{k} + p_{2}^{k} − ỹ_{2}^{k} ∈ N_{U}(ỹ_{2}^{k})  (41)

Now (40) is equivalent to

ỹ_{1}^{k} = (I + \frac{2h^{2}}{Q^{h}} Λ)^{−1} [y_{1}^{k} + p_{1}^{k} − c(x)]  (42)

and (41) is equivalent to

ỹ_{2}^{k} ∈ U and ⟨y_{2}^{k} + p_{2}^{k} − ỹ_{2}^{k}, u − ỹ_{2}^{k}⟩ ≤ 0 ∀u ∈ U  (43)

i.e. if we choose

ỹ_{2}^{k} = P_{U}(y_{2}^{k} + p_{2}^{k})  (44)

then condition (43) is satisfied and consequently condition (41) is satisfied, too. Since A and B are orthogonal complementary linear spaces y^{k+1} = P_{A}(ỹ^{k}) means that y^{k+1} ∈ A and there is a p ∈ B with ỹ^{k} = y^{k+1} + p. If we add the components of this equation we get

ỹ_{1}^{k} + ỹ_{2}^{k} = y_{1}^{k+1} + p_{1} + y_{2}^{k+1} + p_{2}.

Since from p ∈ B follows p_{1} + p_{2} = 0 and from y^{k+1} ∈ A follows y_{1}^{k+1} = y_{2}^{k+1} we get

y_{1}^{k+1} = y_{2}^{k+1} = \frac{1}{2}(ỹ_{1}^{k} + ỹ_{2}^{k}).  (45)

Further we get

p_{1}^{k+1} = −p_{2}^{k+1} = p_{1}^{k} + \frac{1}{2}(ỹ_{2}^{k} − ỹ_{1}^{k}).  (46)

Since y_{2}^{k} = y_{1}^{k} and p_{2}^{k} = −p_{1}^{k} for all k we only need to compute y_{1}^{k}, p_{1}^{k} explicitly and thus equations (42), (44) - (46) yield the following proximal point algorithm for the solution of inclusion (32):

Step 1: Choose arbitrary starting points y_{1}^{1}, p_{1}^{1} ∈ IR^{m}.

Step 2: For given y_{1}^{k}, p_{1}^{k} ∈ IR^{m} compute

ỹ_{1}^{k} = (I + \frac{2h^{2}}{Q^{h}} Λ)^{−1} [y_{1}^{k} + p_{1}^{k} − c(x)] ,
ỹ_{2}^{k} = P_{U}(y_{1}^{k} − p_{1}^{k}),
y_{1}^{k+1} = \frac{1}{2}(ỹ_{1}^{k} + ỹ_{2}^{k}),
p_{1}^{k+1} = p^{k} + \frac{1}{2}(ỹ_{2}^{k} − ỹ_{1}^{k}).

Step 3: Repeat Step 2 until ∥y_{1}^{k+1} − y_{1}^{k}∥^{2} + ∥p_{1}^{k+1} − p_{1}^{k}∥^{2} < ε^{2}.

References

[1] H. Benker, A. Hamel, and C. Tammer. A proximal point algorithm for control approximation problems. Mathematical Methods of Operations Research, 43(3):261–280, 1996.

[2] H. Idrissi, O. Lefebvre, and C. Michelot. A primal-dual algorithm for a constrained fermat-weber problem involving mixed norms. Recherche operationnelle / Operations Research, 22(4):313–330, 1988.

[3] H. Kushner. Weak Convergence Methods and Singularly Perturbed Stochastic Control and Filtering Problems, volume 3 of Systems and Control. Birkhaeuser, Boston, 1990.

[4] H. J. Kushner and P. G. Dupuis. Numerical Methods for Stochastic Control Problems in Continuous Time. Springer, New York, 1992.

[5] C. Michelot and O. Lefebvre. A primal-dual algorithm for the fermat-weber problem involving mixed gauges. Mathematical Programming, 39:319–335, 1987.

[6] R. Rockafellar. Monotone operators and the proximal point algorithm. SIAM J.Contr.Optim., 14:877–898, 1976.

[7] J. Spingarn. Partial inverse of a monotone operator. Applied Mathematics and Optimization, 10(3):247–265, 1983.