Accelerated Point-wise Maximum Approach to Approximate Dynamic Programming

TL;DR

This paper introduces a novel approximate dynamic programming method that iteratively constructs lower bounds on the value function using a point-wise maximum approach, employing a gradient ascent algorithm with convergence guarantees.

Contribution

It proposes a new approach to approximate dynamic programming that maximizes point-wise maximums of value functions, with a gradient ascent method and convergence analysis.

Findings

Computes tighter sub-optimality bounds than existing methods.

Requires less computation time for the same or better bounds.

Provides convergence guarantees for the proposed algorithm.

Abstract

We describe an approximate dynamic programming approach to compute lower bounds on the optimal value function for a discrete time, continuous space, infinite horizon setting. The approach iteratively constructs a family of lower bounding approximate value functions by using the so-called Bellman inequality. The novelty of our approach is that, at each iteration, we aim to compute an approximate value function that maximizes the point-wise maximum taken with the family of approximate value functions computed thus far. This leads to a non-convex objective, and we propose a gradient ascent algorithm to find stationary points by solving a sequence of convex optimization problems. We provide convergence guarantees for our algorithm and an interpretation for how the gradient computation relates to the state relevance weighting parameter appearing in related approximate dynamic programming approaches. We demonstrate through numerical examples that, when compared to existing approaches, the algorithm we propose computes tighter sub-optimality bounds with less computation time.

Figures

Equations

$$V^{\ast }(x):=\underset{\pi \in \Pi }{\inf }$$

$$x_{t+1}=g\left(x_t,u_t,{\xi }_t\right),\forall t\ge 0,$$

$$u_t=\pi (x_t),\forall t\ge 0,$$

$$x_t\in \mathcal{X},u_t\in \mathcal{U},\forall t\ge 0,$$

$$x_0=x.$$

$$V^{\ast}(x)=\underbrace{\inf\limits_{u\in\mathcal{U}}\,\,\overbrace{\Big{\{}\,l(x,u)\,+\,\gamma\,\mathbb{E}\left[V^{\ast}\left(g(x,u,\xi)\right)\right]\,\Big{\}}}^{\left(\mathcal{T}_{u}V^{\ast}\right)(x,u)}}_{\left(\mathcal{T}V^{\ast}\right)(x)}\,,\hskip 2.84544pt\forall\,x\!\in\!\mathcal{X}\,.$$

$${\pi }^{\ast }(x)=\left\{\underset{u\in \mathcal{U}}{\mathrm{argmin}}l(x,u)+\gamma \mathbb{E}\left[V^{\ast }\left(g(x,u,\xi )\right)\right]\right\}.$$

$$\underset{V}{\max }$$

$$V\in \mathcal{F}\left(\mathcal{X}\right)$$

$$V(x)\le \left({\mathcal{T}}_{u}V\right)(x,u),\forall x\in \mathcal{X},u\in \mathcal{U}$$

$$\underset{V_1,\dots ,V_J}{\max }$$

$$V_j\in \mathcal{F}\left(\mathcal{X}\right),j=1,\dots ,J$$

$$V_{\mathrm{PWM}}(x)\le \left({\mathcal{T}}_u{V}_{\mathrm{PWM}}\right)(x,u),\forall x\in \mathcal{X},u\in \mathcal{U}$$

$$V_{\mathrm{PWM}}(x)=\underset{j=1,\dots ,J}{\max }{V}_j(x),\forall x\in \mathcal{X}$$

$$\hat{\mathcal{F}}\left(\mathcal{X}\right)=\left\{{\alpha }^{⊺}\varphi (x)|\alpha \in {\mathbb{R}}^K\right\},\text{with}\varphi (x)=\left[\begin{array}{c}{\varphi }_1(x)\\ ⋮\\ {\varphi }_K(x)\end{array}\right].$$

$${\bar{V}}_{\mathrm{obj}}(x)=\underset{{\bar{\alpha }}_{\mathrm{obj}}\in {\mathcal{A}}_{\mathrm{obj}}}{\max }{\bar{\alpha }}_{\mathrm{obj}}^{⊺}\varphi (x),\forall x\in \mathcal{X},$$

$${\bar{V}}_{\mathrm{con}}(x)=\underset{{\bar{\alpha }}_{\mathrm{con}}\in {\mathcal{A}}_{\mathrm{con}}}{\max }{\bar{\alpha }}_{\mathrm{con}}^{⊺}\varphi (x),\forall x\in \mathcal{X},$$

$$f_{\mathrm{pwm}}\left(\alpha \right)=\frac{1}{N_c}\sum _{i=1}^{N_c}\max \left\{{\alpha }^{⊺}\varphi (x_{c,i}),{\bar{V}}_{\mathrm{obj}}(x_{c,i})\right\},$$

$$\underset{\alpha \in {\mathbb{R}}^K}{\max }$$

$${\alpha }^{⊺}\varphi (x)\le \left({\mathcal{T}}_u{\bar{V}}_{\mathrm{con}}\right)(x,u),\forall x\in \mathcal{X},u\in \mathcal{U}.$$

$$\alpha\in BI\left(\mathcal{A}_{\mathrm{con}}\right)\subseteq\mathbb{R}^{K}\hskip 5.69046pt\Rightarrow\hskip 5.69046pt\alpha\text{ satisfies \eqref{eq:canonical_inner_problem:pwm_ineq}}\,,$$

$${\partial }^{+}{f}_{\mathrm{pwm}}\left(\alpha \right)\ni \frac{1}{N_c}\sum _{i=1}^{N_c}\{\begin{array}{ll}\varphi (x_{c,i})& \text{if}{\alpha }^{⊺}\varphi (x_{c,i})\ge {\bar{V}}_{\mathrm{obj}}(x_{c,i})\\ 0& \text{if}{\alpha }^{⊺}\varphi (x_{c,i})<{\bar{V}}_{\mathrm{obj}}(x_{c,i})\end{array},$$

$$\underset{\alpha \in {\mathbb{R}}^K}{\max }\left\{{\alpha }^{⊺}\varphi (x_{c,i});\alpha \in BI({\mathcal{A}}_{\mathrm{con}})\right\}$$

$$\underset{k\to \infty }{\lim }\left(\underset{\alpha \in BI({\mathcal{A}}_{\mathrm{con}})}{\min }{\left(\alpha -{\alpha }^{(k)}\right)}^{⊺}{d}^{(k)}\right)=0.$$

$${\alpha }^{⊺}\varphi (x)\le \left({\mathcal{T}}_u\left(\max \left\{{\alpha }^{⊺}\varphi (x),{\bar{V}}_{\mathrm{con}}(x)\right\}\right)\right)(x,u),$$

$$\underset{\alpha ,\beta \in {\mathbb{R}}^K}{\max }$$

$$\alpha \in BI({\mathcal{A}}_{\mathrm{con}}),\beta \in BI({\mathcal{A}}_{\mathrm{con}}).$$

$${\partial }^{+}\left(\max \left\{{\alpha }^{⊺}\varphi (x),{\beta }^{⊺}\varphi (x),{\bar{V}}_{\mathrm{obj}}(x)\right\}\right)$$

$$=$$

$$x_{t+1}=x_t-0.5u_t,x_0\sim \mathcal{N}(0,10),$$

$$|u|\le 1,\gamma =0.95,l(x,u)=x^2+0.1u^2,$$

$$\underset{\hat{V}\in \hat{\mathcal{F}}\left(\mathcal{X}\right)}{\mathrm{argmax}}\left\{{\int }_{\mathcal{X}}\hat{V}(x)\nu (\mathrm{d}x);\hat{V}(x)\le \mathcal{T}\hat{V}(x),\forall x\in \mathcal{X}\right\},$$

$$x_{t+1}=A{x}_t+B_u{u}_t,$$

$$l(x,u)+\gamma \mathbb{E}\left[\underset{k=1,\dots ,J}{\max }{V}_k\left(g(x,u,\xi )\right)\right].$$

$$V_j(x)\le l(x,u)+\gamma \underset{k=1,\dots ,J}{\max }\mathbb{E}\left[V_k\left(g(x,u,\xi )\right)\right],$$

$$\forall x\in \mathcal{X},u\in \mathcal{U},j=1,\dots ,J.$$

$$V_j(x)\le l(x,u)+\gamma s_V^2,\forall (x,u,s_V)\in \mathcal{S},$$

$$\mathcal{S}=\left\{x,u,s_V|\begin{array}{c}x\in \mathcal{X},u\in \mathcal{U},s_V\in \mathbb{R},\\ s_V^2\ge \mathbb{E}\left[V_k\left(g(x,u,\xi )\right)\right]\forall k=1,\dots ,J\end{array}\right\},$$

Full text

I Introduction

Many important challenges in science and engineering can be cast in the problem formulation of infinite horizon stochastic optimal control (SOC), from climate control of a building [1] to control of a cell population [2]. The goal of such problems is to find a state feedback policy that minimizes an infinite-horizon discounted cost function. For a general SOC problem instance, the solution, i.e., the optimal policy, is typically characterized by the theory of dynamic programming (DP) [3, 4, 5]. However, in all but a few special cases, solving the SOC directly or applying the DP theory is intractable due to the so-called curse of dimensionality. As such, an extensive body of literature has proposed approximation techniques for computing sub-optimal solutions to SOC problems, ranging from model-free and simulation-based algorithms [6], to model-based approaches [7, 8, 9, 10]. Any technique based on DP theory falls in the category of Approximate DP (ADP), see [11, 12, 13] for an overview. Although the optimal policy is intractable to compute, techniques have been developed to bound the sub-optimality of an approximate policy. These type of bounds provide the designer with valuable information about the potential benefit of synthesizing and evaluating alternative policies. In this paper we consider and propose approaches that provide sub-optimality bounds based on the so-called Linear Programming (LP) approach to ADP [14].

The LP approach to ADP provides sub-optimality bounds by computing approximations that are lower bounds of the so-called value function, i.e., the solution of the Bellman equation for DP [3]. The LP approach parameterizes approximate value functions as a linear combination of fixed basis functions and uses the so-called Bellman inequality to restrict consideration to only those linear combinations that are point-wise lower bounds of the value function. To compute an approximate value function the designer specifies the regions of the state space that are of interest via the state relevance weighting function, and then solves an optimization problem to find a linear combination of basis functions that maximizes the integral with respect to this weighting. The LP approach was first proposed for finite state and input spaces in [15], and equipped with theoretical guarantees in [16]. The authors of [16] also provide a discussion on the importance and difficulty of choosing the state relevance weighting to give the best performance and lower bound. An iterated version of the Bellman inequality was proposed in [17] and used to compute tighter lower bounds, however, the topic of choosing the state relevance weighting is not addressed. The subsequent works [18, 19] avoid the need for a state relevance weighting by focusing on the design of policies rather than providing tighter lower bounds. In [20] the authors use sum-of-squares programming techniques to compute high-order polynomial approximate value functions using the iterated Bellman inequality. The use of high-order polynomials would reduce the difficulty of choosing the state relevance weighting, however, the optimization problem to solve becomes formidable.

Given a family of lower bounding approximate value functions, computed via the LP approach to ADP, taking a point-wise maximum over the family will yield the same or better approximation of the value function. However, to the best of our knowledge, no algorithm exists that explicitly aims to maximize this point-wise maximum combination; attempting to do this directly leads to a non-convex optimization problem. The benefit of a point-wise maximum combination is empirically demonstrated in [17] for a simple example, with the set of state relevance weighting parameters hand-picked using problem-specific insight. In our previous work [21], we proposed a problem formulation with the point-wise maximum combination used in the Bellman inequality. The formulation was used to develop an iterative algorithm for computing lower bounding approximate value functions, however, the quality of the approximation, comparable with that of [17], still relies on the designer choosing a sequence of state relevance weightings. The algorithm proposed in [22] also uses the point-wise maximum combination in the Bellman inequality, and the authors propose an algorithm that computes the sequence of state relevance weightings based on simulating the evolution of the system in a so-called forward pass. They consider a finite horizon setting and it is not clear how to extend the algorithm to an infinite horizon setting. A variety of other ADP algorithms compute lower bounds using theoretical tools different from the Bellman inequality, for example [23, 24, 25], each with its advantages and disadvantages, and none of which are similar to the algorithms we propose.

In this paper, we propose a formulation that explicitly aims to maximize a point-wise maximum of lower bounding approximate value functions, and we use this to develop novel algorithms for computing sub-optimality bounds. In particular, the contributions of the paper are:

• •

We introduce the point-wise maximum formulation of DP and prove that it is equivalent to the LP approach under standard assumptions for SOC problems.

• •

We propose a gradient ascent algorithm for finding approximate solutions to the point-wise maximum formulation, and prove that it converges to stationary points. We provide an interpretation of the gradient ascent steps as an algorithmic choice of the state relevance weighting for the non-convex problem of maximizing the point-wise maximum objective.

• •

We propose an algorithm for computing initial conditions for the gradient ascent and prove that it converges in finite iterations for any tolerance. This algorithm is required because, for non-convex problems, the quality of a gradient ascent solution is influenced by the initial condition.

In support of the contributions, we provide numerical results to demonstrate the sub-optimality bounds achieved and computation time required. Section II presents the point-wise maximum DP formulation. Section III introduces the approximation methods and our proposed algorithms. Section IV demonstrates the performance through numerical examples.

II Dynamic Programming (DP) Formulation

II-A Stochastic Optimal Control Formulation and Assumptions

We consider discrete time, infinite horizon, discounted cost, stochastic optimal control problems over continuous state and action spaces. The state of the system at time $t$ is denoted by $\smash{x_{t}\!\in\!\mathcal{X}\subseteq\mathbb{R}^{n_{x}}}$. The system state is influenced by the control decisions $\smash{u_{t}\!\in\!\mathcal{U}\subseteq\mathbb{R}^{n_{u}}}$, and by the stochastic exogenous disturbance $\smash{\xi_{t}\!\in\!\Xi\subseteq\mathbb{R}^{n_{\xi}}}$. In this setting, the states evolves according to the function $\smash{g:\mathcal{X}\!\times\!\mathcal{U}\!\times\!\Xi\rightarrow\mathcal{X}}$ as $\smash{x_{t+1}=g\left(x_{t},u_{t},\xi_{t}\right)}$, incurring the stage cost $\gamma^{t}\,l\left(x_{t},u_{t}\right)$ at each time step, where $\smash{\gamma\!\in\!\left[0,1\right)}$ is the discount factor. By $\Pi$ we denote the set of all feasible deterministic Markov policies, defined as $\smash{\left\{{\pi(\cdot)}\middle|\pi(x)\!\in\!\mathcal{U},\,\forall x\!\in\!\mathcal{X}\right\}}$. The goal is to find a policy $\smash{x_{t}\mapsto\pi\left(x_{t}\right)}$ that minimizes the cumulative cost over an infinite horizon, with initial condition $x\!\in\!\mathcal{X}$,

[TABLE]

The function $\smash{V^{\ast}:\mathcal{X}\rightarrow\mathbb{R}}$ is the value function that represents the optimal cost-to-go from any state of the system if the optimal control policy is played.

To ensure that the problem is well posed we work in that same setting as [5, §6.3], specifically under [5, Assumption 4.2.1(a)] that the stage cost is lower semi-continuous, non-negative, and inf-compact, and also under [5, Assumptions 4.2.1(b), 4.2.2]. The assumptions ensure that from the class of time-varying stochastic policies, the minimum is attained by a stationary deterministic policy, see [5, Theorem 4.2.3]. Finally, $\mathcal{F}\!\left(\mathcal{X}\!\times\!\mathcal{U}\right)$ and $\mathcal{F}\!\left(\mathcal{X}\right)$ are defined as the vector spaces of bounded, real-valued, Borel-measurable functions on $\mathcal{X}\!\times\!\mathcal{U}$ and $\mathcal{X}$ respectively, where [5, Definition 6.3.2, 6.3.4] provides the definitions of boundedness.

II-B Linear Programming (LP) Formulation of DP

Solving the stochastic optimal control problem is equivalent to finding $V^{\ast}$ as the solution of the Bellman equation [3],

[TABLE]

$\mathcal{T}$ is known as the Bellman operator, and the $\mathcal{T}_{u}$ operator represents the cost of making decision $u$ now and then playing optimally from the next time step forward. The optimal policy can be defined using $V^{\ast}$ by,

[TABLE]

The existence of a $V^{\ast}$ and $\pi^{\ast}$ that are Borel measurable and attain the infimum is ensured by [5, Assumptions 4.2.1(a), 4.2.1(b), 4.2.2]. If $\nu(\cdot)$ is a finite measure on $\mathcal{X}$ that assigns positive mass to all open subsets of $\mathcal{X}$, then it can be shown that the solutions of the following linear program,

[TABLE]

satisfy (2) for $\nu$-almost all ($\nu$-a.a.) $\smash{x\!\in\!\mathcal{X}}$, see [5, §6.3]. Constraint (4c) is referred to as the Bellman Inequality. A key feature of the LP formulation is that any choice of $\nu(\cdot)$ that places mass over the whole state space $\mathcal{X}$ leads (4) to recover a solution of the stochastic optimal control problem.

II-C Point-wise Maximum Formulation of DP

Following [21], we introduce additional decision variables and use a point-wise maximum of value functions in the objective and the Bellman inequality constraint,

[TABLE]

where $\smash{J\!\in\!\mathbb{N}}$ specifies the number of value function decision variables. We refer to problem (5) as the point-wise maximum formulation, and the key difference from [21] is the use of the point-wise maximum $V_{\mathrm{PWM}}$ in the objective (5a). The following lemma establish some important properties of (5).

Lemma II.1

Problems (4) and (5) are equivalent in the sense that there exist mappings between the feasible solutions and the optimal solutions of the two problems. Moreover, objective (5a) is jointly convex in the decision variables $V_{j}$, $\smash{j\!=\!1,\dots,J}$.

Proof:

Under the assumptions and definitions of Section II-A, one can easily see there is a mapping between feasible solutions since the space $\mathcal{F}\!\left(\mathcal{X}\right)$ is closed under the maximum operation, i.e., $V_{j}\!=\!V$ for all $x\!\in\!\mathcal{X}$, $\smash{j\!=\!1,\dots,J}$ in one direction, and $V\!=\!V_{\mathrm{PWM}}$ for all $x\!\in\!\mathcal{X}$ in the other direction. This gives equivalent objective value by construction, and thus $\int V^{\ast}d\nu$ is the optimal value for both (4) and (5).

The function $V_{\mathrm{PWM}}$ is convex in $V_{j}$, $\smash{j\!=\!1,\dots,J}$ by definition of the $\max$ function over a finite number of elements. Thus (5a) is convex as integration is a linear operation, see for example [26, Lemma 2.1]. ∎

The reason for the point-wise maximum formulation becomes apparent when we consider approximating the solution to (4) and (5) by restricting the space of the decision variables to subspaces of $\mathcal{F}\!\left(\mathcal{X}\right)$. To gain some initial insight, observe that the Bellman inequality constraint (5c) implies that feasible decisions for (5) will be point-wise under-estimators of $V^{\ast}$. Thus a point-wise maximum is a natural way to combine a family of feasible but sub-optimal decisions.

Computing a solution of problem (5) poses the following difficulties

1. (D1)

$\mathcal{F}\!\left(\mathcal{X}\right)$ is an infinite dimensional space; 2. (D2)

Objective (5a) involves a multidimensional integral over $\mathcal{X}$; 3. (D3)

The $\mathcal{T}_{u}$-operator involves a multidimensional dimensional integral over $\Xi$; 4. (D4)

Constraint (5c) involves an infinite number of constraints; 5. (D5)

Constraint (5c) is non-convex in the decision variables; 6. (D6)

The objective (5a) involves the maximization of a convex function;

Difficulties (D1-D4) apply also to problem (4) and a variety of approaches have been proposed to address them, see for example [27, 28, 29, 30, 31]. In Section III we take inspiration from previous approaches to propose an approximation algorithm that additionally overcomes difficulties (D5-D6).

III Approximate Dynamic Programming (ADP)

This section proposes an algorithm for computing an approximate value function that is feasible for problem (5) at every iteration and analyzes its convergence.

III-A Approaches Adopted for Difficulties (D1), (D3) and (D4)

To overcome difficulty (D1), as suggested in [14], we restrict the value functions candidates to the span of a finite family of Borel-measurable basis functions $\smash{\phi_{k}\!:\!\mathcal{X}\!\rightarrow\!\mathbb{R}}$, $\smash{k\!=\!1,\dots,K}$. We parameterize the restricted function space as,

[TABLE]

The benefit of this parameterization is that it is linear in the $\alpha$ parameter. Each approximate value function $\hat{V}_{j}$ is parameterized by its own vector that we denote $\alpha_{j}$, i.e., $\smash{\hat{V}_{j}(x)\!=\!\alpha_{j}^{\intercal}\,\phi(x)}$ for all $x\!\in\!\mathcal{X}$. For the numerical examples in Section IV we use the space of polynomial functions up to a certain degree by choosing the $\phi_{k}$ to be each of the monomials up to that degree.

To overcome difficulty (D3) we first use Jensen’s inequality to switch the order of expectation and maximisation in the $\mathcal{T}_{u}V_{\mathrm{PWM}}$ term, thus providing a sufficient condition for constraint (5c). We then require that for each basis function $\smash{\mathbb{E}\left[\phi\left(g(x,u,\xi)\right)\right]}$ has an analytic expression. In the case of polynomial basis functions and polynomial dynamics, this requires knowledge of the moments of the distribution of $\xi$ up to the maximum degree of $\xi$ in $\phi\left(g(x,u,\xi)\right)$. If the required moments are not analytically available, then the Monte Carlo sampling can be used to approximate them, and, as the distribution is stationary, this only needs to be computed once.

To overcome difficulty (D4), a variety of convex sufficient conditions techniques are proposed in the literature for approximating (5c) with a finite number of constraints, for example [17, 20, 32, 33]. The applicable reformulation depends on the problem data and basis functions, and the algorithm we propose in the sequel applies for all such convex inner approximations. For example, when all problem data is polynomial and polynomial basis functions are used, then constraint (5c) can be inner approximated using the sum-of-squares (SOS) S-procedure [20].

III-B Proposed Approach for Difficulties (D2), (D5) and (D6)

The inclusion of the point-wise maximum value function in the objective (5a) is pivotal in the algorithm we propose, however, it precludes the use of previous approaches for evaluating the integral in the objective. To overcome difficulty (D2) we replace $\nu$ by a finitely supported measure denoted $c$. Specifically, we choose $c$ as a finite sum of $N_{c}$ Dirac pulses located at $\smash{\{x_{c,i}\}_{i=1}^{N_{c}}\subset\mathcal{X}}$. This violates the hypothesis for equivalence between (2) and (4), but reduces the multidimensional integral in (5a) to a sum over the locations of the Dirac pulses.

For clarity of presentation, we consider now an auxiliary problem that highlights our proposed approach for overcoming difficulties (D5) and (D6). We consider two families of functions defined by two finite sets $\smash{\mathcal{A}_{\mathrm{obj}},\mathcal{A}_{\mathrm{con}}\!\subset\!\mathbb{R}^{K}}$. The first family, parameterized by $\smash{\bar{\alpha}_{\mathrm{obj}}\!\in\!\mathcal{A}_{\mathrm{obj}}}$ is used in a point-wise maximum objective, while the second, parameterized by $\smash{\bar{\alpha}_{\mathrm{con}}\!\in\!\mathcal{A}_{\mathrm{con}}}$, is used in a point-wise maximum constraint. We then define,

[TABLE]

and assume that $\mathcal{A}_{\mathrm{obj}}$ and $\mathcal{A}_{\mathrm{con}}$ have been selected such that $\smash{\bar{V}_{\mathrm{obj}}\!\leq\!V^{\ast}}$ and $\smash{\bar{V}_{\mathrm{con}}\!\leq\!V^{\ast}}$ for all $\smash{x\!\in\!\mathcal{X}}$. Additionally, we introduce $\smash{f_{\mathrm{pwm}}:\mathbb{R}^{K}\rightarrow\mathbb{R}}$ as the point-wise maximum objective function when adding an additional function $\smash{\alpha^{\intercal}\phi(x)\!\in\!\hat{\mathcal{F}}\!\left(\mathcal{X}\right)}$ to the function $\bar{V}_{\mathrm{obj}}$, i.e.,

[TABLE]

where $x_{c,i}$ are the points selected to overcome difficulty (D2). The auxiliary problem for maximizing $f_{\mathrm{pwm}}$ in the presence of $\bar{V}_{\mathrm{obj}}$ and $\bar{V}_{\mathrm{con}}$ is,

[TABLE]

With $\bar{V}_{\mathrm{con}}$ as fixed parameters in (9) the constraint is convex in the decision variable $\alpha$, thus overcoming difficulty (D5). Moreover, if $\bar{V}_{\mathrm{con}}$ satisfies the Bellman inequality, then (9b) implies that $\max\left\{\alpha^{\intercal}\phi,\bar{V}_{\mathrm{con}}\right\}$ also satisfies the Bellman inequality. The steps to show convexity of constraint (9b), first presented in [21], are provided in Appendix A for completeness. In Section III-D we present the proposed algorithm for iteratively adding elements to $\mathcal{A}_{\mathrm{con}}$ in a greedy fashion, where the difference compared to [21] is the choice of objective weighting. To simplify the presentation, we introduce the notation,

[TABLE]

to represent the convexified Bellman inequality constraint (9b).

Despite the convexified constraint, (9) is still a non-convex problem due to (D6). In general, problem (9) will have multiple distinct local maxima and stationary points. The convexity of the objective means that given element of the sub-differential, constructed at a particular point in the decision variable space $\alpha$, it parameterizes a hyperplane that is a point-wise lower-bound on the objective function. Thus we propose to iteratively maximize along sub-gradient directions to overcome difficulty (D6), and in Section III-C we introduce the algorithm and its convergence properties.

III-C First-order method for the point-wise maximum objective

We propose Algorithm 1 to improve $f_{\mathrm{pwm}}$ from a given feasible initial condition $\alpha^{(0)}$ using only first-order information of the objective function. The objective $f_{\mathrm{pwm}}$ is in general non-smooth as it is a maximum of functions, thus in line 4 we use the upper sub-differential for selecting gradient ascent directions, denoted as $\partial^{+}f_{\mathrm{pwm}}$ and defined in Appendix B-A; an upper sub-differential is considered because (9) is a maximization problem. Given a non-zero element from the upper sub-differential, i.e., an upper sub-gradient, in line 8 we update the decision variable by maximizing along the sub-gradient direction within the feasible region. The algorithm terminates when the change in objective value between two subsequent iterations is less than a pre-specified tolerance.

To compute an element from the upper sub-differential of $f_{\mathrm{pwm}}$, to be used in line 4 of Algorithm 1, we introduce the following assumption on the basis functions.

Assumption III.1

The basis functions in the set $\smash{\left\{\phi_{k}\right\}_{k=1}^{K}}$ are continuous for all $x\!\in\!\mathcal{X}$ and include the constant function.

Without loss of generality we take $\smash{\phi_{1}(x)\,=\,1}$ for all $x\!\in\!\mathcal{X}$. For a general choice of basis functions $\phi_{k}$ it is difficult to characterize the upper sub-differential set at non-smooth points. Instead, we work with a particular element from the upper sub-differential of $f_{\mathrm{pwm}}$ that is readily computable at $\alpha$ under Assumption III.1,

[TABLE]

where the term inside the sum is an element from the upper sub-differential of $\smash{\max\left\{\alpha^{\intercal}\phi(x),\bar{V}_{\mathrm{obj}}(x)\right\}}$. Given the element $\smash{\bar{d}\in\partial^{+}f_{\mathrm{pwm}}}\left(\bar{\alpha}\right)$ computed as per (10) at a point $\smash{\bar{\alpha}\!\in\!\mathbb{R}^{K}}$, the hyperplane $\smash{\left(\alpha-\bar{\alpha}\right)^{\intercal}\bar{d}+f\left(\bar{\alpha}\right)}$ is a supporting hyperplane of the convex function $f_{\mathrm{pwm}}$. However, we note that for maximization of a convex function not all supporting hyperplanes are in the upper sub-differential. In the proof of Theorem III.3 we show that (10) is indeed an element of the upper sub-differential.

Algorithm 1 can be seen as a method that iteratively adjusts the objective of line 8 along sub-gradient directions of problem (9). To ensure that an element from the argmax can always be computed in line 8, we introduce the following assumption on the choice of basis functions and inner approximation set.

Assumption III.2

The basis function set $\smash{\left\{\phi_{k}\right\}_{k=1}^{K}}$, Bellman inequality inner approximation set $BI(\mathcal{A}_{\mathrm{con}})$, and problem data are such that the following optimization problem,

[TABLE]

attains its maximum for all $\smash{\{x_{c,i}\}_{i=1}^{N_{c}}}$.

This assumption is not overly restrictive as it can be ensured for any general problem instance by placing an upper bound on a norm of $\alpha$, see [34] for example. For a particular problem instance the assumption can be verified by, for example, showing the existence of a strictly feasible point in the dual of (11) [35, Theorem 3.1] [36, Corollary 30.5.2].

Theorem III.3

Under Assumptions III.1 and III.2, for any initial condition $\smash{\alpha^{(0)}\!\in\!BI(\mathcal{A}_{\mathrm{con}})}$ and any $\smash{\epsilon\!>\!0}$, Algorithm 1 generates a non-decreasing sequence $\smash{f_{\mathrm{pwm}}\left(\alpha^{(k)}\right)}$ and terminates after a finite number of iterations. With $\smash{\epsilon\!=\!0}$, the sequence $\smash{f_{\mathrm{pwm}}\left(\alpha^{(k)}\right)}$, converges to a finite value, and the sequences $\smash{\alpha^{(k)}}$, $\smash{d^{(k)}}$, satisfy the following necessary optimality condition for $\smash{\max\nolimits_{\alpha\in\mathbb{R}^{K}}\left\{f_{\mathrm{pwm}}(\alpha)\hskip 0.28436pt;\hskip 1.42271pt\alpha\!\in\!BI(\mathcal{A}_{\mathrm{con}})\right\}}$ in the limit,

[TABLE]

Proof:

see Appendix B-D. ∎

Theorem III.3 guarantees that if the initial condition $\alpha^{(0)}$ strictly improves on $f_{\mathrm{pwm}}(0)$, then $f_{\mathrm{pwm}}(\alpha^{(i)})$ returned also strictly improves on $f_{\mathrm{pwm}}(0)$. The convergence in finite iterations ensures that the algorithm is practical to implement, and the limiting behaviour suggests that in the best case the $\alpha^{(i)}$ returned could be close to a local maxima. Although Theorem III.3 provides no insight into the rate of convergence, the numerical examples in Section IV demonstrate that significant improvement in $f_{\mathrm{pwm}}$ can be achieved with only a handful of iterations.

Algorithm 1 is a so-called Minorize Maximize algorithm for maximizing a convex function, and we now contrast with generic algorithms that exist in the literature for this same purpose. In the case where $f_{\mathrm{pwm}}$ is differentiable, then line 4 of Algorithm 1 becomes $\smash{d^{(k)}\leftarrow\nabla f_{\mathrm{pwm}}(\alpha^{(k)})}$ and is a special-case of the so-called convex-concave procedure introduced in [37], and for which convergence guarantees are given in [38, Theorem 4]. Algorithms applicable for non-smooth problems like (9) are presented together with convergence guarantees in [39, Theorem 3] and [40, Proposition 1]. Applying the algorithm from [39] or [40] to problem (9) would require using the lower sub-differential of $f_{\mathrm{pwm}}$ in line 4 of Algorithm 1. For a non-smooth convex function the lower sub-differential contains the upper sub-differential, thus allowing more flexibility on line 4 of Algorithm 1. However, [39] and [40] use the lower sub-differential also for defining necessary optimality conditions. This means that, compared to Algorithm 1, the algorithms from [39] and [40] may have additional points in their convergence set that are not local maxima of the non-smooth convex maximization problem.

III-D Point-wise Maximum ADP Algorithm

In this section we propose Algorithm 2, which iteratively updates the value function estimates used in the objective and constraints of problem (9), i.e., $\bar{V}_{\mathrm{obj}}$ and $\bar{V}_{\mathrm{con}}$. At each iteration of lines 7–12, a candidate approximate value function $\alpha^{(0)}$ is generated by solving (11) with $x_{c,i}$ as one of the Dirac pulse locations from $c(\cdot)$. Algorithm 1 refines this candidate before it is added to the collections $\mathcal{A}_{\mathrm{obj}}$ and $\mathcal{A}_{\mathrm{con}}$. This process of generating, refining, and adding is repeated for all $x_{c,i}$, $\smash{i=1,\dots,N_{c}}$. The algorithm terminates when the improvement in $\bar{V}_{\mathrm{obj}}$ is below some pre-specified threshold. The following theorem formalises the convergence properties of Algorithm 2.

Theorem III.4

For any sets $\mathcal{A}_{\mathrm{obj}}$ and $\mathcal{A}_{\mathrm{con}}$ such that $\bar{V}_{\mathrm{obj}}$ and $\bar{V}_{\mathrm{con}}$ are point-wise under-estimators of $V^{\ast}$, and for any $\smash{\epsilon_{\mathrm{IN}},\,\epsilon_{\mathrm{OUT}}>0}$, Algorithm 2 terminates after a finite number of iterations.

Proof:

By Theorem III.3 we have that line 9 of Algorithm 2 terminates after finite iterations for all $\smash{\epsilon_{\mathrm{IN}}\!>\!0}$. The sequence $f^{(m)}$ in non-decreasing by definition as a point-wise maximum of functions and because elements are never removed from the set $\mathcal{A}_{\mathrm{obj}}$. The same reasoning as Appendix B-D establishes that $\smash{\max\nolimits_{\left\{\bar{\alpha}_{\mathrm{obj}}\in\mathcal{A}_{\mathrm{obj}}\right\}}\hskip 2.84544pt\bar{\alpha}_{\mathrm{obj}}^{\intercal}\phi(x)}$ is bounded above for all $\smash{\{x_{c,i}\}_{i=1}^{N_{c}}}$ at all iterations of Algorithm 2. Hence $f^{(m)}$ is bounded above and is thus a convergent sequence. Therefore, for all $\smash{\epsilon_{\mathrm{OUT}}\!>\!0}$ there must exist an $\smash{m\!\geq\!1}$ such that the condition on line 15 triggers. ∎

The convergence of Algorithm 2 is guaranteed even without the refinement steps of Algorithm 1. However, our numerical results in Section IV show that without refinement convergence tends to be much slower, and that significant improvements are achieved with only a few iterations of Algorithm 1. Hence, the gradient-based motivation and theoretical guarantees of Theorem III.3 suggests that performing the refinement steps of Algorithm 1 is beneficial.

The objective $\alpha^{\intercal}\phi(x_{c,i})$ in line 8 of Algorithm 2 is chosen so that the $\alpha^{(0)}$ passed to Algorithm 1 has a non-zero sub-gradient $d^{(0)}$ (line 4 of Algorithm 1). To see this, note that the sub-gradient in (10) is non-zero if $\alpha^{\intercal}\phi(x_{c,i})$ weakly dominates $\bar{V}_{\mathrm{obj}}(x_{c,i})$ for at least one $\smash{i=1,\dots,N_{c}}$. Thus, by Assumption III.2, line 8 of Algorithm 2 computes an $\alpha^{(0)}$ that weakly dominates $\bar{V}_{\mathrm{obj}}$ at the chosen point $x_{c,i}$ if such a solution exists in the feasible set $\alpha\!\in\!BI(\mathcal{A}_{\mathrm{con}})$. Different objectives for line 8 of Algorithm 2 can be considered and still enjoy the convergence guarantee of Theorem III.4. However, this would introduce a tuning parameter and empirical testing has shown no benefit when hand-tuning the objective.

III-E Discussion and extensions

Considering the motivating problem (5) with objective $\int V_{\mathrm{PWM}}d\nu$, the obvious choice for $c$ is to draw samples from $\nu$, and to choose $\nu$ as the initial state distribution. However, sampling $c$ in different ways may improve the objective $\int V_{\mathrm{PWM}}d\nu$, and for this reason $c$ is commonly referred to as the state relevance weighting [16]. The sub-gradient computed by (10) effectively sub-samples the points from $c$ where the current approximate value function dominates the fixed $\bar{V}_{\mathrm{obj}}$. Thus Algorithm 1 can be seen a method for automatically choosing the state relevance weighting parameter to maximize $f_{\mathrm{pwm}}$, i.e., the surrogate for $\int V_{\mathrm{PWM}}d\nu$.

If the goal is to optimize the on-line performance of the greedy policy, it is again likely that difference choices of $\nu$, and hence different samples for $c$, lead to differing on-line performance. Motivated by the performance bounds provided in [16], a reasonable choice is to place Algorithm 2 inside another iteration that updates $\nu$ as the discounted occupancy measure for the current greedy policy, computed empirically by simulating the system evolution using Monte Carlo sampling.

For real-time applications where the greedy policy must be computed very fast, it is necessary that the cardinality of $\mathcal{A}_{\mathrm{obj}}$ is small, and perhaps even a singleton. For examples with linear dynamics, quadratic stage costs, polytopic spaces, and using the space of quadratics for $\hat{\mathcal{F}}\!\left(\mathcal{X}\right)$, then the greedy policy is a Quadratically Constrained Quadratic Program (QCQP), with the number of quadratic constraints equal to the cardinality of $\mathcal{A}_{\mathrm{obj}}$. In such examples, a low cardinality of $\mathcal{A}_{\mathrm{obj}}$ has clear benefits from an on-line compution perspective. In these cases it is beneficial to run Algorithm 2 twice. First, Algorithm 2 is run for as long as practical to achieve a good under-estimate of $V^{\ast}$, with a simple initialization, for example $\smash{\mathcal{A}_{\mathrm{obj}}\!=\!\mathcal{A}_{\mathrm{con}}\!=\!\{0\}}$. Second, Algorithm 2 is run for as many iterations as the desired cardinality of $\mathcal{A}_{\mathrm{obj}}$, with $\mathcal{A}_{\mathrm{con}}$ initialized as the under-estimate resulting from the first run.

The approximate value function computed by Algorithm 2 can be used off-line to certify the empirical performance of alternative policies that do not use the approximate value function. In this case Algorithm 2 is run for as long as practical, then the chosen policy is simulated from a particular initial state, $\hat{x}$, for a time horizon such that $\gamma^{t}$ has decayed sufficiently. The approximate value function evaluated at the initial state is a lower bound on $V^{\ast}(\hat{x})$ and thus provides a bound on the sub-optimality of the policy, and hence indicates the potential benefit of considering further alternatives.

In [21] the value function decision variable was also included in the right-hand-side of constraint (9b), i.e.,

[TABLE]

for all $x\!\in\!\mathcal{X}$ and $u\!\in\!\mathcal{U}$. This results in a bi-linear term in the constraint, and in that work the authors suggest gridding the multiplier of the bi-linear term. We do not consider this extension in the numerical examples because it adds significant computation time and empirically it provides little or no benefit for the examples considered.

Algorithm 1 can be extended to fit multiple new lower bounding functions at the same time. To exemplify, consider the case of adding two new lower bounding functions. The non-convex optimization problem then becomes,

[TABLE]

We construct an element from the upper subdifferential in a similar fashion,

[TABLE]

As the constraints are separable we see that once the subdifferential element is computed, then line 8 of Algorithm 1 can be solved in parallel for $\alpha$ and $\beta$, differing only in the objective vector.

IV Numerical Examples

In the following numerical example we consider problems with linear dynamics, convex quadratic stage costs, hyper-cube constraints on the input space, and use convex quadratics for the restricted function space $\hat{\mathcal{F}}\!\left(\mathcal{X}\right)$. In Appendix C we provide the definition of $\hat{\mathcal{F}}\!\left(\mathcal{X}\right)$, formulate line 8 of Algorithm 1 so that it can be passed to a standard solver, and verify that Assumption III.1 holds. Assumption III.2 was observed to hold empirically, in that the solver returned a finite, optimal solution at each iteration.

IV-A Algorithm insight on 1-dimensional example

To provide visual insight into how Algorithms 1 and 2 adapt the linear combination of basis functions, we use the simple 1-dimension example from [17], i.e., with $\smash{n_{x}\!=\!n_{u}\!=\!1}$. The dynamics, costs, constraints, and initial state distribution are given by,

[TABLE]

and we use the space of univariate quadratics as $\hat{\mathcal{F}}\!\left(\mathcal{X}\right)$. We initialise $\bar{V}_{\mathrm{obj}}$ and $\bar{V}_{\mathrm{con}}$ with the solution of the LP approach using a single Bellman inequality constraint, i.e.,

[TABLE]

with $\nu$ chosen as the initial state distribution. The solution of (13) represents the lower bound proposed in [16]. In this setting (13) is a convex semi-definite optimization program, with the objective requiring the first and second moments of $\nu$ and the constraint reformulated as a linear matrix inequality with respect to the quadratic coefficient decision variables, see [17, §6].

Choosing $c$ as $\smash{N_{c}\!=\!10^{6}}$ samples from the initial state distribution, Figure 1 shows the first iteration of Algorithm 2, with sub-figure (b) as zoomed view of sub-figure (a). Sub-figure (a) shows that the candidate approximate value function generated by line 8 of Algorithm 2 (blue) is a significantly worse lower bound compared to that resulting from four refinement iterations of Algorithm 1 (red), when integrated with respect to the $N_{c}$ samples (grey histogram).

Figure 1(b) shows that at the $x_{c,i}$ used on line 8 of Algorithm 2 (blue dot), the generated value function (blue line) strictly improves on $\bar{V}_{\mathrm{obj}}$ (dotted black line). All of the refinement iterations (red) trade-off a decrease at this $x_{c,i}$ for a significant increase at the other samples. For this example and at this iteration, the significant increase in the point-wise maximum objective value $f_{\mathrm{pwm}}$ is gained in regions away from the origin.

Figure 2 (top left) shows the lower bound and online performance integrated with respect to the $\smash{N_{c}\!=\!10^{6}}$ samples from the initial state distribution. This shows that 1000 iterations of Algorithm 2 combined with the refinement steps of Algorithm 1 (red line) allow the clipped-LQR controller to be certified as within 1.5% of the optimal. Without the refinement steps (blue line) the sub-optimality bound is 4.5%, even after $10^{4}$ iterations. Figure 2 (bottom left) shows the number of iterations of Algorithm 1 performed at each iteration of Algorithm 2. Together with Figure 2 (top left), this shows that the handful of Algorithm 1 iterations improves the bound (red line) with an order of magnitude fewer iterations than Algorithm 2 without refinement (blue line). In fact, computations performed for $10^{5}$ iterations empirically suggest that Algorithm 2 without refinement will not give better than a 4.5% sub-optimality bound within a practical number of iterations, for this example.

IV-B Higher Dimensional Linear-Quadratic Problems

We consider again the an input constrained linear-quadratic system, this time with dimension $n_{x}\!=\!10$, and $n_{u}\!=\!3$. The system dynamics take the form,

[TABLE]

where $A$ and $B_{u}$ are matrices of compatible size, and the quadratic stage cost is, $\smash{l(x,u)\!=\!x^{\intercal}I_{n_{x}}x+u^{\intercal}I_{n_{u}}u}$, where $I_{n}$ denotes an identity matrix of size $n$, and we use discount factor $\smash{\gamma\!=\!0.99}$. The initial state is normally distributed as $x_{0}\!\sim\!\mathcal{N}\left(0,\Sigma_{\nu}\right)$ with $\Sigma_{\nu}\!=\!9I_{n_{x}}$. The $A$ and $B_{u}$ matrices are randomly generated, with the $A$ matrix scaled to be marginally stable, i.e., a spectral radius equal to 1.

Figure 2 (top right) shows the lower bound achieved by running Algorithm 2 without (blue) and with (red) the refinement iterations of Algorithm 1, with the number of refinement iterations shown in Figure 2 (bottom right). To demonstrate the benefit of our proposed Algorithm relative to previous work, the sets $\mathcal{A}_{\mathrm{obj}}$ and $\mathcal{A}_{\mathrm{con}}$ are initialised with the solution of (13), with $\nu$ chosen as the initial distribution, i.e., $\mathcal{N}\left(0,\Sigma_{\nu}\right)$. Thus, the value for iteration 1 of the blue and red lines in Figure 2 (top right), approximately 291, is the lower bound achieved by the method proposed in [16]. The key feature of the result is that, although Algorithm 2 without refinement is guaranteed to converge, the number of iterations required to reach a reasonable lower bound is significant. Algorithm 2 with refinement, on the other hand, achieves a significantly better lower bound with orders of magnitude fewer iterations.

The improvement in the lower bound achieved with the refinement steps of Algorithm 1 is only meaningful if it significantly tightens the online performance bound for a particular policy. For this example, an MPC policy with a time horizon of $10$ achieves an online performance of $870$, shown by the dotted black line on Figure 2 (top right). Thus the refinement steps of Algorithm 1 (red) certify this policy to be within 11% of the optimal, while the bound without the refinement steps (blue) provides only a 70% sub-optimality certificate.

As another point of comparison, we show in Figure 3 (top) the lower bound achieved by two alternative methods (green), and then we use the output of these methods to initialize our proposed algorithm (blue and red lines). For the iterated Bellman inequality method proposed in [19] (solid green), we use $100$ Bellman inequality iterations, achieving a lower bound of $656$, and we note that more iterations did not improve the bound. This method computes $100$ approximate value functions from the solution of one optimization problem, hence the computation time in Figure 3 (bottom) is constant for the first $100$ iterations, while the lower bound increases because we add the $100$ approximate value functions sequentially to the point-wise maximum approximation. For implementing the hand-tuning method suggested in [21] (dotted green) we run Algorithm 2 without refinement and with the modification that on line 8 we manually select a different objective at each iteration. We performed this for a range of options and show in Figure 3 the option achieving the tightest bound. For both of the comparisons, we then use the approximate value functions computed to initialize the sets $\mathcal{A}_{\mathrm{obj}}$ and $\mathcal{A}_{\mathrm{con}}$ and again run Algorithm 2 without (blue) and with (red) the refinement steps of Algorithm 1. In both cases, only a handful of iterations of Algorithm 1 are required to achieve the improved lower bound of $783$, an $11\%$ and $19\%$ improvement respectively.

The combination of results in Figure 2 (top right) and Figure 3 (top) suggests that the lower bound achieved by Algorithm 2 with refinement (red lines) is not overly sensitive to the initialization of $\mathcal{A}_{\mathrm{obj}}$ and $\mathcal{A}_{\mathrm{con}}$, provided that a reasonable number of iterations are performed. The results also suggest that Algorithm 2 without refinement (blue lines) requires an impractically large number of iterations to achieve comparatively small improvements of the lower bound. Figure 3 (bottom) shows the cumulative computation time required for computing the approximate value functions, i.e., solving line 8 of Algorithm 1 and line 8 of Algorithm 2. This indicates that the hand-tuning method suggested in [21] can achieve a tighter lower bound with less computation time compared to the method suggested in [19]. However, both methods have a limit beyond which extra computation did not improve the lower bound. We note that the computation times for the results in Figure 2 fall between the lines shown in Figure 3 (bottom), and are not shown for the sake of clarity. The results on this higher dimensional example empirically support our claim that sub-gradients of the point-wise maximum objective function are effective for computing tighter lower bounds.

V Conclusion

We proposed an algorithm that computes a family of lower bounding approximate value functions in an iterative fashion, with the choice of initial state distribution as the only parameter to be selected by the designer. We motivate our algorithm by considering the non-convex objective of maximizing the point-wise maximum of lower bounding value functions, and use sub-gradient information to find (potentially) sub-optimal solutions. Testing our algorithm on linear-quadratic examples, we demonstrated a significant tightening of the lower bound compared to existing methods, achieved with a modest or negligible increase in the computation time.

As future work, we will investigate adaptations of the proposed algorithm that are tailored to policy performance. This is a more challenging setting because the computation restriction are more stringent for evaluation of a policy. Moreover, counter-examples can readily be constructed where an approximate value function that provides a relatively tight lower bound leads to a greedy policy with relatively poor online performance.

A weakness of the proposed method is that by sampling from the initial state distribution we forgo the direct theoretical connection to the Bellman equation. An interesting direction for extending our algorithm is to consider methods that directly maximize the integral with respect to the initial state distribution, for example stochastic gradient ascent. This would also address the open question about whether the lower bound quality is sensitive to the choice of samples in our proposed algorithm.

Appendix A Reformulation of Point-wise Maximum Inequality

This appendix summarises our previous work [21] in the context of this paper.

A-A Jensen’s inequality and epigraph reformulation

The point-wise maximum constraint (5c) is equivalent to $J$ separate constraints of the form,

[TABLE]

As $\max(\cdot)$ is a convex function, by Jensen’s inequality a sufficient condition for constraint (5c) is,

[TABLE]

An exact epigraph reformulation can now be applied [41, Theorem 1] with the epigraph variable denoted $s_{V}$, i.e.,

[TABLE]

for $j=1,\dots,J$, where the set $\mathcal{S}$ is defined as,

[TABLE]

We choose to square the epigraph variable $s_{V}$ without loss of generality because [5, Assumptions 4.2.1(b)] implies that $V^{\ast}$ is non-negative.

A-B S-procedure reformulation

The S-procedure [42] is used to obtain a sufficient condition for (15). Applying the S-procedure to the relevant part of $\mathcal{S}$ leads to,

[TABLE]

with $\lambda_{k}\!\in\!\mathbb{R}_{+}$ as the non-negative decision variables introduced by the S-procedure. Reformulation (16) still suffers from difficulty (D5): there will be $J$ bilinear terms of the form $\smash{\lambda_{k}\,\mathbb{E}\left[V_{k}\left(g(x,u,\xi)\right)\right]}$ in each of the $J$ constraints.

The following implications summarize the approximation steps described,

[TABLE]

In words, this reformulation is sufficient in the sense that if a family of functions $\smash{V_{1},\dots,V_{J}}$ satisfies (16) then it also satisfies (5c) (but not necessarily the other way around). An equivalent or tighter approximation can be found by allowing the S-procedure multipliers to depend on the state and input, i.e., $\smash{\lambda_{k}:\mathcal{X}\!\times\!\mathcal{U}\rightarrow\mathbb{R}_{+}}$. This would require the introduction of a restricted function space on $\smash{\left(\mathcal{X}\!\times\!\mathcal{U}\right)}$, denoted $\hat{\mathcal{F}}\!\left(\mathcal{X}\!\times\!\mathcal{U}\right)$, and defined similar to $\hat{\mathcal{F}}\!\left(\mathcal{X}\right)$ in (6).

A-C Overcoming difficulty (D5) for (9)

The auxiliary problem (9) introduced in Section III-B has a form similar to (14) except that the only constraint included is the one with the decision variable on the left side of the inequality. Letting $V_{1}$ in (14) correspond to $\alpha^{\intercal}\phi$ in (9) and applying reformulation (16) we get the following sufficient condition for (9b),

[TABLE]

for all $x\!\in\!\mathcal{X},\,u\!\in\!\mathcal{U},\,s_{V}\!\in\!\mathbb{R}$, where the multipliers $\lambda_{\bar{\alpha}_{\mathrm{con}}}$ are additional non-negative decision variable for each element of $\mathcal{A}_{\mathrm{con}}$. As the $\bar{\alpha}_{\mathrm{con}}$ are fixed parameters in problem (9), it is clear that this reformulation is linear in the decision variables $\alpha$ and $\lambda_{\bar{\alpha}_{\mathrm{con}}}$. When the problem data and basis functions are polynomial, the infinite constraints are reformulated in the usual way, see [17, Appendix A] for example, the result is a single Linear Matrix Inequality (LMI) constraint.

Appendix B Properties of the Inner Problem of Section III-C

All the material in this appendix is formulated for a minimization optimization objective, chosen to make the results readily comparable with existing optimization literature. Problem (9) and Algorithm 1 are readily converted to minimization problems by taking the negative of the objective.

B-A Differentiability definitions

We provide for completeness the definitions of the regular and general sub-differential as taken from [43, §7.D, §8.A]. The definition of the sub-differential commonly used for convex optimization problems is special case of the regular sub-differential defined below [43, Proposition 8.12], required here because (9a) is non-convex when cast as a minimization problem. We require additionally the general sub-differential definition because (9a) is non-smooth.

Given a function $\smash{f:\mathbb{R}^{n}\rightarrow\mathbb{R}\cup\{-\infty,\infty\}}$, a vector $\smash{d\!\in\!\mathbb{R}^{n}}$ is a regular lower subgradient of $f$ at the point $\smash{x\!\in\!\mathbb{R}^{n}}$ if the following one-sided limit condition holds,

[TABLE]

The regular lower subdifferential of $f$ at $x$, denoted $\smash{\hat{\partial}f(x)}$ is the set of regular lower subgradients of $f$ at $x$. A vector $\smash{d\!\in\!\mathbb{R}^{n}}$ is a general lower subgradient of $f$ at the point $x$ if there exists sequences $\smash{x^{(i)}\!\stackrel{{\scriptstyle f}}{{\rightarrow}}\!x}$ and $\smash{d^{(i)}\!\rightarrow\!d}$ with $\smash{d^{(i)}\!\in\!\hat{\partial}f(x^{(i)})}$, where the notation $\stackrel{{\scriptstyle f}}{{\rightarrow}}$ stands for $f$-attentive, defined as,

[TABLE]

The general lower subdifferential of $f$ at $x$, denoted $\smash{\partial f(x)}$ is the set of general lower subgradients of $f$ at $x$. At a point $x$ where $f$ is finite, the set $\smash{\partial f(x)}$ and $\smash{\hat{\partial}f(x)}$ are closed, with $\smash{\hat{\partial}f(x)}$ convex and $\smash{\hat{\partial}f(x)\!\subseteq\!\partial f(x)}$. The function $f$ is subdifferentially regular at a point $x$ if $\smash{\hat{\partial}f(x)\!=\!\partial f(x)}$. These definitions and properties correspond to [43, Definition 8.3, Theorem 8.6, Definition 7.25]. Note that if $f$ is differentiable at $x$, then $\smash{\hat{\partial}f(x)\!=\!\{\nabla f(x)\}}$, i.e., a singleton, and if additionally $f$ is smooth on a neighbourhood of $x$, then $\smash{\partial f(x)\!=\!\{\nabla f(x)\}}$ also. For the standard definitions of the gradient $\nabla f(x)$ of a function $f$ at a differentiable point $x$, the reader is referred to [44, §B.5]. The regular and general upper subdifferential are computed as $\smash{-\hat{\partial}\left(-f\right)}$ and $\smash{-\partial\left(-f\right)}$, and denoted ${\hat{\partial}^{+}f}$ and ${\partial^{+}f}$, respectively.

B-B Necessary condition for local optimality

A function $\smash{f:\mathbb{R}^{n}\rightarrow\mathbb{R}}$ is proper, for a minimisation objective, if $\smash{f(x)\!<\!+\infty}$ for at least one $\smash{x\!\in\!\mathbb{R}^{n}}$, and $\smash{f(x)\!>\!-\infty}$ for all $\smash{x\!\in\!\mathbb{R}^{n}}$. Consider the minimization of a proper, lower-semi-continuous function $\smash{f:\mathbb{R}^{n}\rightarrow\mathbb{R}}$ over a closed set $\mathcal{C}\subseteq\mathbb{R}^{n}$, i.e., $\smash{\min_{x\in\mathcal{C}}f(x)}$. As per [43, Theorem 8.15], a necessary condition for the local optimality of a point $\overline{x}\in\mathcal{C}$ is:

[TABLE]

where $N_{\mathcal{C}}$ is the general normal cone of the set $\mathcal{C}$ at the point $\overline{x}$, see [43, Definition 6.3]. If in addition $\mathcal{C}$ is a convex set, then this condition is equivalent to the existence of a $\smash{d\!\in\!\partial f(\overline{x})}$ satisfying

[TABLE]

see [43, Theorem 6.9]. Note further that if $f$ is convex then these conditions are necessary and sufficient for $\overline{x}$ to be globally optimal.

A stationary point of the optimisation problem $\smash{\min_{x\in\mathcal{C}}f(x)}$ is one satisfying $\smash{0\!\in\!\partial\left(f(\overline{x})+\delta_{\mathcal{C}}(\overline{x})\right)}$, where $\delta_{\mathcal{C}}$ is the indicator function of the set $\mathcal{C}$. All stationary points satisfy (17) as,

[TABLE]

If $\mathcal{C}$ is convex, then the inclusion becomes an equality at a point $\overline{x}$ where $f$ is sub-differentially regular, [43, Corollary 10.9].

B-C Proof of convergence for a more general problem statement

To streamline the proof of Theorem III.3, we consider now a more a general problem statement, and then in Appendix B-D below we show that problem (9) and Algorithm 1 has this form. Given a proper, lower-semi-continuous, concave function $\smash{f:\mathbb{R}^{n}\rightarrow\mathbb{R}}$ and a convex constraint set $\mathcal{C}\subseteq\mathbb{R}^{n}$ such that $f$ is bounded below on $\mathcal{C}$, we consider the optimization problem,

[TABLE]

We show that Algorithm 3 finds points that satisfy (18); note that in line 4 we use the lower sub-differential because (19) is a minimization problem.

Theorem B.1

For any initial condition $\smash{x^{(0)}\!\in\!\mathcal{C}}$ and any $\smash{\epsilon\!>\!0}$, Algorithm 3 generates a non-decreasing sequence $\smash{f\left(x^{(k)}\right)}$ and terminates after a finite number of iterations. With $\smash{\epsilon\!=\!0}$, and assuming that the $\arg\min$ on line 8 is always attained, the sequence $\smash{f\left(x^{(k)}\right)}$, converges to a finite value, and the sequences $\smash{x^{(k)}}$, $\smash{d^{(k)}}$, satisfy condition (18) in the following sense,

[TABLE]

Proof:

The sub-differential gives majorizing functions

We first show that given any element of the general lower subdifferential, $\smash{d^{(k)}\!\in\!\partial f\left(x^{(k)}\right)}$, the surrogate function,

[TABLE]

is a point-wise upper-bound of the concave function $f$. This is trivial for a differentiable point $\smash{x^{(k)}}$ as we have that the gradient is the only element of both the general lower and upper subdifferential of $f$ at $\smash{x^{(k)}}$ and hence $\smash{s_{k}\left(x\right)}$ is a global upper-bound of $f$. Pathological functions where the gradient is not an element of general subdifferential at a differentiable point are excluded by virtue of the $f$ being concave.

At a non-differentiable point $\smash{x^{(k)}}$, as $f$ is concave, the regular lower subdifferential is empty at this point. Thus the general lower subdifferential is defined by the limits along all sequences of differentiable points leading to $\smash{x^{(k)}}$. As the regular lower and upper subdifferential are equal at all points along any such sequence we have that,

[TABLE]

Thus is remains to show that $\partial\left(-f\right)(x^{(k)})$ contains only supporting hyperplanes of the hypograph of $f$ at $\smash{x^{(k)}}$. As $-f$ is convex, we have by [43, Proposition 8.12] that the general lower subdifferential of $-f$ is,

[TABLE]

Thus we have shown that the surrogate function $\smash{s_{k}\left(x\right)}$ is a point-wise upper-bound of the concave function $f$ at any point $\smash{x^{(k)}}$ that it is constructed.

Termination in finite iterations

By definition, the minimization problem on line 8 returns $\smash{x^{(k+1)}}$ satisfying the optimality condition,

[TABLE]

Combining the properties of the surrogate function $s_{k}$ with the definition of line 8 as a minimization problem, we have that,

[TABLE]

with $\smash{x^{(k)},x^{(k+1)}\!\in\!\mathcal{C}}$ ensured by the constraints of line 8. The equality is by (20), the first inequality is by definition of the minimization on line 8, and the final inequality is by the fact that the surrogate is a point-wise upper-bound.

By the assumption that $f$ is bounded below on $\mathcal{C}$, the sequences $\smash{f\left(x^{(k)}\right)}$ and $\smash{s_{k}\left(x^{(k)}\right)}$, for $\smash{k\!\geq\!0}$, are convergent, hence Cauchy. Therefore, for all $\smash{\epsilon\!>\!0}$ there must exist a $\smash{k\!\geq\!1}$ such that the condition on line 11 triggers.

Convergence to necessary conditions for optimality

For $\smash{\epsilon\!=\!0}$ we have from the argument above that the sequences $\smash{f\left(x^{(k)}\right)}$ and $\smash{s_{k}\left(x^{(k)}\right)}$ converge to a finite value. To show that the sequence $\smash{x^{(k)}}$ satisfies condition (18) in the limit, we need to show that,

[TABLE]

To show this it is sufficient to show that the sequence $\smash{x^{(k)}}$ converges to an optimal point of $\smash{\min\nolimits_{x\in\mathcal{C}}s_{k}(x)}$, i.e., we show that sequences $\smash{x^{(k)}}$, $\smash{d^{(k)}}$, satisfy,

[TABLE]

The $\min$ here is attained by the assumption in the theorem statement that line 8 of Algorithm 3 attains at every iteration. To show that the limit in (23) exists and equals zero, we first consider for the sake of contradiction that the sequences $\smash{x^{(k)}}$, $\smash{d^{(k)}}$ satisfy,

[TABLE]

By definition of the $\liminf$, for every $\smash{k\!\geq\!0}$ there exists a $\smash{j\!\geq\!k}$ for which,

[TABLE]

By definition of line 8 as a minimization problem we have for this pair $k,j$ that,

[TABLE]

Repeating this argument starting from $j+1$, we readily establish that,

[TABLE]

which contradicts the previous conclusion that the sequence $\smash{s_{k}\left(x^{(k)}\right)}$ converges to a finite value. Moreover we have that,

[TABLE]

because $\smash{x^{(k)}\!\in\!\mathcal{C}}$, for $\smash{k\!\geq\!0}$. Thus, by contradiction we have shown that,

[TABLE]

and hence the limit in (23) exists and equals zero. ∎

Note that if line 5 of Algorithm 3 triggers, then the subgradient is zero and condition (18) is satisfied. In this case the $x^{(k)}$ returned is a global maximizer of the concave function $f$. Note also that for a positive $\epsilon$, if the condition on line 11 triggers with $\smash{f\left(x^{(k)}\right)=f\left(x^{(k-1)}\right)}$, then $x^{(k-1)}$ satisfies (18). To show this, first note that by (20) and (22) we have,

[TABLE]

From this we substitute $\smash{x^{(k)\intercal}d^{(k-1)}=x^{(k-1)\intercal}d^{(k-1)}}$ into the optimality condition (21) that $\smash{x^{(k)}}$ satisfies, and we get that the $\smash{d^{(k-1)}\!\in\!\partial f\left(x^{(k-1)}\right)}$ from line 4 of Algorithm 3 satisfies condition (18) at $\smash{x^{(k-1)}}$.

B-D Proof of convergence for Algorithm 1

Proof:

We show that the objective function $f_{\mathrm{pwm}}$ and the convex constraint $\smash{\alpha\!\in\!BI(\mathcal{A}_{\mathrm{con}})}$ satisfy the assumptions of Theorem B.1. Casting (9) as a minimization problem, the objective is,

[TABLE]

The two elements of the $\max$ are linear in the decision variable $\alpha$, and thus the objective is concave in $\alpha$.

We now show that $-f_{\mathrm{pwm}}$ is bounded below on the constraint set. The assumption that $\bar{V}_{\mathrm{con}}$ is a point-wise lower bound of $V^{\ast}$ means that for any $\alpha$ satisfying constraint (9b), the function $\smash{\alpha^{\intercal}\phi(x)}$ is also a point-wise lower bound of $\smash{V^{\ast}(x)}$ for all $x\!\in\!\mathcal{X}$. This is ensured by the Bellman operator $\mathcal{T}$ being monotone and $\gamma$-contractive. Moreover, under [5, Assumptions 4.2.1(a), 4.2.1(b), 4.2.2] we have that $\smash{V^{\ast}(x)}$ is finite for all $x\!\in\!\mathcal{X}$. As $\bar{V}_{\mathrm{obj}}$ is also a point-wise lower bound of $V^{\ast}$, we have that all elements of the sum in $f_{\mathrm{pwm}}$ are bounded above, and hence $-f_{\mathrm{pwm}}$ is bounded below for all $\smash{\alpha\!\in\!BI(\mathcal{A}_{\mathrm{con}})}$.

Next we show that equation (10) correctly computes an element of the general upper subdifferential of $f_{\mathrm{pwm}}$. For an $\alpha$ where $f_{\mathrm{pwm}}$ is differentiable, we have that $\smash{\alpha^{\intercal}\phi(x_{c,i})\!\neq\!\bar{V}_{\mathrm{obj}}(x_{c,i})}$ for all $\smash{i=1,\dots,N_{c}}$, and thus equation (10) computes the gradient at this point. The objective function $f_{\mathrm{pwm}}$ is non-differentiable for an $\alpha$ where $\smash{\alpha^{\intercal}\phi(x_{c,i})\!=\!\bar{V}_{\mathrm{obj}}(x_{c,i})}$ for at least one point $\smash{i=1,\dots,N_{c}}$. Letting $\mathcal{I}_{<}(\alpha)$, $\mathcal{I}_{=}(\alpha)$, and $\mathcal{I}_{>}(\alpha)$ denote the indices $\smash{i=1,\dots,N_{c}}$ where $\smash{\alpha^{\intercal}\phi(x_{c,i})}$ is respectively less than, equal, and greater than $\smash{\bar{V}_{\mathrm{obj}}(x_{c,i})}$, we define $\delta_{\min}(\alpha)$ as,

[TABLE]

Recall that under Assumption III.1, $\phi_{1}$ is taken to be the constant function and let $e_{1}$ denote a vector with $1$ as the first element and zero otherwise. Thus for all $\smash{\delta\!\in\!(0,\delta_{\min})}$ we have that $f_{\mathrm{pwm}}$ is differentiable at $\smash{\left(\alpha+\delta\,e_{1}\right)}$ with gradient given by equation (10). For any sequence $\delta\rightarrow 0$, the sequence $\smash{\left(\alpha+\delta\,e_{1}\right)}$ is $f_{\mathrm{pwm}}$-attentive, i.e., $\smash{f_{\mathrm{pwm}}\left(\alpha+\delta\,e_{1}\right)\rightarrow f_{\mathrm{pwm}}\left(\alpha\right)}$ by continuity of $f_{\mathrm{pwm}}$. As the gradient is the same for all $\smash{\delta\!\in\!(0,\delta_{\min})}$, equation (10) correctly computes an element of the general upper subdifferential of $f_{\mathrm{pwm}}$ at $\alpha$.

Finally, we need to show that the maximum on line 8 of Algorithm 1 is always attained. First note that the objective coefficient vector on line 8 of Algorithm 1 is given by,

[TABLE]

By Assumption III.2 we have that,

[TABLE]

attains its maximum for all $\smash{x_{c,i}}$, $\smash{i\!=\!1,\dots,N_{c}}$, and denote $f_{i}^{\ast}$ as the optimal value. Thus the hyperplanes $\smash{\alpha^{\intercal}\phi(x_{c,i})\!\leq\!f_{i}^{\ast}}$ are all supporting hyperplanes of the convex constraint set $BI(\mathcal{A}_{\mathrm{con}})$. The following finite dimensional linear program relaxation of line 8 of Algorithm 1 also attains its maximum,

[TABLE]

To show this, first observe that (25) is feasible and bounded above by $\sum\nolimits_{\left(\mathcal{I}_{>}\,\cup\,\mathcal{I}_{=}\right)}\smash{f_{i}^{\ast}}$, and thus by [36, Corollary 27.3.2] problem (25) attains its maximum. Finally, by [36, Corollary 27.3.3] we have that attainment for (25) implies attainment for line 8 of Algorithm 1.

We have shown that the assumptions of Theorem III.3 satisfy also the assumptions of Theorem B.1 and hence the claims follow from Theorem B.1. ∎

Appendix C Problem setting for Section IV

C-A Quadratic Basis Functions

The space of quadratic functions is parameterized by a constant offset $\smash{s\in\mathbb{R}}$, a linear co-efficient $\smash{p\in\mathbb{R}^{n_{x}}}$, and a quadratic coefficient as a symmetric matrix $\smash{P\in\mathbb{S}^{n_{x}}}$. Thus we express the restricted function space as,

[TABLE]

Thus the $\alpha$ is the stacked vector of $s$, $p$, and the unique elements of $P$, and the basis functions $\phi$ are the monomials of $x$ up to degree two, which clearly satisfy Assumption III.1. Similar to Section III-A we use a subscript on $s$, $p$, and $P$ to label the approximate value function they correspond to, for example, $\smash{\hat{V}_{j}\!\in\!\hat{\mathcal{F}}\!\left(\mathcal{X}\right)}$ is equivalent to $\smash{\hat{V}_{j}(x)\!=\!x^{\intercal}P_{j}x+p_{j}^{\intercal}x+s_{j}}$ for all $x\!\in\!\mathcal{X}$. This space of convex quadratic functions is considered by restricting matrix $P$ to be positive semi-definite.

C-B Formulating line 8 of Algorithm 1 for commercial solver

See Section IV for the definitions of $A$, $B_{u}$, and $B_{\xi}$ as the linear dynamics, and Appendix C-A for the specification of the quadratic basis functions. We introduce $\smash{\underline{u}_{i},\overline{u}_{i}\in\mathbb{R}}$, $\smash{i\!=\!1,\dots,n_{u}}$, with $\smash{\underline{u}_{i}<\overline{u}_{i}}$, to denote the lower and upper bounds that describe each coordinate of the $\smash{\mathcal{U}\subseteq\mathbb{R}^{n_{u}}}$ space. The quadratic stage cost is condensed into the matrix $\smash{L\in\mathbb{R}^{(n_{x}+n_{u}+1)\times(n_{x}+n_{u}+1)}}$ that takes the form $\smash{l(x,u)=[x^{\intercal},u^{\intercal},1]\,L\,[x^{\intercal},u^{\intercal},1]^{\intercal}}$. The notation $\smash{\mathrm{diag}\left(\cdot\right)}$ places the vector argument on the diagonal of an otherwise zero matrix, and $e_{i}$ is the standard basis column vector with $1$ in the $i^{\mathrm{th}}$ element and zeros elsewhere, with the dimension clear from context. We overload the notation $\hat{V}$ and introduce the notation $\hat{\mathbb{V}}$ as the following matrices,

[TABLE]

where $\star$ indicates that the matrix is symmetric. Again, any subscript $\hat{V}_{(\cdot)}$, $\mathbb{V}_{(\cdot)}$ also applies to $s$, $p$, and $P$. Both matrices are symmetric with dimension $\smash{(n_{x}+n_{u}+1)}$.

Using this notation, the point-wise maximum Bellman inequality (9b), repeated here for convenience,

[TABLE]

is sufficiently reformulated as the following LMI:

[TABLE]

The $s$, $p$, and $P$ in $\hat{V}$ are decision variables, as well as the $\smash{\lambda_{i}\in\mathbb{R}_{+}}$ and $\smash{\lambda_{\bar{\alpha}}\in\mathbb{R}_{+}}$, with everything else as fixed problem data. The $\smash{\lambda_{i}}$ are the auxiliary variables introduced when using the S-procedure to reformulate the for all $\smash{u\in\mathcal{U}}$ part of the constraint in Appendix A-C, while the $\smash{\lambda_{\bar{\alpha}}}$ are the auxiliary variables described in Appendix A-C. The objective function on line 8 of Algorithm 1 is linear in the decision variables, and when computed as per line 4 of Algorithm 1 it requires computation of the first and second moments of the $x_{c,i}$ for the indices, $\smash{i=1,\dots,N_{c}}$, where the approximate value function under consideration dominates $\bar{V}_{\mathrm{con}}$. Letting $\mu_{c}$ and $\Sigma_{c}$ denote the first and second moments respectively, the problem on line 8 of Algorithm 1 becomes

[TABLE]

where $\mathrm{tr}\left(\cdot\right)$ denotes the trace of a square matrix. Note that the constraint $\smash{P\,\succeq\,0}$ can be added to restrict to the space of convex quadratic functions.