Closed-Loop Control of a Piezo-Fluidic Amplifier
Nicholls, Chris;Bacic, Marko
2020-05-06 00:00:00
Closed-loop control of a piezo-fluidic amplifier ∗ † Chris Nicholls and Marko Bacic Oxford Thermo-fluids Institute, University of Oxford, Oxford, OX2 0ES, United Kingdom Fluidic valves based on the Coandă effect are increasingly being considered for use in aero- dynamic flow control applications. A limiting factor is their variation in switching time, which often precludes their use. The purpose of this paper is to demonstrate the closed-loop control of a recently developed, novel piezo-fluidic valve that reduces response time uncertainty at the expense of operating bandwidth. Use is made of the fact that a fluidic jet responds to a piezo tone by deflecting away from its steady state position. A control signal used to vary this deflec- tion is amplitude modulated onto the piezo tone. Using only a pressure measurement from one of the device output channels, an output-based LQG regulator was designed to follow a desired −1 reference deflection, achieving control of a 90 ms jet. Finally, the controller’s performance in terms of disturbance rejection and response time predictability is demonstrated. Nomenclature Miscellaneous α Feed-forward gain (-) F Augmented disturbance input vector (-) aug δx(t) State deviation from equilibrium (-) f Dynamic response initial frequency (Hz) δx(t) State deviation estimate (-) f (t) Dynamic response step signal (Hz) δy(t) Output deviations from reference (-) f (t) Modulating tone frequency (Hz) δy (t) Noisy measurement of δy (-) f (t) Carrier tone frequency (Hz) f Sampling frequency (Hz) δy(t) Output deviation estimate (-) G(ω) Jet dynamic system transfer function (-) δu(t) Control input deviation from equilibrium (-) Φ (ω) Output-input cross-spectral density (dB/Hz) G(ω) Empirical transfer function estimate (-) yu Φ (ω) Input power spectral density (dB/Hz) g (t) Audio amplifier signal (V) uu amp ψ( f ) Quasi-steady jet response (-) g (t) Carrier tone signal (-) ψ ( f ) ψ( f ) at design flow rate (-) g (t) Modulation signal / feed-back term (-) D m γ Chirp sweep rate (Hz/s) g (t) Control input (-) ξ Initial chirp frequency (Hz) H(x) Limiting function (saturation limits) (-) ρ Density of air (kg/m ) H (z) Transfer function fitted to plant (-) plant 2 ′ K Anti wind-up gain (-) σ Sensor noise variance for δy (t) (-) AW K System DC gain (-) A State matrix (-) DC K LQR control gain vector (-) A Augmented state matrix (-) aug LQR mÛ Inlet mass flow rate (slpm) b Inlet nozzle width (mm) n Number of time series (-) d Inlet nozzle height (mm) N Number of segments (-) B Input vector (-) ΔP Pressure difference across jet (Pa) B Augmented input vector (-) aug r(t) Controller reference (-) C Output vector (-) R Jet radius of attachment (mm) C Augmented output matrix (-) aug T Sampling time (s) C(z) Controller transfer function (-) s S(z) Sensitivity transfer function (-) e(t) State estimation error (-) e (t) Augmented state estimation error (-) u(t) Control input (-) aug F(x) System nonlinearity (-) u Control input equilibrium (-) F Disturbance input vector (-) u (t) Feed-forward term in control input (-) FF PhD Student, Oxford Thermo-fluids Institute, University of Oxford, christopher.nicholls@eng.ox.ac.uk. University Research Lecturer, Oxford Thermo-fluids Institute, University of Oxford. Member AIAA. On secondment from Rolls-Royce Plc. arXiv:2005.02968v1 [eess.SY] 6 May 2020 v(t) Sensor noise (-) ETFE Empirical transfer function estimate V Sensor noise covariance matrix (-) FFT Fast Fourier Transform W Process noise variance (-) FPGA Field programmable gate array w(t) Process noise signal (-) GM Gain margin x(t) State vector (-) IMC Internal model control x (t) Augmented state vector (-) LQR Linear quadratic regulator aug x State vector equilibrium (-) LQG Linear quadratic gaussian X Unforced DC output (-) LTR Loop transfer recovery unforced th y(t) i system output (-) LUT Look-up table y¯(t) Ensemble averaged output time series (-) MC 1 Measurement connection 1 MC 2 Measurement connection 2 z(t) Integrator state (-) NMP Non-minimum phase z (t) Noisy measurement of z(t) (-) zˆ Integrator state estimate (-) PM Phase margin PT A Pressure transducer A Abbreviations PT B Pressure transducer B AM Amplitude modulation RMS Root mean square BM Burst modulation ZNMF Zero-net mass flux I. Introduction This paper demonstrates closed-loop control of a novel piezo-fluidic amplifier that may be used in high-speed flow control applications [1–3]. Although we present here the example problem of a fluidic diverter (see Fig. 1) being used as an amplifier for a piezo actuator, the technique described in the paper can be used for control of any device based on fluidic jets. Our work is motivated by the challenge of controlling fluid flows in aerospace applications [4, 5] that offers potential for fuel burn reduction through airframe [6] as well as propulsion [2, 7] performance improvements. A major challenge to the practical use of active flow control concepts [8] is the reliability, the authority, and relatively high bandwidth (> 100 Hz) required of the actuators for control of high-Reynolds number flows often encountered in aerospace applications. While actuation on O(100Hz) is achievable with mechanical components like solenoid valves in benign low temperature environments, the same is not true for high-temperature applications such as jet engines [7]. While zero-net mass flux (ZNMF) actuators like piezoelectric buzzers may be practicable from a reliability perspective, they have limited authority in comparison with mechanical devices [8]. One way of addressing the bandwidth and authority requirements is through the use of passive fluidic oscillators [9]. These have been utilised for a range of flow control applications from improving the effectiveness of a vertical tail [10, 11] to noise control [12]. Another means to achieve both high bandwidth and high authority is by amplifying the effect of plasma or piezo actuators by switching fluidic diverters [13–16]. However, unlike conventional flow control valves, neither of these approaches allows for controlled modulation of the output flow according to a desired reference trajectory. Piezo-fluidic amplifiers [14–16] rely on the fundamental physics of jet dynamics when subjected to sound injection in a transverse direction. The purpose of this paper is to demonstrate the effective closed-loop control of such dynamics thereby enabling continuous output flow modulation according to a desired output trajectory. The physical understanding, modelling, and control of the nonlinear dynamics of fluids subject to external excita- tion is an important research topic. Much of the published work in the field consists of numerical studies [17–19], and more practical work typically comprises either open-loop control demonstrations [20–22] or extremum-seeking schemes [23–26]. There are comparatively few examples of feedback controllers based on dynamic models [27–30]. In [22], a square jet was excited with piezoelectric actuators, which were driven with an amplitude modulated signal about their resonant carrier frequency. This study demonstrated that the jet demodulated the excitation signal and responded at the modulation frequency. This result has been used to control jets in several studies since [24, 29] and is also used in the research presented here. For example, in [24], an adaptive closed-loop control scheme is used to explore the optimal AM or BM (burst modulated) forcing frequency of a synthetic jet actuator, which injected flow into the boundary layer of a NACA airfoil to promote reattachment of the separated flow. A strain gauge was used to determine the lift and drag forces on the airfoil, which were used to assess the degree of separation in the cost function for the adaptive control algorithm. This approach yielded a doubling of the lift-to-drag ratio. In [27], the drag over a bluff body was reduced by excitation with synthetic jet actuators on the trailing edge to suppress the natural vortex shedding frequency. A first order model with a static nonlinear map and delay was used for synthesizing a robust controller. The controller was able to track a reference pressure coefficient whilst the flow Reynolds number varied significantly, highlighting the benefits of using feedback. A second control scheme based on an extremum-seeking method made use of an extended Kalman filter to estimate amplitude, frequency, and phase of the fluctuations in the flow. This controller was found to achieve the same recovery of pressure and reduction of drag for half the actuation energy as the first scheme. In [29], a jet issuing from a nozzle terminated with a wide-angle diffuser was encouraged to attach to the diffuser by thrust vectoring using a synthetic jet actuator. The main jet responded to the modulation signal of an AM-driven synthetic jet. System identification was used to fit a second order dynamic model and a controller was designed using an internal model control (IMC) scheme. The controlled jet responded up to 30 - 50 Hz, achieving an order of magnitude higher bandwidth than conventional thrust vectoring mechanisms. In this paper, we employ the same piezo actuation method as in [14] where a highly reliable piezo buzzer is used to deflect the jet inside a Coandă diverter causing the device to switch its state. However, the response of the device is limited in Mair et al. [14–16] by the application of open-loop control only. The use of open-loop control has two main disadvantages: (i) no disturbance rejection, and (ii) significant variation of the switching time. The former makes the system susceptible to upstream and downstream pressure changes whereas the latter precludes the use of the device Outlet A Inlet Piezo transducer A Pressure transducer A Pressure transducer B Outlet B Piezo transducer B Fig. 1 Fluidic amplifier used in this paper. 1.6 2.1 1.6 40.0 3.4 22.0 4.6 Fig. 2 Device dimensions. Unless otherwise indicated, the units are mm. The depth of the fluid path is d = 4.8 mm. in applications that require strong synchronisation. In this paper we propose the use of closed-loop control to tackle both of these deficiencies by utilising pressure measurements from the total pressure tappings in the output channels as feedback signals (see Fig. 1) and by relying on the demodulating properties of jet dynamics as first demonstrated by Wiltse et al. [22]. Section II reviews the piezo-fluidic concept. Section III describes the experimental set-up and carries out system identification of the jet dynamics inside the device. Section IV describes the design of the output LQG controller. Finally, Section V demonstrates experimentally the effectiveness of the closed-loop scheme in controlling the jet. II. Piezo-Fluidic Concept The device used is shown in Fig. 1, with detailed dimensions given in Fig. 2. The inlet channel has a rectangular cross-section with width b = 1.6 mm by height d = 4.8 mm. Fluidic amplifiers operate using the Coandă effect, which is the tendency of a jet to attach itself to a nearby surface [31]. When supplied with a pressurised fluid, a jet issues into the interaction region of the device from the nozzle orifice and entrains the surrounding stationary fluid, which lowers the surrounding pressure. Any asymmetry causes the jet to bend to the side with slightly lower pressure, which acts to increase the pressure difference across the jet by confining the side towards which the jet bends and relieving the opposite side. The jet eventually strikes the wall to enclose a low-pressure separation bubble, which sucks in fluid from the jet that has insufficient total pressure to continue downstream [32], counteracting the pressure-reducing effect of the entrainment and stabilising the bubble pressure. A steady state is reached once there is equilibrium between the entrainment and recirculation flows - this is known as the Coandă effect [33], which gives all wall-attachment devices their stability and supports a pressure difference across the jet. This results in a bistable device, with flow attaching to one of the two walls and exiting through the corresponding outlet channel. Classically, such a diverter can be switched by injection or extraction of the transverse mass flow, causing the main jet to move past the splitter (which is set back and blunt in this device) and attach to the opposite wall. Recently, however, Mair et al [14, 15] demonstrated that zero-net-mass-flux piezoelectric transducers can be used to switch a similar device using open-loop control at a characteristic frequency. The mechanism causing the jet deflection and subsequent switching depends on which side of the device the acoustic signal from the piezoelectric transducer (the excitation) is applied from relative to which channel the jet is attached to. It is applied either from the same side or the opposite side to the channel to which the jet is attached. In [14] it is shown that, when exciting from the unattached side, the deflection is caused by exciting the shear layer roll-up mode of the jet, or one of its subharmonics. This is the most unstable natural mode, and has been shown in several studies to be at a nondimensional frequency of St = 0.012, e.g. [34, 35]. However, Mair et al. [14] note that switching is possible across a broad range of frequencies when the flow rate is sufficiently low. Exciting the jet shear layer results in a steady deflection of the jet so that a greater portion of it exits via the unattached side outlet. To understand this, consider the steady jet deflection equation, which is easily derived by considering the radial acceleration of a curved jet [36] mÛ 1 = ΔP, (1) ρbd R where mÛ is the jet mass flow, b and d are the width and height of the nozzle cross-section from which the jet emerges respectively, R is the radius of attachment of the jet, and ΔP is the pressure difference across the jet. The radius of attachment, R, results from assuming that the jet centreline follows the arc of a circle between the nozzle orifice and the attachment point (where it strikes the wall) [37]. Therefore, the smaller the value of R, the tighter the attachment and the greater pressure difference across the jet, ΔP, is required to supply the centripetal acceleration to maintain the arc shape. This equation describes the strength of the Coanda effect - smaller values of ΔP and larger values of R for a given mass flow rate, mÛ , indicate that the Coanda effect is weaker and the jet less firmly attached to the wall. Exciting the shear layer promotes vortex production, which increases the entrainment flow on both sides of the jet, but more so on the unattached side [14]. This results in a biased reduction in pressure on either side of the jet, such that ΔP decreases in magnitude. Therefore the radius of attachment, R, increases, the Coandă effect is weakened, and the jet is deflected away from the wall. In the limiting case that R → ∞, i.e. a straight jet, the flow would be divided equally between the two outlet channels by the splitter. Hence, increasing R by acoustically ex- citing the unattached side shear layer results in a portion of the jet exiting the fluidic device via the unattached side outlet. In this paper, we propose the idea of using a carrier tone to which amplitude modulation is applied as a control signal. Therefore, our control signal u(t) will amplitude modulate the carrier tone, u (t) = u(t)sin(ω t). The flow piezo c −1 rate used in this paper was 40 lpm, corresponding to a mean inlet channel velocity of 90 ms , an inlet-to-outlet pressure ratio of 1.1, and a Reynolds number based on the hydraulic diameter of 2.2 × 10 . At this flow rate, it is not possible to make the jet switch, although large deflections of the jet are possible, and a large portion of the flow can be directed out of the unattached side channel. Hence, we aim to provide a second mode of operation for the device. This mode operates at higher jet speeds relative to conventional operation so that we obtain a faster response, but due to limited piezo amplitudes does not lead to a full switch. This results in lower effective gain (from piezo amplitude to total pressure in the unattached-side outlet channel) as only part of the jet is directed out of the unattached side outlet. However, since the jet is never fully detached, the slow dynamics of detachment and reattachment are avoided, details of which can be found in Epstein [38]. This, in combination with higher jet speeds ensures, a faster device response, leading to higher effective bandwidth. The output used is the total pressure in the unattached side channel, as measured by a total pressure tapping in the centre of the channel. While it would be more accurate to use, for example, a hot-wire anemometer at the unattached side outlet orifice, this measurement strategy would not be sufficiently robust to be used in a real application. III. System Identification III.A. Experimental Set-up The experimental set-up is shown in Fig. 3. The FPGA (field programmable gate array) used is the National In- struments (NI) cRIO-9035, with the NI 9205 Analogue Input (AI) card and NI 9263 Analogue Output (AO) card. The piezo used is the Kingstate 108 dB Panel Mount Continuous External Piezo Buzzer, the audio amplifier is the Kemo 40 W M034N, and the pressure transducer is the First Sensor 10 mbar HCXPM005D6V (henceforth referred to as PT A), which has a response time of 100 μs, and the Kulite XCQ-series pressure transducer (PT B). PT A was used in all cases except when stated otherwise. A measurement connection (MC) was required to connect PT A to the total pressure tapping on the device. This consisted of a short length of 1.65 mm Scanivalve tubing to connect the tapping to an adapter, which comprised a short length of 1.6 mm Scanivalve soldered to a thicker brass cylinder with the same internal diameter. A second piece of Scanivalve tubing of appropriate (larger) diameter connected the adapter to the pressure 100 PSIG FMA-2612A air supply flow controller Amplifier (Kemo M034N) DAC (NI 9263) Piezo buzzer NI cRIO 9035 Pressure ADC transducer (NI 9205) (HCXPM010D6) Fig. 3 Experimental set-up transducer. A diagram of the measurement connection is shown in Fig 4. There were two versions of this measurement connection: MC 1 and MC 2. The latter version simply had the shortest possible lengths of each section of the connec- tion, thus reducing its filling time and increasing its bandwidth. The First Sensor pressure transducer has a response time of 100 μs. However, the limiting factor in the measurement frequency response is the filling associated with the measurement connection. Another pressure transducer, a Honeywell SDX series device, was also used but showed no improvement in the frequency response for the same reason. It is the measurement connection that causes the roll-off. The dynamics that we have studied in the present work are those of the bulk jet rather than those associated with shear layer instabilities. This would usually inform the choice of sampling rate, however in our case it was necessary to Scani-tubing Pressure 1 Scani- t v Pressure transducer membrane Brass with hole machine Fig. 4 Measurement connection diagram F(u(t)) y(t) u(t) G(s) ( !"#$%&’)*+ ,-./023456789 ) Fig. 5 Approximation of plant: Hammerstein model generate input waveforms, and it was convenient to match the controller loop rate to the sampling rate. The frequency of the input signals required was O(1kHz), and to synthesize these signals accurately we required a loop rate at least one order of magnitude greater. The loop and sampling rates were therefore set to f = 50 kHz in all cases. An analogue, first order, 25 kHz (the Nyquist frequency) anti-aliasing filter was applied to the pressure transducer measurements before sampling by the FPGA, and the flow controller used was the Omega FMA-2612A. The frequency resolution of signals obtained varied from 30 Hz for the ETFEs in Section III.D to 0.02 Hz for the quasi-steady jet response obtained from simulations, also in Section III.D. For stationary signals, a standard frequency-domain averaging method was used, where sampled time series were split into segments, their FFTs were calculated then the collection of FFTs were averaged [39]. The number of segments (N) used determined the frequency resolution ( f /N), which varied depending on the signal processing task and was chosen to keep the signal-to-noise ratio above 10 dB. III.B. Linearity There are several sources of nonlinearity. The pressure tapping measures the pressure across a central section of the unattached side channel. This is in the shear layer of the deflected jet, which does not have a linear velocity profile, so that as the deflection varies the pressure measurement varies nonlinearly. The piezo-amplifier system is another source of nonlinearity. Finally, the jet response to excitation is faster than the speed of the jet’s movement back to its natural, unexcited steady state. The system model can be approximated with a Hammerstein model, shown in Fig. 5, where the system nonlinearities have been incorporated into the static nonlinearity, F(x). The function F(x) was determined at several flow rates by driving the audio amplifier at 2.75 kHz with ampli- tude linearly increasing from 0 to 10 V . The resulting characterisations are shown in Fig. 6 with corresponding fitted pp rational functions. Input amplitudes greater than ∼ 0.8 V are omitted because the resulting deflection was not strictly monotonically increasing. The curves in Fig. 6f are the scaled, inverted, fitted functions that were implemented in look-up tables for the dynamic system identification experiments in Section III.D in order to preserve linearity. These curves demonstrate that the function F(x) is relatively insensitive to flow rate, suggesting that any feedback controller will not be limited in tracking reference jet positions away from the design flow rate because of a variation in the system nonlinearity. 80 100 90 90 80 80 70 70 60 60 50 50 40 40 30 30 20 20 10 10 0 -10 0 -10 0 0.2 0.4 0.6 0.8 0 0.2 0.4 0.6 0.8 0 0.2 0.4 0.6 0.8 Piezo-amplifier voltage (V) Piezo-amplifier voltage (V) Piezo-amplifier voltage (V) (a) 30 lpm (b) 35 lpm (c) 40 lpm - design case 80 70 0.8 30 lpm 70 35 lpm 0.7 40 lpm 45 lpm 0.6 50 lpm 0.5 0.4 0.3 0.2 0.1 -10 0 0 0.2 0.4 0.6 0.8 0 0.2 0.4 0.6 0.8 0 0.2 0.4 0.6 0.8 Piezo-amplifier voltage (V) Piezo-amplifier voltage (V) Pre-scaled control input (V) (f) Scaled, inverted functions (d) 45 lpm (e) 50 lpm Fig. 6 System nonlinearities, F(x), at several flow rates: ensemble averaged data (N = 2, blue) and fitted rational functions (red), f = 2.75 kHz; Fig. 6f shows the scaled, inverted functions for each flow rate for look-up table implementation. An experiment was conducted in order to determine the significance of the contributions to the nonlinearity curves in Fig. 6, which represent the overall input-output nonlinearity at each flow rate. A 2.75 kHz tone was used to drive the piezo-amplifier system, linearly increasing in amplitude from 0 to 0.8 V over 50 seconds, and the resulting sound pp pressure level was measured by PT B (the Kulite) in one of the outlet channels. The time series was split into 400 sets of 6250 samples, and the RMS was taken of each set, resulting in a 400-sample curve. These data were scaled up to the amplitude of the 40 lpm input-output nonlinearity rational function so as to draw a comparison, and both of these curves are shown in Fig. 7. The RMS full scale error (i.e. error relative to the maximum value) of input nonlinearity (red) relative to the overall nonlinearity (blue) is 4.2%. This justifies our use of a Hammerstein model structure, which assumes that the system nonlinearity is entirely at the input to the system (Fig. 5). III.C. Frequency Sensitivity of the Jet to Perturbation When the piezo is driven by a single frequency tone in the range 1.3 to 4.8 kHz, the jet is deflected, resulting in a steady state increase in the unattached side total pressure. To identify frequencies at which the piezo-jet system is most responsive, the piezo on the unattached side was driven through the audio amplifier with a chirp input signal with an Total pressure in unattached side Total pressure in unattached side channel due to jet deflection (Pa) channel due to jet deflection (Pa) Total pressure in unattached side Total pressure in unattached side channel due to jet deflection (Pa) channel due to jet deflection (Pa) Total pressure in unattached side Compensated piezo-amplifier voltage (V) channel due to jet deflection (Pa) 90 90 80 80 70 70 60 60 50 50 40 40 30 30 20 20 10 10 0 0 -10 -10 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 Piezo-amplifier voltage (V) Fig. 7 Input-output nonlinearity rational function at 40 lpm (solid blue) and scaled input nonlinearity as measured by PT B (Kulite) RMS (dash-dot red) amplitude of 70 mV from 1.3 to 4.8 kHz over 100 s, i.e. pp h i u(t) = Asin 2π ξt + t , (2) −1 where γ = 35 Hzs , and ξ = 1.3 kHz. The signal mean was sampled to determine the degree of deflection. Fig. 8 shows the deflection curves. It should be noted that this curve shows deflections away from the natural bias, which is itself a function of flow rate. III.D. Jet Dynamics To identify the jet’s dynamic behaviour, we carried out system identifications where the deflection was varied dynamically. The carrier signal at 2.75 kHz, g (t) = sin(2π f t) (where f = 2.75 kHz), is amplitude modulated by c c c another signal, g (t). The resulting jet is deflected dynamically and the deflection follows the shape of g (t). As an m m example, if g (t) = Asin(2π f t), where f is low enough to avoid exciting the jet dynamics, the signal driving the m m m audio amplifier is given by −1 g (t) = g (t)g (t) = sin(2π f t)F {(Asin(2π f t) + B)}, (3) amp c m c m Total pressure in unattached side channel due to jet deflection (Pa) Scaled PT B signal RMS (-) 90 30 lpm 35 lpm 40 lpm 45 lpm 50 lpm -10 -20 1000 1500 2000 2500 3000 3500 4000 4500 5000 Piezo tone frequency (Hz) Fig. 8 Static jet deflection vs perturbation tone frequency at several flow rates: 30 lpm (blue, solid), 35 lpm (red, dash-dot), 40 lpm (yellow, dash), 45 lpm (purple, dot), 50 lpm (green, bold solid) −1 where F (x) is the inverted system nonlinearity, i.e. the relevant curve in Fig. 6f, and the deflection varies −1 between that which would result from applying the constant excitations g (t) = F {(B − A)} sin(2π f t) and amp c −1 g (t) = F {(B + A)} sin(2π f t) at a frequency of f . The system considered is between the modulating signal amp c m g (t) and the total pressure in the unattached side channel. The amplitude of the offset of the carrier signal, B in the example above, was chosen to be 0.3 V because it is pp in the middle of the range of possible input amplitudes, and the modulation signal amplitude, A above, was also 0.3 −1 V . The FPGA produces the signal g (t) = g (t)F {g (t) + B}. The signal g (t) is a sinusoid which changes in pp amp c m m frequency from 10 to 1960 Hz in 30 Hz steps over 200 s, so that g (t) = Asin (2πt ( f + f (t))), (4) m 0 1 where f (t) increments by 30 Hz periodically. We chose a stepped sinusoid input rather than a chirp signal because of the high noise levels in the device. A chirp signal spreads the input energy over a broad, continuous range of frequencies, and the jet response at each frequency was found to be dominated by its random fluctuations rather than the response to the input excitation. While the stepped sinusoid does not result in a continuous range of frequencies, the response it produces dominates at a discrete set of frequencies, thus allowing accurate magnitude and phase calculations. This was Total pressure measured in unattached side channel due to deflection (Pa) -5 -10 -15 -20 -25 -30 -35 -40 -45 -50 1 2 3 10 10 10 Frequency (Hz) Fig. 9 Open-loop Bode magnitude plot from ETFE at 40 lpm and f = 2.75 kHz. ETFE from MC 1 (blue, solid) and from MC 2 (red, dash-dot). done at several flow rates around the design case (40 lpm) in order to indicate the sensitivity of the plant to variations −1 in the inlet flow rate, with the relevant nonlinearity compensator, F (x), from Fig. 6f, implemented in a LUT. The empirical transfer function estimate (ETFE) is defined as Φ (ω) yu G(ω) = , (5) Φ (ω) uu where Φ (ω) and Φ (ω) are the cross-spectral density of the output and the input and the power spectral density of yu uu the input respectively. Note that the input is taken as u(t) = g (t). Initially, the ETFE captured at 40 lpm using PT A and measurement connection 1 (MC 1) gave the blue response in Fig 9. It was found that the roll-off in the response, at ∼ 150 Hz, was a result of the dynamics of the tube connecting the pressure tapping to the pressure transducer (MC 1), despite the faster response time of the transducer itself (100 μs). As such, the connection between the pressure tapping and the transducer was redesigned (MC 2), and the experiment was repeated. The corresponding ETFE is shown in red in Fig. 9. Additionally, the same system identification was carried out at several flow rates in order to indicate the sensitivity of the plant to this parameter. The ETFEs of these data at each flow rate are shown in Fig. 10. The true jet deflection system roll-off was measured to be slightly higher by repeating the system identification experiments with PT B, which requires no measurement connection and has a higher bandwidth limitation. However, PT B suffers from poor temperature compensation which causes the mean signal to vary, making it unsuitable for use in a tracking Magnitude (dB) 0 0 0 -20 -20 -20 -40 -40 -40 1 2 3 1 2 3 1 2 3 10 10 10 10 10 10 10 10 10 Frequency (Hz) Frequency (Hz) Frequency (Hz) 0 0 0 -500 -500 -500 -1000 -1000 -1000 1 2 3 1 2 3 1 2 3 10 10 10 10 10 10 10 10 10 Frequency (Hz) Frequency (Hz) Frequency (Hz) (a) 30 lpm (b) 35 lpm (c) 40 lpm - design case 0 0 -20 -20 -40 -40 1 2 3 1 2 3 10 10 10 10 10 10 Frequency (Hz) Frequency (Hz) 0 0 -500 -500 -1000 -1000 1 2 3 1 2 3 10 10 10 10 10 10 Frequency (Hz) Frequency (Hz) (d) 45 lpm (e) 50 lpm Fig. 10 Open-loop Bode plots from ETFE with MC 2, f = 2.75 kHz. problem. Therefore, the small reduction in bandwidth was considered acceptable and PT A was used for the purposes of control. A transfer function model was fitted to G(ω) for the design flow rate, 40 lpm, which is shown in Fig. 11. The transfer function was found using MATLAB’s System Identification Toolbox, and is given by −6 −1 4.86×10 z −36 H = z . (6) plant −1 −2 −3 −4 1 − 3.88z + 5.67z − 3.68z + 0.899z The estimated input-output delay is 36 time steps, which is 0.72 ms and corresponds to a transport delay in the device. This is the time taken for the jet to travel from the inlet orifice, where upon it is acted by the piezo, to the pitot probe in the unattached side outlet channel. This distance is approximately 50 mm, which gives an time-averaged jet speed of 69 -1 -1 ms between these points. Based on the orifice mean jet velocity of around 90 ms , this seems to be a reasonable value. The curves in Fig. 10 demonstrate a significant variation in the plant dynamics with inlet flow rate away from the 40 lpm design case. In practice, it is possible that the device could be supplied with an unsteady pressure ratio, leading to varying and potentially unknown mass flow rates. The use of closed-loop control makes the system robust to these variations to some degree, which will be tested later. Phase (deg) Magnitude (dB) Phase (deg) Magnitude (dB) Phase (deg) Magnitude (dB) Phase (deg) Magnitude (dB) Phase (deg) Magnitude (dB) 0 -10 -20 -30 -40 1 2 3 10 10 10 Frequency (Hz) -500 -1000 1 2 3 10 10 10 Frequency (Hz) Fig. 11 Open-loop Bode plot at 40 lpm: ETFE (blue, solid) and fitted transfer function (red, dash-dot). Phase (deg) Magnitude (dB) We now consider the effect of the static deflection curves on the dynamic response. It is well known that the fre- quency content of an amplitude modulated signal is given by the sum and difference frequencies of the carrier and modulation signals, sin(2π f t) (Asin(2π f t) + B) = [cos(2π ( f − f ) t) − cos(2π ( f + f ) t)] + Bsin(2π f t). (7) c m c m c m c Wiltse and Glezer [22] demonstrated that jets demodulate amplitude modulated acoustic signals and the bulk jet response is seen at the modulation frequency if the amplitude is high enough. This can be expressed mathematically by multiplying the signal by the carrier tone and then low-pass filtering. In practice, the tones on the right-hand side in (7) have different magnitudes from one another according to the shape of the relevant static deflection curve in Fig. 8. The effective magnitude response resulting from this variation alone (i.e. ignoring any jet or sensor dynamics) at a given frequency f is therefore a combination of the magnitude of the static deflection curve in Fig. 8 at f − f and f + f , m c m c m along with the DC offset in jet position in response to the term Bsin(2π f t) in (7). This effective magnitude response, which has no corresponding phase response, is referred to as the quasi-steady jet behaviour. If we fix f = 2750 Hz, this can be written as a single-input function, which we name ψ( f ). The overall system dynamic response to the input (7), given that the amplitude nonlinearity, F(x), has been compensated for, is then given by y(t) = ψ( f ) |G(2π f )| sin (2π f t + φ (G(2π f )) ), (8) m m m m where G is the true jet dynamic transfer function, |G| is its magnitude response and φ(G) is its phase response. In the frequency domain, without compensating for ψ( f ), the ETFE captured from system identification experiments is G(2π f ) = ψ( f )G(2π f ). (9) To determine the quasi-steady jet behaviour, ψ( f ), a simulation was created where two tones, initially at 2750 Hz and -1 out of phase, linearly increased and decreased in frequency respectively at a rate of 28 Hzs over 50 seconds. This -1 corresponds to two copies of the signal in (2), with ξ = 2750 Hz, γ = ±28 Hzs . The magnitude of these tones were assigned by reference to a look-up table, which interpolated the values of the relevant static deflection curve in Fig. 8. The tones were then summed and the signal demodulated by multiplying by the 2750 Hz carrier tone. A block diagram of the simulation is shown in Fig. 12. The power spectral density of the resulting signal, ψ( f ), is plotted up to 1350 Hz in Fig. 13 at all of the flow rates considered. The curve for 40 lpm is referred to as ψ ( f ). These responses in Fig. 13 justify our choice of carrier tone at D Create ramps to indicate current frequency for referencing LUTs Assign amplitudes according to static reponse LUT (static jet response) Create chirps acsending & decsending from carrier tone Demodulation Output Create carrier tone for demodulation LUT (static jet response) Fig. 12 Block diagram of simulation used to determine quasi-steady jet response, ψ( f ), from static jet deflection response (Fig. 8). 30 lpm 35 lpm 40 lpm - (f ) D m 45 lpm 50 lpm -2 -4 -6 -8 1 2 3 10 10 10 Frequency (Hz) Fig. 13 Quasi-steady state jet behaviour, ψ( f ), at several flow rates Magnitude (dB) 2750 Hz, since the resulting responses are independent of flow rate to within 1 dB, and the roll-off is higher than that of the dynamic response. To consider the effect of these quasi-steady curves on the dynamic response, we repeated the dynamic response identification at the design flow rate (40 lpm) with the amplitude of the modulating tone set to invert −1 the shape of the 40 lpm quasi-steady response curve in Fig. 13. This was implemented in a look-up table (ψ { f }), which was referenced by the frequency of the modulating tone being produced, such that the audio amplifier signal was given by −1 −1 g (t) = sin(2π f t)F ψ { f + f (t)} sin (2πt ( f + f (t))) + B . (10) amp c D 0 1 0 1 While it would be more precise to set the amplitudes of the sum and difference tones resulting from the product of the modulating and carrier tones according to the inverse of the 40 lpm static deflection curve in Fig. 8, as in the simulation −1 described above, this was not possible in practice due to the static nonlinearity compensator, F (x), in (10), which makes it impossible to evaluate the first term on the right-hand side of (10) analytically. A dynamic system identification was conducted at 40 lpm using (10), with the carrier offset set to 0.25 V (B in pp (10)). The expected result was predicted by numerically inverting ψ ( f ) and applying it to the original 40 lpm ETFE D m −1 in Fig. 11, i.e. rearranging (9) to give G = ψ G . The ETFE resulting from the experiment, the predicted result −1 given by ψ G , and the original ETFE for 40 lpm from Fig. 11 are shown in Fig 14. There is good agreement between the predicted and measured magnitude responses, and the phase response is the same as the original ETFE data, G, as expected. This demonstrates that (9) is a reasonable model for the effect of the quasi-steady jet response. −1 However, it was not possible to implement ψ { f } in practice in a real-time controller because it is defined in the frequency domain and has no corresponding phase response, such that a linear filter could not be used as a model. Additionally, we ignored the 36 T transport delay in the plant model (6) for the purposes of linear control, which is justified in Section IV. IV. Controller Design and Simulation We adopted an LQR control strategy, and therefore required an observer to estimate the system states. Converting (6) to a state-space system, x(t + T) = Ax(t) + Bu(t) (11) y(t) = Cx(t) + Du(t), 0 -10 -20 -30 -40 -50 1 2 3 10 10 10 Frequency (Hz) -500 -1000 1 2 3 10 10 10 Frequency (Hz) Fig. 14 ETFEs at 40 lpm: original curve, G (blue, cross), predicted magnitude response from numerical inversion of quasi-steady −1 jet response, ψ G (red, star) and the result of experimentally inverting ψ ( f ) (yellow, square). D m Phase (deg) Magnitude (dB) gave the matrices 3.88 -2.83 1.84 -0.899 2.00 0 0 0 A = ; D = 0; 0 1.00 0 0 0 0 0.500 0 (12) 0.00391 B = ; C = . 0.00297 0 0 0 The equilibrium of the system is shifted by subtracting the steady state values of the states, input and output once the output has been driven to a reference signal, r. These are x , u , and r respectively. New variables are introduced to 0 0 describe variation from the steady state values, namely δx(t) = x(t) − x , δu(t) = u(t) − u , and δy(t) = y(t) − r, so 0 0 that the system dynamics take the form δx(t + T) = Aδx(t) + Bδu(t) (13) δy(t) = Cδx(t). To eliminate steady state error, the plant is augmented with an integrator state defined by (14) z(t + T) = z(t) + δy(t), such that when process and sensor noise is added, the dynamics become δx(t + T) A 0 δx(t) B F = + δu(t) + w(t) z(t + T) C 1 z(t) 0 0 (15) δy (t) C 0 δx(t) = + v(t), z (t) 0 1 z(t) ′ ′ where δy (t) and z (t) are the noisy measurements of δy(t) and z(t). The regulator minimised the cost function given by i=∞ T 2 J = δx Qδx + Rδu , (16) i i i=0 where δx = δx(t + iT) and δu = δu(t + iT). The tuned LQR cost parameters are i i 1 0 0 0 0 0 0 0 0 0 Q = , R = 10 , (17) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 10 where Q is the state-error cost matrix and R is the cost of control actions. Q was chosen as per standard output-weighted LQR, i.e. C C, and the integrator cost was tuned manually. The states were not measured directly so a Kalman filter was used for state estimation. The observer dynamics take the form δxˆ (t + T) A 0 δxˆ(t) B δy (t) − δyˆ(t) = + δu(t) + K (18) zˆ(t + T) C 1 zˆ(t) 0 z (t) − zˆ(t) δxˆ(t) δyˆ(t) = , (19) C 0 zˆ(t) ˆ ˆ where δx, δy and zˆ are the estimates of the state deviations, the output error and the integrator respectively. The augmented system matrices are substituted as follows A 0 B C 0 A = , B = , C = , aug aug aug C 1 0 0 1 (20) F δx δx − δxˆ F = , x = , e = , aug aug aug 0 z z − zˆ leading to closed-loop dynamics given by x (t + T) A − B K B K x (t) F 0 w(t) aug aug aug LQR aug LQR aug aug = + , (21) e (t + T) 0 A − K C e (t) F −K v(t) aug aug f aug aug aug f where e is the state estimation error. aug There is no automatic guarantee of robustness despite the individual gain and phase margin guarantees of the Kalman ^_ Increasing confidence Chosen controller in the measurements LQR - GM = 47.2, PM = 90.4 W/V = 1e-5, GM = 61.1, PM = 89.4 W/V = 1e-3, GM = 41.8, PM = 84.9 Z[\ W/V = 1e-1, GM = 26.8, PM = 75.1 W/V = 1e1, GM = 19.1, PM = 80.4 W/V = 1e3, GM = 22.6, PM = 83.7 W/V = 1e5, GM = 25.8, PM = 86.2 WXY :;>AEH 10<=?@BCFG IJKLMNOPQRTUV -1 0 1 2 3 10 10 10 10 10 Frequency (Hz) Fig. 15 LTR procedure demonstrated with sensitivity functions for a variety of process-to-sensor noise ratios. Units for gain and phase margins in legend are dB and degrees respectively. LQR curve is dash-dot blue. filter and LQR separately in continuous-time [40]. Furthermore, in the practical application considered here, discrete- time LQR has typically inferior stability margin properties [41]. It is well known that the loop-transfer recovery (LTR) methodology presented by Doyle and Stein in [42] for continuous-time systems allows the robustness properties of the full state feedback case to be recovered. It is also well understood that there lies a trade off between recovering these robustness margins and the system noise performance in continuous-time. In discrete-time it may not be possible to recover fully the LQR loop transfer properties [43]. The value of the sensitivity transfer function across the frequency space, S(z) = , (22) 1 + C(z)H (z) plant determines a controller’s ability to reject disturbances, follow the reference signal as well its noise performance [44]. The sensitivity transfer functions of the system with several process-to-sensor noise ratios are shown in Fig. 15 2 2 to demonstrate the LTR procedure. The integrator state measurement noise variance is fixed at 100σ , where σ is ′ ′ y y the sensor noise variance for the measurement δy (t). Fig. 15 shows that the LQG curves approach the LQR curve as the process-to-sensor noise ratio is increased. However, practical controller design is not as simple as choosing Phase (deg) Magnitude (dB) -30 -40 -50 -60 -70 -80 -90 -100 -110 1 2 3 4 10 10 10 10 Frequency (Hz) Fig. 16 Jet spectrum: PSD of total pressure in unattached side channel at 40 lpm, no excitation. the system found by allowing → ∞, which would likely suffer from poor noise performance. In order to select a value of , the system noise spectrum must be considered. A time series was captured by sampling the unattached side channel pressure tapping whilst the jet was unexcited. The turbulence is technically an output disturbance but cannot be rejected by control action - the acoustic signals only control the jet position. Therefore, the turbulence can be thought of as sensor noise for control purposes. The power spectral density (PSD) of the time series is shown in Fig. 16. As can be seen in Fig. 16, the roll-off is caused in part by the measurement system, as described in section III.D, and the frequency content is in the same bandwidth as the system. This means that choosing a larger causes the controller not only to reject disturbances that can be controlled (e.g. fluctuations in the inlet mass flow) and follow the reference more effectively, but also to react to the jet turbulence. Therefore a compromise must be made to give both acceptable disturbance rejection and noise performance. The value of chosen is 1e1, with a gain margin of 19.1 dB and a phase margin of 80.4 . The tuned LQR costs and Kalman filter variances that give these results are 1e-1 0 W = 1, V = (23) 0 1e1 th where W and V are the process and sensor noise variance and covariance respectively. Note that the 5 state in the nd system is the integrator state, and that the 2 dimension in the covariance matrix V is the variance of the integrator state measurement. The LQG controller magnitude response is shown in Fig. 17. The cross-over frequency for this Magnitude (dB) 50 -50 -100 -1 0 1 2 3 10 10 10 10 10 Frequency (Hz) -200 -400 -1 0 1 2 3 10 10 10 10 10 Frequency (Hz) Fig. 17 LQG transfer function Bode plot Phase (deg) Magnitude (dB) _ Plant Demodula‘on Hammerstein model LPF Fig. 18 Closed-loop system block diagram controller is around 50 Hz, while the transport delay (0.72 ms) is a period corresponding to ∼ 1400 Hz. Therefore, since the closed-loop bandwidth is two orders of magnitude below this, our earlier decision not to include the input-output delay in the controller design process is justified. V. Implementation and Results The LQG controller simulated in Sec. IV was implemented in the FPGA and the system block diagram is shown in Fig. 18. As shown in the diagram, a feed-forward term was added to help increase the speed of the initial rise of the step response, giving the control law δx(t) u(t) = −K + u (t), (24) LQR FF zˆ(t) where u (t) is the feed-forward term. The DC level measured by the pressure transducer due to the unforced response FF is X , K is the DC gain of the model, and α is a parameter used to vary the contribution from the feed-forward unforced DC term when tuning the step response. The implementation of the control law is therefore (r − X ) unforced g (t) = g (t) + α . (25) DC (r−X ) unforced It can be seen by comparison of (24) and (25) that g (t) is the feedback term and u (t) = α . Fig. 18 also m FF DC ′ ′ shows that g (t) was limited to between 0 and 0.8 V, denoted by H g (t) . This was necessary because the linear m m −1 model fails to predict the plant behaviour outside this region. If F (g (t)) > 0.8 V, the deflection does not increase because of the lack of strict monotonicity of the system characterisation in this range (see Sec. III.B). The model −1 indicates that values of F (g (t)) < 0 V will cause the jet to deflect away from the centre of the device, which is not m |} u (t) no abcdefghijklm y(t) yW Fig. 19 Anti wind-up scheme the case; the deflection is caused, in steady state, by g (t), whose phase does not affect the direction of deflection. The −1 audio amplifier signal is therefore g (t) = F (H {g (t)}) g (t). amp m c The combination of an integrator state and input saturation limits necessitates an anti wind-up scheme to avoid an unstable integrator state. The details of this scheme are shown in Fig. 19, where the value of the anti wind-up gain, K , was set by manual tuning, and took a value of 1e-5. This was not included in the overall closed-loop block AW diagram in Fig. 18 for brevity. V.A. Step response The ensemble averaged (n = 50) response of this system to a reference step from 0 to 74 Pa deflection in the unattached side channel was recorded in closed-loop, with the inlet flow at 40 lpm. As a comparison, an ensemble averaged (n = 50) full switch was recorded at a lower flow rate, such that the pressure in the initially unattached side channel was the same as the deflection in the closed-loop case once the jet had switched, i.e. 74 Pa. To achieve this, the flow rate was set to 15 lpm, and the amplitude of the piezo tone was set to 0.8 V. The ensemble-averaged time series were then filtered with a notch filter at 2.75 kHz to reduce the direct coupling between the piezo tone and the pressure transducer. A broad notch was also applied at 650 Hz to reduce the jet background noise in order to judge when the deflection first reached its steady state value. The resulting step responses are shown in Fig. 20a, while Fig. 20b shows an ensemble averaged control signal (n = 50), g (t) (see Fig. 18). In terms of rise time, the closed-loop response at 1.1 ms is thirty times faster than the open-loop full switch at 34 ms. However, this improvement comes from the faster jet dynamics at higher velocities and from using a feed-forward term in the control law rather than the feedback term. Nevertheless, using feed-forward control alone does not compensate for model uncertainty or disturbances such as mass flow variations. V.B. Disturbance rejection - mass flow variation The controller’s ability to reject disturbances was tested in two ways. First the inlet mass flow was varied. This was done in the steady state since the response time of the FMA-2612A was insufficient to test the closed-loop bandwidth 150 0.5 0.4 0.3 0.2 0.1 -50 -0.1 -100 -0.2 0.45 0.5 0.55 0.45 0.5 0.55 Time (s) Time (s) (a) Ensemble averaged & filtered step responses (n = 50): (b) Control signal (ensemble averaged, n = 50), g (t) in Fig. closed-loop at 40 lpm (blue) and open-loop full switch at 15 18. lpm (red). Reference (dash-dot black). Fig. 20 System step response: Open- and closed-loop system responses and closed-loop control input. of the controller. The flow rate was set to values between 19 and 54 lpm, which was the range of flow rates that the controller was able to reject whilst tracking a reference of 50 Pa deflection from the natural bias state. The mean control inputs required to maintain these deflections are shown in Fig. 21. The figure shows how the DC gain of the jet deflection system varies over different flow rates, with a broad maximum DC gain at the minimum of the curve, i.e. between 28 and 35 lpm. V.C. Disturbance rejection - input disturbance The response of the controller to input disturbances was tested at the design flow rate (40 lpm) by removing the feed-forward term whilst the controller was tracking a constant reference deflection of 50 Pa. The pressure in the unattached side channel due to the deflection and the control input signal are shown in Fig. 22. The signal processing for these data were as follows: as for the step responses in Section V.A concerning the step response, a 2.75 kHz high-Q factor notch as well as a broad notch at 650 Hz were applied to ensemble averaged time series (n = 50). This was repeated 6 times and the resulting signals were also ensemble averaged. This additional ensemble averaging step was required due to the higher variance of the open-loop response times (relative to the mean). The input signal was ensemble averaged without any filtering (n = 350). Fig. 22 indicates that the closed-loop response time to disturbances is around 28 ms, whereas the case with no feedback term has a faster response time (1.4 ms) but is unable to reject the input disturbance (as expected). It is interesting to consider the response times of the open- and closed-loop cases over several flow rates in order to evaluate the robustness of the input disturbance rejection properties to varying operating Total pressure in unattached side channel due to jet deflection (Pa) Control input (V) 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 15 20 25 30 35 40 45 50 55 Flow rate (lpm) Fig. 21 Mean control input required to maintain 50 Pa jet deflection over several flow rates. 0.2 -0.2 -0.4 -0.6 -0.8 -20 -1 0.45 0.5 0.55 0.6 Time (s) Fig. 22 Controller response to input disturbances. Pressure in unattached side channel due to deflection: closed-loop response (blue) and open-loop response (yellow). Control signal (g (t)) showing step change in u (t) at t = 0.5s (red). Signal processing FF scheme described in text. Total pressure in unattached side Control input required to sustain channel due to jet deflection (Pa) 50 Pa jet deflection (V) Control input, g '(t) (V) m 1.4 1.3 Open loop response 1.2ƒ§¤'“«‹›fifl–