A Dual Ensemble Kalman Filter Approach to Robust Control of Nonlinear Systems: An Application to Partial Differential Equations

TL;DR

This paper introduces a dual ensemble Kalman filter-based robust control method for nonlinear PDEs, combining optimal learning control with Lyapunov-based robustification, demonstrated on heat and Burgers equations.

Contribution

It develops a novel dual EnKF approach integrated with Lyapunov redesign for robust control of nonlinear PDE discretizations.

Findings

Effective stabilization of nonlinear PDEs demonstrated in simulations.

Dual EnKF enhances robustness against disturbances.

Method applicable to various nonlinear PDE control problems.

Abstract

This paper considers the problem of data-driven robust control design for nonlinear systems, for instance, obtained when discretizing nonlinear partial differential equations (PDEs). A robust learning control approach is developed for nonlinear affine in control systems based on Lyapunov redesign technique. The robust control is developed as a sum of an optimal learning control which stabilizes the system in absence of disturbances, and an additive Lyapunov-based robustification term which handles the effects of disturbances. The dual ensemble Kalman filter (dual EnKF) algorithm is utilized in the optimal control design methodology. A simulation study is done on the heat equation and Burgers partial differential equation.

Equations

$$\dot{x}(t)=a(x(t))+b(x(t))u(t)+d(t,x),x(0)=x_0$$

$$\dot{x}(t)=Ax(t)+Bu(t)+d(t,x)$$

$$\underset{u(\cdot )}{\min }$$

$$\text{system }\text{eqref}eq:lin-control-system\text{ with zero disturbance, that is, }d\equiv 0$$

$$-\dot{P}=A^{\top }P+PA-PB{R}^{-1}{B}^{\top }P+Q,P_T=G$$

$$0=A^{\top }\bar{P}+\bar{P}A-\bar{P}B{R}^{-1}{B}^{\top }\bar{P}+Q.$$

$$u_d:=-\frac{\lambda (t,x)}{|\bar{P}x|}{B}^{†}\bar{P}x,B^{†}:=(B^{\top }B{)}^{-1}{B}^{\top },$$

$$\dot{V}(x)$$

$$\le |d|\cdot |\bar{P}x|-\lambda |\bar{P}x|$$

$$\le (|d|-\lambda )|\bar{P}x|<0.$$

$$\mathcal{H}(x,u)$$

$$=({\bar{P}}^{(N)}x{)}^{\top }\mathcal{S}(x,u)+\frac{1}{2}\mathcal{L}(x,u)$$

$$\underset{u(\cdot )}{\min }$$

$$\text{system }\text{eqref}eq:control-system\text{ with zero disturbance, that is, }d\equiv 0$$

$$\varphi (s,x)$$

$$\text{system }\text{eqref}eq:control-system\text{ with }{X}_s=x\text{ and }d\equiv 0$$

$$u_d:=-\frac{\lambda }{|\nabla V|}{b}^{†}\nabla V,b^{†}:=(b^{\top }b{)}^{-1}{b}^{\top }.$$

$$\dot{V}(x)$$

$$\displaystyle\quad+b(x(t))u_{d}(t)+d(t,x)\bigg{)}$$

$$\le (|d|\cdot |\nabla V(x)|-\lambda |\nabla V(x)|)$$

$$\le (|d|-\lambda )|\nabla V(x)|<0.$$

$$\mathcal{H}(x,u)$$

$$\frac{\partial z}{\partial t}(t,y)+\mathcal{F}(z(t,y))=\omega (t,y),$$

$$\underset{t\to \infty }{\lim }\Vert z(t,\cdot ){\Vert }_{L^2}:=\underset{t\to \infty }{\lim }{\int }_0^L|z(t,y){|}^2\mathrm{d}y=0.$$

$$\frac{\mathrm{d}{z}_p}{\mathrm{d}t}+F(z_p)=BU$$

$$\frac{\partial z}{\partial t}(t,y)-\nu \frac{{\partial }^{2}z}{\partial y^2}(t,y)=\omega (t,y).$$

$$\frac{\partial z}{\partial t}(t,y)+z(t,y)\frac{\partial z}{\partial y}(t,y)-\nu \frac{{\partial }^{2}z}{\partial y^2}(t,y)=\omega (t,y).$$

$$\dot{x}=Ax+Bu$$

$$z(0,x)$$

$$\alpha$$

$$\mathrm{d}{Y}_t^i$$

$$Y_T^i$$

$$S_t^{(N)}:=\frac{1}{N}\sum _{j=1}^N(Y_t^j-{\hat{n}}_t^{(N)})(Y_t^j-{\hat{n}}_t^{(N)}{)}^{\top }.$$

$$\mathrm{d}{Y}_t^i$$

$$+\left(\sum _{j=1}^N(Y_t^j-n_t^{(N)})(c(Y_t^j)-{\hat{c}}_t^{(N)}{)}^{\top }\right)\frac{(c(z)+{\hat{c}}^{(N)})}{N-1}$$

$$Y_T^i$$

Full text

I Introduction

In this paper, we are primarily interested in the robust control of nonlinear affine in control systems modelled as

[TABLE]

where $x(t)\in\mathbb{R}^{n}$, $u(t)\in\mathbb{R}^{m}$ and $d(t,x)\in\mathbb{R}^{l}$ for all $(t,x)$, and $d$ is regarded as an unknown disturbance bounded in norm by a known real valued function $\lambda(t,x)$. The emphasis is on systems where we do not have explicit access to $a$ and $b$, but we can access trajectories of the system, either generated by a simulator or collected from real-life experiments. For simplicity, we focus in the remaining of the paper on the case of simulated trajectories. Furthermore, we target the application of our algorithms to the control of PDEs, when their discretized model (1) is available as a simulator.

Indeed, PDE control is challenging because a PDE is by nature infinite dimensional and may have strong nonlinearities. The reduce-then-design approach is a well researched control methoology, which involves discretization of the PDE followed by dimensionality reduction, which yields a model amenable to the application of standard model based control approaches [20, 14, 2, 25, 28]. There is a recent body of work using data driven methods for building more accurate reduced order models of the PDE [3, 11, 16]. While more accurate reduced order models are beneficial to the performance of the controller, their complexity makes them computationally challenging to implement.

In this work, we focus more on the design of the controller with emphasis on obtaining the control policy with only simulator for the model, without access to the actual model parameters. See for example [6, 5, 8, 12] for some other recent efforts in obtaining data-driven controllers using a simulator of the PDE. In particular, we work on the stabilization problem, which aims to drive the PDE state to zero, by posing it as an optimal control problem. The first step is to discretize the PDE in space, to yield a high-dimensional nonlinear system of ordinary differential equations (ODEs). We consider two different cases for the simulator availability:

a linear reduced order simulator, obtained using Dynamic Mode Decomposition with control (DMDc) [24, 32] is available. 2. 2.

an extension to the case where a simulator for the full high-dimensional nonlinear discretized model of the PDE is available.

We build on previous work in this area [32] by taking a robust learning control approach. We construct the robust control as a sum of two parts: an optimal learning controller, which we expect to produce stabilization in the absence of disturbances, and an additional term that builds on the optimal control term and is based on the Lyapunov redesign approach [18] to suppress the effect caused by disturbances. The optimal control is approximated using the dual ensemble Kalman filter (dual EnKF) algorithm [15], which simulates multiple interacting copies of the system.

The ensemble Kalman filter (EnKF) [30, 26] has historically been an algorithm used for filtering, and is a key numerical method especially for high-dimensional systems, for example in weather prediction [7] (see [4] for more references). It features the design of an interacting particle system to sample from the posterior density of a filtering problem to provide a state estimate. The dual EnKF algorithm, inspired from the EnKF, computes the optimal control by converting the control problem into the problem of sampling from an appropriate probability density, through the log-transform duality between optimal control and filtering [10]. Such an approach for solving optimal control by posing it as a sampling problem is well studied (see [27, 17, 29, 13, 21]). The novelty of the dual EnKF lies in the design of an interacting particle system, inspired from the EnKF, to solve the sampling problem. The dual EnKF-based control exhibits a distinct advantage in its performance on systems of high dimension – a trait inherited from the EnKF. It performs almost two orders of magnitude faster when compared with policy gradient type approaches [23, 9] as was computationally demonstrated in [15].

The paper is organized as follows. In Section II we introduce and solve the problem for the case when the system is linear time invariant (LTI) to clarify the main ideas. Then we present the extension to the nonlinear affine in control case in Section III. Finally, in Section IV we present an application to the stabilization of PDEs, with the heat equation and Burgers’ equation as examples.

Notation: $|\cdot|$ denotes the Euclidean norm, ⊤ denotes the transpose of a matrix, and $\mathcal{I}$ refers to the identity matrix.

II Solution for linear reduced order simulator

Consider the simplified case when (1) is a linear time invariant (LTI) system:

[TABLE]

where as earlier, $x(t)\in\mathbb{R}^{n}$, $u(t)\in\mathbb{R}^{m}$ and $d(t,x)\in\mathbb{R}^{l}$ for all $(t,x)$, and $d$ is regarded as a disturbance. From here on we may suppress the $t$ argument to improve readability. We construct the robust control as a sum of a stabilizing control and an additional robustification term. We first present a recipe to obtain the stailizing control, and then to obtain the robust control, and lastly a simulator based method to obtain the two.

II-A Stabilizing control

Consider, for any arbitrary but fixed $T>0$, the optimal control problem

[TABLE]

where $C\in\mathbb{R}^{n_{1}\times n}$ (for arbitrary $n_{1}>0$), $G\in\mathbb{R}^{n\times n}$ and $R\in\mathbb{R}^{m\times m}$, and define $Q:=C^{\top}C$.

Assumption 1

We make the following assumptions about the structure of the optimal control problem (3)

1. (i)

$(A,B)$* is controllable and $(A,C)$ is observable* 2. (ii)

$R,G\succ 0$**

Then there exists a positive definite solution $\{P(t):t\in[0,T]\}$ to the differential Riccati equation (DRE) [19, Chapter 3]:

[TABLE]

which converges to $\bar{P}\succ 0$ as $T\to\infty$ and $\bar{P}$ solves the algebraic Riccati equation (ARE)

[TABLE]

Moreover, the control $u=-\bar{K}x$ with $\bar{K}:=R^{-1}B^{\top}\bar{P}$ makes the system asymptotically stable [19, Theorem 3.7].

II-B Robust control

Suppose that for the system (2) in the case of no disturbance, there exists a stabilizing control $u=-\bar{K}x$ and a strictly positive definite Lyapunov function $V(x)=\frac{1}{2}x^{\top}\bar{P}x$ with $\dot{V}=-x^{\top}\bar{P}(A-B\bar{K})x\leq 0$ along the controlled trajectories with zero disturbance. We let $\bar{K}$ be the gain obtained from the optimal control demonstrated previously and $\bar{P}$ be the solution of the ARE.

We are interested in the idea of disturbance rejection using Lyapunov redesign [18, Chapter 14]. The idea is to design a robust control $u_{d}$ such that the controller $u=-\bar{K}x+u_{d}$ makes the system (1) asymptotically stable in the presence of disturbance. To that end, we make the following assumption, and then present a design methodology for $u_{d}$.

Assumption 2

1. (i)

The rank of $B$ is $n$. 2. (ii)

There exists a known $\lambda$ such that $0\leq|d(t,x)|<\lambda(t,x)<\infty$ for each $(t,x)\in[0,\infty)\times\mathbb{R}^{n}$.

Remark 1

Assumption 2-(i) is required to motivate the theoretical derivation for the linear system, but we relax it in our implementations for PDE control.

Consider the controller $u=-\bar{K}x+u_{d}$ with

[TABLE]

where $B^{\dagger}$ denotes the Moore-Penrose pseudoinverse. The following result demonstrates the effectiveness of the robust control.

Proposition 1

Using the control $u=-\bar{K}x+u_{d}$ renders the system (1) asymptotically stable.

Proof:

We follow the method in [18, Chapter 14.2], using the same Lyapunov function $V(x)=\frac{1}{2}x^{\top}\bar{P}x$ as before. The quadratic form $V$ is strictly positive definite (by assumption) hence is a valid Lyapunov function. Taking derivative along system trajectories,

[TABLE]

For the first inequality we recall that $-x^{\top}\bar{P}(A-B\bar{K})x\leq 0$ and use Cauchy-Schwarz inequality. Next by properties of the Moore-Penrose pseudoinverse, $BB^{\dagger}\bar{P}x$ is the orthogonal projection of $\bar{P}x$ onto the column span of $B$. Hence, under Assumption 2-(i), we have $BB^{\dagger}\bar{P}x=\bar{P}x$ therefore $Bu_{d}=-\frac{\lambda}{|\bar{P}x|}\bar{P}x$. ∎

Remark 2

When implementing $u_{d}$ we add a regularizing parameter to avoid division by zero, which makes the system asymptotically stable till it enters a ball around the origin, i.e., practical stability.

II-C Data-driven (Simulator-based) implementation

In this section, we present a method to implement the control $u=-\bar{K}x+u_{d}$ with only access to a disturbance-free system simulator of (1). We use the dual ensemble Kalman filter (dual EnKF) algorithm [15] for the same. To use the algorithm, we need the following assumptions:

Assumption 3

We have access to the following:

1. (i)

Knowledge of optimal control matrices $Q,R,G$. 2. (ii)

Simulator of the dynamical system with no disturbance, that is, we have access to function evaluations of $\mathcal{S}(x,u):=Ax+Bu$. Moreover, we assume we can run the simulator backward in time, that is, to find a trajectory of the system by specifying the terminal condition.

Remark 3

With access to a perfect simulator, one may exactly find the model matrices $A$ and $B$ in $n+m$ evaluations of the simulator (set $u=0$ and evaluate the simulator at the basis vectors of $\mathbb{R}^{n}$ to find $A$ and equivalently for $B$). However, the utility of the simulator is revealed when we consider the nonlinear case, especially for high-dimensional systems like partial differential equations, where estimating the state dynamics in this manner is not possible.

The emphasis is on obtaining the controller in a model-free way. The following three steps are done (which are elaborated upon after listing them):

Find an approximation to $\bar{P}$, the solution of the ARE, using the dual EnKF algorithm [15]. See Appendix A-A for details 2. 2.

Find an approximation $\bar{u}^{(N)}$ for $\bar{u}:=-\bar{K}x$ using [15, Algorithm 2] (recalled in Appendix A-C) 3. 3.

Find an approximation $u_{d}^{(N)}$ for $u_{d}$ using Algorithm 1

Step 1: Using the simulator, we find $\bar{P}^{(N)}$, an approximation to the solution of the ARE $\bar{P}$ by running the dual EnKF algorithm of [15], which simulates $N$ copies of the dynamical system along with a mean-field coupling term to approximate the solution of the ARE. The algorithm is provided in Appendix A-A.

Step 2: To help evaluate the optimal control in a model-free way, define the Hamiltonian,

[TABLE]

The Hamiltonian is constructed so that it can be evaluated using function calls of $\mathcal{S}$. For a fixed $x$, it is a quadratic function of $u$ with the unique minima at the optimal control. Therefore, $\bar{u}=\operatorname*{arg\,min}_{u}\mathcal{H}(x,u)$ and the minimization can be carried out using gradient estimation as shown in [15, Algorithm 2] (recalled in Appendix A-C) or using zero-order methods, such as [1].

Step 3: Similarly, to find $u_{d}$ we solve $u_{d}=-\lambda\operatorname*{arg\,min}_{u}|Bu-\frac{\bar{P}x}{|\bar{P}x|}|$ in Algorithm 1. If $B$ is known, one may directly use the Moore-Penrose pseudo inverse. If $B$ is unknown, one may use zero order optimization methods [1] where $Bu=\mathcal{S}(0,u)$ can be obtained using only access to simulator or one may estimate $B$ as specified in Remark 3 and use the pseudo inverse. Using zero order optimization is preferred over estimating $B$ in cases when the simulator is very high dimensional.

III Extension for nonlinear discretized model of the PDE

Let us consider now the main result of this paper, dealing with the case of the nonlinear affine in control system (1). We will extend the robust control design methodology presented in Section II to (1). Similar to the previous section, we will obtain the robust control as a sum of a stabilizing control and a robustification term, and then present a simulator based methodology to obtain both terms.

III-A Stabilizing control

In an effort to find a stabilizing control $\bar{u}$, we consider, for any arbitrary but fixed $T>0$, the optimal control problem

[TABLE]

where $c,\mathcal{G}$ are non-negative real valued function, and $R\in\mathbb{R}^{m\times m}$ is symmetric and strictly positive definite. The value function for the problem is defined as the cost-to-go,

[TABLE]

and it satisfies the Hamilton Jacobi Bellman partial differential equation. The optimal control is computed as $\bar{u}(t,x)=-R^{-1}b(x)^{\top}\nabla\phi(t,x)$ [22].

III-B Robust control

Assumption 4

1. (i)

There exists a control law $\bar{u}$ which makes (1) with zero disturbance asymptotically stable. Moreover, there exists a strictly positive definite Lyapunov function $V$ with $\dot{V}<0$ along the controlled trajectories with zero disturbance. 2. (ii)

The rank of $b(x)$ is $n$ for all $x\in\mathbb{R}^{n}$. 3. (iii)

There exists a known $\lambda$ such that $0\leq|d(t,x)|<\lambda(t,x)<\infty$ for each $(t,x)\in[0,\infty)\times\mathbb{R}^{n}$.

Remark 4

Assumption 4-(ii) is required for the theoretical result, but will be relaxed in the PDE control implementation.

Similar to the linear case, consider the controller $u=\bar{u}+u_{d}$ with

[TABLE]

Proposition 2

Using the control $u=\bar{u}+u_{d}$ renders the system (1) asymptotically stable.

Proof:

We follow the method in [18, Chapter 14.2]. Taking the derivative of $V$ along system trajectories,

[TABLE]

For the first equality we recall Assumption 4 and Cauchy-Schwarz inequality. Moreover, under Assumption 4-(ii), $bu_{d}=-\frac{\lambda}{|\nabla V(x)|}\nabla V(x)$. ∎

III-C Data-driven (Simulator-based) implementation

To approximate $\bar{u}$ and $u_{d}$ we use an approach similar in spirit to Section II. First, we use the nonlinear dual EnKF algorithm [15] to approximate the gradient of value function $\nabla\phi^{(N)}(x)$. Implementation details can be found in Appendix A-B. Then we define the nonlinear counterpart of the Hamiltonian

[TABLE]

where the simulator is now nonlinear, that is, $\mathcal{S}(x,u)=a(x)+b(x)u$. The Hamiltonian can again be evaluated using function calls of the simulator. Moreover, it is quadratic in the control, and can be minimized using [15, Algorithm 2] (recalled in Appendix A-C), or zero order optimization approaches [1]. Similar to the linear control case, to find $u_{d}$ we use $u_{d}=-\lambda\operatorname*{arg\,min}_{u}|b(x)u-\frac{\nabla V(x)}{|\nabla V(x)|}|$ as given in Algorithm 1 (replacing $\bar{P}^{(N)}x$ by $\nabla V^{(N)}(x)$) where if $b$ is not known, the optimization can be done by zero-order approaches [1] or by estimating $b$ similar to Remark 3.

IV Application to forced nonlinear PDEs

We consider PDEs of the form

[TABLE]

where $\mathcal{F}$ is a differential operator that specifies the structure of the PDE, $z$ is the state of the PDE and $\omega$ is the external input. The functions $z,\omega:\mathbb{R}_{+}\times[0,L]\to\mathbb{R}$. Mathematically, the goal of stabilization is to make

[TABLE]

For numerical implementation, we consider a time interval of $[0,T]$ and discretize the PDE in space on a grid of $p$ uniformly-spaced points in $[0,L]$ to obtain $z_{p}(t)\in\mathbb{R}^{p}$ for each $t\geq 0$, so that the PDE is reduced to a set of $p$ ODEs. We choose $m$ basis functions $\{\chi_{j}\}_{j=1}^{m}$ to discretize the control as $\omega(x,t)=\sum_{j=1}^{m}\chi_{j}(x)U_{j}(t)$ where $\chi_{j}:[0,L]\to\mathbb{R}$ is the indicator function of $[\frac{j-1}{m},\frac{j}{m}]$ and $U_{j}:[0,T]\to\mathbb{R}$. Then $U:=(U_{1},U_{2},\ldots,U_{m})$ is interpreted as the control. The discretized PDE has nonlinear affine in control form,

[TABLE]

where $F$ approximates the derivatives and $B$ obtained by discretizing $\{\chi_{j}\}_{j=1}^{m}$. To study the effect of disturbances, we let $U(t)=u(t)+d(t)$ where $u(t)$ is the control action applied to the system and $Bd(t)$ is the effective disturbance acting on the system. We consider here disturbances added directly to the control, which is a special case of the theory presented earlier. In the simulations, we consider only time varying disturbances, so we denote disturbance as $d(t)$. We use the Controlgym library in Python [31] for numerical implementation. Now we first recall the heat equation and Burgers equations, give some information about the implementation of both PDEs, then go on to discuss the simulation results.

IV-A Heat equation

The heat equation is given by

[TABLE]

The value of $\nu=0.002$ is used. Since the PDE is linear, upon discretization, (8) reduces to an LTI system, hence we use the robust control method design discussed in Section II to implement the robust control. The dimension of the control used is $m=10$, while the dimension of the state is $n=100$, thus relaxing Assumption 2-(i) in the implementation. Other simulations details can be obtained in Appendix B-A.

IV-B Burgers’ equation

Recall the Burgers’ equation,

[TABLE]

We present two sets of results to demonstrate the performance of the robust control algorithm developed in the paper when applied to the Burgers equation. In the first set of results, we calculate the robust control using a simulator of the linear reduced order model of the PDE, and in the second set of results, the robust control is calculated using a simulator of the full nonlinear discretized PDE. In both results, the obtained control is tested by applying it on the full nonlinear discretized PDE. We study both sets of results for two values of viscosity, $\nu=0.02,0.002$.

Control design using linear reduced model: We perform a model reduction on the Burgers’ equation to obtain a linear system using the Dynamic Mode Decomposition for control (DMDc) algorithm [32]. The DMDc yields a projection matrix $\Phi\in\mathbb{R}^{n\times p}$ with orthogonal rows, where $n$ is the dimension of the reduced state (usually $n\ll p$), and an LTI system

[TABLE]

where $x\in\mathbb{R}^{n}$ is the reduced state and satisfies the relation $x=\Phi z_{p}$, and $u\in\mathbb{R}^{m}$ where as described, $n$ and $m$ are chosen by the user. The discretized PDE state can be reconstructed as $z_{p}=\Phi^{\top}x$. DMDc yields a discrete-time system, but we transform it into its continuous-time equivalent, since we design control in the continuous time domain.

Using a simulator for the obtained linear system (9), we get an approximation to the optimal stabilizing LQ controller, with the linear EnKF methodology described in Section II. That is, we choose optimal control weights $Q,R,G$ (as described in (3)) and find $\bar{P}^{(N)}$ for the LTI system (9). The stabilizing control is calculated as $u=-R^{-1}B^{\top}\bar{P}x$ where $x=\Phi z_{p}$. We use the robust control methodology described in Section II to design $u_{d}$.

Control design using full nonlinear model: Using a simulator for the full nonlinear PDE (8), we obtain a stabilizing control $\bar{u}$ and robust controller $u_{d}$ with the nonlinear dual EnKF methodology discussed in Section III. We choose $c(x):=|x|^{2}$ and $\mathcal{G}(x)=|x|^{2}$ in (5). The dimension of the state used is $p=128$ and the dimension of the control is $m=10$, thus relaxing Assumption 4-(ii) in the implementation.

IV-C Discussion of results

We study the effect of two types of disturbances, sinusoidal with $d(t):=d_{0}\sin(t)$ and constant $d(t):=d_{0}$. The $\lambda$ function is chosen as a constant function. The obtained control is applied to the full discretized nonlinear PDE (8). Both PDEs are initialized for 100 randomly sampled iid initial conditions (inspired from the previous work [32]):

[TABLE]

The 100 trajectories are simulated for the case of zero control (that means $\bar{u}=u_{d}=0$), and then with the robust control given by various values of $\lambda$ (the stabilizing control $\bar{u}$ is obtained from the dual EnKF). These simulations are repeated for various values of $\lambda$, $d_{0}$, both types of disturbances mentioned, and the values of $\nu$ mentioned.

Recall that the control objective (7) in finite time is to drive $\|z(T,\cdot)\|_{L^{2}}$ as close to zero as possible. To illustrate the performance of the controllers we make two plots for both PDEs. The first plot is a heat map depicting the mean value (over the 100 simulations) of $\frac{\|z(T,\cdot)\|_{L^{2}}}{\|z(0,\cdot)\|_{L^{2}}}$ – the ratio is plotted to highlight the order of magnitude by which $\|z(T,\cdot)\|_{L^{2}}$ has reduced relative to its initial value $\|z(0,\cdot)\|_{L^{2}}$. Moreover, in a second plot, we plot the mean and variance (over the 100 simulations) of $\|z(t,\cdot)\|_{L^{2}}$ as a function of $t$ for $d_{0}=0.1$ and $\lambda=0.2$. The trajectory is generated using three different control policies – “uncontrolled” trajectories are when the control is zero, “optimal controlled” trajectories have $\lambda=0$ (that is, only the optimal stabilizing control is applied and robust control is zero) and “robust controlled” trajectories are with $\lambda=0.2$. The plots for heat equation, Burgers with DMDc and Burgers with full nonlinear model are found in Figures 1, 2, and 3, respectively. Other simulations details can be obtained in Appendix B-B.

For both PDEs, we observe that in the presence of disturbances, the robust control works better than the optimal control (control with $\lambda=0$) in stabilizing the PDE. For the Burgers’ equation in particular, the trajectories with robust control exhibit an order of magnitude lower value of $\|z(T,\cdot)\|_{L^{2}}$ compared to those with only the optimal control. Additionally, for the Burgers’ equation, we also observe that the optimal control obtained using the full nonlinear model is significantly more effective than using the reduced-order DMDc model – the controlled state settles much faster and closer to zero with the former than the latter. However, while the transient performance of the robust control is much better using the full nonlinear model, the settling performance is similar using both DMDc and the full nonlinear model. Thus the robust control term also compensates for model mismatch between DMDc and the full nonlinear model. The conclusions presented are true for all values of $\lambda$ chosen, both types of disturbances, and both values of viscosity.

V Future work

Some avenues for future work are: (i) to extend this approach to different PDEs such as the Allen-Cahn equation or the Korteweg deVries equation, (ii) to design a controller using information from sensors, thus turning a full state feedback problem into a partially observed problem (iii) to extend to the case when $d$ is Gaussian white noise.

Appendix A Dual EnKF algorithm

In this algorithm, we simulate over the time horizon $[0,T]$ an ensemble of $N$ particles $\{Y_{t}^{i}:1\leq i\leq N,\quad 0\leq t\leq T\}$ where the evolution equation for each particle is given by an Ito stochastic differential equation (SDE).

A-A Linear system

The SDE [15, Section 2] is given by

[TABLE]

${\eta}^{i}$ are an i.i.d Brownian motions with covariance $R^{-1}$ and $\hat{n}^{(N)}_{t}:=\left(\frac{1}{N}\sum_{j=1}^{N}Y_{t}^{j}\right)$, and $L_{t}^{(N)}:=S_{t}^{(N)}C^{\top}$ with

[TABLE]

Finally, $\bar{P}^{(N)}:=(S^{(N)}_{0})^{-1}$.

A-B Nonlinear system

The SDE [15, Section 3] for the particle system is

[TABLE]

and $\hat{c}^{(N)}_{t}:=N^{-1}\sum_{i=1}^{N}c(Y^{i}_{t})$. Finally, $\nabla\phi^{(N)}(x):=(S^{(N)}_{0})^{-1}x$, where $S^{(N)}_{t}$ is the empirical covariance of $\{Y_{t}^{i}\}$, defined same as in the linear case.

A-C *Algorithm for Hamiltonian minimization *

Algorithm 2 minimizes Hamiltonian and calculates $\bar{u}$ for Section II. To minimize Hamiltonian in Section III, in Algorithm 2 make the following two changes: replace $\bar{P}^{(N)}$ by $\nabla\phi^{(N)}$ and use the Hamiltonian defined in (6).

Appendix B Simulation details

B-A Heat equation

Simulation parameters are as follows. The simulation time $T=0.1$ with simulation time step = 0.001. The number of states of the discretized PDE (8) is $p=100$. The number of control basis functions is $m=8$, and they are $\chi_{j}$ is the indicator functions of $[\frac{i}{10},\frac{(i+1)}{10}]$. The regularization parameter for robust control is $r=0.002$. The matrices $Q=\mathcal{I}$, $G=\mathcal{I}$, $R=\mathcal{I}$ for both, the full nonlinear control and control using DMDc model. The number of dual EnKF particles is chosen as $N=10000$. The controlgym library [31] is used for PDE simulation.

B-B Burgers equation

Simulation parameters are as follows. The simulation time $T=3$ with simulation time step = 0.001. The number of states of the discretized PDE (8) is $p=128$. The number of control basis functions is $m=10$, and they are $\chi_{j}$ is the indicator functions of $[\frac{i}{10},\frac{(i+1)}{10}]$. The regularization parameter for robust control is $r=0.002$. The number of states in the DMDc model (9) is $n=10$. The matrices $Q=\mathcal{I}$, $G=\mathcal{I}$, $R=0.1\mathcal{I}$ for both, the full nonlinear control and control using DMDc model. The number of dual EnKF particles is chosen as $N=1000$. The controlgym library [31] is used for PDE simulation.