Robust feedback control of nonlinear PDEs by numerical approximation of high-dimensional Hamilton-Jacobi-Isaacs equations

TL;DR

This paper introduces a numerical method for designing robust feedback controllers for high-dimensional nonlinear PDEs using Hamilton-Jacobi-Isaacs equations, enabling stabilization and disturbance rejection.

Contribution

It develops a separable polynomial approximation approach to solve high-dimensional Hamilton-Jacobi-Isaacs PDEs for robust control of nonlinear PDEs.

Findings

Effective stabilization of systems up to 12 dimensions

Robust controllers achieve optimal stabilization and disturbance rejection

Applicable to nonlinear parabolic PDEs with uncertainties

Abstract

We propose an approach for the synthesis of robust and optimal feedback controllers for nonlinear PDEs. Our approach considers the approximation of infinite-dimensional control systems by a pseudospectral collocation method, leading to high-dimensional nonlinear dynamics. For the reduced-order model, we construct a robust feedback control based on the $\cH_{\infty}$ control method, which requires the solution of an associated high-dimensional Hamilton-Jacobi-Isaacs nonlinear PDE. The dimensionality of the Isaacs PDE is tackled by means of a separable representation of the control system, and a polynomial approximation ansatz for the corresponding value function. Our method proves to be effective for the robust stabilization of nonlinear dynamics up to dimension $d\approx 12$. We assess the robustness and optimality features of our design over a class of nonlinear parabolic PDEs, including nonlinear advection and reaction terms. The proposed design yields a feedback controller achieving optimal stabilization and disturbance rejection properties, along with providing a modelling framework for the robust control of PDEs under parametric uncertainties.

Tables

Table 1. Table 1: Cardinality of the set of d − limit-from 𝑑 d- dimensional monomials with total degree smaller or equal to M 𝑀 M . The dependence on d 𝑑 d is combinatorial, mitigating the curse of dimensionality for low-degree polynomials.

	Full monomial basis
$d$\$M$	2	4	6	8
6	27	209	923	3002
8	44	494	3002	12869
10	65	1000	8007	43757
12	90	1819	18563	125969
14	119	3059	38759	319769

Table 2. Table 2: 1D polynomial approximation for the HJI equation with nonquadratic running cost. The number n 𝑛 n denotes the polynomial degree of the approximated value function with a monomial and Legendre basis elements. Errors are shown for the value function V n ​ ( x ) subscript 𝑉 𝑛 𝑥 V_{n}(x) and for the optimal feedback u n ​ ( x ) subscript 𝑢 𝑛 𝑥 u_{n}(x) . The number of iterations includes the total number of subiterations related to the disturbance update.

	Monomial basis			Legendre basis
$n$	error $V_n$	error $u_n$	iterations	error $V_n$	error $u_n$	iterations
2	1.4128	0.5919	71	1.4035	0.5846	71
4	0.3643	0.1468	79	0.8184	0.3172	71
6	0.0202	9.80$\times {10}^{-4}$	75	0.0224	0.0101	75
8	6.39$\times {10}^{-4}$	3.39$\times {10}^{-4}$	75	6.97$\times {10}^{-4}$	3.55$\times {10}^{-4}$	75
10	1.43$\times {10}^{-5}$	8.01$\times {10}^{-6}$	76	1.45$\times {10}^{-5}$	8.23$\times {10}^{-6}$	77

Table 3. Table 3: Disturbance support ω 1 subscript 𝜔 1 \omega_{1} , control support ω 2 subscript 𝜔 2 \omega_{2} , and γ ∗ superscript 𝛾 \gamma^{*} for R = 0.1 𝑅 0.1 R=0.1 and P = 1 𝑃 1 P=1 .

	support of disturbance (${\omega }_1$)	support of control (${\omega }_2$)	${\gamma }^{\ast }$
case 1	$(0.5,0.8)$	$(-1,1)$	0.3937
case 2	$(-1,-0.5)$	$(0.5,0.8)$	0.8087
case 3	$(-1,1)$	$(0.5,0.8)$	3.1281

Table 4. Table 4: HJI and HJB feedback laws under sinusoidal disturbances, Burgers’ equation with exponential source term, 𝒥 ( u , w = 0 ; X 0 ) = ∫ 0 ∞ ( ∥ X ( t ) ∥ 2 + ∥ u ∥ R 2 ) d t \mathcal{J}(u,w=0;X_{0})=\int_{0}^{\infty}(\|X(t)\|^{2}+\left\lVert u\right\rVert^{2}_{R})\,dt .

$w\left(t\right)$	${\parallel X\parallel }^2$ with $u_{\gamma }$	${\parallel X\parallel }^2$ with $u_{\infty }$	$𝒥(u_{\gamma },0;X_0)$	$𝒥(u_{\infty },0;X_0)$
$0.1sin\left(10t\right)$	2.371	2.666	2.377	2.669
$10sin\left(10t\right)$	3.891	4.278	4.184	4.605
$0.1sin\left(40t\right)$	2.361	2.659	2.369	2.661
$10sin\left(40t\right)$	2.676	3.022	2.837	3.149

Table 5. Table 5: CPU-assembly corresponds to the amount of time spent in offline assembly of the different terms of the Galerkin residual equation ( 21 ). CPU-iterative refers to the amount of time spent inside Algorithm   2 .

Test	Source term	CPU-assembly	CPU-iterative
1	(Without source term)	$1.161\times {10}^4\left[s\right]$	$3.5078\times {10}^4\left[s\right]$
1	(With source term)	$2.1508\times {10}^4\left[s\right]$	$1.1808\times {10}^5\left[s\right]$
2		$3.0262\times {10}^3\left[s\right]$	$4.8975\times {10}^4\left[s\right]$

Equations

$$\underset{u(\cdot )\in \mathcal{U}}{\min }\mathcal{J}(u(\cdot );x):=\underset{0}{\overset{\infty }{\int }}\ell (y(t))+\Vert u(t){\Vert }_R^{2}dt$$

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

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

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

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

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

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

$$\underset{0}{\overset{T}{\int }}\ell (x)+\Vert u(t){\Vert }_R^{2}dt\le {\gamma }^2\underset{0}{\overset{T}{\int }}{\Vert w(t)\Vert }_P^{2}dt,$$

$$\underset{u\in U}{\min }\underset{w\in W}{\max }\{DV(x{)}^t(f(x)+g(x)u+h(x)w)+\ell (x)+\Vert u{\Vert }_R^2-{\gamma }^2\Vert w{\Vert }_P^2\}=0.$$

$$u_{\gamma }(x)$$

$$w_{\gamma }(x)$$

$$D{V}_{\gamma }(x{)}^{t}f(x)+\frac{1}{4}D{V}_{\gamma }(x{)}^{t}Q(x)D{V}_{\gamma }(x)+l(x)=0,$$

$$Q(x)=\frac{1}{2{\gamma }^2}h(x)P^{-1}h(x{)}^t-g(x)R^{-1}g(x{)}^t.$$

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

$${\partial }_{t}X(\xi ,t)=\mathcal{L}X(\xi ,t)+\mathcal{N}(X(\xi ,t),\nabla X(\xi ,t))+\mathcal{G}(\xi ,X)u(t)+\mathcal{H}(\xi ,X)w(t),$$

$$X(\xi ,t)\approx \sum _{i=1}^d{X}_i(t){\Psi }_i(\xi )$$

$$\dot{X}(t)=AX(t)+N(X(t))+G(X)u(t)+H(X)w(t),$$

$$D{V}_{\gamma }(x{)}^{t}f(x)+\frac{1}{4}D{V}_{\gamma }(x{)}^t\left(\frac{1}{2{\gamma }^2}h(x)P^{-1}h(x{)}^t-g(x)R^{-1}g(x{)}^t\right)D{V}_{\gamma }(x)+\ell (x)=0,$$

$$\underset{u\in U}{\min }\underset{w\in W}{\max }\{DV(x{)}^t(f(x)+g(x)u+h(x)w)+\ell (x)+\Vert u{\Vert }_R^2-{\gamma }^2\Vert w{\Vert }_P^2\}=0.$$

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

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

$$\displaystyle DV^{(i)}(x)^{t}\big{(}f+g(x)u^{(i)}(x)\big{)}+\ell(x)+\|u^{(i)}(x)\|^{2}_{R}=0\,,$$

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

$$\displaystyle DV^{(i,j)}_{\gamma}(x)^{t}\big{(}f+gu^{(i)}_{\gamma}+hw^{(i,j)}_{\gamma}\big{)}+\ell(x)+\|u^{(i)}_{\gamma}\|_{R}^{2}-\gamma^{2}\|w^{(i,j)}_{\gamma}\|^{2}_{P}=0\,,$$

$$w_{\gamma }^{(i,j+1)}(x)=\frac{1}{2{\gamma }^2}{P}^{-1}h(x{)}^{t}D{V}_{\gamma }^{(i,j)}(x),$$

$$u_{\gamma }^{(i+1)}(x)=-\frac{1}{2}{R}^{-1}g(x{)}^{t}D{V}_{\gamma }^{(i,\infty )}(x).$$

$$\displaystyle DV(x)^{t}\big{(}f(x)+g(x)u(x)+h(x)w(x)\big{)}+\ell(x)+\|u(x)\|_{R}^{2}-\gamma^{2}\|w(x)\|^{2}_{P}=0\,.$$

$$\mathcal{G}(V_{\gamma },D{V}_{\gamma };u,w)=0,V_{\gamma }(0)=0,\text{ where }$$

$$\mathcal{G}(q,p;u,w)=p^t(f+gu+hw)+\ell +\Vert u{\Vert }_R^2-{\gamma }^2\Vert w{\Vert }_P^2$$

$$V(x)=\sum _{j=1}^{\infty }{c}_j{\varphi }_j(x),\text{with}{\varphi }_j(0)=0\forall j,$$

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

Full text

{AMS}

49J20, 49L20, 49N35, 93B52, 93B36

1 Introduction

Robustness is a feature of paramount importance in control systems design. That is, given a dynamical system onto which we act through an external control signal, we expect this signal to achieve asymptotic stabilization and to simultaneously compensate the effects of the different disturbances arising in the control loop. Perturbations may appear in the form of noisy measurements, exogenous disturbances, as well as parametric and structural uncertainties. A first step towards robustness in the control action is the synthesis of feedback laws, i.e. a control expressed as a mapping from the current state of the system to the space of possible control actions. Feedback controls, in contrast to open-loop controls which are expressed solely as a function of time, exhibit enhanced stability properties. The design of feedback maps is a central topic in control theory, and has been thoroughly studied in its many facets, spanning ad-hoc and optimization-based approaches for linear as well and nonlinear dynamics (we refer the reader to [9] for a comprehensive overview).

In this paper we are particularly concerned with the design of optimal feedback controllers for nonlinear dynamics with enhanced robustness properties. The natural course of action is the use of Dynamic Programming techniques which establish a functional relation to be satisfied by the value function of the underlying control problem. In this context, once the value function associated to the control has been solved, an optimal feedback map is trivially obtained as a generalized gradient of the value function. This approach leads to a characterization of the value function as a viscosity solution of a nonlinear Hamilton-Jacobi-Bellman (HJB) or Hamilton-Jacobi-Isaacs (HJI) type Partial Differential Equations (PDEs), depending on whether additional robustness properties are required. In the control engineering literature these approaches are known as the $\mathcal{H}_{2}$ control synthesis whenever an optimal feedback is sought through the solution of a HJB equation, and as $\mathcal{H}_{\infty}$ control design when an optimal robust feedback is computed upon the HJI equation or a related variation. The interplay between Hamilton-Jacobi-Bellman and Isaacs’ PDEs and $\mathcal{H}_{2}$ and $\mathcal{H}_{\infty}$ optimal control theory have been thoroughly discussed in [23, 41, 5, 45]. We specifically point the reader to [43, 44] for a study of the link between Isaacs’ equation, differential game theory, and the computation of robust feedback controls via nonlinear $\mathcal{H}_{\infty}$ synthesis. We recall that in the case of linear dynamical systems with no constraints on the control action and a quadratic cost, the Dynamic Programming approach reduces to the computation of Algebraic Riccati Equations, a topic which has been extensively discussed in the mathematical and engineering literature, see for instance [10]. Here, our focus is on nonlinear dynamical systems, and therefore our discussion will be centered around the numerical approximation of HJB and HJI equations (HJ PDEs in short).

The numerical approximation of Hamilton-Jacobi type equations is a topic dating back to the seminal work by Crandall and Lions [18], including a wide range of discretization strategies such as finite differences [21], finite element methods [30], level-set methods [36], domain decomposition algorithms [19, 27], filtered schemes [13] and most notably, semi-Lagrangian schemes [2, 25]. In particular, the numerical approximation of the Hamilton-Jacobi-Isaacs equation has been studied in [7, 6, 15, 24, 20]. For a comprehensive overview of the different approximation techniques available for HJ PDEs we refer the reader to [26]. The aforementioned techniques have proven to overcome the difficulties associated to the nonlinear character of the HJ PDEs. However, since they are grid-based discretization schemes they suffer from the curse of dimensionality and become overwhelmingly expensive for problems with dimension higher than 5. This is particularly challenging in the context of nonlinear optimal control, as the dimension of the associated HJ PDE is determined by the dimension of the state space of the control problem. A partial remedy to this problem is the coupling of model reduction techniques for the dynamics with a grid-based discretization of the HJ equation, an alternative which has been applied in [1, 4, 32, 34]. Another alternative is the use of approximate dynamic programming techniques in the context of reinforcement learning, see e.g. [8]. However, the design of numerical methods for the solution of high-dimensional HJ PDEs remains a daunting task. Along this direction, some encouraging results have been obtained over the last years in connection with the use of sparse grids [11, 28], causality-free methods [33, 22, 47], machine learning techniques [39, 38], tensor calculus [46], graph-tree structures [3], Taylor series expansions [17], and polynomial approximation [31]. In this latter work, we develop a numerical scheme based on a high-dimensional polynomial ansatz for the value function coupled with a Newton-type (policy iteration) method for the solution of the Galerkin residual equation associated to the HJB equation. This method allows to solve HJB PDEs up to dimension 14, becoming an effective tool for the control of dynamical systems arising from the spectral semi-discretization of nonlinear parabolic PDEs.

In this paper, our focus is on the construction of feedback laws for nonlinear parabolic PDEs which are robust with respect uncertainties and structural perturbations. These feedback laws are obtained on the basis of approximations to the Isaacs equation. We are also particulary interested in differences between the HJB and HJI based feedback controls. As a first step towards this goal, we consider a semi-discretization in space of the infinite-dimensional control system by means of a pseudospectral collocation method [37, 35]. While this approach dramatically reduces the dimensionality of the resulting state space, we are still left with a nonlinear dynamical system of prohibitively large for grid-based schemes. Therefore, a fundamental step in our work is the construction of a numerical scheme for the approximation of the high-dimensional Isaacs equation. For this, we extend our previous work [31] to the HJI PDE by considering a high-dimensional polynomial approximation for the value function, together with the use of a bilevel policy iteration algorithm for the nonlinear Galerkin residual equation. The final outcome of our method is an optimal feedback map, constructed by the solution of the HJI equation, with enhanced robustness properties with respect to additive disturbances in the control loop and/or parametric uncertainties in the dynamics. The method is successfully applied to the control of nonlinear parabolic equations such as the viscous Burgers’ equation and the Newell-Whithead equation under different uncertainty scenarios. To our knowledge, this is the first attempt to design optimal robust feedback controls for nonlinear PDEs by directly approximating the high-dimensional HJI equation rather than relying on the application of the classical linear $\mathcal{H}_{\infty}$ method applied over the linearized dynamics.

The rest of the paper is structured as follows. In Section 2 we state the robust optimal feedback control problem for nonlinear dynamics and the associated Hamilton-Jacobi-Isaacs equation. In Section 3 we introduce an approximation scheme for the HJI PDE, and Section 4 is devoted to specific aspects of its implementation for high-dimensional dynamics. Section 5 concludes our work with a thorough computational assessment of the proposed methodology over a class of nonlinear parabolic PDEs.

2 Preliminaries

In the following, we discuss the synthesis of optimal feedback controllers for nonlinear dynamics based on the dynamic programming principle and more precisely, on the solution of static Hamilton-Jacobi-Bellman and Hamilton-Jacobi-Isaacs equations. In the control literature these methods are referred to as the $\mathcal{H}_{2}/\mathcal{H}_{\infty}$ control syntheses, and we briefly revisit them.

2.1 Nonlinear $\mathcal{H}_{2}$ control and the Hamilton-Jacobi-Bellman equation

We consider the following infinite horizon optimal control problem:

[TABLE]

subject to the nonlinear dynamical constraint

[TABLE]

where we denote the state $y(t)=(y_{1}(t),\ldots,y_{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(y)\geq 0$, and the quadratic control penalization term $\|u\|_{R}^{2}=u^{t}Ru$ with $R\in\mathbb{R}^{m\times m},\,R>0$. Furthermore, we assume the running cost $l(y)$, the system dynamics $f(y):\mathbb{R}^{d}\rightarrow\mathbb{R}^{d}$, and $g(y):\mathbb{R}^{d}\rightarrow\mathbb{R}^{d\times m}$ to be $\mathcal{C}^{1}(\mathbb{R}^{d})$, $f(0)=0$ and $\ell(0)=0$. Our focus is therefore on asymptotic stabilization to the origin. By the application of the Dynamic Programming Principle, 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}$. We restrict our analysis to the unconstrained case $U\equiv\mathbb{R}^{m}$ where the minimizer $u^{*}$ of (2) is given explicitly by

[TABLE]

By inserting (3) into (2), we obtain the HJB equation

[TABLE]

to be understood in the classical sense.

Remark 1**.**

*Under the ansatz $\ell(x)=x^{t}Qx$, $f(x)=Ax$ and $V(x)=x^{t}\Pi x$ with with $Q,A,$ and $\Pi\in\mathbb{R}^{d\times d}$, equation (4) becomes the Algebraic Riccati Equation associated to infinite horizon linear-quadratic control. The motivation behind our work is to depart from this framework by considering a nonlinear representation of $f(x)$ and a high-order polynomial expansion for $V(x)$, in order to reduce the curse of dimensionality which arises when attempting to solve (4) over high-dimensional spaces. *

It is always worth to stress that by solving the HJB equation (4), the optimal feedback map $u^{*}(x)$ is obtained as a by-product via (3). For a real-time realization of the optimal control signal $u^{*}(t)$ in (1), the optimal feedback map is stored and evaluated online at the current state of the dynamics $y(t)$.

2.2 Nonlinear $\mathcal{H}_{\infty}$ control and the Hamilton-Jacobi-Isaacs equation

A variation of the $\mathcal{H}_{2}$ control synthesis is generated by considering the system dynamics

[TABLE]

where an additional disturbance signal $w(\cdot)\in\mathcal{W}$, with $\mathcal{W}=\{w(t):\,\mathbb{R}_{+}\rightarrow W\subset\mathbb{R}^{p}\}$ enters the system through $h(y):\mathbb{R}^{d}\rightarrow\mathbb{R}^{d\times p}$. We assume that $y=0$ is an equilibrium of the system for $u=w=0$. In this case, the $\mathcal{H}_{\infty}$ control goal is to achieve both internal stability of the closed-loop dynamics and disturbance attenuation. This translates into considering the design of a feedback law $u=u(x)$ such that for given $\gamma>0$, and for all $T\geq 0$ and $w\in\mathcal{L}_{2}(0,T)$

[TABLE]

where $\left\lVert w(t)\right\rVert^{2}_{P}=w^{t}Pw$ with $P\in\mathbb{R}^{p\times p},\,P>0$.

We say that system (5) has ${\mathcal{L}}_{2}$-gain smaller or equal than $\gamma$, if (6) holds. Subsequently we are interested in finding a minimum value $\gamma^{*}$, such that for all $\gamma\geq\gamma^{*}$ it is possible to find an asymptotically stabilizing feedback law. The calculation of $\gamma^{*}$, the so-called $\mathcal{H}_{\infty}$-norm of the system, is a challenging problem in its own right. Here we shall obtain an estimate by means of a bisection algorithm [16].

Proceeding analogously as in the $\mathcal{H}_{@}$ control synthesis, we follow the application of the Dynamic Programming Principle. For a given $\gamma\geq\gamma^{*}$, the $\mathcal{H}_{\infty}$ synthesis is based on the solution of the Hamilton-Jacobi-Isaacs equation

[TABLE]

Similarly as in the Hamilton-Jacobi equation (2), given an unconstrained disturbance $\mathcal{W}=\mathbb{R}^{p}$ we synthesize a stabilizing feedback law $u_{\gamma}$ and its associated disturbance $w_{\gamma}$ by setting

[TABLE]

where $V_{\gamma}(x)$ solves the Hamilton-Jacobi-Isaacs equation

[TABLE]

with

[TABLE]

For further details see e.g Van Der Schaft [40]. When there is no confusion, we denote $V_{\gamma}(x)$ by $V$. Note that if the disturbance attenuation is neglected by taking $\gamma\to\infty$, (10) becomes the HJB equation related to $\mathcal{H}_{2}$ control.

2.3 Separability assumptions and the control of semidiscretized PDEs

In the following sections, we will discuss the approximation of the aforedescribed HJB and HJI equations linked to $\mathcal{H}_{2}$ and $\mathcal{H}_{\infty}$ nonlinear control synthesis, respectively. In particular, we will focus on the approximate solution of such PDEs when the dimension of the state space of the dynamics to be controlled is large ($d>10$). As it will become clear later, our numerical scheme is based on the use of globally defined polynomials for approximating the value function $V(x)$ along with a Galerkin residual equation, and we will be faced against the computation of high-dimensional integrals. In order to mitigate the computational difficulties associated with the calculation of high-dimensional integrals, the following assumption will be fundamental in the construction of our control algorithm.

Assumption 1**.**

The free dynamics $f(x):\mathbb{R}^{d}\to\mathbb{R}^{d},f(x):=(f_{1}(x),\ldots,f_{d}(x))^{t}$ are separable 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. We shall also assume a similar separable structure for $g(x)$, $h(x)$, and $\ell(x)$. *

Let us note that while this assumption includes a large class of nonlinear optimal control problems, it does not cover dynamics such as those arising in particle systems with interactions governed by distances between agents, e.g. $f(x)=f(\|x_{i}-x_{j}\|)$ (see [14] and references therein for control-related examples).

Here, we are interested in the application of nonlinear optimal control methods for the robust stabilization of nonlinear parabolic PDEs of the form

[TABLE]

where $\mathcal{L}$ is a linear elliptic operator, $\mathcal{N}$ is a nonlinearity, while $\mathcal{G}$ and $\mathcal{H}$ are control and disturbance operators, respectively. At this point two natural questions arise: how can we apply finite-dimensional control techniques to infinite-dimensional dynamics, and whether it is plausible to assume the separability of the dynamics. The answer to both questions rely on the semi-discretization of the dynamics in space and the use of the so-called method of lines [42]. In fact, under the Galerkin ansatz

[TABLE]

with $\{\Psi_{i}\}_{i=1}^{\infty}$ a suitable set of basis functions in space such as local spline approximants, spectral elements or global orthogonal polynomials, the infinite-dimensional dynamics are approximated by a $d$-dimensional dynamical system of the form

[TABLE]

for $X(t)=(X_{1}(t),\ldots,X_{d}(t))$ and where every term is corresponds to a separable expression. The linear operator $A$ is naturally written in separable form, and for instance, in the case of a polynomial nonlinearity $\mathcal{N}(X)=p(x)$ the use of collocation points in space leads to a finite-dimensional realization of the nonlinearity consisting on fully decoupled pointwise evaluations of the polynomial at each collocation point. While the use of finite elements/splines also leads to separable representations, in this work we opt for the use of spectral/pseudospectral collocation methods as the dimension $d$ of the finite-dimensional dynamics must be kept moderate. Once a suitable separable finite-dimensional state space representation of the system has been obtained, both $\mathcal{H}_{2}$ and $\mathcal{H}_{\infty}$ control synthesis methods can be applied, amounting to the solution of HJB/HJI equations over a $d$-dimensional domain. In the following sections we shall focus on the construction of high-dimensional HJI equations arising in $\mathcal{H}_{\infty}$ control. The details concerning the construction of an analogous method for HJB equation in $\mathcal{H}_{2}$ control can be found in [31].

3 Iterative solution of the Hamilton-Jacobi-Isaacs equation

This section addresses the construction of a numerical scheme for the solution of the HJI equation

[TABLE]

which corresponds to the unconstrained version $\mathcal{U}=\mathbb{R}^{m},\mathcal{W}=\mathbb{R}^{p}$ of the Isaacs’ equation

[TABLE]

Both versions of the HJI equation share the same computational difficulties: they correspond to nonlinear PDEs which need to be solved over a high-dimensional space. At first glance, a natural idea to address the quadratic nonlinearity on $DV$ arising in (15) is to utilize a Newton iteration over the continuous PDE or its discretization. In fact, in the context of the Hamilton-Jacobi-Bellman equation

[TABLE]

the application of Newton’s method is equivalent to the well-known policy iteration or Howards’ algorithm [12] for

[TABLE]

which we recall in Algorithm 1.

This algorithm takes as an input an asymptotically stabilizing controller for the nonlinear dynamics (1), and iterates in two steps: a policy evaluation step in which the current feedback $u^{(i)}(x)$ is inserted in the HJB equation (18) to solve for $V^{(i)}(x)$ avoiding the nonlinearity associated to the minimization, and a policy update step in which the formula $u^{(i+1)}(x)=-\frac{1}{2}R^{-1}g(x)^{t}DV^{(i)}(x)$ updates the control. Under the assumption that the initial controller $u^{0}(x)$ is an asymptotically stabilizing feedback, it has been proven in [7] that Algorithm 1 generates a sequence converging to the optimal feedback control $u^{(\infty)}(x)$, and its associated value function $V^{(\infty)}(x)$ which solves (18). This iterative procedure can be extended in a straightforward manner to HJI equations of the form (16) by nesting inside every policy evaluation step, an inner iterative loop to find the optimal adversarial perturbation $w^{(i,\infty)}$ associated to $u^{(i)}$. A version of the policy iteration procedure for HJI equations is presented in Algorithm 2 below.

The convergence of policy iteration-type algorithms for differential games has been discussed in [6] and [7] , among many others. Under the same assumption on the asymptotic stability of the initial guess $u_{\gamma}^{(0)}(x)$ for $w=0$, the algorithm has been shown to converge towards an optimal stabilizing feedback $u_{\gamma}^{(\infty)}(x)$ with an associated perturbation $w_{\gamma}^{(\infty,\infty)}$, and a value function $V_{\gamma}^{(\infty,\infty)}$, which is the unique stabilizing solution of the HJI equation (15).

Remark 2**.**

*Both Algorithm 1 and Algorithm 2 assume the existence of a globally stabilizing controller for initialization. It is often the case that only a locally stabilizing controller is available, i.e. there exists a bounded domain $\Omega\subset\mathbb{R}^{d}$ where the initial conditions are stabilized with this feedback law. In that case the convergence of the algorithms is only valid inside $\Omega$, and this must be taken into account when defining a domain of interest for deriving a numerical approximation. *

Remark 3**.**

Algorithm 2* is valid for $\gamma\geq\gamma^{*}$. In order to estimate $\gamma^{*}$ for a given nonlinear problem, we proceed similarly as in the linear $\mathcal{H}_{\infty}$ problem [16]. We perform a bisection method over $\gamma$ until finding the smallest value for $\gamma$ for which Algorithm 2 converges to a stabilizing feedback law. *

Given an asymptotically stabilizing initialization of Algorithm 2, and having fixed $u(x)$ and $w(x)$ at a policy evaluation step (19), we must solve the linear PDE for $V(x)$

[TABLE]

In the following section, we discuss the construction of a numerical scheme for the solution of this equation over high-dimensional domains by means of a Galerkin formulation with global polynomial basis functions.

4 Solving the linear equations

In the previous section we described how the solution to (7) can be approximated iteratively by solving at each iteration level a high-dimensional linear equation of the form

[TABLE]

for the unknown $V_{\gamma}=V_{\gamma}(x)$, given $u,w,f,g$ and $\ell$ all functions of $x$.

4.1 Galerkin Approximation

To derive a Galerkin approximation to (20) let $\{\phi_{j}\}_{j=1}^{\infty}\in C^{\infty}(\Omega,\mathbb{R})$ be a complete set of basis functions in $L^{2}(\Omega,\mathbb{R})$ so that

[TABLE]

and introduce $V_{n}(x)$ as the approximation of $V_{\gamma}$, of the form

[TABLE]

where $\Phi_{n}:=(\phi_{1}(x),\ldots,\phi_{n}(x))$ and $\mathbf{c}=(c_{1},\ldots,c_{n})^{t}$. The coefficients $c_{j}$ are obtained by imposing the Galerkin residual system

[TABLE]

We now compute the different terms involved in the approximation of (21) once we apply the Galerkin ansatz to $V_{\gamma}^{(i,j)}$ appearing in (19). The value function $V_{\gamma}^{(i,j)}$ is approximated with the expansion

[TABLE]

For each iteration the values of $w^{(i,j)}$ and $u^{(i)}$ are determined according to

[TABLE]

where $V^{(i-1,\infty)}_{n}$ denotes the value function obtained in the last policy evaluation step of the previous $i-$iteration. We next specify the additive terms of the Galerkin residual equation corresponding to (19) with $V_{\gamma}^{(i,j)}$ replaced by $V_{n}^{(i,j)}$, or alternatively (21) with $u=u^{(i)}$ and $w=w^{(i,j)}.$ These expressions are in part already available in [31], but for the sake of coherence we provide all of them here. Below we shall write $\mathbf{c}^{(i,j)}$ for $(c^{(i,j)}_{1},\dots,c^{(i,j)}_{n})^{t}$ ($i,j=0,\ldots,\infty$). Also note, that superscripts $i,j$ refer to the iterations in the loops of Algorithm 2, whereas subscripts $i,j$ represent running indices of the Galerkin discretization:

1. 1)

$\boldsymbol{\langle DV_{n}^{t}f,\phi_{i}\rangle_{L^{2}(\Omega)}}$: The expansion $V_{n}$ leads to $DV_{n}^{t}f=\sum_{j=1}^{n}c_{j}D\phi_{j}^{t}f,$ and hence

[TABLE] 2. 2)

$\boldsymbol{\langle DV_{n}^{t}hw,\phi_{i}\rangle_{L^{2}(\Omega)}}$: We use relation (23) to obtain

[TABLE]

such that

[TABLE] 3. 3)

$\boldsymbol{\langle DV_{n}^{t}gu,\phi_{i}\rangle_{L^{2}(\Omega)}}$: By inserting the expansion (23), we arrive at

[TABLE]

such that

[TABLE] 4. 4)

$\boldsymbol{\langle\ell,\phi_{i}\rangle_{L^{2}(\Omega)}}$: Assuming $\ell(x)=x^{t}Qx$, leads to

[TABLE] 5. 5)

$\boldsymbol{\langle\gamma^{2}\left\lVert w\right\rVert^{2}_{P},\phi_{i}\rangle_{L^{2}(\Omega)}}$ : Note that $\gamma^{2}\left\lVert w\right\rVert^{2}_{P}=\frac{P^{-1}}{4\gamma^{2}}DV_{n}^{(i,j-1)t}hh^{t}DV_{n}^{(i,j-1)},$ hence $\langle\gamma^{2}\left\lVert w\right\rVert^{2}_{P},\phi_{i}\rangle_{L^{2}(\Omega)}=(\mathbf{c}^{(i,j-1)})^{t}\mathbf{W}_{(i,\bullet)}\mathbf{c}^{(i,j-1)},$ where $\mathbf{W}\in\mathbb{R}^{n\times n\times n}$ is given by

[TABLE] 6. 6)

$\boldsymbol{\langle\left\lVert u\right\rVert^{2}_{R},\phi_{i}\rangle_{L^{2}(\Omega)}}$ : Since $\left\lVert u\right\rVert^{2}_{R}\equiv u^{t}Ru=\frac{1}{4}R^{-1}DV_{n}^{(i-1,\infty)t}gg^{t}DV_{n}^{(i-1,\infty)},$ this leads to $\langle\left\lVert u\right\rVert^{2}_{R},\phi_{i}\rangle_{L^{2}(\Omega)}=(\mathbf{c}^{(i-1,\infty)})^{t}\mathbf{U}_{(i,\bullet)}\mathbf{c}^{(i-1,\infty)},$ where $\mathbf{U}\in\mathbb{R}^{n\times n\times n}$ is given by

[TABLE]

After discretization, the discretized generalized Hamilton Jacobi Isaacs equation (21) reduces to the linear system for $\mathbf{c}$

[TABLE]

where $\mathbf{b}$ is given by the expansion of -\Big{(}\ell+\left\lVert u\right\rVert^{2}_{R}-\gamma^{2}\left\lVert w\right\rVert^{2}_{P}\Big{)} ( terms 4 , 5 and 6 in the list above).

4.2 Separable approximation

In order to tackle the computational difficulties associated to the curse of dimensionality and the calculation of high-dimensional integrals, we resort to 1 concerning the separability of the different terms appearing in the control system. The multi-dimensional basis functions $\Phi_{n}:=(\phi_{1}(x),\ldots,\phi_{n}(x))$ used in the expansion of $V_{n}$ are generated by products of one-dimensional elements chosen from a polynomial basis $\varphi_{M}:\mathbb{R}\to\mathbb{R}^{M}$. Here we choose the monomial basis $\varphi_{M}=(1,x,\ldots,x^{M})^{T}$ with $M\in\mathbb{N}$ the degree of the polynomial, however it is also possible to use orthogonal polynomials. The multidimensional basis is then generated as a subset of the $d$-dimensional tensor product of one-dimensional bases such that the maximum total degree is fixed,

[TABLE]

A canonical basis element $\phi_{i}$ in $\Phi_{n}$ is then given by $\phi_{i}=\Pi_{j=1}^{d}x_{j}^{\nu_{j}}$ , with $\sum_{j=1}^{d}\leq M$. The cardinality of the set of monomials in $d$ dimensions with a total degree not greater than $M$ is given by

[TABLE]

which partially mitigates the curse of dimensionality affecting mesh-based discretizations where the dependence of the total number of degrees of freedom is exponential with respect to $d$. Table 1 shows the cardinality of $\Phi_{n}$ for different values of interest for $M$ and $d$.

Having generated a set of separable multidmensional basis function, we proceed as in the previous section, obtaining the summands in (21) . The integration is carried out over the hyperrectangle $\Omega=\Omega_{1}\times\ldots\times\Omega_{d}$. In [31], some of the terms are given in details. For completeness, we repeat them here.

1. 1)

$\boldsymbol{\langle DV_{n}^{t}f,\phi_{i}\rangle_{L^{2}(\Omega)}}$: In this case, using the expansion of $V_{n}$, we require the computation of

[TABLE]

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

[TABLE]

where

[TABLE] 2. 2)

$\boldsymbol{\langle DV_{n}^{t}hw,\phi_{i}\rangle_{L^{2}(\Omega)}}$: This term demands to calculate $\langle h^{t}D\phi_{k}D\phi_{j}^{t}h,\phi_{i}\rangle_{L^{2}(\Omega)}.$ Using the separable expansion, we can write $h_{i}(x)=\sum_{j=1}^{n_{h}}\prod_{k=1}^{d}H_{(i,j,k)}(x_{k}),$ where $H(x):\mathbb{R}^{d}\rightarrow\mathbb{R}^{d\times n_{h}\times d}$ is a tensor valued function. Hence

[TABLE]

with

[TABLE]

Therefore, we obtain

[TABLE] 3. 3)

$\boldsymbol{\langle DV_{n}^{t}gu,\phi_{i}\rangle_{L^{2}(\Omega)}}$: Similar to 2). 4. 4)

$\boldsymbol{\langle\ell(x),\phi_{i}\rangle_{L^{2}(\Omega)}}$: Here, the expression involves

[TABLE]

where $\langle x_{j}x_{k},\phi_{i}\rangle_{L^{2}(\Omega)}$ can be expanded as

[TABLE] 5. 5)

$\boldsymbol{\langle\left\lVert w\right\rVert^{2}_{P},\phi_{i}\rangle_{L^{2}(\Omega)}}$ : This term requires the computation of the inner product

[TABLE]

which is expanded with the help of $\boldsymbol{\langle DV_{n}^{t}hw,\phi_{i}\rangle_{L^{2}(\Omega)}}$ as

[TABLE] 6. 6)

$\boldsymbol{\langle\left\lVert u\right\rVert^{2}_{R},\phi_{i}\rangle_{L^{2}(\Omega)}}$: This term is similar to 5).

5 Numerical experiments

In this section we present numerical experiments illustrating that a combination of spectral techniques for the semi-discretization of a PDE control system, together with a polynomial approximation for the value function solving the associated HJI equation can mitigate the curse of dimensionality. We apply this methodology for the robust feedback design for nonlinear parabolic PDEs.

5.1 Convergence of the polynomial approximation

We begin by assessing the convergence of the polynomial approximation for the HJI equation in a 1D test, setting

[TABLE]

with $R=P=1$, and $\gamma=2$. In this case the exact solution of the HJI equation is given by

[TABLE]

We initialize Algorithm 2 with $u^{0}=-0.01x$ and a threshold value $\epsilon=1\times 10^{-6}$. The relative error for an $n-$degree approximation of $V(x)$, denoted by $V_{n}(x)$, is defined as

[TABLE]

with $\Omega=[-1,1]$. The total number of iterations and errors for different polynomial degree approximations are shown in Table 2 and Fig. 1. We observe a consistent error decay as the approximation degree is increased, while the number of iterations of the policy iteration algorithm remains independent. The use of a monomial basis for approximating the value function yields similar results as using a set of orthogonal Legendre polynomials.

5.2 Robust control of nonlinear parabolic PDEs

In the following tests, we study the design of robust feedback controllers for a class of nonlinear parabolic PDEs including nonlinear advection and reaction terms such as the Burgers or Allen-Cahn equations. In order to embed the infinite-dimensional PDE dynamics into our finite-dimensional robust control framework, we semi-discretize the dynamics in space by means of a pseudospectral collocation method, leading to nonlinear dynamical systems with as many dimensions as degrees of freedom in the semi-discretization.

Test 1: a viscous Burgers’ equation with forcing term

We study the robust optimal feedback stabilization of a control system related to a Burgers’ equation with a nonlinear reaction term. The control problem is stated as

[TABLE]

where $(\xi,t)\in\mathcal{I}\times\mathbb{R}^{+}$ with $\mathcal{I}=(-1,1)$ and $\mathcal{U}=\mathcal{W}=L^{2}([0,+\infty))$. Here $\chi_{\omega_{1}}$ and $\chi_{\omega_{2}}$ are indicator functions for the action of the disturbance and the control supported over subsets $\omega_{1}$ and $\omega_{2}$, respectively. The additional source term $1.5X(\xi,t)e^{-0.1X(\xi,t)}$ prevents the origin to be origin asymptotically stable in the uncontrolled/unperturbed case $u=w=0$. In our numerical tests the viscosity parameter is set $\sigma=0.2$.

In order to reduce the infinite dimensional state equation $\Sigma_{\infty}$ to a finite dimensional dynamical system, we apply a pseudospectral collocation method in space based on Chebyshev polynomials (see e.g. [35, 37], Chapter 4) which leads to a state space representation of the form

[TABLE]

where the discrete state $X(t)=(X_{0}(t),\ldots,X_{d}(t))^{t}\in\mathbb{R}^{d+1}$ corresponds to the approximation of $X(\xi,t)$ at $d$ collocation points, $X_{i}(t)=X(\xi_{i},t)$ with $\xi_{i}=-cos(\pi i/(d+1))$, $i=1,\ldots,d$, while $X(t)\circ DX(t)$ and $X\circ e^{-0.1X}$ denote coordinatewise operations. The matrices $A\in\mathcal{M}^{d\times d}$, $D\in\mathcal{M}^{d\times d}$, $B\in\mathbb{R}^{d}$ and $C\in\mathbb{R}^{d}$ are finite dimensional approximations of the Laplacian, gradient, control and disturbance operators, respectively. The $L^{2}$ norm of the state penalization of the running cost is approximated accordingly. Note that such a state-space representation is consistent with 1, as for each $i=1,\ldots,d$ the state equation reduces to

[TABLE]

where $A_{i,\cdot}$ denotes the $i$-th row of the matrix $A$. In this case, the dynamics have a separability degree $n_{f}=2d+1$. We shall solve the HJI equation

[TABLE]

with $f(x),h(x),g(x)$ and $\ell(x)$ generated with the aforedescribed semi-discretization (25) . The value function $V(x)$ is approximated with $d=12$ over $\Omega=(-2,2)^{12}$ with a monomial basis of total degree up to 4.

In the following numerical tests we will study the effect of the location of domains $\omega_{1}$ and $\omega_{2}$ on the solution of the HJI problem and in particular on the resulting $\gamma^{*}$. We further investigate numerically the behaviour of the HJI control in the case $\gamma$ is chosen differently from the optimal $\gamma^{*}$, and we compare the HJB and HJI controllers in terms of their robustness against additive noise in the state equation.

Effect of $\omega_{1}$ and $\omega_{2}$. Here we are interested in the effect of the control and disturbance supports on the solution of the HJI equation. In particular, we concentrate on the consequences in terms of the $\mathcal{H}_{\infty}$-norm $\gamma^{*}$, the smallest value for which the HJI equation leads to an asymptotically stabilizing feedback control law. This value is computed by a bisection algorithm discarding values of $\gamma$ where Algorithm 2 fails to converge. For this test, we consider the dynamics (25) without the exponential source term. Table 3 reports three cases varying $\omega_{1}$ and $\omega_{2}$. For the largest ratio of size of control to size of disturbance, we obtain the smallest $\gamma^{*}$, see case 1. Case 3 arises from case 1 by increasing the support of the disturbance while decreasing the support of the control, leading to a significant increase of $\gamma^{*}$. In case 2 we decrease the disturbance support compared to case 3, in such a way that $w_{1}\cap w_{2}=\emptyset$, which results in decrease of $\gamma^{*}$. The three cases are consistent with the expectation that as the ratio $|\omega_{2}|/|\omega_{1}|$ increases, stabilization and noise cancellation capabilities are improved resulting in a smaller $\gamma^{*}$ value.

Effect of choosing $\gamma$. For the following test, $\omega_{1}$ and $\omega_{2}$ are chosen as in case 1, above. We report that similar qualitative results were obtained in the other cases. We begin by setting $\gamma=\gamma^{*}$, results for the initial condition $X(\xi,0)=sign(\xi)$ are shown in Fig. 2 (a) for $u$ and $w$ given by the feedback formulas

[TABLE]

For the results in Fig. 2(b) the HJI equation is solved for an increasing sequence of values for $\gamma$. We observe that the corresponding norms of the trajectories $\left\lVert X(\xi,t)\right\rVert^{2}_{L^{2}(\mathcal{I})}$ converge as $\gamma$ increases to the graph for $\left\lVert X(\xi,t)\right\rVert^{2}_{L^{2}(\mathcal{I})}$ which corresponds to the HJB synthesis, equivalent to $\gamma=+\infty$. We report that for $\gamma^{*}$ the graph of $\left\lVert X(\xi,t)\right\rVert^{2}_{L^{2}(\mathcal{I})}$ also tends to 0 as $t\to\infty$, but at a much smaller rate.

Comparison between HJI and HJB feedback laws in the presence of noise. We compare the behaviour of the (25) in the case that the control is computed by means of the HJI eqn. with optimal $\gamma^{*}$, and the HJB feedback law corresponding $\gamma=+\infty$. Different types of noise signals $w$ are tested. Here we include the exponential source term in the Burgers’ equation, we set $R=0.01$ as the weight in the running cost for the control, and $\omega_{1}=\omega_{2}=(-0.8,0.5)$.

With the additional source term, the uncontrolled-unperturbed system is unstable. This is highlighted in Fig. 3 (a), where the norm of the state of the uncontrolled system with initial condition $X(\xi,0)=sign(\xi)$ is compared against the state norms produced by the HJB and HJI feedback laws with $w=0$. We note that both feedback laws generate an asymptotically stable closed loop system.

Next we compare the HJI/HJB synthesis capabilities to reject additive noise. We compute $V_{\gamma^{*}}$ and $V_{\infty}$ to generate the associated feedback control laws, which are evaluated for different choices of the perturbation signal $w(t)$. For a sinusoidal noise of the form $w=\eta\sin(\omega t)$, the total state contributions $\left\lVert X\right\rVert^{2}=\int_{0}^{\infty}\|X(t)\|^{2}\,dt$ as well as the total costs for HJI and HJB closed loops, $\mathcal{J}(u_{\gamma},w=0;X_{0})$ and $\mathcal{J}(u_{\infty},w=0;X_{0})$ respectively, are presented in Table 4 for $X(\xi,0)=sign(\xi)$.

We note that the state cost given by the HJI synthesis is consistently smaller than associated to the HJB feedback control law. This can also be confirmed in Fig. 3 (b). We proceed with other choices of the perturbation signals given by:

• •

$w_{1}(t)=\kappa w_{\gamma}+\eta sin(\omega t)$

• •

$w_{2}(t)=\kappa u_{\gamma}+\eta sin(\omega t)$

• •

$w_{3}(t)=\kappa u_{\infty}+\eta sin(\omega t)$,

where $\kappa$ and $\eta$ are parameters, and $w_{\gamma}$, $u_{\gamma}$ and $u_{\infty}$ are the HJI disturbance, the HJI control, and the HJB control signals respectively, the latter two obtained by simulating first the unperturbed dynamics.

The results (a)-(b) in the first row of Fig. 4, are obtained for $w_{1}$, for the values of $\kappa$ and $\eta$ indicated in the figures. Unlike the results with sinusoidal perturbation only shown in Table 4, where there is no significant difference between the HJI and the HJB closed loops, here the HJI synthesis asymptotically stabilizes the dynamics to the origin, whereas the HJB feedback law fails. Also, from Fig. 4 (c) and (d) we observe that in the presence of the perturbation $w_{2}(t)$, the HJI feedback controls achieve stabilization to the origin, while the HJB feedback fails in this respect. In Fig. 4 (e) and (f) we present results for the disturbance $w_{3}(t)$ with different values of $\kappa$. For $\kappa=-0.8$ the HJB synthesis not only fails to stabilize but generates a finite time blow up, while the HJI feedback is again successful. For $\kappa=0.8$ both control strategies reject the disturbance, with the HJI law having a more effective transient behaviour.

Test 2: the Degenerate Zeldovich equation

In this second test we consider the stabilization of the following diffusion-reaction model with Neumann boundary conditions arising in combustion theory [29]

[TABLE]

with $(\xi,t)\text{ in}\;\mathcal{I}\times\mathbb{R}^{+}$, $\mathcal{I}=(-1,1)$, $\sigma=0.5$, $\omega_{1}=\omega_{2}=(-0.8,0.5)$. In this case $X\equiv 0$ is an unstable and $X\equiv 1$ is a stable steady state for the uncontrolled state equation. The dynamics are semi-dscretized in space with the same collocation method as in Test 1. The reduced state space is chosen to be $\Omega=(-2,2)^{12}$, and the basis functions for the value function are monomials up to degree 4. We take $R=0.01$, $P=1$. Then the $\mathcal{H}_{\infty}$-norm is found to be $\gamma^{*}=0.125$.

In Fig. 5 (a), we depict the norm of the uncontrolled solution and compare it to the HJI and HJB controlled solutions. We observe the two feedback laws generating a similar response in the absence of noise. We also compared HJI and HJB with sinusoidal perturbations and found only marginal differences between the two cases, with a tendency or HJI to be superior in the case of small frequencies and vice versa for higher frequencies.

As in Test 1, we continue to compare HJI to HJB for non-sinusoidal disturbances, and we first consider the piecewise constant signal given by

[TABLE]

For this case we see from Fig. 5 (ii) that in the transient phase the HJI synthesis is superior to the HJB closed-loop.

Now we modify the disturbance term in the dynamics to observe the corresponding effect close to the worst case scenario and choose $w=0.9\,w_{\gamma^{*}}+\eta sin(\omega t)$. We can see from Fig. 6 that the HJB feedback control is not able to attenuate this disturbance level. The HJB feedback controller fails to steer the trajectories towards the unstable equilibrium, and the closed-loop dynamics remain attracted to $X=1$. In the presence of the same disturbance, the HJI controller succeeds in directing the state to the desired equilibrium steers the state to the origin. The choice of signals $w_{2}$ and $w_{3}$ as in Test 1 leads to a similar type of behaviour.

Test 3: Stabilization under model uncertainty

We conclude by studying the case of state-dependent perturbations. Alternatively, we can interpret this setting as modelling parametric uncertainties in the dynamics as for instance, uncertainties in the diffusion and/or reaction coefficients. We first investigate the dynamics

[TABLE]

with $(\xi,t)\text{ in}\;\mathcal{I}\times\mathbb{R}^{+}$, $\mathcal{I}=(-1,1)$, $\omega_{2}=(-0.8,0.5)$ and $\kappa$ to be specified below. As in the previous tests, the computational domain is $\Omega=(-2,2)^{12}$.

We compare between the HJB and HJI synthesis for trajectories generated by choosing the constant perturbation $w=1$ in the above model. Note that for $w=1$ the above reaction-diffusion model corresponds to the Newell-Whitehead /Allen-Cahn equations

[TABLE]

The HJB synthesis does not incorporate the destabilizing term $X(\xi,t)$ . We study the behaviour of the closed loop under additional perturbations of the type $\eta sin(\omega t)$.

In Fig. 7 we close the loop with HJI and HJB based control, where in the three consecutive rows we increase the magnitude of the initial condition by considering the cases $\kappa=0.1$, $\kappa=0.5$, and $\kappa=1$, respectively.

For all the three cases the HJI controller stabilizes the system towards the origin and the three corresponding controls tend to [math] as time increases. The HJB feedback law, which does not anticipate the modelling error in the reaction term, fails to stabilize this system to [math].

Finally, we consider a model uncertainty for a cubic reaction term in

[TABLE]

In this case a robust control design is critical, as the cubic reaction term can generate a finite time blow-up of the dynamics. We observe that for the case $w=1$, the HJB control fails to stabilize the dynamics whilst the HJI closed-loop is asymptotically stable. This is documented in Fig. 8

Remarks on CPU times

Table 5 provides CPU-assembly and CPU-iterative times for solving the HJI equation related to Tests 1 and 2. Tests were run on a muti-core architecture 8x Intel Xeon E7-4870 with 2,4Ghz, 1 TB of RAM. MATLAB pseudoparallelization distributes the tasks among 12 workers.

Concluding remarks

We have proposed a numerical scheme for the design of robust optimal stabilizing feedback controllers, and assessed its capabilities in the context of stabilization of nonlinear parabolic PDEs. The overall technique consists of two main steps: i) the projection of the infinite-dimensional nonlinear parabolic dynamics onto a finite-dimensional relevant subspace by means of a pseudospectral collocation method and ii) the synthesis of a robust control by solving a Hamilton-Jacobi-Isaacs equation, for which we have resorted to a Galerkin-type method using a global polynomial ansatz for the value function, together with iterative techniques for the solution of the nonlinear HJI. Overall, the proposed methodology proves to be successful in achieving the fundamental features of the $\mathcal{H}_{\infty}$ control design, that is, to generate an optimal stabilization of the closed-loop with enhanced disturbance mitigation properties. We have numerically assessed that our synthesis scheme offers a more robust control law than the classical HJB/$\mathcal{H}_{2}$ feedback. Furthermore, the proposed scheme mitigates the curse of dimensionality, allowing to approximate the $\mathcal{H}_{\infty}$ synthesis for nonlinear dynamics of high dimensions. The developed framework also opens the door to the study of nonlinear feedback control of PDEs under uncertainty.

Acknowledgments

The research of this paper was supported by the ERC advanced grant 668998 (OCLOC) under the EU’s H2020 research program and the Imperial College London Research Fellowship program.