Point-actuated feedback control of multidimensional interfaces

We consider the application of feedback control strategies with point actuators to stabilise desired interface shapes. We take a multidimensional Kuramoto--Sivashinsky equation as a test case; this equation arises in the study of thin liquid films, exhibiting a wide range of dynamics in different parameter regimes, including unbounded growth and full spatiotemporal chaos. In the case of limited observability, we utilise a proportional control strategy where forcing at a point depends only on the local observation. We find that point-actuated controls may inhibit unbounded growth of a solution, if they are sufficient in number and in strength, and can exponentially stabilise the desired state. We investigate actuator arrangements, and find that the equidistant case is optimal, with heavy penalties for poorly arranged actuators. We additionally consider the problem of synchronising two chaotic solutions using proportional controls. In the case when the full interface is observable, we construct feedback gain matrices using the linearised dynamics. Such controls improve on the proportional case, and are applied to stabilise non-trivial steady and travelling wave solutions.


Introduction
The study of evolving interfaces is at the core of many areas of applied mathematics, ranging from crystal growth problems in chemistry (Langer, 1980;Kobayashi, 1993;Pimpinelli & Villain, 1999) to the study of flame front propagation in combustion theory (Michelson & Sivashinsky, 1977;Sivashinsky, 1977Sivashinsky, , 1980. In fluid dynamics, there is interest in the evolution of interfaces separating immiscible phases, e.g. in thin film flows on inclined planes (Homsy, 1974;Michelson & Sivashinsky, 1980) or jet and droplet breakup problems (Papageorgiou, 1995;Eggers, 1997). Starting with complicated multiphase systems comprising numerous coupled equations, it may be possible to employ reducedorder modelling techniques to isolate the evolution of the interface alone; this is particularly desirable when information about the bulk dynamics (away from the interface) is not of interest. It is often challenging to isolate the interfacial dynamics while still retaining all desired physical effects, and in many cases it is found that the obtained low-dimensional models only replicate the true dynamics well in restricted parameter regimes. This pitfall is balanced by the relative simplicity of the interface evolution equations along with a large decrease in computational complexity for numerical simulations.
It may be useful to control interfacial dynamics in order to optimize an industrial process. Controlling jet and droplet breakup is of interest for both inkjet printing (Martin et al., 2008) and the fabrication of uniform emulsions (Christopher & Anna, 2007). Miller et al. (2010) developed a microfluidic device incorporating a feedback controller, which enabled the production of a tunable emulsion. Cooling and coating processes involving thin film flows also arise in microfluidic applications; for cooling, waviness of fluid interfaces is desirable as it improves heat transfer (Miyara, 1999;Serifi et al., 2004), whereas for coating, a flat interface is sought. These thin film flows may be controlled with air/liquid actuators, external electric or magnetic fields, surfactants and substrate coating or topography, among others. For crystal growth processes, controls may be introduced in the form of heat sources to modify the rate of growth (as utilized in experiments by Kokh et al., 2005). The stochastic Kardar-Parisi-Zhang equation is used as a continuum model for the interfacial dynamics of crystal growth via the artificial deposition process known as molecular-beam epitaxy (Pimpinelli & Villain, 1999). Controlled changes in temperature influence the strength of the nonlinearity in the Kardar-Parisi-Zhang equation; this is reflected in a variation of the lateral growth of the interface, influencing its statistical properties. Timedelayed feedback control of this system was considered by Block et al. (2007). For most real-world problems, interfaces are described in terms of two spatial variables. For problems where variations in one direction are negligible, such as the growth of a flat crystal, simplification of the interface problem to one spatial dimension may be viable, but it is important that the controls utilized preserve this property. It may be the case that such a simplification overlooks important instabilities or mechanisms that are only observed from the full 3D formulation of the original multiphase problem, e.g. Rayleigh-Taylor instabilities or electrostatically induced instabilities in liquid films (Tomlin et al., 2017). This paper investigates feedback control strategies for multidimensional evolving interfaces, taking as a test case the 2D Kuramoto-Sivashinsky equation (KSE) η t + ηη x + (1 − κ)η xx − κη yy + Δ 2 η = ζ . (1.1) This nonlinear PDE in two space dimensions, x, y and time, t, may be derived to describe the weakly nonlinear evolution of small-amplitude, long-wave perturbations of gravity-driven thin liquid films on flat substrates. The surface η(x, y, t) represents a perturbation of the exact flat film solution, and κ is a parameter that we discuss below. The control ζ(x, y, t) arises physically as same-fluid blowing and suction through apertures on the substrate surface. Through reduced-order modelling, boundary controls at the level of the full physical system, comprising the injection and extraction of mass from the bulk flow, are manifest as internal controls acting on the fluid-gas interface. In the present work, controls are applied using a finite set of point-actuators, which are commonly the most physically realizable in applications. The point-actuation is represented mathematically with Dirac delta functions; smoothed alternatives have been utilized throughout the literature. Point-actuators have also been employed in the dynamic control of 2D deformable membrane mirrors for adaptive optics (see the review of Ruppel, 2012, and the references therein). Equation (1.1) is an extension to two spatial dimensions of the 1D KSE with controls, which, for ζ = 0, is the paradigmatic model for a class of active-dissipative nonlinear evolution equations. Usually, (1.2) is supplemented with periodic boundary conditions on the interval [0, L]. As L is increased beyond 2π (at which point the first Fourier mode destabilizes), the dynamics cascades through steady and travelling wave, time-periodic and quasi-periodic attractors, to full spatiotemporal chaos. The latter is characterized by the exponentially fast separation of solutions with neighbouring initial conditions under the time evolution of the PDE. There is a wealth of literature exploring the nature of the chaotic attractor for (1.2)-see Kevrekidis et al. (1990); Papageorgiou & Smyrlis (1991); Smyrlis & Papageorgiou (1996). There are a number of studies considering the control of (1.2) with point actuators: Christofides (1998); Armaou & Christofides (2000); Lunasin & Titi (2017); Gomes et al. (2017);-the latter work motivates much of the present paper. We supplement the 2D KSE (1.1) with periodic boundary conditions on the rectangular domain Q = [0, L 1 ] × [0, L 2 ]. Different dynamical regimes are found by varying κ. For κ ≥ 1, the flat film solution is stable-we do not consider this case in the present paper. For κ < 0 (which physically corresponds to a film on the underside of a substrate), solutions can grow without bound due to an unsaturated linear (Rayleigh-Taylor) instability in the transverse modes. For 0 ≤ κ < 1, chaotic dynamics are observed as investigated by Akrivis et al. (2016) and Tomlin et al. (2018) in the case of κ = 0 (corresponding to a vertical film set-up). The 2D KSE (1.1), albeit a deterministic equation, provides a natural and challenging test case for the control strategies considered in this work. Efficient and convergent numerical schemes allow us to study the control of (1.1) on large domains with many unstable modes and solutions exhibiting full spatiotemporal chaos; the majority of existing numerical studies consider parameter regimes for model problems not far from the onset of instability (one or two unstable modes).
The two feedback (closed-loop) control strategies we consider are proportional control and feedback control with full state observations ('full-state feedback control'). These strategies are polar opposites in terms of the assumed observability of the interface and knowledge of the governing dynamics. Unsurprisingly, we find that more information (observations/knowledge of governing equation) results in a much improved control performance, but both methods are effective. A primary interest is in finding the most favourable arrangement of the point-actuators, which is independent of initial conditions. We emphasize that we do not study the (open-loop) optimal control problem for (1.1), whereby the locations and time-varying strengths of the actuators, and the resulting solution, minimize a given cost functional; this was studied for the 1D KSE (1.2) by Gomes et al. (2017). For such a direction, optimization over short time intervals gives actuator positions that are heavily dependent on initial conditions, and optimizing over large time intervals is numerically expensive due to the adjoint methods involved. The optimal control problem for (1.1) was studied in Tomlin et al. (2019) under the assumption that any forcing profile can be realized.
Proportional controls are the most simplistic form of feedback control. In this study, each actuator is co-located with a point-observer, and applies a forcing proportional to the difference between the current observation of the interface and the corresponding observation of a chosen desired state (e.g. a steady or travelling wave solution). No information about the dynamical system is required or even beneficial since the actuation at a particular point does not utilize observations from other spatial locations. This study can alternatively be viewed as continuous-time data assimilation for the 2D KSE, as from pointobservations of a given solution we are able to reconstruct it in full (Azouani & Titi, 2014;Lunasin & Titi, 2017). We note that upstream (phase-shifted) observers were found to improve the control of thin film models in Thompson et al. (2016a). We perform a number of numerical experiments, investigating the control of exponentially growing or chaotic interfaces and different actuator arrangements. For the latter, we perform proportional control with a range of actuator arrangements, and compare them based on the decay of the solution to the desired state, finding that an equidistant spacing is the best choice.
Full-state feedback control is a more advanced closed-loop strategy, which assumes both the ability to observe the full interface and knowledge of the governing system. The linearization of the dynamics is used to create a function (the feedback gain matrix) that maps the observation of the interfacial state at an instant in time to the controls required to provide linear stability of the desired state. We present two methodologies for full-state feedback control to non-trivial interface shapes based on studies of the 1D KSE (1.2) by Al Jamal & Morris (2018) and Gomes et al. (2017). The former ensures exponential stabilization through a rigorous analytical result, whereas the latter is much more feasible numerically if the desired state is non-trivial. The capabilities of full-state feedback control in this multidimensional setting are tested with comparisons against the proportional control results. There is a long list of extensions and hybridizations of these control strategies, which we do not consider in the current work, such as dynamical observers and time-delayed/phase-shifted observers. However, the present study considers methods that are readily extendable to more complicated systems and experiments.
The current paper is organized as follows: Section 2 introduces the control problem for the 2D KSE (1.1) following a brief derivation. We then provide the analytical setting of the equation, and discuss the numerical methods utilized and the various actuator arrangements considered. Sections 3 and 4 contain the studies of the proportional and full-state feedback control strategies, respectively. The concluding remarks are given in Section 5.

Multidimensional KSE with point-actuated controls
We begin with an expedited derivation of the 2D KSE (1.1). Consider a thin liquid film on a substrate inclined some non-zero angle θ to the horizontal-see the schematic in Fig. 1 for an overlying liquid film, i.e. θ ∈ (0, π/2). As shown in the figure, the x-coordinate is in the streamise direction, y in the transverse direction (cross-flow and parallel to the substrate) and z is in the direction perpendicular to the substrate surface. The substrate and interface are located at z = 0 and z = h(x, y, t), respectively. We allow hanging liquid films, for which θ ∈ (π/2, π) and the z-component of gravity has a destabilizing effect on the interface. The bulk fluid flow in 0 < z < h(x, y, t) is governed by the Navier-Stokes equations, which we give following a long-wave (lubrication) rescaling in a non-dimensional form as Here, the velocity field is (u, v, w), the pressure is p and is a small parameter that measures the ratio of the average film thickness to typical interfacial wavelengths. In our non-dimensional variables, the undisturbed fluid film has unit thickness, and the Reynolds number measuring the ratio of inertial to viscous forces, Re, is assumed to be O(1) and is as defined in Tomlin et al. (2019). We control the system with fluid-mass actuators at the substrate surface, z = 0. These are modelled by altering the impermeability condition. Specifically, we have the usual no-slip condition and a modified impermeability condition at the substrate, For the flow past a fluid-filled porous substrate, Beavers & Joseph (1967) proposed a tangential slip boundary condition where the slip length scales with the typical pore size. As discussed by Thompson et al. (2016b), condition (2.2) is valid for actuator apertures that are small in comparison to the film thickness (so that the Beavers-Joseph slip is negligible). The spread of actuator apertures relative to the long interfacial waves is also important; in the present work, actuators are spaced comparatively to the interfacial deformations (as indicated in Fig. 1), and so the function f is given as a sum of localized forcing profiles. For the optimal control study in Tomlin et al. (2019), it is assumed that the actuators are closer together, so that any forcing function may be realized.
To leading order, the interfacial stress balances simplify under the lubrication scaling to where Δ ≡ ∂ 2 x + ∂ 2 y , and C is a capillary number that measures the ratio of viscous to surface tension forces (for the definition of C, see Tomlin et al., 2019). We assume that C ∼ 2 so that surface tension effects appear in the leading order pressure. The full formulation of the film flow problem is completed with the kinematic condition, In order to obtain a reduced-order system governing the controlled interface dynamics alone, we apply a systematic asymptotics procedure-this was carried out in the 2D case by Thompson et al. (2016b), and in 3D (giving 2D interfacial dynamics) by Tomlin et al. (2017) with an electric field'instead of fluid-mass blowing and suction controls. Combining (2.1c) with the above normal stress balance, we obtain the fluid pressure in the liquid layer (2.5) The flow field may be obtained from (2.1a),(2.1b),(2.1d) with the substrate and interfacial boundary conditions as Substituting the leading order terms into the kinematic condition (2.4) gives the degenerate evolution equation h t + 2h 2 h x = f . This equation is regularized at the next order, for which we must obtain the velocity corrections (u 1 , v 1 , w 1 )-these corrections involve inertial terms, as well as gravity and surface tension contributions through variations of the pressure at the interface. The resulting Benney evolution equation, with errors of O( 2 ), is where ∇ = (∂ x , ∂ y ) with the pressure given in (2.5). Ignoring transverse variations, this model reduces to the Benney equation with controls considered by Thompson et al. (2016a,b). Notice that f is not only on the right hand side of (2.7) but also appears in a nonlinear term in the streamwise direction at O( ). In the derivation of the above equation it is assumed that f t = O(1), otherwise a term involving f t is promoted from the O( 2 ) error to the O( ) terms of (2.7)-this restriction still permits feedback controls that vary on the same timescale as the interfacial dynamics. The weakly nonlinear dynamics is obtained by substituting h = 1 + η and f = 4 ζ into (2.7), resulting in equation (7) in Tomlin et al. (2019). By performing the change of variables given in equation (8) in Tomlin et al. (2019), but omitting the Galilean component of the x-rescaling, we arrive at For κ > 0 we have overlying film flows, a vertical film flow for κ = 0 and hanging flows when κ < 0. The value κ = 1 corresponds to taking the critical Reynolds number, Re c = 5 cot θ/4, for an overlying flow. As κ is decreased from 1, the flat film solution first becomes unstable to long waves (Re > Re c ); for κ > 1 we have subcritical Reynolds number flows (Re < Re c ), where all initial conditions decay to zero uniformly in the absence of controls. The advection parameter χ measures the speed of waves travelling downstream in the static 'lab' frame of reference. Recalling that C ∼ 2 , we have χ ∼ −1 , which is large. If the set of actuators is invariant to shifts in the streamwise direction (not true for the present study), it is viable to employ the Galilean transformation x → x + χ t (substitute x = x + χ t into (2.8) and drop the bar) to remove the advective term, resulting in equation (1.1). The advection term provides no additional complexity to the problem, and the fluid dynamicist's perspective alone merits the study of its influence on the controlled dynamics. We found that increasing χ improved the controlled dynamics (i.e. increased decay rates) until the effect plateaus; this can be rationalized by the fact that, for large χ , a single actuator influences the interface along lines of constant y. The case of χ = 0 is seemingly the most challenging, and so we continue without the linear advective term; it is noteworthy that previous studies of the 1D KSE with point-actuated controls also ignored advection effects. We conjecture that the reintroduction of advection to the problem would allow actuators to be located more sparsely in the streamwise direction than for χ = 0. We supplement (1.1) with periodic boundary conditions on the rectangle Q = [0, L 1 ] × [0, L 2 ], in which case the spectrum of solutions is restricted to wavenumbersk = (k 1 ,k 2 ) wherẽ (2.10) Denoting x = (x, y), we may write η and ζ in terms of Fourier series as (2.11) Equation (1.1) may thus be replaced by an equivalent infinite system of ODEs for the Fourier coefficients η k and ζ k , where η −k and ζ −k are the complex conjugates of η k and ζ k , respectively, since both the solution and control are real-valued. The effect of the parameter κ can be seen more clearly from this system of ordinary differential equations (ODEs), and the different dynamical regimes are given in Table 1. The unbounded growth for κ < 0 is due to a linear (Rayleigh-Taylor) instability in transverse modes (assuming L 2 is sufficiently large), which is not saturated by the nonlinearity-this can be seen from the ODE system (2.12) where the nonlinear (summation) term has no contribution ifk 1 = 0. Further details of the dynamical regimes of equations (1.1),(2.12) can be found in Tomlin et al. (2019); in particular, the vertical falling film case has been considered extensively both in analytical and numerical studies (Akrivis et al., 2016;Nepomnyashchy, 1974a,b;Pinto, 1999Pinto, , 2001Tomlin et al., 2018). We focus our attention on stabilizing the dynamics in the unstable regimes, κ < 1. In contrast, for κ ≥ 1, the relevant control problem is in destabilizing the flat interface. It is evident from (2.12) that the solution mean, η 0 , is preserved by the dynamics if the controls are zero mean (i.e. ζ 0 (t) ≡ 0), otherwise it will drift. We denote the desired state by η; this will usually be the trivial zero solution, but we also consider travelling waves and interfaces exhibiting full spatiotemporal chaos. Although not always necessary, we assume throughout that η is an exact solution (stable or unstable) of the uncontrolled equation; the full-state feedback control methodology requires this assumption, and furthermore, convergence to non-solutions usually requires a dense set of actuators. We consider point-actuated controls, with N ctrl actuators located at {x j } N ctrl j=1 ⊂ Q. The actuator function and time-dependent control corresponding to the location x j are denoted by b j (x) and φ j (t), respectively. In the thin film setting, these controls correspond to blowing or suction applied through holes in the substrate surface as depicted in Fig. 1. In general, the control may be written as (2.13) The Fourier coefficients of ζ may thus be expressed as a linear combination of the b j k depending on the controls φ j . Many studies of point-actuated controls employ smoothed actuators that are centred at the given actuator locations, e.g. Gaussians, or functions that can be obtained from rescalings and translations of exp((cos x − 1)w −2 ) as used by Thompson et al. (2016a) in their study of the control of long-wave (Benney, weighted-residual) models related to the KSE in 1D. Note that the latter approximation converges to a 2π -periodic extension of the usual Dirac delta δ(x) as w → 0. We do not make such an analytic approximation, and consider the Dirac delta actuators in two space dimensions (2.14) The division by |Q| = L 1 L 2 in this expression ensures that the corresponding distribution is dimensionless, so that for v ∈ H 2 , the Sobolev space of periodic functions with spatial derivatives up to second order in L 2 (note that δ is in the dual space H −2 and H 2 ⊂ C 0 in 2D). It follows (by taking v = 1) that the spatial integral of b j over Q is well-defined and is unity; additionally, any truncations of b j that include the zero mode have unit spatial integral. For numerical experiments, we truncate the Fourier series of the control at the same refinement as for the solution. In this way, the two limits of improving the resolutions of the solution and control are taken together. We note that the forcing is highly singular, and we require many Fourier modes for good spatial convergence in numerical simulations. We define the spatial L 2 -inner product and corresponding norm as where v k and w k are the Fourier coefficients of v and w, respectively. The scaling in (2.16) gives meaning to the L 2 -norm as a measure of interfacial energy density. The concerns of the current work are primarily numerical; however, we make some brief remarks on the analytical aspects of the problem. In 2D, we have δ ∈ H −2 , and thus for φ j ∈ L 2 (0, T) we conclude that ζ ∈ L 2 (0, T; H −2 ), i.e. the H −2 -norm of ζ is L 2 -in-time. This is the minimal regularity of forcing needed for existence and uniqueness of solutions to (1.1) in L 2 (0, T; H 2 ) ∩ C 0 ([0, T]; L 2 ), assuming an initial condition η 0 ∈ L 2 . This is due to standard results for parabolic equations-see Ch. III §3 of Temam (2001) or Ch. 9.4 of Robinson (2001) for the corresponding results for the Navier-Stokes equations, with extension to KS-type equations following  Temam (1997). We have also considered the possibility of obtaining analytical estimates for the number of controls and control strength sufficient for exponential stabilization of (1.1) in the proportional control case, as done for similar problems in 1D by Azouani & Titi (2014) and Lunasin & Titi (2017). We observed for the 1D KSE (1.2) that the analytical result appears far from optimal, and the extension to multiple spatial dimensions is not trivial.

Numerical methods and data analysis
For our numerical study of (1.1) on Q-periodic domains, we utilize backwards differentiation formula (BDF) methods for the time discretization and spectral methods in space. The BDF schemes belong to the family of implicit-explicit methods constructed by Akrivis & Crouzeix (2004) for a class of nonlinear parabolic equations-see the appendix of Akrivis et al. (2009) for the first-to sixth-order schemes. They considered evolution equations of the form where A is a positive definite, self-adjoint linear operator, and B is a nonlinear operator that satisfies a local Lipschitz condition. It was shown that these numerical schemes are efficient, convergent and unconditionally stable. The applicability of these schemes for (1.1) without controls was shown in Akrivis & Smyrlis (2011) and a convergence study was performed in Akrivis et al. (2016) for the choice of κ = 0. It was observed that the BDF schemes of order three to six (which are not unconditionally stable) achieved convergence to machine accuracy as soon as the time-step was small enough for the stability of the scheme. These BDF methods were also utilized in Tomlin et al. (2017) for a non-local problem, and have been employed for both the 1D and 2D optimal control problems for the KSE in Gomes et al. (2017) and Tomlin et al. (2019), respectively. We predominantly utilized the fourth-order BDF scheme, and performed tests with the other schemes for validation. For us, with the addition of the forcing, the operators in (2.17) are defined as where c is chosen to ensure that A is positive definite. Although the forcing ζ is very singular, being a summation of Dirac delta functions, B still satisfies the required Lipschitz condition if the controls φ j themselves are Lipschitz in the state η. This can be seen from the following calculation where we assume that the controls consist of one point actuator at for the control schemes, where μ is a Lipschitz constant. This is not true in general for the proportional controls with point-observations (since L 2 ⊂ C 0 ), but holds trivially for our full-state feedback controls since it is assumed that the control is a linear function of the Fourier coefficients of the state-see the details in Sections 3 and 4. However, since our controlled solutions remain sufficiently regular, this Lipschitz bound is not an issue; a convergence study follows in the next section. We discretize the spatial domain Q with 2M equidistant points in the streamwise x-direction, and 2N equidistant points in the transverse y-direction, producing a grid of 2M × 2N spatial points, and a corresponding frequency truncation in Fourier space that resolves modes with wavenumbers |k 1 | ≤ M − 1 and |k 2 | ≤ N − 1. The BDF method is then applied to the 2D Fast Fourier Transform (FFT) of the discretized interface, and since (η 2 ) x /2 = ηη x , the nonlinearity may be calculated using the 2D FFT of η 2 . With (2.13) and (2.14), it is clear that actuators do not need to be centered at computational grid points; additionally, having computed the Fourier coefficients of the solution at a given time, observations at any point location in Q may be obtained to spectral accuracy using (2.11). However, taking observations at grid points reduces the computational cost of each time-step.
In many parts we take random initial conditions for our numerical simulations; for these we use where the coefficients a k and b k are random numbers from the interval (−0.05, 0.05). All initial conditions considered have zero spatial mean (the mean is a conserved quantity for the uncontrolled system).
To measure control performance, we track the time-dependent costs where C 1 measures the cost of the solution deviation from the desired state, and C 2 measures the cost of the controls. We emphasize that our aim is not to minimize the above costs, as in optimal control. Both costs are spatially dimensionless, which is appropriate for our study on spatially periodic domains-the costs over one period are the same as the costs over any number of periods considered together. Note that C 1 and C 2 are not necessarily equivalent (in the analytical sense); however, they yield the same behaviour in our numerical simulations as will be seen below.

Actuator arrangements.
In this paper, we consider a variety of actuator (and observer) arrangements. Gomes et al. (2017) studied the optimal actuator placement for the 1D KSE (1.2). For a given initial condition and finite time interval, they obtained the optimal (in terms of some cost functional) open-loop control with contributions from a finite set of optimally placed point actuators. The control locations were heavily dependent on the initial condition, desired state and time interval over which controls were applied, thus the results do not imply that any particular arrangement of actuators performed well across a broad range of initial conditions. In this work, we do not seek results that are initial condition dependent or optimized for a cost functional, instead seeking actuator arrangements that give the best control performance in general. We utilize the following families of actuator arrangements in our numerical simulations, examples of which are shown in Fig. 2: (a) Equidistant. Actuators are equally spaced in the x-and y-directions with separations d 1 and d 2 , respectively, forming a rectangular lattice (or a square lattice if d 1 = d 2 ). In order to comply with the periodic boundary conditions, both L 1 /d 1 and L 2 /d 1 must be positive integers.
(b) Perturbed equidistant. Beginning with an equidistant arrangement, we randomly perturb each actuator slightly (so that one actuator remains in each d 1 × d 2 rectangular region). The random shifts in x and y are sampled from a normal distribution with zero mean.
(c) Random. The locations of the actuators are obtained by sampling coordinates from uniform distributions, i.e. x j ∼ Unif(0, L 1 ), y j ∼ Unif(0, L 2 ). Repeated locations are discarded.
(d) Quasirandom. Also known as low-discrepancy sequences, quasirandom sequences are commonly used to sample space more evenly than uniform distributions. We use the 2,3-Halton sequence (Halton, 1960) to obtain quasirandom sequences of points in [0, 1] × [0, 1], which are appropriately rescaled to yield actuator locations in Q.
Equidistant actuators are commonly considered in studies of 1D control problems; however, such arrangements are unsuccessful when the actuators are located at zeros of unstable eigenfunctions. Obviously there are actuator arrangements of interest not in the above classes, which we do not study here, e.g. parallelogram lattices or hexagonal patterns. The success of a control strategy for a particular actuator arrangement is equation and boundary condition dependent. The arrangements we consider appeared most natural for a study on rectangular periodic domains-for a problem where periodicity is enforced on a hexagonal domain, a hexagonal actuator arrangement would be the most appropriate.
We measure three areas to quantify the spacing of the various arrangements-examples of each are shown as shaded regions in Fig. 2(b-d). The first, A 1 , is defined as the area of the largest circle centred at an actuator that contains no other actuators, shown in panel (b) for the case of a perturbed equidistant arrangement. If all actuators are placed in groups of more than one, then A 1 will be small, with an infimum of zero found in the limit of each actuator approaching another. The supremum of (L 2 1 +L 2 2 )π/4 corresponds to moving one actuator the maximum distance away from a cluster of all the other actuators. The value of A 1 is thus not entirely informative, but we find that its deviation from the value for the equidistant case, A E 1 = π · min{d 1 , d 2 } 2 , is a more appropriate metric for our study. The quantity A 2 is defined via the Voronoi tessellation, shown in panel (c) with dashed lines, which separates the plane into the sets of points that are closest in Euclidean distance to each actuator. We define A 2 to be the area of the largest Voronoi cell, with the minimum value attained for the equidistant case, A E 2 = d 1 d 2 . Lastly, A 3 , shown in Fig. 2(d), is defined as the area of the largest circle that can be inscribed in the plane without containing any actuators-this is the solution of the well-known largest empty circle problem (Toussaint, 1983). This last area is minimized for equidistant actuator arrangements with A E 3 = π(d 2 1 + d 2 2 )/4, and in contrast to A 1 and A 2 , quantifies the size of the gaps between the actuators.

Proportional control
In this section, we consider proportional point-actuated controls. We assume that both observation and actuation occur at the same set of locations in the periodic domain Q. The strength and orientation of actuation at a point are based only on the observation of the local interface height at that instant in time.
There is no communication from observers at other actuator locations, and no a priori knowledge of the governing equation is utilized; this is the most basic level of feedback control.
To motivate our study, we allow actuation and observation at every location in space, and consider controls of the form ζ = −αη for a strength α ≥ 0 to stabilize the trivial zero solution, i.e. η = 0. The dispersion relation for the controlled KSE (1.1) is then where s is the linear growth rate andk is the scaled wavenumber vector for k ∈ Z 2 given by (2.10).
Multiplying (1.1) by η and integrating by parts, we find that the energy of the solution (given by the L 2 -norm) evolves according to the global energy equation If s < 0 for all arguments, the L 2 -norm of η decays exponentially for all choices of L 1 and L 2 ; it follows from (3.1) that this is achieved for α > α c (κ), where For our setting with periodic boundary conditions, taking α > α c is sufficient to exponentially stabilize the zero state, but a sharp stabilization condition may be obtained in terms of L 1 and L 2 . For a non-trivial desired state, this control generalizes to ζ = −α(η − η), and condition (3.3) becomes more complicated since the system must be linearized about a non-trivial state.
To explicitly obtain a critical actuation strength we used knowledge of the governing equation; if a PDE is well-posed in the very weak sense that its linear part is dissipative at small scales, then there exists an α c such that for all α > α c , the linear controlled system is stable. For an unknown system, this may be found experimentally by simply increasing the actuation strength α. It is not necessarily true that this will yield full nonlinear stability of the given system; however, we expect that many interfacial problems with weak nonlinear interactions may respond well to this form of control.
It may not be viable to actuate at every spatial location, or even observe the entire interface. Therefore, for the remainder of the section, we consider point-actuated controls (2.13) with where α ≥ 0 is a strength parameter. The above is often referred to as 'pinning control' in the literature (Grigoriev et al., 1997). We note that the solution average is not necessarily preserved by this choice of forcing, but successful controls ensure that solutions do not deviate significantly from zero mean (since η has zero mean).

Convergence study
Since convergence of the BDF scheme is not guaranteed analytically for point-actuated proportional controls, and we are utilizing particularly singular actuator functions, it is appropriate to perform a convergence study. For this, we take L 1 = L 2 = 21 and κ = 0.25 (corresponding to an overlying liquid film). We use N ctrl = 49 actuators/observers in a quasirandom arrangement, located according to the 2,3-Halton sequence-see Fig. 2(d). Here, we fix the initial condition to be (3.5) and let the system evolve until t = 1 without controls. For t ∈ [1, 2], we actuate proportionally with strength α = 150 to drive the solution to η = 0. We give the value of the L 2 -norm at time t = 2, C 1 (2), for a range of spatial discretizations M, N, and time discretization Δt in Table 2.
Recall that we are employing a fourth-order BDF scheme that is not unconditionally stable, unlike the first-and second-order BDF schemes. In a convergence study of the unforced equation, Akrivis et al. (2016) observed, for a particular set of parameters, that once the higher-order BDF schemes (third-to sixth-order) are stable, the solution converges to machine accuracy-even for the worst case in Table 2, we found that C 1 (1) = 3.2652720 was accurate to eight significant figures. Analogously, with the addition of point actuators, we find that once stability of the fourth-order scheme has been achieved (the scheme is not stable for Δt = 2 × 10 −4 ), then variations of C 1 (2) due to changes in Δt are on the order of 10 −6 . Variations due to the refinement of the spatial discretization are larger, but it is evident that all results are accurate to at least two significant figures. For the unforced equation, numerical experiments showed that the spectrum of solutions decays exponentially, indicating analyticity (Tomlin et al., 2018).  It was also observed that the spectrum decays faster in the transversek 2 modes than the streamwisẽ k 1 modes; this is expected given the asymmetry of the linear part of (1.1). In numerical simulations on a square periodic domain, the Fourier mode truncation N is not required to be as large as M for good accuracy. This is no longer true with the introduction of Dirac delta forcing-the decay of the Fourier spectrum becomes more symmetric in wavenumber space, placing similar restrictions on the constants M and N to obtain good accuracy. Thus, in Table 2, we only consider M = N. The spectrum at t = 2 from the most well-resolved simulation is shown in Fig. 3. Analyticity (exponential decay of the spectrum) is lost due to the Dirac delta forcing, and we find numerically that |η k | ∼ |k| −4 . Thus, Δ 2 η has an approximately constant spectrum, |k| 4 |η k | ∼ O(1), which balances the constant spectrum of the Dirac delta forcing.

Controlling unbounded exponential growth
In this subsection, we investigate the efficacy of proportional point-actuated controls in suppressing the unbounded growth of solutions to the 2D KSE (1.1) with κ = −0.5 < 0. This choice of κ physically corresponds to a hanging fluid film, with linear Rayleigh-Taylor instabilities in the transverse modes. Beyond the weakly nonlinear regime that is modelled by (1.1), this instability gives rise to spanwise rivulet structures in fluid films-see the experiments of Charogiannis et al. (2018). Additionally, there are the usual linear instabilities in the streamwise and mixed modes present for 0 ≤ κ < 1. Optimal controls involving the transverse modes alone were applied to (1.1) in Tomlin et al. (2019), revealing windows of chaotic and travelling wave attractors for the streamwise and mixed modes. The success of the proportional controls will be measured by two objectives; the first being the successful suppression of transverse growth resulting in bounded solutions, and the second being the stronger property of exponential stabilization of the desired state. We set η = 0 for simplicity, and consider a square domain with L 1 = L 2 = 18, for which 30 Fourier modes are linearly unstable in total, including the transverse (0, ±1)and (0, ±2)-modes. For our numerical simulations, we take the initial condition (3.5), which has contributions from all unstable transverse modes. Controls are applied from the initial time, rather than letting the transverse instabilities develop further. For this subsection, we use random actuator arrangements that are constructed recursively as follows. Actuator locations are obtained by sampling coordinates from a uniform distribution, and are ordered in a list so that taking N ctrl = p corresponds to 'switching on' the first p control actuators/observers in the list. This arrangement and ordering is fixed in this subsection, across all simulations for different values of α. Fixing the initial condition and having consistency in how the actuators are arranged and the order in which they are 'switched on' allows us to draw robust conclusions. The success or failure of the controls is analysed in the α-N ctrl plane. Figure 4 shows the costs C 1 , C 2 , and solution contours for α = 55 and N ctrl = 100, 110. The actuator locations are superimposed on the solution contours in panels (b,d); the actuators represented by black circles with white edges are 'switched on' for N ctrl = 110, but not N ctrl = 100. Panels (a,b) show the plots for N ctrl = 100; the costs plateau at a finite value and cellular humps are visible in regions of the interface where no actuators are active. For the N ctrl = 110 case in panels (c,d), the zero solution is exponentially stabilized. One of the actuators that is 'switched on' for N ctrl = 110 is located at the site of a cellular hump which forms for the N ctrl = 100 simulation, emphasizing the importance of actuator placement.
In order to determine how the control success depends on the parameters α and N ctrl , numerical experiments were carried out with N ctrl ranging from 0 to 130, and α ∈ [0, 150]. Figure 5(a) provides the numerical results over a section of the α-N ctrl plane, where the parameter space is seen to be divided into three regions depending on the behaviour of the costs C 1 and C 2 . Unsurprisingly, for sufficiently small values of α and N ctrl , the growth of the transverse modes is not saturated. For α ≥ 1.5 and N ctrl ≥ 18, approximately, bounded controlled solutions emerge as seen in Fig. 4(b). Exponential stabilization is obtained for much larger values of the control parameters; we find that for α = 24, the solution is bounded and non-zero for N ctrl = 108, but find exponential stabilization of the zero solution for N ctrl = 109. The regime boundaries on the right of Fig. 5(a) continue at a constant value of N ctrl for the strengths α ∈ [40, 150] (not shown in the figure). In Fig. 5(b), the exponential decay rate λ of the costs C 1 and C 2 (the decay rate of these two quantities is approximately the same) is plotted against α for the arrangement with N ctrl = 110. It is evident that λ is a monotonically increasing function of α, but that there is a maximal exponential decay rate associated with each actuator arrangement. For Fig. 5(b), the maximal value of λ = 0.09 is predicted by curve fitting in the limit of large α.
The sharp critical values in Fig. 5(a) indicate that there is no trade-off between the quantity or placement of actuators and control strength; the same level of control may not necessarily be achieved for smaller α and increased N ctrl (or vice versa). Pinning the solution at a set of spatial locations with large α will not ensure exponential stabilization of the desired state if there are sufficiently large regions of Q in which no actuators are located, since cellular humps form as in Fig. 4(b). Moreover, many actuators cannot stabilize an unstable state if α is too small. Costs and solution contours for α = 55 and N ctrl = 100, 110. Panels (a,c) display the evolution of the costs C 1 and C 2 for the cases of N ctrl = 100 and N ctrl = 110, respectively. The uncontrolled case with exponentially growing C 1 is included for reference. Panels (b,d) display contours of the respective solution profiles at t = 200. The actuator locations are superimposed on the solution contours in panels (b,d)-the white circles with black edges denote the point actuators used in both cases, and the black circles with white edges denote the remaining 10, which are 'switched on' for N ctrl = 110.

Comparison of actuator arrangements
In this subsection, we test different actuator arrangements with the aim of understanding the relationship between actuator spacing and control performance. We vary the domain size |Q|, while keeping the ratio |Q|/N ctrl constant, i.e. we have a fixed number of actuators per unit area, and set α = 150 (this is chosen to be large so that the control strength is not a limiting factor-see the previous subsection). A finite energy density is observed for solutions of the uncontrolled system with κ = 0 on large domains (Tomlin et al., 2018), and thus it is expected that the number of actuators required per unit area for successful control should be constant in the limit of large periodicities, L 1 , L 2 → ∞.
We set κ = 0.25, which yields bounded solution dynamics, restrict to square domains with L 1 = L 2 = L and take N ctrl to be a square number so that we may compare the results of random and quasirandom arrangements with results for equidistant and perturbed equidistant actuator arrangements Downloaded from https://academic.oup.com/imamat/advance-article-abstract/doi/10.1093/imamat/hxz031/5678719 by University of Warwick user on 18 December 2019  (with d 1 = d 2 ). We additionally retain the choice of η = 0. The domain dimensions and numbers of actuators are chosen so that |Q|/N ctrl = 9, i.e. one control actuator per nine units of domain area, and are summarized in Table 3. For each case in Table 3, we perform simulations for the unique equidistant and quasirandom (using the 2, 3-Halton sequence) actuator arrangements, one perturbed equidistant and five different random actuator arrangements, constituting a total of eight arrangements for each of the nine choices of (L, N ctrl ).
For each simulation, we take a random initial condition with zero spatial average, containing sufficiently many unstable low modes. In contrast to the previous subsection, we allow the system to evolve without controls until it reaches the global attractor (200 time units suffices for the above cases), and then apply proportional point-actuated controls. For our choice of α = 150, exponential stabilization of the zero solution is observed for all eight actuator arrangements with L = 21. At this domain size, we observe the emergence of bimodal states in the absence of controls, with the onset of chaos for slightly larger L.
The results of the numerical experiments are shown in Fig. 6. The decay rate of the costs, λ, is plotted against the areas |A 1 − A E 1 |, A 2 and A 3 in panels (a-c), respectively. The equidistant cases performed the best, closely followed by the perturbed equidistant arrangements, both with λ ≈ 0.4. The quasirandom arrangements, which are much more regularly spaced than the random ones, all gave exponentially decaying costs with λ ≈ 0.2. The solutions controlled with random actuator arrangements either reached a non-trivial steady state (for which we assign λ = 0 in the figure), or decayed to zero with rate no greater than 0.1. Although not obvious from Fig. 6, we found that the proportion of random arrangements that resulted in exponential stability of the flat state decreased as L increased, with no successes for L = 45. Fig. 6. Performance of various actuator arrangements. The decay rate λ of the costs C 1 and C 2 is plotted against the areas |A 1 − A E 1 |, A 2 and A 3 . In cases where exponential decay is not observed, the value of λ = 0 is assigned. Markers: -Equidistant; ♦ -Perturbed equidistant; -Quasirandom; -Random. The markers are shaded according to the values of L in Table 3, with the lightest corresponding to L = 21 and the darkest to L = 45. For A 2 and A 3 , the data for L = 39, 42, 45 are used to construct threshold models (details given in the text), denoted by the dotted lines in panels (b,c). From Fig. 6(a), we see that the maximum attainable decay rate is a monotonically decreasing function of |A 1 − A E 1 | (the points are bounded above by a monotonically decreasing curve). Actuator arrangements with |A 1 − A E 1 | 10 perform poorly, with exponential decay rates below 0.1. For A 2 and A 3 , we fit a threshold model using the data for L = 39, 42, 45 of the form λ = a j (A j − A c j ) 2 for A j ≤ A c j , and λ = 0 (no exponential decay) for A j > A c j ; the constants a j and A c j are computed using least squares fitting and optimization. The threshold model is not appropriate for fitting with |A 1 − A E 1 |, since the data-points do not clearly lie on a curve, whereas it is evidently appropriate for A 2 and A 3 . We obtain λ = 2.4 × 10 −3 (A 2 − 22.4) 2 for A 2 ≤ 22.4 and λ = 1.5 × 10 −4 (A 3 − 67.2) 2 for A 3 ≤ 67.2 (with λ = 0 otherwise). The least squares errors between the data-points and the threshold models are approximately the same across panels (b,c). These threshold models appear to be equally suitable predictors for the success/failure of controls, providing reasonable estimates of the decay rate based on values of A 2 and/or A 3 for a given actuator arrangement. We expect that similar measures of spacing would also be useful for the point-actuated control of other physical systems. In particular, such measures may be helpful in problems where geometrical constraints are placed on the arrangement of actuators/observers.
For equidistant arrangements, exponential stability was achieved for all cases in this numerical experiment. However, these fail when the actuators lie at the zeros of unstable eigenfunctions (in this section, the spacing between actuators was finer than the shortest unstable wavelength of the system). Although the results are not presented here, we found that in such situations, slightly shifting the pointactuator locations, i.e. using a perturbed equidistant actuator arrangement, does not prevent failure of the controls.

Non-trivial desired states: synchronization of chaotic dynamics
The proportional control methodology is very robust in the sense that it allows for any choice of desired state η, even non-solutions. For full convergence to an arbitrary non-solution, a dense set of control actuators is required, whereas if η is a solution of the uncontrolled equation, then it can usually be stabilized with a finite number of actuators. Non-trivial steady and travelling waves are popular target states for stabilization in the control literature; however, in this section, we show that even more complicated situations can be tackled, and employ proportional controls to stabilize a given chaotic orbit of the 2D KSE (1.1).
The problem of synchronizing a pair of chaotic solution orbits of the 1D KSE (1.2) was considered by Junge et al. (1999) and Tasev et al. (2000). The authors used actuators of non-zero width, and performed proportional control using local spatial averages of η − η over the actuator regions-in the limit as the actuator/averaging width becomes small, this converges to (3.4) with the Dirac delta actuators. They found that the number of actuators required for synchronization scaled with the size of the periodic domain, consistent with the results of Gomes et al. (2017) and our own observations in the previous subsection. It was also observed that the equidistant arrangement was close to optimal-they suggest that the discrepancy may be due to the imposed rigid boundary conditions. Their study was extended to a generalized KSE by Basnarkov et al. (2014) with the addition of a third order dispersion term, γ η xxx ; they discussed the synchronization of solutions to systems with different values of the dispersion strength γ .
For our numerical experiments, we take κ = 0.25 and L 1 = L 2 = 45 (52 unstable modes in total), and consider both the quasirandom and equidistant arrangements with N ctrl = 225 from the previous subsection. We take the initial condition η(x, 0) = η 0 (x) given in (3.5), and the desired state is the orbit starting from η(x, 0) = 2η 0 (x). Proportional controls (recall (3.4) for the case of a non-trivial desired state) are applied with α = 5 from time t = 550, with the numerical simulation running until 600 time units. Figure 7 shows the results for the quasirandom and equidistant arrangements in panels (a,b) and (c,d), respectively. The costs in panel (a) decay exponentially after controls are applied (although the rate is not constant), but the decay rate visible in panel (c) for the equidistant case is much larger (note the differing limits of the Cost-axes in these panels). The projection of the orbits onto a 3D phase space tells a similar story, with the orbit of the controlled solution tracking that of η much more closely in panel (d) than (b). In the former case, the orbits very rapidly become indistinguishable after controls are applied. See Movie 1 in the supplementary material for a movie of the synchronizing chaotic interfaces using the quasirandom actuator arrangement for times t ∈ [500, 600].

Feedback control with full state observations
In this section we study feedback control strategies with observation of the entire interface, which, along with the knowledge of the linearized dynamics, is employed to construct controls. This is in contrast to the previous section, where actuation at a point was governed by the local film thickness alone. We note that, as will be elaborated below, this control methodology only requires observation of a finite set of Fourier modes of the interface-this is advantageous for the numerical implementation.
In order to apply the following theory, the assumption that η is an exact solution of (1.1) with ζ = 0 is required. Convergence to non-solutions can only be achieved transiently with our control approach (see the discussion in Thompson et al., 2016a). The difference w = η − η evolves according to the We define A to be the diagonal matrix with and let B denote the control actuator matrix with entries B k,j = b j k . We assume that the controls φ j depend linearly on the Fourier coefficients of w through a constant matrix K (to be chosen) by the relation where w is the vector containing the Fourier coefficients of w. This is known as a linear state feedback control law, and K is the feedback gain matrix. The control ζ has Fourier coefficients We may then rewrite (4.1) in terms of matrices as where ν is the vector with entries being the Fourier coefficients of −ww x , and J is the matrix such that Jw gives the Fourier coefficients of −(ηw) x . The entries of J may be computed as J k,l = −ik 1 η k−l . The task now is to choose the feedback gain matrix K so that this nonlinear infinite-dimensional dynamical system is stable about w = 0. It is well known that the characterizations of stability of linear infinite-dimensional dynamical systems are more intricate than the finite-dimensional case (Zabczyk, 2009), and further complications arise when nonlinear terms are added. We now describe two methodologies that are applicable to our problem, both having been considered for the 1D KSE (we note that a comparison was not performed in 1D). Both follow an early-lumping approach, in which the full infinite-dimensional system is first truncated to a finite-dimensional one, for which controls are constructed (Morris & Levine, 2010). Such controllers are susceptible to 'spillover effects', in which the controls based on a reduced system excite unmodelled modes, leading to instability (Balas, 1978). 'Controller spillover' can be avoided for our problem if the reduced system involves both the unstable modes and modes that are weakly damped (eigenvalue shifts due to controls can make these unstable), and if the feedback gains are not too large.
Al Jamal & Morris (2018) showed for the 1D KSE (1.2) that the full infinite-dimensional nonlinear system may be controlled by ensuring the stability of a suitable truncation of the linearized system. Linearizing our problem about the desired state η, i.e. linearizing (4.5) about w = 0 (or (4.1) about w = 0), gives (4.5) with ν = 0. We let A n denote the truncation of the matrix operator A to the Fourier modes with |k 1 |, |k 2 | ≤ n, and take similar notations for the other matrices and vectors. The analysis of Al Jamal & Morris (2018) can be lifted to our 2D setting, and their Theorem 5.1 may be recast as: Theorem 4.1 (Theorem 5.1 in Al Jamal & Morris, 2018) Consider the sequence of approximations of the linear matrix problem (4.5) with ν = 0 defined by the Galerkin truncation in Fourier space onto the modes with |k 1 |, |k 2 | ≤ n, w n t = J n + A n + B n K n w n . (4.6) Assume that there exists a convergent sequence of controllers K n that stabilize the problem (4.6), and such that the limit K exponentially stabilizes (4.5) with ν = 0. Then for sufficiently large n, the feedback gain matrix K n exponentially stabilizes the full nonlinear problem (4.5).
To realize this numerically, we take a sufficiently large truncation n and compute a (possibly time varying) matrix K n such that (4.6) is stable-this is then extended with zeros to obtain a viable choice of K. We expand upon how such a K n is chosen later, but note that it must possess certain symmetries to ensure that the numerical solution remains real-valued. If η is a steady state, the matrix J n is constant in time (J n is zero if η = 0). In this case, K n is constant in time and only needs to be computed once, giving a static feedback law. If η is an exact travelling wave solution with period τ , then K n needs to be prescribed at each time-step in [0, τ ] as J n varies, resulting in a dynamic feedback law. The matrix K n should vary continuously in time; it is computationally excessive and unfeasible to compute K n based on the current value of J n at each time-step. Constructing feedback controllers for (4.6) becomes even less feasible for choices of η with more complicated dynamics, such as quasi-periodicity or chaos.
A different approach was considered by Gomes et al. (2017) in their study of the 1D problem. The authors linearize about the zero solution regardless of the desired state η (this methodology coincides with the former one for η = 0). The matrix K n is chosen to stabilize the linear system w n t = A n + B n K n w n , (4.7) under a further restriction that averts the problem caused by disregarding J n , as described next. For λ > 0, we define λ-stability for a complex matrix C to be the property that for any complex vector z, (4.8) where † denotes the conjugate transpose. If C is a normal matrix, satisfying C † C = CC † , condition (4.8) is equivalent to the real parts of all eigenvalues of C being bounded above by −λ. Following the same argument as Gomes et al. (2017), in addition to stabilizing (4.7), the truncated feedback gain matrix K n must be chosen such that the bracketed term in (4.7) is λ-stable for Then, if the truncation is suitably large, the full nonlinear system (4.5) is exponentially stabilized with decay rate given by the left hand side of (4.9). In practice, we replace the λ-stability requirement with the weaker condition of having the real parts of all eigenvalues bounded above by −λ, as proved successful in Gomes et al. (2017)-this is not equivalent since the bracketed term in (4.7) is non-normal, and exponential stabilization is no longer guaranteed analytically. Although transient growth is observed (as is possible for non-normal systems), we find that this constraint on the eigenvalues is sufficient to ensure long-time exponential stability in our numerical simulations. Moreover, we find that the decay rate surpasses the left hand side of (4.9) significantly. We construct the feedback gain matrix K n via an exact pole placement approach to yield desired eigenvalues for the linear systems (4.6),(4.7). For this we use the Matlab function place, which was developed by Kautsky et al. (1985). The procedure is difficult to carry out in the complex Fourier mode framework, since in this basis, the entries of K n may be complex yet must satisfy symmetry requirements so that the resulting forcing is real-valued. Thus we translate the problem into real-valued, trigonometric basis functions by applying linear transformations to the matrices J n , A n and B n . Given the matrices J n , A n , B n (in the new basis), and a vector of desired eigenvalues, place computes a suitable real matrix K n such that the linearized system possesses the desired eigenvalues. The matrix K n is then transformed back into the original complex Fourier mode basis for use in numerical simulations. The solution given by the exact pole placement algorithm place is robust to changes in the input matrices (Kautsky et al., 1985), an attractive property for the present study. It is useful for implementing the approach of Gomes et al. (2017) since the true linearized system (4.6) and that which is stabilized (4.7) differ by the matrix J n . Additionally, a robust controller that is constructed for the 2D KSE may be applicable to more complicated thin film flow models, such as the Benney equation (2.7).
The precise details of how the desired eigenvalues are chosen is given in the following text, yet in broad terms, we choose the same set of eigenvalues as in the uncontrolled system, but replacing the real parts of those beyond a given threshold with a negative value. It is important to note that, with this procedure of choosing the desired spectrum, the matrices K n do not form a convergent sequence as described in Theorem 4.1. For large values of n, the method breaks down, with K n having large entries ('large gains'), resulting in a heavily sensitive problem that is vulnerable to spillover effects and/or numerical instability. Optimized regional or partial pole placement algorithms may be more appropriate for this situation (and would be required for systems with many more unstable modes), which guarantee that the spectrum is located in a desired subset of the left half of the complex plane, and also minimize the entries of the feedback gain matrix, such as the algorithm of Datta & Chakraborty (2014); this would also reduce the cost of the applied controls. The investigation of other pole placement algorithms in beyond the scope of the present work.
We now present numerical experiments testing the effectiveness and applicability of full-state feedback controls in three different situations.

Stabilizing the trivial flat solution
We first apply full-state feedback controls to stabilize the flat film solution, taking η = 0. In this case, the methodologies in Al Jamal & Morris (2018) and Gomes et al. (2017) agree. We take parameters κ = 0.25, L 1 = L 2 = 42 and N ctrl = 196 with one of the random actuator arrangements previously used for numerical simulations in subsection 3.3. With proportional controls and α = 150, exponential stability of the zero solution was not achieved for this particular actuator arrangement. As in subsection 3.3, controls are applied after 200 time units, and the (random field) initial condition is the same as used there also. We construct a controller based on the linearized system with a truncation of n = 19, which importantly covers the linearly unstable modes-note that since J n is a matrix of zeroes, the eigenvalues of the uncontrolled system are real, being the diagonal entries of A given in (4.2). The feedback gain matrix K n is computed using the Matlab function place so that the eigenvalues of the linearized system (4.7) are real and at most −0.1. More precisely, the eigenvalues of the linearized system are unchanged if they are less than −0.1, with the rest replaced by −0.1. The evolution of the costs is shown in Fig. 8(a), and it can be seen that exponential stabilization is achieved with decay rate 0.1 (correct to 7 decimal places) for the full-state feedback controls. This is because −0.1 was chosen to be the largest eigenvalue of the truncated linearized system, improving on the decay rates obtained using random actuator arrangements and proportional controls in subsection 3.3. We find that C 2 (200) (the initial control cost) is an order of magnitude smaller for full-state feedback control than for proportional control. In Fig. 8(b) we plot the actuation strengths applied at the same two actuators, |φ 1 | and |φ 2 |, for both the proportional and full-state feedback control strategies. The actuation strengths for the full-state feedback controls (shown with solid lines) all decay exponentially after an initial phase of transient growth, but the proportional controls (dotted lines) 'overcontrol' the system (indicated by the singularities in the log-linear plot corresponding to a sign change).
The full-state feedback controls are thus shown to perform well for the case of η = 0, even with actuator arrangements that are far from optimal. The costs decayed exponentially at the prescribed rate, and convergence to the desired state was achieved to machine precision. In the next two subsections, we investigate the stabilization of non-trivial desired states.

Stabilizing a non-trivial steady state
In this subsection, we consider the stabilization of a non-trivial steady solution of (1.1) with full-state feedback controls; the approaches of Al Jamal & Morris (2018) and Gomes et al. (2017) differ in this case, yet are both computationally feasible. For this numerical experiment, we take the parameters κ = −0.5, L 1 = L 2 = 18, and use N ctrl = 100 randomly located actuators (the arrangement used for the numerical experiments in subsection 3.2, see Fig. 4). We take initial condition (3.5), and choose η to be a 1D steady state, which is constant in y; the profile is plotted in Fig. 9(a). Such 1D steady states may be computed with ease using the continuation and bifurcation software AUTO-07P (Doedel et al., 2007), for example. This exact solution of the uncontrolled 2D KSE is unstable to transverse perturbations, and is even unstable to streamwise perturbations; chaos is already prevalent in the 1D KSE with analogous parameters.
Downloaded from https://academic.oup.com/imamat/advance-article-abstract/doi/10.1093/imamat/hxz031/5678719 by University of Warwick user on 18 December 2019 Guided by the analytical result of Gomes et al. (2017), we choose the desired eigenvalues for the exact pole placement algorithm to satisfy (4.9); however, we do not ensure λ-stability. For the steady state shown in Fig. 9(a), the infimum of the x-derivative is computed to be approximately −6.762, and so we choose the eigenvalues of the controlled system (4.7) to be bounded above by −3.5. For the construction of the feedback gain matrix, we use a mode truncation of n = 9. Since the number of repeated eigenvalues cannot be more than the number of controls (place cannot assign poles with multiplicity greater than the rank of the matrix B n ), the desired spectrum for the controlled problem is chosen similarly as in the previous subsection, but with eigenvalues larger than −3.5 replaced with −(3.5 + U) where U ∼ Unif(0, 0.1). Although not required, we chose the desired eigenvalues for our computations using the method of Al Jamal & Morris (2018) in a similar way to ensure a fair comparison, replacing the real parts of the complex eigenvalues of J n + A n , which exceed −3.5 with −(3.5 + U).
The costs corresponding to the methodology of Al Jamal & Morris (2018), shown with thick lines in Fig. 9(b), decay at the expected rate of 3.5 approximately. For the approach of Gomes et al. (2017), shown with thin lines in Fig. 9(b), the costs initially proceed as in the previous case, before switching to a less extreme decay rate, approximately 1.7, at around t = 1. This decay is much greater than the left hand side of (4.9), which is about 0.1. It is apparent that a weaker bound on the spectrum of the controlled linear system may be viable for the approach of Gomes et al. (2017) and still obtain long-time exponential stabilization-an investigation of this is beyond the scope of the present study. Deviations of the solution from the desired state, η − η, are plotted at t = 2 for the two approaches in Fig. 9(c,d); the difference plotted in panel (d) is an order of magnitude smaller than that in panel (c).

Stabilizing a travelling wave solution
We now focus our attention on the stabilization of an exact travelling wave solution. For this, we take the same parameters and initial condition (3.5) as in the previous subsection, but use a random actuator arrangement with N ctrl = 50 (previously used in subsection 3.2). We choose η to be the 1D travelling wave of speed c ≈ 1.356 with cross-section shown in Fig. 10(a); this is computed using AUTO-07P. As was true for the steady state considered above, this exact travelling wave solution is unstable to both streamwise and transverse perturbations.
To construct a controller following the approach of Gomes et al. (2017), only one feedback gain matrix is calculated. The infimum of the x-derivative of the desired wave shown in Fig. 10(a) is approximately −4.802; to satisfy (4.9), the eigenvalues of the linearized system (4.7) are chosen to have real parts smaller than −2.5. We do this in a similar way to the previous numerical experiment, choosing K n so that we retain the same spectrum as in the uncontrolled problem, but with the real eigenvalues larger than −2.5 replaced with −(2.5 + U) where U ∼ Unif(0, 0.1). For the methodology of Al Jamal & Morris (2018), it is not computationally feasible to construct feedback controllers for the linearized system (4.6) at each individual time-step. The pole placement algorithm is not fast enough, nor is it guaranteed that the obtained feedback gain matrices vary continuously in time as J n varies. We proceed instead by constructing a hybrid controller, for which the feedback gain matrix K n linearly interpolates between feedback gain matrices that stabilize (4.6) at μ ∈ N times {0, τ/μ, 2τ/μ, . . . , τ −τ/μ} ⊂ [0, τ ], where τ is the time-period of the travelling wave. Precisely, we take K n (t) = μ−1 m=0 ψ m (t; τ , μ)K n m , (4.10) where the ψ m are τ -periodic tent functions, commonly used for piecewise linear interpolation, chosen so that K n (mτ/μ) = K n m with a linear variation in between. Again, to guarantee a fair comparison between this and the approach of Gomes et al. (2017), the K n m are constructed so that the real parts of the eigenvalues of (4.6) also satisfy (4.9) at the μ times. The success of the hybrid controller relies on the robustness of the pole placement in a similar way to the controller based on the linearized system (4.7), which disregards J n entirely.
For the numerical simulations we present, controls were constructed with a mode truncation of n = 19 (we note that n = 14 also performed well). We used μ = 8 for the hybrid controller-the period of the travelling wave is τ ≈ 13.3, so that τ/μ ≈ 1.66 time units. The costs for the controlled evolutions are plotted in Fig. 10(b), with all exhibiting a brief period of transient growth (due to nonnormality), and then exponential decay at non-constant rates. The effective decay rate attained by the Gomes et al. (2017) controller is approximately 1.1, with the hybrid controller providing faster decay with an effective rate of about 2. The evolution of the interface controlled using the approach of Gomes . Panel (a) shows the profile of the 1D travelling wave that moves in the streamwise direction with speed c ≈ 1.356 (downstream in the direction of increasing x). Panel (b) plots the costs C 1 and C 2 for the applied full-state feedback controls, both using the approach of Gomes et al. (2017) and a hybrid controller with μ = 8 (details given in the text), and also plots C 1 in the uncontrolled case for reference. Panels (c,d) show the controlled interface at times t = 0.75, 2.5, respectively, for which the feedback controller was constructed following the approach of Gomes et al. (2017Gomes et al. ( ). et al. (2017 is shown in Movie 2 in the supplementary material. Fig. 10(c,d) presents snapshots of this movie at t = 0.75 and t = 2.5, respectively.
Although not presented here, feedback controls can be constructed to stabilize chaotic orbits following the method of Gomes et al. (2017); it is not feasible to construct controllers for the latter based on the fully linearized dynamics as in Al Jamal & Morris (2018). However, with the impressive robustness of place, it may be possible to stabilize quasi-time-periodic solutions using the hybrid controller (4.10) where the K n m are constructed using a typical attractor cycle.

Conclusions
In this paper we considered the feedback control of a multidimensional KSE (1.1) using point-actuators. The equation describes the dynamics of a liquid film interface, and yields steady, travelling, time-periodic and quasi-time-periodic waves, as well as chaotic and unbounded solutions depending on the parameter regime. We applied two closed-loop control strategies, proportional control and feedback control with full state observations. For proportional control, we investigated the limitations of the method depending on the strength, number, and placement of the actuators/observers. The controls were able to prevent the unbounded growth of the interface, and exponentially stabilize the flat film solution. We used three measures of actuator spacing, and found that they were strongly correlated with the decay of the controlled system, with a maximum decay rate for equally spaced actuator arrangements. The proportional controls performed well for non-trivial desired states; we were able to synchronize two chaotic orbits of the system. Knowledge of the governing dynamics is not necessary to apply proportional controls, thus they are more easily applicable in experiments. Throughout, we monitored two time-dependent costs that measured the success of controls and the total actuation strength, respectively. If a controller exponentially stabilizes a desired state, then the cost due to actuation also decreases exponentially; if stabilization can be achieved, the cost to maintain it is relatively small. It may even be viable to utilize on/off controls that are inactive when the distance between the solution and the desired state is sufficiently small.
We used feedback control with full state observations to stabilize the zero solution, a steady state and a travelling wave solution; in all test cases, a large number of modes were linearly unstable, with chaos or exponential growth prevalent in the uncontrolled dynamics. Feedback controllers were constructed with a Fourier-Galerkin truncation of the infinite-dimensional system, following the approaches of Al Jamal & Morris (2018) and Gomes et al. (2017), which were originally formulated for the 1D equation (1.2). Feedback gain matrices were computed using the exact pole placement algorithm place, which provided a controlled system whose eigenvalues had negative real parts. Both approaches performed well, with long-time exponential stability obtained for all desired states. For choices of parameters giving many more unstable modes, a different pole placement algorithm is required as the problem becomes too high-dimensional for place.
Although not presented here, the authors also considered controls with dynamical (Luenberger) observers. These supply an estimate of the state from partial measurements (point observations), which is then utilized in the feedback controller. This was studied for the 1D KSE (1.2) by Christofides (1998) and Armaou & Christofides (2000), and for the control of 1D Benney and weighted-residual equations by Thompson et al. (2016a). We found this method to be largely unsuccessful for the situations presented above, and that the results did not improve on those of proportional controls. Due to the highly nonlinear nature of the dynamics, more point observers, or a better state estimator such a sliding mode observer, may be needed to secure exponential stabilization.
The success of controls constructed for interfacial evolution equations when applied to more complicated models (for the same physical system) is unknown. In the thin film scenario, the authors are investigating the possibility of controlling the interface for a flow governed by the full Navier-Stokes equations using controls constructed for long-wave models. Current work by the authors also involves data-driven control strategies for both deterministic and stochastic evolution equations; the aim is to achieve similar levels of success as feedback control with full state observations, without knowledge of the governing equations and only partial observability.