Feedback control of unstable flows: a direct modelling approach using the Eigensystem Realisation Algorithm

Obtaining low-order models for unstable flows in a systematic and computationally tractable manner has been a long-standing challenge. In this study, we show that the Eigensystem Realisation Algorithm (ERA) can be applied directly to unstable flows, and that the resulting models can be used to design robust stabilising feedback controllers. We consider the unstable flow around a D-shaped body, equipped with body-mounted actuators, and sensors located either in the wake or on the base of the body. A linear model is first obtained using approximate balanced truncation. It is then shown that it is straightforward and justified to obtain models for unstable flows by directly applying the ERA to the open-loop impulse response. We show that such models can also be obtained from the response of the nonlinear flow to a small impulse. Using robust control tools, the models are used to design and implement both proportional and $\mathscr{H}_{\infty }$ loop-shaping controllers. The designed controllers were found to be robust enough to stabilise the wake, even from the nonlinear vortex shedding state and in some cases at off-design Reynolds numbers.


Model reduction
In the previous section, we have given a brief overview of the large spectrum of tools that can be used once a linear low-order model of the relationship between the inputs and the outputs of the system has been obtained. In order to construct such a model, two dominant approaches can be identified from previous studies. The first approach is based on the Galerkin projection of known (linearised) flow equations onto a set of modes, in order to reduce the dimension of the system. These modes can hold valuable information, which can help both to understand the flow structure and to design the controller, for instance by guiding actuator-sensor placement. Additionally, they provide a way to reconstruct the full flow field from the reduced-order model (ROM) state. This in turn facilitates the interpretation of the effect of the controller on the flow and enables the use of full-state feedback algorithms (e.g. LQR).
Choosing the nature of the modes onto which the dynamics are projected is the first crucial step of this approach. 'Global modes' (e.g. Åkervik et al. 2007;Henningson & Åkervik 2008;Barbagallo et al. 2011;Ehrenstein, Passaggia & Gallaire 2011) and 'proper orthogonal modes' (e.g. Aubry et al. 1988;Podvin & Lumley 1998;Graham, Peraire & Tang 1999;Ravindran 2000aRavindran ,b, 2002Ravindran , 2006Prabhu, Collis & Chang 2001;Ma & Karniadakis 2002;Noack et al. 2003;Gloerfelt 2008;Siegel et al. 2008; Barbagallo et al. 2009Barbagallo et al. , 2012Tadmor et al. 2010) can yield successful ROMs, but 'balanced modes' have been shown to lead to superior performance in terms of robustness and required model order in a number of studies (e.g. Willcox & Peraire 2002;Ilak & Rowley 2008;Bagheri et al. 2009c; Barbagallo et al. 2009;Dergham et al. 2011). If the system is stable, such modes can be approximately identified at a low computational cost, by using a method (Willcox & Peraire 2002;Rowley 2005) referred to as approximate balanced truncation or balanced proper orthogonal decomposition (BPOD), which relies on snapshots from a set of impulse responses of the forward and associated adjoint system (one of each for SISO systems). This technique assumes the flow is evolving around a stable equilibrium point (base flow) and that the impulse response of the (linearised) forward and adjoint systems can be computed. It has been successful in modelling the dynamics of a large range of flows (e.g. Willcox & Peraire 2002;Rowley 2005;Bagheri et al. 2009a,b;Ahuja & Rowley 2010;Dergham et al. 2011;Semeraro et al. 2011;Tu & Rowley 2012).
One drastic limitation of BPOD is the fact that it is restricted to stable flows. Extensions have therefore been designed to circumvent this issue, but these require either the evaluation of all the unstable global modes of the flow field (Barbagallo et al. 2009;Ahuja & Rowley 2010) or the inversion of large matrices (Dergham et al. 2011), making them often computationally intractable. These restrictions have recently been lifted, as it was shown by Flinois, Morgans & Schmid (2015) that the snapshot-based approach can in fact be used directly to obtain ROMs of unstable flows, as long as the snapshots are chosen appropriately. This is therefore one of the approaches that is considered in this study. An outline of the BPOD method is provided in § A.1.

System identification
The second main approach for obtaining low-order linear models is system identification (for an overview, see Ljung 1987). In this case, only information collected by sensors and chosen actuator inputs is used to obtain the model. These methods are therefore more readily applicable in experiments than the projection approach introduced in § 1.2. Indeed, they do not rely on full-state information or knowledge of the governing equations, and hence do not require a linearised or an adjoint solver. Although no modal or sensitivity information is readily available, input-output models are sufficient for the purposes of controller design.
In a limited number of cases, it is possible to first model each physical phenomenon affecting the system independently. The resulting submodels can then be combined in order to create the full input-output model. For instance, Illingworth et al. (2012) considered the behaviour of the shear layer, the acoustics, the scattering and the receptivity of a cavity flow separately in order to construct an overall model. For most flow set-ups, however, this is not practical, and more general approaches are required. To this end, the relationship between the actuators and the sensors is often assumed to have a general mathematical description with unknown parameters. These parameters are then calibrated by minimising the error between a representative set of input-output data and the model predictions. This can be done either in the frequency domain or in the time domain. In the former case, the data is fitted directly to a chosen transfer function structure (e.g. Kwong & Dowling 1994;Zhu, Dowling & Bray 2005;Dahan et al. 2012). In the latter case, 'prediction-error methods' are often used, whereby the output is considered to be a linear combination of earlier inputs, outputs and a moving average of an error term. In flow control, only some of these components are typically necessary to obtain an accurate model, as for instance in Morgans & Dowling (2004), Zhu et al. (2005), Huang & Kim (2008), Hervé et al. (2012), Juillet, McKeon & Schmid (2014), Roca et al. (2014) and Gautier et al. (2015).
System identification can also be performed using 'subspace identification' methods, which result in an approximation for the system in state-space form directly. Subspace identification has been successfully implemented in some flow control studies (for instance Huang & Kim 2008;Juillet et al. 2013;Guzman Inigo, Sipp & Schmid 2014) and a review of these methods can be found in Qin (2006). Despite the more advanced mathematical framework, one of the advantages of this approach is that it yields noise covariance matrices that can be used in an LQG framework.
Finally, a method that has attracted an increasing amount of attention in recent years is the 'Eigensystem Realisation Algorithm' (ERA) introduced by Juang & Pappa (1985) and subsequently used in numerous flow control studies (e.g. Ma, Ahuja & Rowley 2011;Illingworth et al. 2011Illingworth et al. , 2012Dahan et al. 2012;Belson et al. 2013;Brunton, Rowley & Williams 2013;Dadfar et al. 2013;Brunton, Dawson & Rowley 2014;Illingworth, Naito & Fukagata 2014). It is straightforward to implement and results in state-space models that have recently been shown to be theoretically equivalent to those obtained with BPOD for stable systems (Ma et al. 2011). The focus of this article is the ERA, and an outline of the method is thus included in § A.2. As with most system identification techniques (and BPOD), the ERA was originally only designed for stable systems. One solution is to first stabilise the system with an ad hoc controller, then identify the (stable) closed-loop dynamics, from which the plant's open-loop (unstable) dynamics can finally be extracted. This approach has been successful in some studies (e.g. Morgans & Dowling 2007;Illingworth et al. 2011Illingworth et al. , 2012, but is not always practical, as it can be challenging to find a stabilising controller without a model. Here, we instead use the theoretical equivalence of the ERA and BPOD and the fact that BPOD can be applied to unstable systems to build ERA models directly from the impulse response of unstable systems.

Article structure
In this paper, we build on the findings of Ma et al. (2011) and to show that the ERA is a practical and computationally cheap approach to obtain low-order models of unstable flows that are useful for feedback control. We confirm that such models are in practice very similar to the ones obtained with BPOD and also show that they can be obtained directly from the nonlinear system dynamics, without the need for a linearised or an adjoint solver. As a result, the method may even be applicable in some experimental set-ups.
After introducing the numerical approach and set-up in § 2, we consider the unstable flow over a D-shaped body and compare three different ROMs in § 3: the first is obtained with BPOD ( § 3.2), the second with the ERA ( § 3.3), and the third also with the ERA but based on the flow's nonlinear dynamics ( § 3.4). We show that the three models are very similar, as expected, and then design controllers based on the last set of models using proportional control and H ∞ loop-shaping in § 4. We show that the controllers are able to stabilise the nonlinear flow, even from the vortex shedding limit cycle and in some cases at off-design Reynolds numbers. We compare two sensor configurations, where in one case the sensor is located in the wake and in the other the sensor is body-mounted. Concluding remarks are included in § 5.

Numerical set-up
The work presented in this article is based on the fast multi-grid immersed-boundary fractional step (IBFS) algorithm introduced in Taira & Colonius (2007) and Colonius & Taira (2008)  FIGURE 1. Input/output and flow set-up for the simulations. The actuators (dark grey) are disk-shaped body forces with a diameter of 10 % of the body height, acting in opposite horizontal directions as shown by the arrows. The first sensor measures the vertical velocity of the flow along the symmetry plane, two body heights downstream of the trailing edge (black triangle). The second sensor (light grey) is a distributed sensor, which measures the antisymmetric component of the force acting on the base of the body. and used in many studies (e.g. Taira & Colonius 2007;Colonius & Taira 2008;Ahuja & Rowley 2010;Joe et al. 2010;Joe, Colonius & MacMynowski 2011;Ma et al. 2011;Tu & Rowley 2012;Brunton et al. 2013Brunton et al. , 2014Choi, Colonius & Williams 2015;Flinois & Colonius 2015) so only an overview of the formulation of this code is given here. The reader is referred to the cited studies and the references therein for further details.
The incompressible Navier-Stokes equations in vorticity form are considered and discretised on a uniform Cartesian grid. A set of discrete immersed-boundary forces is defined along the surface of the body and regularised onto the Cartesian grid in order to enforce the no-slip boundary condition on the surface. The flow equations are solved sequentially on a number of nested grids at each time step. Each grid is identical to the one it is nested into, but has half the physical extent and hence is twice as fine. The flow state on the next grid level is used to obtain boundary conditions at the edges of each grid. For the largest grid, a far-field boundary condition is used (zero vorticity). In order to solve for the linearised and adjoint dynamics of the flow, the unstable steady state is identified using selective frequency damping (Åkervik et al. 2006;Jordi, Cotter & Sherwin 2014), and the equations are linearised about this state. We emphasise that no time-averaged flow is used in this study: all linearised and adjoint simulations evolve about an unstable steady base flow, which is a solution of the Navier-Stokes equations. The time-continuous adjoint equations are then derived from the spatially discretised vorticity equations. The linearised and associated adjoint equation solvers are based on the same unstable base flow and both use the same discretisation and time-marching schemes as the nonlinear equations.
The flow around a D-shaped body is used as a test problem in this article. The body is in the shape of a half-ellipse with a blunt vertical base. The length of the body is twice as long as its height (see figure 1). The Reynolds number Re = U ∞ H/ν (where U ∞ is the incoming flow velocity, H is the body height and ν is the kinematic viscosity) is 80, corresponding to an unstable flow. A time step of 0.005 convective time units was used; the finest grid is of dimensions 30H × 10H and contains 1500 × 500 grid points. Three nested grids were used in total. The trailing edge of the body is located at the centre of the finest grid since the 'flow' is advected upstream in the adjoint simulations. Note that, although the body shape is different, the grid used is the same as in Tu & Rowley (2012) and , where the flow over a cylinder was considered in both cases.
In order to measure the fluctuations in the wake, two sensor configurations are considered. The first measures the vertical velocity at a point located along the symmetry plane, two body heights downstream of the trailing edge. The second is a body-mounted, distributed sensor that measures the antisymmetric component of the force on the base, as illustrated in figure 1. The first sensor is used in § 3 to demonstrate the modelling procedure, and both sensors are used in § 4, in order to compare the performance of the resulting closed-loop systems. The flow is actuated using two disk-shaped horizontal body forces of diameter 0.1H and centred 0.15H above and below the trailing edge corners, directly in line with the base. The two forces are of equal magnitude but opposite sign -thus this is a SISO system. The input forces are acting purely in the x direction and the total force on the flow is always zero but with a non-zero antisymmetric component.

Comparison between modelling approaches
In this section, three different approaches are used to obtain ROMs for the input-output dynamics between the actuator and the wake velocity sensor described above. First, the linearised forward and adjoint impulse responses are computed and the projection-free snapshot-based balanced truncation (BPOD) approach is used to obtain a ROM. We use this as our starting point, as this method has been shown  to yield balanced ROMs for unstable systems. Second, the ERA is applied to the linearised impulse response and we show that we obtain a very similar model to the BPOD-based ROM. Finally, we show that another effectively identical model can be obtained for the linearised dynamics about this steady state, simply by considering the early response of the nonlinear system to a small impulsive input. This is a key point, as it shows that it is possible to obtain ERA models for unstable systems with standard (nonlinear) computational fluid dynamics codes and even potentially with experiments.

Unstable steady state and unforced flow
The flow over the D-shaped body considered in this article was also studied computationally by Palei & Seifert (2007), who identified the critical Reynolds number of the unforced flow to be Re c ≈ 63. As the simulations in the present work focus on Re = 80, only one complex conjugate pair of eigenvalues is expected to be unstable, since this is only slightly higher than the critical value (it will become apparent that this is indeed the case in later sections of this article). In order to validate the mesh used for our simulations, the vortex shedding Strouhal number (St = fH/U ∞ , where f is the dimensional frequency) and the root mean square (r.m.s.) value of the lift coefficient fluctuations C L were compared to those found in Palei & Seifert (2007) for several Reynolds numbers. The results are summarised in table 1. The agreement is very good, and this provides a good indication that our grid is sufficiently well resolved at all the considered Reynolds numbers.
The unstable base flow at Re = 80, obtained using selective frequency damping (Åkervik et al. 2006;Jordi et al. 2014), is compared to the unforced fully developed vortex shedding flow field (also at Re = 80) in figure 2. The length of the recirculation bubble in the base flow was found to be x rec = 3.58H and the drag coefficient was found to be C D = 1.138. The mean drag coefficient of the unforced vortex shedding flow was found to be C D = 1.234.

The BPOD model
In order to compute the balanced truncation ROM, the impulse response of the linearised forward and adjoint systems were computed (with the wake velocity  3.2.1. Actuator and sensor placement As discussed for instance in Bagheri et al. (2009c) and Chen & Rowley (2011), a feedback controller can only stabilise an unstable flow if all unstable modes are controllable and observable using the chosen actuators and sensors, respectively. In other words, the spatial location of the actuators and sensors must overlap with the spatial distribution of the adjoint and forward global modes, respectively. Furthermore, the flow structures that are generated with the least amount of input energy by the actuators are given by the leading eigenmodes of the controllability Gramian. It is therefore important to ensure that sensors can measure the state of the system where the magnitude of these modes is large. Conversely, the flow structures that lead to the most energetic output signal correspond to the leading eigenmodes of the observability Gramian. It is therefore important to ensure that actuators can influence the state of the system where the magnitude of these modes is large. Finally, the leading balanced modes inform us about the states that contribute the most to the input-output dynamics of the flow.
In the case of unstable systems, it was shown in  that, using projection-free BPOD (and assuming the unstable modes are observable and controllable), the direction of the unstable balanced modes tends to the direction of the unstable eigenmodes of the Gramians as t ∞ → ∞. Additionally, the subspace spanned by these modes converges to the subspace spanned by the unstable global modes (note that the modes themselves are not identical in general in this case). Physically, this corresponds to the fact that the behaviour of unstable modes is eventually dominant for an uncontrolled linear system, as the influence of any initial transient growth becomes negligible as t ∞ → ∞. In figure 4, we therefore show the first two primal and adjoint balanced modes, identified with a large final simulation time of t ∞ = 100 convective time units. Note that the ERA and most other system identification methods do not usually provide this type of full-state or adjoint information. Figure 4(c,d) confirms that there is an overlap between the chosen actuators and the unstable modes and figure 4(a,b) confirms that the location of the vertical velocity sensor in the near wake will allow it to measure the growth of the unstable modes. Figure 4(a,b) also shows that, qualitatively, the unstable modes are dominant in the wake, downstream of the body. This suggests that the base (rear face) of the body is a justified location for a body-mounted sensor.
However, the leading global modes and Gramian eigenmodes do not provide sufficient information to select actuators and sensors that guarantee good robustness and performance in a closed-loop setting, especially in the presence of disturbances (Chen & Rowley 2011). On the one hand, considering global modes individually (even unstable ones) is not adequate for highly non-normal systems since the input-output dynamics may be strongly dependent on the interaction between several modes (Sipp et al. 2010). On the other hand, Gramian eigenmodes decouple the effect of actuators and sensors: the controllability Gramian is computed without any knowledge of the sensor matrix C, while the observability Gramian is computed without any knowledge of the actuator matrix B. Gramian eigenmodes thus cannot take into account the effect of time delays between actuation and sensing, which can have a large impact on the performance of feedback control (Chen & Rowley 2011). In the case we are considering here, placing the velocity sensor further downstream of the body than x = 2 would provide a better observability of the unstable modes, as indicated by figure 4(a,b). However, in practice, it was found that the design of robust stabilising controllers became increasingly difficult as the sensor was shifted further downstream.
In order to find a compromise between these conflicting requirements, Chen & Rowley (2011) suggested that a 'structural sensitivity' analysis, resulting in the 'wavemaker' region introduced by Giannetti & Luchini (2007), might provide a promising starting point for actuator-sensor placement in a closed-loop setting. This region is identified by taking the pointwise product of the forward and adjoint global modes. Using projection-free BPOD, however, the unstable balanced modes only converge to a subspace spanned by the global modes. We therefore show the average of the two wavemaker regions predicted by the unstable modes in order to estimate the actual wavemaker of the system in figure 5. It is encouraging that the identified region is in good qualitative agreement with the one resulting from the structural sensitivity analysis performed by Giannetti & Luchini (2007) for the circular cylinder at similar Reynolds numbers. Figure 5 shows that the wake sensor is at the centre of the wavemaker region and also predicts that sensors located further downstream than x ≈ 4-5 may not perform well. Figure 5 also shows that the actuators and the force sensor are located as close as possible to the region where the structural sensitivity is high, while remaining 'body-mounted'.

Reduced-order model
As explained in , it is necessary to select an adequate final simulation time t ∞ in order to strike a balance between the convergence of the model and the exponential growth of the unstable modes. time t ∞ increases, while the other singular values converge. HSVs of low dynamical importance are also affected by t ∞ , as described in : as the unstable modes start to dominate the response, information related to these less significant modes is lost. The corresponding singular values then increase as a whole with t ∞ as opposed to converging. For stable systems, there exists an upper bound on the H ∞ norm of the difference between the full-order model G n (where n is the full system order) and the reducedorder model G r (of order r) obtained with balanced truncation (e.g. Rowley 2005): where σ j is the jth HSV of the system. This bound holds for exact balancing procedures but can in practice be observed by using a sufficient number of snapshots in the BPOD approach. In the case of unstable systems, the H ∞ error norm is infinite by definition. An upper bound on the H ∞ error of the system's stable subspace can nevertheless be recovered if a projection approach is used to first identify and project out the unstable subspace. This can be done either by first identifying the unstable modes as in Barbagallo et al. (2009) and Ahuja & Rowley (2010), or by using the extension to the projection-free BPOD approach suggested in , whereby a first set of impulse responses is computed to identify the unstable balanced modes only. This is not performed in the present work. In the present case, it was found by trial and error that an accurate 10th-order ROM could be obtained with only t ∞ = 25 convective time units. In figure 7, we show the impulse response of the model as a thin solid line, superimposed on the full impulse response (thick dashed line). Clearly, the long-term response is well very predicted, even for much longer times than the final simulation time used to identify the model. The transients are also well approximated by the model. Figure 8 is a Bode plot showing the transfer function of the ROM as a thin solid line. It can be used to design controllers using classical or H ∞ loop-shaping methods. This plot is also useful to compare this model with the ROMs obtained in § § 3.3 and 3.4 in the frequency domain.

The ERA model
The construction of ERA models is computationally very efficient as it relies only on the sensor measurement of the impulse response of the forward linearised equations. The procedure for obtaining such models is described in § A.2. As with BPOD, this method was originally designed exclusively for stable systems since the output becomes unbounded as t ∞ → ∞ if the system is unstable. However, as shown by Ma et al. (2011) and summarised in § A.2, the ERA and BPOD are essentially the same algorithm and can thus theoretically yield exactly identical models. Furthermore, as shown in , the snapshot-based BPOD method also leads to balanced ROMs for unstable systems. As a result, the ERA must also lead to balanced ROMs, even for unstable systems.
In order to construct a ROM that is comparable to the BPOD ROM of § 3.2, every 40th sensor measurement and a final simulation time of t ∞ = 50 convective time units were used (this matches the size of the Hankel matrices on which the two models are based, as shown in § § A.1 and A.2). The computational cost for obtaining the ROM is negligible, since we are simply stacking the output signal from the sensor into two matrices and performing a singular value decomposition of a (125 × 125) matrix. Much larger Hankel matrices could thus be used at virtually no additional expense if required.
Using again a 10th-order model, the ERA ROM impulse response and transfer function are compared with the BPOD ones in figures 7 and 8, respectively (dashed lines). Clearly there is a good match between the two models, in both the time and frequency domains. We provide a quantitative measure of how similar the two models are in § 4.
The models are not exactly identical, mainly due to small numerical and algorithmic differences in their implementation. The adjoint simulations are not the exact discrete adjoint of the forward simulations due to the multi-grid solver, which is not self-adjoint despite being used in both solvers, as noted by Ahuja & Rowley (2010). Additionally, the continuous version of BPOD was used in the present work and the forward and adjoint responses were not computed exactly in the same manner: in the adjoint simulation, the initial state was set to z(1) = C † and this state is shown in figure 3(b). On the other hand, in the forward simulation, the B matrix was approximated by a one-step pulse at the first time step of the forward simulation. This was done partly for convenience and partly because the same simulation was used for the BPOD and ERA models, where we attempt to keep the procedure pseudo-experimental. The flow state, just as the pulse is applied, is shown in figure 3(a) in this case.
3.4. ERA model based on nonlinear dynamics 3.4.1. Obtaining the reference model Using the ERA with unstable systems as above allows balanced ROMs to be obtained for unstable systems, without the need for an adjoint solver. However, in many relevant flow control scenarios, the linearised impulse response cannot be computed directly either. This is the case if a linear solver is not available or in experimental set-ups. As a third approach, therefore, we show that ERA models can be readily obtained directly from the nonlinear code: if the impulse is sufficiently small, the early response of the system is approximately linear until disturbances grow enough for nonlinear effects to become significant. It is therefore theoretically possible to obtain an arbitrarily long (approximately) linear response by imposing a correspondingly small impulse. Of course, in any realistic situation, having an excessively small impulse will result in significant issues due to the associated low signal-to-noise ratio of the output signal. In practice, an excessively long impulse response is not desirable anyway for unstable flows, as t ∞ is limited by the exponential growth of the unstable modes.
Nevertheless, if a sensor measurement can be recorded with a large enough portion in the linear regime and an acceptable signal-to-noise ratio, it is straightforward to scale the response back to recover the output that would result from a full impulse applied to the linearised system. The results are shown again in figures 7 and 8 as dot-dashed lines for a ROM obtained with exactly the same procedure as in § 3.3, but scaled by a factor of 400 to retrieve the full-sized signal from the small impulse response of the nonlinear flow. The two ERA models are indistinguishable in both the time domain and the frequency domain. The model obtained in this section was thus used as the nominal or reference model for controller design in § 4.

Effect of nonlinearity
In a linear setting, the effect of letting t ∞ → ∞ on the quality of balanced models was investigated in . In a nonlinear setting, the appearance of a limit cycle (or other nonlinear effects) may impose a more stringent upper bound on t ∞ than the exponential growth of unstable global modes if the applied impulse is large. It is therefore relevant to examine how nonlinearity affects the ERA models. To this end, the system is subjected to a large impulse (200 times greater than the one used to obtain the model above) in order to promote the saturation of the unstable modes. Figure 9(a) shows the evolution of the sensor signal in both the linear and nonlinear simulations.
The main purpose of obtaining reduced-order models in this study is to use them to design stabilising feedback controllers for the flow. The models we obtain are therefore only useful if they are able to accurately represent the behaviour of the unstable subspace of the actual system (this is a necessary but not always a sufficient condition). The influence of the final simulation time on the extent to which this requirement is satisfied is shown in figure 9(b,c): recalling that the considered linearised system has two unstable modes, we analyse how the real and imaginary parts of the poles of the second-order ERA model obtained from the nonlinear system's response (to the large impulse) are affected by t ∞ . For approximately the first 50 convective time units, both the real and imaginary parts of the poles converge to the values obtained with the linearised impulse response. For larger values of t ∞ , the growth rate of the poles tends to zero and the frequencies tend to the limit-cycle frequency. Note that the maximum acceptable t ∞ value depends on the size of the impulse applied to the system.
In figure 10, we compare the transfer function gain of the ERA models from figure 9(b,c) for several t ∞ values to the 'reference' transfer function gain, corresponding to the 10th-order model described above ( and (c), the growth rate and frequency of the converged linearised system poles (dashed line) and of the saturated limit cycle (dot-dashed line) are also shown. Note that the two poles are a complex conjugate pair, so they have opposite frequencies and identical growth rates.
frequency range of the unstable modes. As t ∞ is further increased, the frequency of the peak tends to the limit-cycle frequency and it becomes increasingly sharp as the growth rate tends to zero. As mentioned above, if the system's unstable poles are not accurately identified by the model, it cannot be expected to be of satisfactory quality for controller design. Figure 9 shows that the poles move away from their 'converged' values shortly after the linear and nonlinear responses begin to differ. Unlike in the linear setting, the saturation of unstable modes affects the (dominant) unstable modes first as opposed to the least significant modes, so the loss of accuracy of the models can be expected to be more sudden in this case, even if a higher-order model is used. This therefore suggests that the ERA should only be applied in the approximately linear portion of the impulse response in order to obtain a useful model for controller design purposes. As the maximum acceptable t ∞ value for a given impulse response can be identified directly from figure 9, it is relatively straightforward to enforce this condition.

Feedback control
In this section, we use the (linear) model obtained from the nonlinear impulse response in § 3.4 (converted to continuous time) to design linear controllers using H ∞ loop-shaping ( § 4.1) and proportional control ( § 4.2). We also show how robust control tools can be used to analyse model uncertainty and controller robustness in § 4.3. Although convenient numerically, velocity sensors in the middle of the wake are not always realistic. In § 4.4, we therefore also design controllers using the body-mounted force sensor on the base of the body (as described in § 2) and compare the results obtained with the two sensor configurations in § 4.5.
The block diagram for the control set-up is shown in figure 11. It consists of a linear system, whose input signals u(s) (s is the Laplace variable) and output signals y(s) are related by the model G(s), the system's open-loop transfer function: y(s) = G(s) u(s). In closed loop, the measurements y(s) are used to automate the actuation u(s) using a mathematical law K(s). This controller can be designed for instance to improve the system's tracking performance or disturbance rejection capabilities. In the case of an unstable G(s), the controller can be used to robustly stabilise the system. The full closed-loop transfer function from r(s) to y(s) is . (4.1) The modelling procedure from the previous sections makes it possible to design controllers for unstable systems in a systematic way, without prior flow stabilisation. Note also that the entire procedure only requires a standard nonlinear solver: the linear model is obtained from an approximately linear impulse response, the controller is designed using this model, and it is applied directly to the nonlinear flow. Although the controllers are also applied to the linearised Navier-Stokes equations in § 4.5.1, this is only done to compare their response with theoretical predictions obtained with the ROM, and hence this step is not strictly necessary.

H ∞ loop-shaping
As mentioned in § 1.1, the method chosen here is a two-step procedure. The end goal is to obtain a robust controller that stabilises the system despite noise and/or disturbances, nonlinearities and model inaccuracies. Ideally, we would like to be able to stabilise the system from an off-design state like the unforced vortex shedding wake or at increased Reynolds numbers. H ∞ loop-shaping provides a framework for obtaining this type of robust controller (McFarlane & Glover 1992).
The H ∞ loop-shaping methodology requires the user to first manually 'shape' the gain of the open-loop frequency response or transfer function using dynamic weights W(s) that ensure the resulting response gain has the desired properties. In the second step, an optimisation problem is solved (using MATLAB), which generates a controller K ∞ (s) that modifies mainly the phase of the shaped plant in order to make the closedloop system robust to a large class of disturbances. Finally, the two components are joined together to form the final controller: K(s) = K ∞ (s)W(s). Further details are included in appendix B.

Manual loop-shaping step
For this system, the first iterative step was used to improve the robustness of the final controller as much as possible and to ensure the shaped model was stabilised in a desired manner. Without this step, it was found that H ∞ optimisation often returned stabilising controllers that were themselves unstable. Although this is not necessarily an issue, it should be avoided if possible for practical reasons.
The result of this procedure is shown in figure 12 (where the minus signs are due to the positive feedback sign convention used in figure 11). The thick solid line shows the Nyquist diagram of the open-loop model G(s), scaled by a factor of −21 (for reasons explained below and in § 4.2). The Nyquist diagram is the locus of a transfer function in the complex plane for values of s = iω and for all real frequencies ω. The thin dashed line is the Nyquist diagram of the shaped plant −WG. In this case, W(s) is a first-order lag compensator.
The standard stability criterion for a system with two unstable poles is that, following the Nyquist curve from ω = −∞ to ω = +∞ rad s −1 , the point with coordinates (−1, 0) should be 'encircled' twice in the anticlockwise direction for the controller to stabilise the plant. Figure 12 shows that this is the case for the scaled open-loop system, which shows that there exists a so-called 'gain window', where proportional control stabilises the ROM. We therefore also implement a proportional controller and compare its performance with the robust controller. The design of the proportional controller is discussed in § 4.2. The shaped system also stabilises the ROM in this case, although we emphasise that this is only an intermediate step and that the aim of this initial procedure is to improve the robustness of the final controller iteratively and not to obtain the best possible controller using only classical loop-shaping.

Robustness optimisation step
The aim of the second step is to optimise the shaped system's robustness, as quantified by the generalised stability margin b, defined formally in appendix B. If the controller does not stabilise the plant, then b = 0, and as the robustness increases, b → 1. A minimum value of 0.2 b 0.3 is usually considered to be acceptable. As implied by its name, b is a generalisation of the standard gain and phase margins used in classical control. It gives a measure of both the robust performance and robust stability characteristics of a given (shaped) model and is also applicable to multiple-input-multiple-output systems.
The shaped plant WG was therefore used as an initial condition for the H ∞ loopshaping algorithm. This resulted in a robust 11th-order controller K 0 (s) = K ∞ (s)W(s). High-order controllers are less easily implemented in practice and also less reliable, so it is desirable to keep the order of the controller as low as possible without excessively compromising the closed-loop performance. The K ∞ part of the controller was thus subsequently reduced using MATLAB's balreal and modred commands, resulting in a final eighth-order controller K(s) = K ∞,r (s)W(s). The generalised stability margins for the optimised and reduced controllers are b(WG, K ∞ ) = 0.427 Feedback control of unstable flows using ERA models 59 and b(WG, K ∞,r ) = 0.422, respectively, which shows first that the robustness of the H ∞ controller is satisfactory, and second that the order reduction process only deteriorated this robustness very marginally.

Proportional control
In order to compare the H ∞ loop-shaping approach to the proportional control approach used in several classical studies such as Roussopoulos (1993) and Park, Ladd & Hendricks (1994), we also design a proportional controller: u(s) = K P y(s) if r(s) = 0, where K P is simply a constant. As mentioned above, it was found that there exists a 'gain window', where proportional control can be expected to stabilise the ROM. In a similar manner to Illingworth et al. (2014), we use our ERA model to select an adequate K P value. Here, we choose to maximise robustness, as quantified by b(K P G, 1), in order to show that this stability and performance parameter is useful even for very simple controller design strategies. The optimal value was found to correspond to K P ≈ 21, for which b(K P G, 1) = 0.180. This value is smaller than b(WG, K ∞r ) = 0.422, indicating that the H ∞ loop-shaping controller can be expected to be more robust than the proportional controller. This is further discussed in the next section.

Controller robustness and model uncertainty
In order to obtain the reference model that was used for controller design in § § 4.1 and 4.2, we selected one of three modelling algorithms, a final simulation time of t ∞ = 50 convective time units and a model order of r = 10. In this section, we show that robust control tools can quantify the model uncertainty associated with these choices. A related approach was used by Jones et al. (2015) in a grid refinement study.
A useful tool that naturally fits in the H ∞ loop-shaping framework is the ν-gap metric introduced by Vinnicombe (1993) and denoted by δ ν . It provides a measure of how differently two systems shaped with the same weights W(s) can be expected to behave in closed loop. The ν-gap takes a value between 0 and 1 (δ ν = 0 between identical models) and it is formally defined in appendix B. This metric can be used to provide stability guarantees for uncertain systems. Specifically, if b(WG, K ∞ ) > δ ν (WG, WG i ) for given W and K ∞ , and for a 'perturbed' model G i , then the controller K 0 = K ∞ W also stabilises G i . For any perturbed model where b(WG, K ∞ ) < δ ν (WG, WG i ), the stability of the perturbed closed-loop model can be checked by calculating b(WG i , K ∞ ).
As summarised in table 2, several 'perturbed' open-loop models G i were constructed by varying the model order, final simulation time and modelling method used. For all models, the ν-gap metric corresponding both to the H ∞ loop-shaping weights δ ν (WG, WG i ) and to the proportional controller δ ν (K P G, K P G i ) were calculated and are shown in table 2.
The ν-gap values in table 2 demonstrate the advantage of using H ∞ loop-shaping: although the ν-gap associated with all models is similar for both controllers, b(WG, K ∞ ) > δ ν (WG, WG i ) for all i, whereas this is not the case with the proportional controller. The robust controller is thus guaranteed to stabilise all the models considered here. On the other hand, it was necessary to verify that the proportional controller in fact does stabilise G BPOD by checking that b(K P G BPOD , 1) > 0 since b(K P G, 1) < δ ν (K P G, K P G BPOD ). It was found that b(K P G BPOD , 1) = 0.144. reference model, whereas larger differences appear with lower t ∞ and r values. This suggests that our chosen t ∞ and r values are allowing us to obtain a converged model and resolve most of the important dynamics of the input-output system. Similarly, the distance between the two ERA models is negligible, which confirms that the early impulse response of the nonlinear flow provides adequate information to obtain a useful ERA model. Finally, a non-negligible value of δ ν is obtained with G BPOD . This indicates that the ROMs may be sensitive to the modelling approach in this case and that, with the chosen weights W and K P , some differences can be expected between the responses of the two models in closed loop. Nevertheless, the model G BPOD is still stabilised by both controllers. We emphasise that these tests are all performed a priori, in the sense that no closed-loop simulation is required. They suggest that the full-order flow and the ROMs are likely to behave in a similar manner since the stability and performance of the closed loop are suitably robust to the modelling technique and model parameters. This analysis is equally instructive for controllers designed explicitly using robust control methods such as K(s) and for much simpler controllers such as K P .

Body-mounted force sensor: modelling and controller design
The same procedure was then followed using the body-mounted force sensor introduced in § 2: first a 10th-order ERA model was obtained by scaling the response of the full nonlinear flow field to a small impulse. Note that the two sensor signals were recorded at the same time so it was not necessary to run a new simulation to obtain this model. In fact, the measurements from an arbitrary number of candidate sensors can be stored simultaneously for a given actuator configuration. As a result, there is effectively no added cost associated with the construction and comparison of models obtained with a large number of sensors. An ERA model was therefore also obtained from the fully linear response using this sensor.
A robust stabilising controller was then designed by shaping the Nyquist diagram with a first-order lag filter and a first-order low-pass filter. In this case, we note that, at this intermediate stage in the robust controller design process, the chosen weights W(s) in fact do not stabilise the open-loop model. The robustness of the shaped system WG was then improved using H ∞ loop-shaping. The stability margins corresponding to the resulting 14th-order controller and the eighth-order reduced controller are b(WG, K ∞ ) = 0.535 and b(WG, K ∞,r ) = 0.526, respectively. The corresponding Nyquist diagrams are shown in figure 13. Note that the robust shaped systems K 0 G and KG have a peak at ω ≈ 0.78 rad s −1 so their Nyquist diagrams have a large magnitude close to this frequency. Using the force sensor configuration, it was noticed that the range of stabilising proportional gains was greatly reduced due to the fact that proportional control triggered instabilities at frequencies close to the time-stepping frequency. In order to circumvent this numerical issue, the proportional controller was first coupled with a simple second-order low-pass filter. The filter has a cutoff frequency of ω = 200 rad s −1 , and hence has no effect on the frequencies where physical phenomena take place in the flow. The final controller can therefore still be considered to be effectively a proportional controller. With this added filter, it was possible to identify a 'gain window' where a stability margin of b(K P G, 1) = 0.182 was obtained by using a scaling factor of 555. The model shaped with this controller K P is also shown in figure 13. Table 3 shows the ν-gap corresponding to the robust controller weights W and the proportional controller K P for all considered models for the force sensor. The results are similar to those in table 2, although the ν-gap values are overall slightly higher. As with the velocity sensor, b(WG, K ∞,r ) > δ ν (WG, WG i ) for all G i with the robust controller. With the proportional controller, on the other hand, b(K P G, 1) < δ ν (K P G, K P G t ∞ =25 ). The stability of the perturbed closed-loop model was thus checked and it was found that b(K P G t ∞ =25 , 1) = 0.202 > 0.
This analysis again confirms that the closed loop is not excessively sensitive to the chosen modelling technique and model parameters. We can therefore confidently proceed to implement the controllers in the full-order flow. Name  TABLE 3. Parameters used to generate the considered family of 'perturbed' models G i for the body-mounted force sensor, and values of the corresponding ν-gap metric, calculated for both the robust controller weights and the proportional controller.
4.5. Controller performance Equipped with robust stabilising controllers, three increasingly demanding tests were performed for the four considered control configurations. The results are presented in this section.

Closed-loop response predictions
First, the closed-loop response predicted by each ROM was compared to that of the full linearised system: the flow was forced with an impulse and allowed to evolve for 50 convective time units, thus allowing the unstable modes to develop before the controllers are switched on. In figure 14 the response predicted by the ROM (velocity sensor) is compared to the full linearised system's response in closed loop. With both controllers, the two responses match very closely, even after the controller is switched on. Note also that the proportional controller leads to a slower stabilisation of the flow in this case. Similar conclusions can be drawn from the results obtained with the body-mounted force sensor, shown in figure 15. In this case, the convergence rate of the output signal is lower and it therefore takes longer for the flow to return to its unstable equilibrium state. On the other hand, the initial response of the controller is less aggressive than with the wake velocity sensor. Again, the stabilisation of the flow is faster with K than with K P .

Robustness to nonlinear initial conditions
The robustness of the controllers to nonlinear effects was then tested. After triggering the instability with an impulse, the nonlinear flow was left to evolve to its limit-cycling (vortex shedding) state. The controllers were only switched on when vortex shedding was fully established (after 200 convective time units). We emphasise that neither the models nor the controllers incorporate any information about the limit-cycle dynamics, which do not even necessarily evolve about the same state as the base flow.
However, as shown in figures 16 and 17, the two proportional controllers and the two H ∞ controllers are robust enough to stabilise the flow. Figure 18 shows the evolution of the forces on the body in all four configurations: starting from their base flow values, the instability is triggered and the forces evolve as the flow reaches the limit cycle. Once the controllers have been switched on, the forces quickly return to the base flow values.
In fact, the entire closed-loop flow returns to a state that is close to the low-drag unstable steady state in all cases. Figure 19 shows snapshots of the flow field (controlled using the robust controller associated with the body-mounted force sensor) at the moment when the control is switched on (after t = 200 convective time units), and then as the flow is gradually stabilised. After 500 convective time units (figure 19d), the vortex shedding is effectively fully suppressed and the flow is similar to the base state in figure 2. For the wake velocity sensor, the convergence of the forces takes place at a similar rate with both controllers. As in the linear framework in § 4.5.1, the flow takes longer to converge to its steady state with the body-mounted force sensor configuration. However, in this case, the proportional controller leads to a faster flow stabilisation than the robust controller.
Overall, these results thus show that the designed controllers are robust to significant nonlinear effects for all four configurations.

Robustness to changes in the Reynolds number
The controllers were then challenged further by repeating the test from § 4.5.2 but with larger (and hence off-design) Reynolds number values of Re = 90 and 100. This introduces additional modifications in the dynamics of the flows that the controllers are attempting to stabilise. In particular, the flow becomes more unstable at higher Re values and is hence harder to control. It was found that, using the body-mounted force sensor, neither controller was able to fully stabilise the flow at Re = 90 or 100. With the velocity sensor configuration (and with both of the associated controllers), the flow was fully stabilised from the limit cycle at Re = 90 and 100. Further simulations were run at Re = 110, 130 and 150: the flow was found to be stabilised from the limit cycle at Re = 110 with the robust controller and almost fully with the proportionally controller. At Re = 130, neither controller fully stabilised the flow.
Even in cases where the flow was not fully stabilised, the control was found to result in a positive net energy gain, as measured by the total average power required to drive the body through the fluid (after transients). The instantaneous total power can be expressed as the sum of the power required to overcome the drag D(t) and the power supplied to the fluid by the actuators (see e.g. Frohnapfel, Hasegawa & Quadrio 2012). In non-dimensional form, this can be written as where P total (t) and P control (t) are the dimensional total and control power per unit width, respectively, ρ is the fluid density, A F is the area where the forcing is applied, f is the local dimensional force per unit area and u is the local dimensional velocity. Note that, in practice, an actuator model would be desirable in order to obtain a more realistic estimate of the required input power. In our case, the control power (second term) tends to zero when the flow is fully stabilised. Moreover, in all cases where the control input converges to a finite value, it was found that the time-averaged value of P control after transients was at most of the order of a hundredth of a per cent of P total and usually much less. This negligible value is due to the fact that the flow is almost symmetric at the actuation location while the forcing is purely antisymmetric. As a result, changes in the converged drag coefficient can be directly linked to changes in the total required power. The results are shown in figure 20: clearly in all cases, the required time-averaged converged power is reduced significantly. For the wake velocity sensor, the robust controller leads to a great power reduction than the proportional controller. This general trend seems to be inverted for the body-mounted force sensor.
The results from this test suggest that the velocity sensor configuration is in general more robust to an increase in the Reynolds number than the force sensor configuration, despite the fact that the stability margin is larger with the body-mounted sensor. This implies that the dynamics of the shaped system experience larger changes with the force sensor than with the velocity sensor. As the Reynolds number is increased, Giannetti & Luchini (2007) have shown that the wavemaker region in a cylinder wake at low Reynolds numbers becomes elongated and shifts further downstream. Since the body-mounted sensor is located in a less observable region of the flow (see figure 4) and is already on the upstream edge of the wavemaker region at the design  FIGURE 20. Comparison of the total power required to drive the body through the fluid in the unstable equilibrium flow state, the unforced vortex shedding state and with the two controllers, for cases where the flow is not fully stabilised with (a) the body-mounted force sensor and (b) the wake velocity sensor.
Reynolds number (see figure 5), this could explain why the associated input-output dynamics are more strongly affected by an increase in the Reynolds number.

Further comments
In § 4.3, it was found that H ∞ loop-shaping yields controllers that are expected to be considerably more robust than proportional controllers (as quantified by b). Additionally, the window of stabilising proportional gains eventually disappears for sufficiently large Re and in this case a dynamic controller is necessary to stabilise the flow (e.g. Illingworth et al. 2014). Nevertheless, the results from this section indicate that, in cases where a gain window does exist, choosing the proportional control gain that maximises the generalised stability margin can result in controllers whose robustness and performance are in practice similar to that of H ∞ loop-shaping controllers. The generalised stability margin b therefore appears to be a useful parameter to guide the design of stabilising controllers for fluid flows in general. We note however that b should not be used blindly: although the force sensor robust controller provided in the largest stability margin, this system was not the most robust to an increase in the Reynolds number.
Compared to other studies where proportional and/or classical control was used (e.g. Roussopoulos 1993;Park et al. 1994), the availability of an accurate input-output model enables the use of b to select a proportional controller gain value that provides good robustness. This information is not available when tuning a controller manually.
It is also important to note that the controllers used in this study were designed in part manually. This iterative procedure is different for each actuator-sensor configuration and may be better optimised in this case with one of the sensor configurations than the other. As a result, the comparison between the two sensors cannot be made quantitatively.
Additionally, recall that the purpose of the simulations performed here was to test the robustness of fixed controllers as the Reynolds number takes on larger off-design values. The aim was therefore not to identify the highest value of the Reynolds number for which flow stabilisation is possible, as optimal control techniques would be better suited for such a study (see for instance Lauga & Bewley 2003).

Conclusions
This paper has shown that it is straightforward and justified to obtain linear loworder models for the dynamics of fluid flows about an unstable equilibrium state, by directly applying the ERA to the system's impulse response. The resulting models were shown to be capable of predicting the transient and long-term flow dynamics in open loop and in closed loop. They were also found to be sufficiently accurate to design robust stabilising controllers.
The method was applied to the unstable flow over a D-shaped body, with bodymounted actuators, and two sensor configurations. We first showed that, if projectionfree BPOD is used to obtain the model as in , then the analysis of actuator-sensor placement can be facilitated by extracting relevant information from the balanced modes. We then confirmed that the ERA can also be used at a much reduced cost and without the need for an adjoint solver, in order to obtain theoretically identical models, directly from the unstable system's impulse response. In experiments or if a linearised solver is not available, it was shown that accurate models can also be obtained using the ERA, from the early, approximately linear, part of the impulse response of the full nonlinear flow. Such ERA ROMs, based on the nonlinear flow dynamics, were then used to design stabilising feedback controllers for two sensor configurations: one located in the wake and one mounted on the body surface.
Both proportional controllers and robust controllers based on the H ∞ loop-shaping procedure of McFarlane & Glover (1992) were designed. In both cases, it was found that robust control tools such as the generalised stability margin and the ν-gap metric are particularly well suited to the design of robust stabilising controllers for unstable flows. These enable the analysis of the sensitivity of the closed-loop behaviour to the chosen modelling method and to model parameters such as the model order. They are central to the H ∞ loop-shaping framework and were also found to provide an effective means to optimise the robustness of proportional controllers using the ERA models.
Comparing the closed-loop response of the full-order linearised flow to the predictions made by the models showed that the ERA ROMs were sufficiently accurate to design stabilising controllers for the actual flow and also capable of predicting its closed-loop behaviour accurately. Furthermore, all controllers led to the full stabilisation of the nonlinear flow from the (unmodelled) vortex shedding limit cycle. Finally, it was found that the wake sensor was more robust to an increase in the Reynolds number than the body-mounted sensor.
All the key steps followed throughout this study can be performed with a standard nonlinear flow solver. This is a crucial point, as linearised and adjoint solvers are not available in a large number of relevant flow control scenarios. It also suggests that some experimental set-ups may benefit from this approach. In order to obtain models for unstable systems, previous studies have relied either on expensive computational techniques or on the difficult task of finding a stabilising controller for the flow, without a model for the dynamics. In contrast, the ERA is straightforward to implement and has a low computational cost. Here, we have shown that it can be used directly to construct models of high enough quality to design and implement robust stabilising controllers for unstable flows.
the forward system into a matrix X and state snapshots from the impulse response of the adjoint system into another matrix Z . Assuming that a snapshot is saved every p time steps and that (N + 1) snapshots are saved in each case, the matrices are structured as follows: Here, for k 1 is an (n x × n u ) matrix and b i is the ith column of B, so that A (k−1) b i is the forward impulse response from the ith input. Similarly, z(k) = A †(k−1) C † is an (n x × n y ) matrix and A †(k−1) c † i is the adjoint impulse response from the ith input. The forward (respectively adjoint) impulse response snapshots can be obtained either by letting the state evolve from the initial condition x(1) = b i (z(1) = c † i ) for each column of B (C † ) or by starting the simulation at x(0) = 0 (z(0) = 0) and forcing the system with an impulsive input at k = 0.
Once the snapshots have been collected, the following singular value decomposition (SVD) is performed to obtain the transformation matrices: HereT is an (n x × r) matrix andŜ is an (r × n x ) matrix. These matrices are used to simultaneously balance and truncate the system and thus obtain an rth-order ROM. The elements of the diagonal matrix Σ are the HSVs σ i mentioned above. The singular values in Σ 2 identified in (A 5a) are zero or negligible since more snapshots than there are dynamically significant states are usually stored. The number of truncated states can be chosen based on the HSV distribution. The transformed and truncated system is then given bŷ x 1 (k + 1) =Âx 1 (k) +Bu(k), y(k) =Ĉx 1 (k) ≈ y(k), (A 6) where (Â,B,Ĉ) = (ŜAT ,ŜB, CT ),x 1 ∈ C r is the ROM state andŷ ∈ C n y is the ROM output, which approximates the full system output y. This method is particularly appealing because it is applicable to even very large systems since only matrix multiplications and an SVD of an (N + 1) × (N + 1) matrix are required regardless of the actual dimension of the full system. However, it was originally designed for stable systems only: clearly, as t → ∞ the state vector x → 0 if the system is stable. As a result, if snapshots are stored over a long enough period of time, the matricesT andŜ converge and so does the identified balanced ROM.
If the system is unstable, on the other hand, the state vector becomes unbounded for large t, so it is not immediately obvious that the transformation matrices converge, that they lead to a converged ROM or that this ROM is balanced. For this reason, a number of methods have been developed to circumvent this issue (e.g. Barbagallo et al. 2009;Ahuja & Rowley 2010;Dergham et al. 2011). They require additional computations that render the algorithm significantly more costly and often completely intractable if the system dimension is very large (as for instance in three-dimensional flows).
In a recent study,  showed that the unmodified algorithm does in fact theoretically lead to a fully converged balanced ROM. In practice, finite-precision arithmetic cannot handle both the small but important transient information and the unbounded long-term growth of the response if snapshots are stored until a large final simulation time t ∞ → ∞. As a result, if the final simulation time is too large, the information related to the stable modes is lost. It was therefore shown that a compromise must be found: on the one hand, t ∞ must be large enough to obtain a converged model; and on the other, it must be small enough to retain the information about stable modes. In , accurate models were obtained both for the unstable Ginzburg-Landau equation and for the unstable flow around a cylinder. This approach therefore provides a way of obtaining balanced ROMs even for potentially very large unstable systems, and it is used in § 3.2 of this article.

A.2. The Eigensystem Realisation Algorithm
The ERA was developed by Juang & Pappa (1985). In this method, the impulse response output of the system given in (A 1) is recorded and stored in the following two Hankel matrices:  where y k = CA (k−1) B is an (n y × n u ) matrix. Once these matrices have been formed, we take the SVD of H 1 . Noting that Z † X = H 1 , this results in H 1 = UΣ V † as in (A 5a). The reduced system matrices are then given by A r = Σ −1/2 1 U 1 H 2 V 1 Σ −1/2 1 , B r = first n u columns of Σ 1/2 1 V † 1 , C r = first n y rows of U 1 Σ 1/2 1 .
   (A 9) y u FIGURE 21. Block diagram considered for H ∞ loop-shaping.
In Ma et al. (2011), it was shown that the ROM obtained with both BPOD and the ERA are in fact identical. Indeed Z † AX = H 2 , and hence (A 10) Similarly, using (A 5a) we have Σ 1/2 1 V † 1 = Σ −1/2 1 U † 1 H 1 and since the first n u columns of H 1 are equal to Z † B, we have B r =B. Finally U 1 Σ 1/2 1 = H 1 V 1 Σ −1/2 1 and again since the first n y rows of H 1 are equal to CX , we have C r =Ĉ.
This approach therefore also provides a way of obtaining balanced ROMs, but at an even lower computational cost, and without the need for an adjoint solver as it is only based on sensor measurements from as little as one impulse response simulation. This article focuses mainly on the ERA and models are obtained using this method in § § 3.3 and 3.4.
In order to use the ν-gap within the H ∞ loop-shaping procedure outlined above, one typically first obtains a nominal model G 0 and a family of 'perturbed' models that collectively represent an uncertain system, where the ith model is G i . Then, after designing the weights W as discussed above, the ν-gap is calculated for each shaped perturbed model: δ ν (WG 0 , WG i ).
The shaped nominal model is then used to obtain K = K ∞ W. At this stage, the νgap can be used to provide guarantees regarding the stability of the model family. As described by Vinnicombe (1993), given weights W and a controller K ∞ that result in a stability margin of b(WG 0 , K ∞ ) = β, two scenarios can arise for each G i : (i) δ ν (WG 0 , WG i ) < β, in which case WG i is guaranteed to be stabilised by K ∞ ; (ii) δ ν (WG 0 , WG i ) > β, in which case there exists a controller that stabilises WG 0 with a stability margin of a least β but destabilises WG i .
Note that the shaped models are again used in all cases. The second scenario is problematic since it is possible that our particular choice of K ∞ does not stabilise WG i . In this case, one option is to alter the system/models in order to reduce the associated uncertainty. Alternatively, the weights W can be changed in order to ensure b(WG 0 , K ∞ ) > δ ν (WG 0 , WG i ) for all G i . An important note therefore is that the ν-gap only provides information regarding the similarity of two models in feedback given the weights W and is not a measure of the similarity between the two open-loop models in general. In other words, modifying W is likely to have a significant impact on the calculated ν-gap.