Polynomial approximation of high-dimensional Hamilton-Jacobi-Bellman equations and applications to feedback control of semilinear parabolic PDEs

TL;DR

This paper introduces a numerical method combining pseudospectral collocation and iterative Galerkin approximation to solve high-dimensional Hamilton-Jacobi-Bellman equations for feedback control of semilinear parabolic PDEs, achieving stabilization in systems up to 14 dimensions.

Contribution

It presents a novel computational approach for high-dimensional HJB equations using polynomial approximation and iterative methods, enabling feedback control synthesis for complex PDE systems.

Findings

Successfully approximated HJB equations in systems up to 14 dimensions.

Developed a stable iterative scheme for nonlinear HJB equations.

Demonstrated effectiveness in feedback control of semilinear parabolic PDEs.

Abstract

A procedure for the numerical approximation of high-dimensional Hamilton-Jacobi-Bellman (HJB) equations associated to optimal feedback control problems for semilinear parabolic equations is proposed. Its main ingredients are a pseudospectral collocation approximation of the PDE dynamics, and an iterative method for the nonlinear HJB equation associated to the feedback synthesis. The latter is known as the Successive Galerkin Approximation. It can also be interpreted as Newton iteration for the HJB equation. At every step, the associated linear Generalized HJB equation is approximated via a separable polynomial approximation ansatz. Stabilizing feedback controls are obtained from solutions to the HJB equations for systems of dimension up to fourteen.

Tables

Table 1. Table 1: Number of elements n 𝑛 n in the basis, as a function of the dimension d 𝑑 d and the total polynomial degree M 𝑀 M . The global polynomial approximation partially circumvents the curse of dimensionality, as the dimension of the basis no longer depends exponentially on the dimension, but rather combinatorially.

	Full monomial basis				Even-degree monomials
$d$\$M$	2	4	6	8	2	4	6	8
6	27	209	923	3002	21	147	609	1896
8	44	494	3002	12869	36	366	2082	8517
10	65	1000	8007	43757	55	770	5775	30085
12	90	1819	18563	125969	78	1443	13819	89401
14	119	3059	38759	319769	105	2485	29617	233107

Table 2. Table 2: 1D polynomial approximation of the infinite horizon control problem with nonquadratic running cost. The number n 𝑛 n denotes the total number of basis functions.

	Monomial basis		Legendre basis
$n$(degree)	error	iterations	error	iterations
2	1.1539	53	1.4127	52
4	0.2541	49	0.3643	58
6	0.015	52	0.0206	52
8	5.01$\times {10}^{-4}$	55	6.41$\times {10}^{-4}$	53
10	8.33$\times {10}^{-6}$	55	1.072$\times {10}^{-5}$	55

Table 3. Table 3: CPU times for different tests and dimensions. CPU-assembly corresponds to the amount of time spent in offline assembly of the different terms of the Galerkin residual equation ( 13 ). CPU-iterative refers to the amount of time spent inside Algorithm 2 .

Test	Dimension	CPU-assembly	CPU-iterative (#)
1	10	2.061$\times {10}^3$[s]	4.221$\times {10}^2$[s](32)
1	12	1.945$\times {10}^4$[s]	3.377$\times {10}^3$[s](32)
4	14	1.557$\times {10}^5$[s]	3.102$\times {10}^4$[s](37)

Equations

$$\underset{u(\cdot )\in \mathcal{U}}{\min }\mathcal{J}(u(\cdot ),x_0):=\underset{0}{\overset{\infty }{\int }}\ell (x(t))+\gamma |u(t){|}^{2}dt$$

$$\dot{x}(t)=f(x(t))+g(x)u(t),x(0)=x_0,$$

$$V(x_0)=\underset{u(\cdot )\in \mathcal{U}}{\inf }J(u(\cdot ),x_0)$$

$$\underset{u\in U}{\min }\{DV(x{)}^t(f(x)+g(x)u)+\ell (x)+\gamma |u{|}^2\}=0,V(0)=0,$$

$$u^{\ast }(x)=\underset{u\in U}{\mathrm{argmin}}\{DV(x{)}^t(f(x)+gu)+\ell (x)+\gamma |u{|}^2\}=-\frac{1}{2\gamma }g(x{)}^{t}DV(x).$$

$$DV(x{)}^{t}f(x)-\frac{1}{4\gamma }DV(x{)}^{t}g(x)g(x{)}^{t}DV(x)+\ell (x)=0,$$

$$f_i(x)=\sum _{j=1}^{n_f}\prod _{k=1}^d{\mathcal{F}}_{(i,j,k)}(x_k),$$

$$\underset{u(\cdot )\in L^2([0;+\infty ))}{\min }\mathcal{J}(u(\cdot ,X_0):=\underset{0}{\overset{\infty }{\int }}\Vert X(\cdot ,t){\Vert }_{L^2(\mathcal{I})}^2+|u(t){|}^{2}dt$$

$${\partial }_{t}X(\xi ,t)$$

$${\partial }_{\xi }X(-1,t)$$

$$\dot{X}(t)=AX(t)-X(t{)}^3+Bu(t),$$

$${\dot{X}}_i(t)=A_{i,1}{X}_1(t)+\dots +A_{i,d}{X}_d(t)-X_i(t{)}^3+B_{i}u(t),$$

$$\underset{u\in U}{\min }\{DV(x{)}^t(f(x)+gu)+\ell (x)+\gamma |u{|}^2\}=0,V(0)=0,$$

$$D{V}^i(x{)}^t(f(x)+g{u}^i)+\ell (x)+\gamma |u^i{|}^2=0,V^i(0)=0.$$

$$u^{i+1}(x)=-\frac{1}{2\gamma }{g}^{t}D{V}^i(x),$$

$$\mathcal{G}(DV;u)=$$

$$\mathcal{G}(p,u):=$$

$$\underset{u(\cdot )\in \mathcal{U}}{\min }\mathcal{J}(u(\cdot ),x_0):=\underset{0}{\overset{\infty }{\int }}{e}^{-\lambda t}(\ell (x(t))+\gamma |u(t){|}^2)dt,\lambda >0,$$

$$\lambda V(x)+\underset{u\in U}{\min }\{DV(x{)}^t(f(x)+gu)+\ell (x)+\gamma |u{|}^2\}=0,V(0)=0,$$

$${\mathcal{G}}_{\lambda }(V,DV;u)=$$

$${\mathcal{G}}_{\lambda }(q,p,u):=$$

$$\lambda V(x)+\underset{u\in U}{\min }\{DV(x{)}^t(f(x)+gu)+\ell (x)+\gamma |u{|}^2\}=0,$$

$$u^0$$

$$\lambda$$

$${\mathcal{G}}_{\lambda }(V,DV;u)=0,V(0)=0.$$

$$V_n(x)=\sum _{j=1}^n{c}_j{\varphi }_j(x)\equiv {\Phi }_n\mathbf{c},$$

$$⟨{\mathcal{G}}_{\lambda }(V_n,D{V}_n;u),{\varphi }_i{⟩}_{L^2(\Omega )}=0,\forall {\varphi }_i\in {\Phi }_n.$$

$$u(x)=-\frac{1}{2}{\gamma }^{-1}{g}^{t}D{V}_n^0(x),$$

$$V_n^0(x)=\sum _{j=1}^n{c}_j^0{\varphi }_j(x).$$

$$⟨\lambda V_n+D{V}_n^t(f(x)+gu)+\ell (x)+\gamma |u{|}^2,{\varphi }_i{⟩}_{L^2(\Omega )}=0,\forall {\varphi }_i\in {\Phi }_n.$$

$$⟨\lambda V_n,{\varphi }_i{⟩}_{L^2(\Omega )}={\mathbf{M}}_{(i,\bullet )}\mathbf{c},\mathbf{M}\in {\mathbb{R}}^{n\times n},{\mathbf{M}}_{(i,j)}=\lambda ⟨{\varphi }_i,{\varphi }_j{⟩}_{L^2(\Omega )}.$$

$$D{V}_n^{t}f=\sum _{j=1}^n{c}_{j}D{\varphi }_j^{t}f,$$

Full text

{AMS}

49J20, 49LXX, 49MXX

1 Introduction

Optimal feedback controls for evolutionary control systems are of significant practical importance. Differently from open-loop optimal controls, they do not rely on knowledge of the initial condition and they can achieve design objectives, as for instance stabilisation, also in the presence of perturbations. Furthermore, the online synthesis of feedback control can be implemented in a real-time setting. It is well-known that their construction relies on special Hamilton-Jacobi-Bellman (HJB) equations, see for instance [4, 21]. The solution of the HJB equation is the value function associated to the optimal control problem, and its gradient is used to construct the optimal feedback control. In the very special, but important case of a linear control system with quadratic cost without constraints on the control or the state variables, the HJB equation reduces to a Riccati equation which has received a tremendous amount of attention, both for the cases when the control system is related to ordinary or to partial differential equations. Otherwise one has to deal with the HJB equation which is a partial differential equation whose spatial dimension is that of the control system. Thus optimal feedback control for partial differential equations leads to HJB equations in infinite dimensions [20]. After semi-discretization in space of the controlled partial differential equation (PDE), the HJB equation is posed in a space of dimension corresponding to the spatial discretization of the PDE [22]. For standard finite element or finite difference discretizations this leads to high-dimensional HJB equations. This is one of the instances which is referred to as the curse of dimensionality [9].

Many attempts to tackle the difficulties posed for numerically solving the HJB equations arising in optimal control have been made in the past or are currently being investigated. We refer, for instance, to [21], which mainly focuses on semi-Lagrangian schemes, and further references given there. A related approach to numerical optimal feedback control of PDEs is to semi-discretize the dynamics and to add a model order reduction step, either with Balanced Truncation or Proper Orthogonal Decomposition, in order to reduce the dimension of the dynamics to a number that is tractable for grid-based, semi-Lagrangian schemes. This approach has been successfully explored, for instance, in [1, 32, 35] and references therein. It strongly relies on a trustworthy representation of the dynamics via low-dimensional manifolds. Such a low-dimensional representation may deteriorate when nonlinear and/or advection effects are relevant. Thus, it is important to strive for techniques, or combinations of techniques, which allow to solve higher dimensional problems.

Another direction of research evolves around generalizing the Riccati-based approach to allow for nonlinearities in the state equation. One such technique is termed state-dependent Riccati equation [19]. Here the coefficients in the ’ordinary’ Riccati equation are functions of the state rather than constants as in the case of linear state equations. Another approach realizes the fact that the Riccati equation can be interpreted as the equation satisfied by the first term arising in the power series expansion of the value function, and attempts to improve by realizing also higher order terms in the expansion. These methods are succinctly explained in [7].

Yet another technique which has received a considerable amount of attention is termed Successive Galerkin Approximation. Roughly speaking, the nonlinear HJB equation associated to the continuous-time optimal control problem is solved by means of a Newton method. At each iteration, the control law is fixed. This leads to a Generalized Hamilton-Jacobi equation (GHJB) which is linear. The iteration is closed by an update of the control law based on the gradient of the value function. This method was intensively investigated in [5, 6], see also [7], and the references given in these citations. It is worth mentioning that the discrete-time counterpart of this method corresponds to the well-known policy iteration or Howards’ algorithm [31, 12, 2].

The numerical examples in [5, 6, 7] do not go beyond dimension five, and most, if not all, of the published numerical results for nonlinear HJB equations do not exceed dimension eight [11, 24, 26]. An alternative sparse grid approach for high-dimensional approximation of HJB equations based on open-loop optimal control has been presented in [33], with tests up to dimension six. Numerical methods relying on tensor calculus have been shown to perform well in high-dimensional settings where the associated HJB equation is a linear PDE [41]. A key feature of these works is the use of sparse tensor products either for the construction of the basis or for the representation of the solution. This idea constitutes a cornerstone of high-dimensional approximation [10, 18], and their applicability ranges from sparse grid approximations [15], to polynomial chaos expansion [40, 29] and uncertainty quantification [28].

In the present paper, to solve optimal control problems for certain classes of semilinear parabolic equations we shall proceed as follows. To accommodate the curse of dimensionality, the discretization of the PDE is based on a pseudospectral collocation method, allowing a higher degree of accuracy with relatively few collocation points. To solve the resulting HJB we utilize a Newton method based on the GHJB equation as described above. Next, the discretization of the GHJB equation is addressed through a Galerkin approximation with polynomial, globally supported, ansatz functions. While this mitigates the curse of dimensionality in terms of removing the mesh structure, it leads to high-dimensional integrals. We therefore resort to separable representations for the system dynamics and for the basis set of the polynomial approximation. The separability assumption reduces the computation of the Galerkin residual equation to products of one-dimensional integrals. The combination of these procedures allowed us to solve HJB equations related to nonlinear control systems up to dimension fourteen by means of basic parallelization tools. The successful use of the Newton procedure requires to provide a feasibly initialization, i.e. a sub-optimal, stabilizing control. Since we do not consider constraints, this is not restrictive for finite horizon problem, but can be challenging for infinite horizon problems, and specifically for the stabilization problems which are considered in the present paper. In this respect we developed a continuation procedure based on the use of a discount factor. Specifically, we consider a nested iterative procedure: within the outer loop the value of a positive discount factor is driven to zero, within the inner loop the HJB equation is solved approximately for a fixed discount factor. With this approach, which, is summarized in Algorithms 1 and 2 below, we managed to solve optimal feedback stabilization problems for semilinear parabolic equations with different stability behavior of the desired steady state.

Let us give a brief outline of the paper. Section 2 sets the stage and provides the discussion of a special case to facilitate the understanding of the following material. In Section 3 the solution process of the HJB equation is detailed. In Section 4 we provide the formulas which are needed to numerically realize the discretized HJB equation after a separable basis has been chosen. Numerical experiments are documented in Section 5. There we can also find comparisons to suboptimal feedback strategies based on Riccati and asymptotic expansion techniques.

2 Infinite horizon optimal feedback control

We consider the following undiscounted infinite horizon optimal control problem:

[TABLE]

subject to the nonlinear dynamical constraint

[TABLE]

where we denote the state $x(t)=(x_{1}(t),\ldots,x_{d}(t))^{t}\in\mathbb{R}^{d}$, the control $u(\cdot)\in\mathcal{U}$, with $\mathcal{U}=\{u(t):\,\mathbb{R}_{+}\rightarrow U\subset\mathbb{R}^{m}\}$, the state running cost $\ell(x)>0$, and the control penalization $\gamma>0$. Furthermore, we assume the running cost and the system dynamics $f(x):\mathbb{R}^{d}\rightarrow\mathbb{R}^{d}$ and $g(x):\mathbb{R}^{d}\rightarrow\mathbb{R}^{d\times m}$ to be $\mathcal{C}^{1}(\mathbb{R}^{d})$. Throughout it is assumed that $f(0)=0$ and $\ell(0)=0$. Our focus is therefore asymptotic stabilization to the origin.

It is well-known that the optimal value function

[TABLE]

characterizing the solution of this infinite horizon control problem is the unique viscosity solution of the Hamilton-Jacobi-Bellman equation

[TABLE]

with $DV(x)=(\partial_{{x}_{1}}V,\ldots,\partial_{{x}_{d}}V)^{t}$. Here we follow the convention of dropping the subscript of $x_{0}$. We study this equation in the unconstrained case, i.e., $U\equiv\mathbb{R}^{m}$, where the explicit minimizer $u^{*}$ of (1) is given by

[TABLE]

note that by inserting this expression for the optimal control in (1), we obtain the equivalent HJB equation

[TABLE]

which under further assumptions can be simplified to the Riccati equation associated to linear-quadratic infinite horizon optimal feedback control.

The methodology we present in this work is applicable to systems fitting the aforedescribed setting, although for the sake of simplicity we restrict the presentation by the following choices:

• (i)

the control $u(t)$ is a scalar variable, i.e. $m=1$.

• (ii)

the running cost $\ell(x)$ is quadratic, i.e. $x^{T}Qx$, with $Q$ positive-definite,

• (iii)

the control term $g(x)\equiv g$ is a constant vector in $\mathbb{R}^{d}$.

At this point, our setting differs from the linear-quadratic case as it allows nonlinear dynamics. For the numerical scheme that we develop the following assumption is crucial:

Assumption 1**.**

The free dynamics $f(x):\mathbb{R}^{d}\to\mathbb{R}^{d},f(x):=(f_{1}(x),\ldots,f_{d}(x))^{t}$ is a sum of separable functions in every coordinate $f_{i}(x)$

[TABLE]

*where $\mathcal{F}(x):\mathbb{R}^{d}\rightarrow\mathbb{R}^{d\times n_{f}\times d}$ is a tensor-valued function. In the case $g=g(x)$, then we shall also assume a similar separable structure for $g(x)$. *

Separated representations are a fundamental tool for mitigating the curse of dimensionality, often leading to algorithms that scale linearly in $d$. Its computational efficiency depends on the number of summands or separation rank ($n_{f}$ above). In this work the we shall assume the separated representation is exact. This can be readily checked as we consider dynamical systems arising from the application of the method of lines to semilinear parabolic PDEs. This translates into working over an ODE control system of the form $\dot{x}(t)=Ax(t)+N(x(t))+u(t)\,,$ where $A$ is a linear operator, and $N(x(t))$ is a polynomial source term. Even though such a setting is quite general and covers a wide class of meaningful problems, it excludes an important set of agent-based control models where the governing dynamics depend on a metric interaction between states, i.e. $f(x)=f(\|x_{i}-x_{j}\|)$ (see [13] and references therein for control-related examples). Nonetheless, it is possible to address the problem of finding an approximate best separated representation of fixed rank $n_{f}$, but this procedure and its error analysis is beyond the scope of the present work. The interested reader can find in [10] a thorough presentation of this topic, with a concrete application to linear HJB equations in [41].

Under the framework provided by Assumption 1, the methodology can be directly applied to multidimensional control signals, non-quadratic state costs, and state dependent $g(x)$ which corresponds to bilinear control systems. The additional computational cost of addressing multidimensional control signals scales linearly with the dimension of $u$, whereas the computational burden associated to a bilinear control system will depend on the separability degree of $g(x)$. Non-quadratic, separable state costs can incorporated at a negligible computational cost with straightforward modifications of our setting.

2.1 Towards optimal feedback control of semilinear parabolic equations

In the following, we illustrate how the presented framework sets the grounds for a computational approach for approximate optimal feedback controllers for nonlinear PDEs. We consider the following optimal stabilization problem:

[TABLE]

subject to the semilinear parabolic equation

[TABLE]

In this case, the scalar control acts through the indicator function $\chi_{\omega}(\xi)$, with $\omega\subset\mathcal{I}$. At the abstract level, this corresponds to an infinite-dimensional optimal control problem. A first step towards the application of the proposed framework is the space discretization of the system dynamics, leading to finite-dimensional state space representation. The use of the pseudospectral collocation methods for parabolic equations has been studied in [37, 39], and leads to a state space representation of the form

[TABLE]

where the discrete state $X(t)=(X_{1}(t),\ldots,X_{d}(t))^{t}\in\mathbb{R}^{d}$ corresponds to the approximation of $X(\xi,t)$ at $d$ collocation points $\xi_{i}=-cos(\pi i/d)$, $i=1,\ldots,d,$ and $X^{3}$ is the coordinatewise power. The matrices $A\in\mathcal{M}^{d\times d}$ and $B\in\mathbb{R}^{d}$ are finite-dimensional approximations of the Laplacian and control operators, respectively. Such a discretization of the dynamics directly fulfills the separability required in Assumption 1, as the i-th equation of the dynamics reads

[TABLE]

with a separability degree $n_{f}=d+1$. It is very important to note that semidiscretization in space of a wide class of time-dependent PDEs will lead to finite-dimensional state space representations of this type, thus the applicability of the presented framework is only limited by the dimensionality of the associated HJB equation. This motivates the choice of a pseudospectral collocation method for the discretization, as it is possible to obtain a meaningful representation of the dynamics with considerably fewer degrees of freedom than classical low-order schemes. However, if pseudospectral collocation is not a suitable discretization method for the dynamics, model reduction procedures such as balanced truncation, proper orthogonal decomposition, or reduced basis techniques shall also lead to separable state-space representations. Once the finite-dimensional state state space representation is obtained, we proceed to approximate the solution of the associated HJB equation (1), leading to the optimal feedback controller (2).

We now present a preview of the numerical results of the proposed approach. Further details of the numerical scheme will be developed in the forthcoming sections. The system dynamics in (5), are approximated in 12 collocation points (14 with b.c.s’), and therefore our approximation scheme seeks for a solution of a 12-dimensional HJB equation, which allows the computation of online optimal feedback controllers. We compare our HJB-based controller (HJB) to the linear-quadratic controller (LQR) obtained by linearization of the system dynamics, and to an approximation method for the HJB equation based on power series expansion (PSE) [25, 42]. In Figure 1 we observe the basic features of the dynamics and the control schemes. The uncontrolled system dynamics (diffusion+dissipative source term) are stable, but stabilization is extremely slow. The control algorithms considerably reduce the transient phase. However, the control signals are different, and the HJB-based controller generates a feedback control with reduced overall cost (4). Observe that at the beginning of the time horizon even the signs of the LQR-, PSE-, and HJB-based controls differ.

3 Approximate iterative solution of HJB equations

In this section, we construct a numerical scheme for the approximation of the HJB equation

[TABLE]

where $U\equiv\mathbb{R}$. We recall two additional features in this equation which render the application of classical approximation techniques difficult: the absence of a variational formulation, and the minimization with respect to the control variable $u$, which makes the HJB equation fully nonlinear. The simplest numerical approach to these problems is the use of monotone, grid-based discretizations (finite differences, semi-Lagrangian), in conjunction with a fixed point iteration for the value function $V,$ which typically depends on the use of a discount factor. The so-called “value iteration” procedure was first presented by Bellman in [8], and although it has become a standard solution method for low-dimensional HJB equations, it suffers from three major drawbacks. First, the grid-based character of the scheme makes it inapplicable for high-dimensional dynamics, as the total number of degrees of freedom scales exponentially with respect to the dimension of the dynamical system. This corresponds to the most classical statement of the so-called curse of dimensionality. Second, the contractive mapping includes a minimization procedure which needs to be solved for every grid point at every iteration. Third, the Lipschitz constant of the contractive mapping goes to 1 when the discretization parameter goes to 0, becoming extremely slow for fine-mesh solutions. In order to circumvent these limitations, we develop a numerical scheme combining an iteration on the control variable rather than the value function, together with a polynomial expansion for the value function to mitigate the computational burden associated to mesh-based schemes.

3.1 Successive approximation of HJB equations

In the following, we revisit the method presented in [5, 6], which is referred as Successive Approximation Algorithm. We begin by defining the set of admissible controls.

Definition 3.1** (Admissible control).**

*We say that a feedback mapping $u:=u(x)$ is admissible on $\Omega\subset\mathbb{R}^{d}$, denoted as $u\in\mathcal{A}(\Omega)$, if $u(x)\in\mathcal{C}(\Omega)$, $u(0)=0$, and $\mathcal{J}(u(x(\cdot)),x_{0})<\infty$ for all $x_{0}\in\Omega$. *

Starting from an admissible initial guess $u^{0}(x)$, the Successive Approximation Algorithm (Algorithm 1 below) generates the pair $(V^{*},u^{*})$ which solves equation (6).

Algorithm 1 corresponds to a Newton method for solving equation (6), and in the linear-quadratic setting it is equivalent to the Newton-Kleinmann iteration for solving the Riccati equation. It can be also directly identified with the policy iteration algorithm for HJB equations (see [2] and references therein), although in this context the usual setting includes a discount factor which relaxes the admissibility assumption, as well as discrete-time dynamics. Consequently, it is applied to a Bellman equation with no continuous gradient. In both cases, the core ingredient of the algorithm is to generate a decreasing sequence of values $V^{i}$ by solving an associated sequence of linear problems. In our case this translates into solving, for a given $u(x)$ at each iteration, the Generalized Hamilton-Jacobi-Bellman (GHJB) equation

[TABLE]

The following result from [5] summarizes relevant properties of the GHJB equation.

Proposition 3.2**.**

If $\Omega$ is a compact subset of $\mathbb{R}^{d}$, $f(x)$ is Lipschitz continuous on $\Omega$ and $f(0)=0$, $l(x)\geq 0$ is strictly increasing in $\Omega$, $\gamma>0$, and $u\in\mathcal{A}(\Omega)$, then:

There exists a unique $V(x)\in C^{1}(\Omega)$ satisfying (8). 2. 2.

$V(x)$* is a Lyapunov function of the controlled system.* 3. 3.

$V(x)=\mathcal{J}(u,x)$, for all $x\in\Omega$. 4. 4.

The update $u^{+}(x):=-\frac{1}{2\gamma}g^{t}DV(x)$ satisfies $u^{+}\in\mathcal{A}(\Omega)$. 5. 5.

If $V^{+}$ satisfies $\mathcal{G}(DV^{+};u^{+})=0$, then $V^{+}\leq V$ for all $x\in\Omega$.

3.2 A continuation procedure

A critical aspect of the Successive Approximation Algorithm 1 is its initialization, which requires the existence of an admissible control $u^{0}(x)$ which in view of (4) means that it asymptotically stabilizes all the initial conditions in $\Omega$. For asymptotically stable dynamics, this is trivially satisfied by $u^{0}(x)=0$. For more general cases, the computation of stabilizing feedback controllers is a challenging task. A partial answer is to consider the stabilizing feedback associated to the linearized system dynamics. However, this feedback is only locally stabilizing, and therefore the identification of a suitable domain $\Omega$ where this control law is admissible becomes relevant. For low dimensional dynamics, this has been studied in the context of Zubov’s method in [17]. An alternative solution that we propose is to consider a discounted infinite horizon control problem

[TABLE]

where the inclusion of the discount factor $\lambda$ relaxes the admissibility condition. Recently, in [23, 38], the link between discounted optimal control and asymptotic stabilization has been discussed, and under certain conditions, the discounted control problem can generate optimal controls that are also admissible for the undiscounted problem. We recall that the associated HJB equation for the infinite horizon optimal control problem is given by

[TABLE]

and the associated GHJB reads

[TABLE]

We consequently modify the Successive Approximation Algorithm in order to embed it within a path-following iteration with respect to the discount factor:

For a sufficiently large $\lambda$, this algorithm can be initialized with $u_{\lambda}^{0}=0$. Continued reduction of the discount factor using hotstart every time when (11) is called with a reduced $\lambda$-value, leads to an approximate solution of equation (6).

3.3 Spectral element approximation of the GHJB equation

So far we have discussed the iterative aspects of a computational method for solving HJB equations. We now address the numerical approximation of the GHJB equation.

[TABLE]

For this purpose, we consider an expansion $V_{n}(x)$ of the form

[TABLE]

where $\Phi_{n}:=(\phi_{1}(x),\ldots,\phi_{n}(x))$, with $\phi_{j}\in C^{\infty}(\Omega,\mathbb{R})$ belonging to a complete set of basis functions in $L^{2}(\Omega,\mathbb{R})$, and $\mathbf{c}=(c_{1},\ldots,c_{n})^{t}$. In particular, we shall often generate $\Phi_{n}$ from a multidimensional monomial basis as illustrated in Figure 2, which directly satisfies the boundary condition $V_{n}(0)=0$. The coefficients $c_{j}$ are obtained by imposing the Galerkin residual equation

[TABLE]

Remark 3.3**.**

*The convergence of $V_{n}$ has been studied thoroughly in [5]. It follows a power series argument, and requires conditions for uniform convergence of pointwise convergent series, in order to guarantee that $u_{n}:=-\frac{1}{2}\gamma^{-1}g^{t}DV_{n}(x)\in\mathcal{A}(\Omega)$ for $n$ sufficiently large. In our particular case, we further assume that the dynamics $(f,g)$ are polynomial (as illustrated in Section 2.1). Therefore, under the assumptions of Theorem 26 in [5], by choosing a multidimensional monomial basis (of degree $\geq 2$) and an admissible control $u^{0}\in\mathcal{A}(\Omega)$, it can be established that, $\forall\epsilon>0$, $\exists K$ such that for $n>K$, $\|V-V_{n}\|_{L^{2}(\Omega)}<\epsilon$, and $u_{n}(x)\in\mathcal{A}(\Omega)$. *

We now focus on the different terms involved in the approximation of the GHJB equation. Since this equation is meant to be solved within the iterative loop described in the previous section, we assume that $u(x)$ can be expressed in the form

[TABLE]

where $V^{0}(x)$ corresponds to the value function of the previous iteration, approximated with the expansion

[TABLE]

Below we shall write $\mathbf{c}^{0}$ for $(c^{0}_{1},\dots,c^{0}_{n})^{t}$. We proceed by expanding case by case the different terms of the Galerkin residual equation

[TABLE]

1. 1)

$\boldsymbol{\langle\lambda V_{n},\phi_{i}\rangle_{L^{2}(\Omega)}}$: it is directly verifiable that

[TABLE] 2. 2)

$\boldsymbol{\langle DV_{n}^{t}f,\phi_{i}\rangle_{L^{2}(\Omega)}}$: by inserting the expansion we obtain

[TABLE]

and therefore

[TABLE] 3. 3)

$\boldsymbol{\langle DV_{n}^{t}gu,\phi_{i}\rangle_{L^{2}(\Omega)}}$: the relation (14) leads to

[TABLE]

such that

[TABLE] 4. 4)

$\boldsymbol{\langle l(x),\phi_{i}\rangle_{L^{2}(\Omega)}}$: we further assume that

[TABLE] 5. 5)

$\boldsymbol{\langle\gamma|u|^{2},\phi_{i}\rangle_{L^{2}(\Omega)}}$ : note that

[TABLE]

leading to

[TABLE]

$\mathbf{U}\in\mathbb{R}^{n\times n\times n}$ is given by

[TABLE]

After discretization, the GHJB (13) reduces to a parameter-dependent linear system for $\mathbf{c}$

[TABLE]

where $\mathbf{b}$ is given by the expansion of $l(x)+\gamma|u|^{2}$ ( terms 4) and 5) in the list above).

4 Computation of integrals via separable expansions

Under Assumption 1 concerning the separability of the free dynamics $f$, and with the construction of a separable set of basis functions by taking the tensor product of one-dimensional basis functions as shown in Figure 2, the calculation of the $d$-dimensional inner products of the Galerkin residual equation of the previous section is reduced to the product of one-dimensional integrals. In the following, we provide further details of this procedure.

4.1 Generation of a multi-dimensional basis

The multi-dimensional basis functions

$\Phi_{n}:=(\phi_{1}(x),\ldots,\phi_{n}(x))$ for the expansion of $V_{n}$ are generated as follows. We start by choosing a polynomial degree $M\in\mathbb{N}$, and a one-dimensional polynomial basis $\varphi_{M}:\mathbb{R}\to\mathbb{R}^{M}$. For the sake of simplicity, we consider the monomial basis $\varphi_{M}=(1,x,\ldots,x^{M})^{T}$, but the same ideas apply for other basis, such as orthogonal polynomials. The multidimensional basis is generated as a subset of the $d$-dimensional tensor product of one-dimensional basis, such that

[TABLE]

i.e., we construct a full multidimensional tensorial basis and then we remove elements according to the approximation degree $M$. The elimination step is fundamental and is twofold. If no elimination is performed, the cardinality of $\Phi_{n}$ would be $M^{d}$, and again one would face the curse of dimensionality that also affects grid-based schemes. By reducing the set to multdimensional monomials of degree at most $M$, the cardinality $n$ of the set $\Phi_{n}$ is given by

[TABLE]

which replaces the exponential dependence on $d$ by a combinatorial one. This formula is evaluated in Table 1 for different values of interest for $M$ and $d$. By considering globally defined polynomial basis functions, the dependence on the dimension is replaced by the combinatorial expression (16). The dimensional reduction of the basis is particularly significant for low order polynomial approximation (up to degree 6). A second justification for the way in which we generate the basis set has a control-theoretical inspiration. A well-known result in optimal feedback control is that if the dynamics are linear, and the running cost is quadratic, the value function associated to the infinite horizon control problem (in the unconstrained case and other technical assumptions) is a quadratic form, i.e. is of the form $V(x)=x^{t}\Pi x$, which fits precisely the elements generated for $\Phi_{n}$ with a monomial basis when $M=2$ and linear elements are eliminated. Therefore, our basis can be interpreted as a controlled increment, accounting for the nonlinear dynamics, of the basis required to recover the solution of the control problem associated to the linearized dynamics around the equilibrium point.

Remark 4.1**.**

Theorem 7.1 in [6] states parity conditions to reduce the polynomial basis $\Phi_{n}$. Under the assumptions $l(x)=x^{t}Qx$, and $g\in\mathbb{R}^{d}$, if

• i)

$\Omega$* is a symmetric rectangle around the origin, i.e., $\Omega=[-l_{1},l_{1}]\times\ldots\times[-l_{d},l_{d}]\,,$*

• ii)

the free dynamics are odd-symmetric on $\Omega$, i.e. $f(-x)=-f(x)$, for all $x\in\Omega\,,$

*then $V_{n}(x)$ is an even-symmetric function, i.e., $V_{n}(-x)=V_{n}(x)$, and therefore odd-degree monomials are excluded from the basis. A direct corollary is that in the linear quadratic case, where the linear dynamics are trivially odd-symmetric, $V(x)$ is a quadratic form. *

Finally, for the calculation presented in the following, it is important to note that due to the construction procedure, the basis elements directly admit a separable representation

[TABLE]

where each component $\phi_{i}^{j}(x)\in\varphi_{M}$.

4.2 High-dimensional integration

We begin by recalling that

[TABLE]

where $\mathcal{F}(x):\mathbb{R}^{d}\rightarrow\mathbb{R}^{d\times n_{f}\times d}$ is a tensor-valued function, and that $g\in\mathbb{R}^{d}$.

As in the previous section, we proceed term by term, to obtain the summands in (15). The integration is carried over the hyperrectangle $\Omega=\Omega_{1}\times\ldots\times\Omega_{d}$.

1. 1)

$\boldsymbol{\langle\lambda V_{n},\phi_{i}\rangle_{L^{2}(\Omega)}}$: this term is directly assembled from the calculation of

[TABLE] 2. 2)

$\boldsymbol{\langle DV_{n}^{t}f,\phi_{i}\rangle_{L^{2}(\Omega)}}$: This term involves the calculation of

[TABLE]

which is expanded by using the separable structure of the free dynamics

[TABLE]

where

[TABLE] 3. 3)

$\boldsymbol{\langle DV_{n}^{t}gu,\phi_{i}\rangle_{L^{2}(\Omega)}}$: In this case, we need to work on the expression

[TABLE]

which is obtained directly from the computations for $\boldsymbol{\langle\gamma|u|^{2},\phi_{i}\rangle_{L^{2}(\Omega)}}$ in 5) below. 4. 4)

$\boldsymbol{\langle l(x),\phi_{i}\rangle_{L^{2}(\Omega)}}$:

[TABLE]

where with a similar argument as in the previous term we expand

[TABLE] 5. 5)

$\boldsymbol{\langle\gamma|u|^{2},\phi_{i}\rangle_{L^{2}(\Omega)}}$ : This term requires the computation of the inner product

[TABLE]

with $\tilde{\mathbf{U}}\in\mathbb{R}^{n\times n\times n\times d\times d}$ given by

[TABLE]

By using the separable representation of the basis functions

[TABLE]

we expand the inner product

[TABLE]

Initialization

The first iteration, with a stabilizing initial guess $u^{0}$, requires special attention. If it is obtained via a Riccati-type argument, then initialization follows directly from (14). Otherwise we shall relax this requirement, and only assume that the initial stabilizing controller is given in separable form

[TABLE]

In this case, we must recompute the term:

• •

${\boldsymbol{\langle\gamma|u^{0}|^{2},\phi_{i}\rangle_{L^{2}(\Omega)}}}$

[TABLE]

As for the term $\langle DV_{n}^{t}gu^{0},\phi_{i}\rangle_{L^{2}(\Omega)}$, which needs to be computed differently in the first iteration, we can proceed in the same way as for $\langle DV_{n}^{t}f,\phi_{i}\rangle_{L^{2}(\Omega)}$, since both $gu^{0}$ and $f$ have the same separable structure, it just takes to assign $f_{i}=g_{i}u^{0}$.

4.3 Computational complexity and implementation

Among the expressions developed in the previous subsection, the overall computational burden is governed by the approximation of

[TABLE]

which requires the assembly of the 5-dimensional tensor $\tilde{\mathbf{U}}\in\mathbb{R}^{n\times n\times n\times d\times d}$. Each entry of this tensor is a $d$-dimensional inner product, which under the aforementioned separability assumptions is computed as the product of $d$, one-dimensional integrals. Thus, the total amount of one-dimensional integrals required for the proposed implementation is $O(n^{3}d^{3})$. A positive aspect of our approach is that the assembly of tensors like $\tilde{\mathbf{U}}$ falls within the category of embarrassingly parallelizable computations, so the CPU time scales down almost directly with respect to the number of available cores. Furthermore, $\tilde{\mathbf{U}}$ can be entirely computed in an offline phase, before entering the iterative loops in Algorithms 1 and 2. However, for values of interest of $n$ and $d$, such as $d>10$ and $n=4$, Table 1 indicates that $n^{3}d^{3}$ is indeed a very large number. A rough estimate of the CPU time required for the assembly of $\tilde{\mathbf{U}}$ is given by

[TABLE]

where $t_{1d}$ corresponds to the time required for the computation of a one-dimensional integral. Therefore, it is fundamental for an efficient implementation to reduce $t_{1d}$ to a bare minimum. From closer inspection of the expression (5)), we observe that all the terms can be identified as elements of the tensors $\mathcal{M},\mathcal{K}\in\mathbb{R}^{M\times M\times M}$

[TABLE]

Both $\mathcal{M}$ and $\mathcal{K}$ can be computed exactly with a Computer Algebra System, or approximated under suitable quadrature rules. We follow this latter approach, implementing an 8-point Gauss-Legendre quadrature rule. After having computed $\mathcal{M}$ and $\mathcal{K}$, the assembly of (5)) reduces to $d$ calls to properly indexed elements of these tensors. This approach requires a careful bookkeeping of the separable components of each multidimensional basis function $\phi_{i}$. In this way, an entry of $\tilde{\mathbf{U}}$ takes of the order of $10^{-7}$ seconds and the overall CPU time is kept within hours for problems of dimension up to 12.

5 Computational implementation and numerical tests

5.1 Convergence of the polynomial approximation

We assess the convergence of the polynomial approximation in a 1D test, with

[TABLE]

such that the exact solution of equation (3) is given by

[TABLE]

We implement the path-following version (Algorithm 2), starting with $u^{0}=0$, $\lambda=1$ and a threshold value $\epsilon=1\times 10^{-6}$, a parameter $\beta=0.5$, and an internal tolerance $tol=10^{-8}$. The relative error for Table 2 is defined as

[TABLE]

and number of iterations for different polynomial degree approximations are shown in Table 2 and Figure 3.

5.2 Optimal feedback control of semilinear parabolic equations

Similarly as in Section 2.1, we consider the following optimal control problem

[TABLE]

subject to the semilinear dynamics

[TABLE]

where the linear operator $\mathcal{L}$ is of the form $\mathcal{L}:=\sigma\partial_{\xi\xi}X(\xi,t)+rX(\xi,t)$ with $r\in\mathbb{R}$, and $\mathcal{N}$ is a nonlinear operator such that $\mathcal{N}(0,0)=0$. The scalar control acts through the indicator function $\chi_{\omega}(\xi)$, with $\omega\subset\mathcal{I}$ The system is closed under suitable boundary conditions. We choose $\mathcal{I}=(-1,1)$, $\omega=(-0.5,-0.2)$, $\sigma=0.2$, and $\gamma=0.1$. The nonlinearity covers both advective, Burgers’-type, and polynomial source terms. In order to generate a low-dimensional state space representation of the dynamics, we resort to a pseudospectral collocation method with Chebyshev polynomials as in [37] (for further details we also refer to [39, p. 107]. By considering $d$ collocation points $\xi_{i}=-cos(\pi i/d)\,,i=1,\ldots,d$, the continuous state $X(\xi,t)$ is discretized into $X(t)=(X_{1}(t),\ldots,X_{d}(t))^{t}\in\mathbb{R}^{d}$, where $X_{i}(t)=X(\xi_{i},t)$. The semilinar PDE dynamics are thus approximated by the $d-$dimensional nonlinear ODE system

[TABLE]

where the operators $(A,N,B)$ correspond to the finite-dimensional realization of $(\mathcal{L},\mathcal{N},\chi_{\omega}(\xi))$ through the Chebyshev pseudospectral method. Therefore, the number of collocation points governs the dimension of the resulting nonlinear ODE system (21), and consequently determines the dimension of the domain $\Omega$ where the associated HJB equation is solved. In the following, Tests 1-3 are computed in 14 collocation points, which after including boundary conditions lead to a 12 dimensional domain $\Omega$ for the HJB equation. Test 4 is solved in 14 dimensions. The high-dimensional solver was implemented in MATLAB, parallelizing the tensors assembly, and tests were run on a muti-core architecture 8x Intel Xeon E7-4870 with 2,4Ghz, 1 TB of RAM. The MATLAB pseudoparallelization distributes the tasks among 20 workers. Representative performance details are shown in Table 3. The assembly of high-dimensional tensor that enter the iterative algorithm accounts for over 80% of the total CPU time. This percentage increases when Algorithm 1 is implemented for asymptotically stable dynamics, as it requires a much lower number of iterations. Note that much of the work done during the assembly phase is independent of the dynamics (see for instance (5))), and therefore can be re-used in latter problems, mitigating the overall computational burden.

We now turn to the specification of parameters for the solution of the HJB equation. We set $\Omega=(-2,2)^{d}$, and consider a monomial basis up to order 4 as described in Section 4. Depending on the dynamics of every example, we will neglect odd-degree basis functions as in Remark 4.1. All the integrals are approximated with an 8 point Gauss-Legendre quadrature rule. Whenever system dynamics are stable at the origin, the value function is obtained from the undiscounted Algorithm 1, initialized with $u^{0}=0$. When the dynamics are unstable over $\Omega$, we implement Algorithm 2, with $\lambda=1$, $\epsilon=10^{-6}$, and $\beta=0.9$. The initializing controller is given by the solution of the associated linear-quadratic optimal feedback, as described below. For both implementations, the tolerance of the algorithm is set to $tol=10^{-8}$. In the following tests, we compare the HJB-based feedback control with respect to the uncontrolled dynamics ($u=0$), the linear-quadratic optimal feedback (LQR), and the power series expansion type of controller (PSE). We briefly describe these controllers. The well-known LQR feedback controller corresponds to the HJB synthesis applied over the linearized system around the origin

[TABLE]

and results in the optimal feedback control law given by

[TABLE]

where $\Pi\in R^{d\times d}$ is the unique self-adjoint, positive-definite solution of the algebraic Riccati equation

[TABLE]

and $X^{t}QX$ corresponds to the finite-dimensional approximation of $\|X(\xi,t)\|_{L^{2}(\mathcal{I})}$. Once this controller has been computed, the high-order PSE feedback is obtained as

[TABLE]

where $N_{l}(X)$ corresponds to the lowest order term of the nonlinearity $N(X)$.Variations of such feedback laws have been discussed in previous publications, see eg. [16] and references therein. For the Burgers’ equation it was observed numerically in [42] that this suboptimal nonlinear controller leads to an increased closed-loop stability region with respect to the LQR feedback applied for the linearized dynamics.

Test 1: Viscous Burgers’-like equation

In this first test we address nonlinear optimal stabilization of advective-reactive phenomena, by considering a 1D Burgers’-like model with $(\xi,t)\in\mathcal{I}\times\mathbb{R}^{+}$ given by

[TABLE]

The feedback stabilization of Burgers’ equation (without the exponential source term) has been thoroughly studied in different contexts, including the work of [16], and the recent work [34]. Since our interest is the study of optimal stabilization, we consider an additional source term $1.5X(\xi,t)e^{-0.1X(\xi,t)}$ such that the origin is not asymptotically stable. This can be appreciated in the numerical results shown in Figure 4. For this model, we consider a reduced-order state space representation of 12 states, solving a HJB equation over $\Omega=(-2,2)^{12}$. The value function is approximated with a monomial basis including both even and odd-degree polynomials up to degree 4. In Figure 4 we can compare the uncontrolled solution to the LQR- and HJB-controlled solutions, where the LQR decay is significantly slower that the one of the HJB synthesis. The HJB controller stabilizes at a higher speed, which is reflected both in the plots and in the total costs. The HJB controller obtains a reduction of approximately 18% with respect to the LQR cost. More importantly, the control signals differ in sign, magnitude, and speed. Such a behavior illustrates the nonlinear character of both the control problem and the feedback law.

Test 2: Diffusion with unstable reaction term

We now turn our attention to a diffusion equation with nonlinearity $\mathcal{N}(X)=X^{3}$ (the case with the reversed inequality sign in front of the cubic term was already treated in Subsection 2.1),

[TABLE]

We close the system with Neumann boundary conditions. The origin $X(\xi,t)\equiv 0$ is an unstable equilibrium of the uncontrolled dynamics. Any other initial condition is unstable with finite time blow-up. In this case, feedback controls can only provide local stabilization, and the purpose of this numerical test is to show that HJB-based synthesis leads to an increased closed-loop asymptotic stability region when compared to LQR, and PSE controllers. For this purpose, we compute feedback controls with the LQR, PSE and HJB approaches, for initial conditions of the form $X_{0}(\xi)=\delta(\xi-1)^{2}(\xi+1)^{2}$, with $\delta\in\mathbb{R}^{+}$. The HJB feedback is computed with Algorithm 2 initialized with a nonlinear feedback control law provided by the PSE approach. The test is carried out over $\Omega=(-2,2)^{12}$, and the value function is approximated with monomial basis elements of degree 2 and 4. Numerical results are presented in Figure 5, for $\delta=2$ and for a series of increased values of $\delta$ in Table 4. As the magnitude of the initial condition grows, the locally stabilizing LQR and PSE controllers are not able to prevent the finite blow-up of the dynamics. This eventually also happens for the HJB feedback, but at a much larger value of $\delta$ (we report the last value $\delta=4$ until which the HJB control stabilizes the dynamics).

Test 3: Newell-Whitehead equation

. The diffusion-reaction equation

[TABLE]

corresponds to a particular case of the so-called Schlögl model, whose feedback stabilization has been studied in [14, 27]. This is a special case of a bistable system with $\pm 1$ as stable and [math] as unstable equilibria. Here we use in an essential manner that we consider Neumann boundary conditions. For Dirichlet conditions the only equilibrium is the origin. Such systems arise for instance in Rayleigh-Benard convection and describe excitable systems such as neurons or axons. As in the previous example, the reduced state-space is chosen as $\Omega=(-2,2)^{12}$, and the basis elements for the HJB approach are even degree monomials of degree 2 and 4. Numerical results for the different controllers are shown in Figure 6. While all the feedback laws effectively stabilize the initial condition $X_{0}(\xi)=cos(2\pi\xi)cos(\pi\xi)+2$ to the origin, the HJB feedback has the smallest overall cost $\mathcal{J}(u,X)$. As in Test 1, it can be observed that the three feedback strategies have a considerably different transient behavior. Note that the LQR controller, which neglects the effect of the nonlinearity $\mathcal{N}(X)=-X^{3}$, has an increased control magnitude with respect to the nonlinear controllers which are able to account the dissipative effect of the nonlinearity.

For the sake of completeness, we also consider this test case with a switch of the sign of nonlinearity, i.e., $\mathcal{N}(X)=-X^{3}$. This case is more demanding than Test 2, as now the linear part is $\sigma\partial_{\xi\xi}X+X$. However, the performance of the controllers is similar as in Test 2, and the results are summarized in Table 5. Again, the HJB feedback law has an increased closed-loop stability region compared to the LQR and PSE controllers.

Test 4: Degenerate Zeldovich equation

In this last test case, we consider the model given by

[TABLE]

This equation, which arises for instance in combustion theory, has $X\equiv 1$ as stable and $X\equiv 0$ as unstable equilibria. For this case, we increase the dimension of the HJB domain to 14, i.e., $\Omega=(-2,2)^{14}$, and the basis functions are monomials of odd and even degree up to 4. Numerical results are shown in Figure 7, where it can be seen that the HJB controller yields the smaller overall cost $\mathcal{J}(u,X)$. Note that the PSE controller for this case has a diminished performance as compared even to the LQR controller. This can be explained by the fact that the PSE controller only takes into account the lowest order nonlinearity, in this case $\mathcal{N}_{l}(X)=X^{2}$, neglecting the cubic term. This is a well-known drawback of this controller, and therefore justifies the need of more complex synthesis methods for nonlinear feedback design, such as the proposed HJB approach.

Concluding remarks

A systematic technique for the computational approximation of HJB equations in optimal control problems related to semilinear parabolic equations was presented. To partially circumvent the curse of dimensionality, the dynamics of the parabolic equation are approximated by a pseudospectral collocation method, and the generalized HJB equation is approximated by separable multi-dimensional basis functions of a given order. The numerical results show that the feedback controls obtained by the proposed methodology differ and improve upon applying Riccati approaches to the linearized equations. The generalized HJB approach has been addressed in earlier publications, reporting on numerical results with lower dimensions than here and in part restrained enthusiasm about the numerical performance, possibly due to the lack of a systematic initialization procedure. For the class of problems considered in this paper the results were consistently better than Riccati approaches. The use of the discount factor path-following technique as proposed in Algorithm 2 is essential for stabilizing to unstable equilibria.

Acknowledgments

The authors gratefully acknowledge support by the ERC advanced grant 668998 (OCLOC) under the EU’s H2020 research program. Finally, D. Kalise wishes to dedicate this paper to the memory of Alexander Vasiliev.