程序代写代做代考 go html graph Hive algorithm 2014/15

2014/15
Development of an Online Identification Algorithm to Stabilize Combustion Systems
Final Report: MP4079/B427
Submitted by: Calum Clarke – N1403292E
Supervisor: A/P Zhao Dan
School of Mechanical & Aerospace Engineering
A Final Year Project Report

Project MP4079/B427 | Calum Alexander Clarke
Table of Contents
Chapter 1 : Introduction………………………………………………………………………………………………….. 5 Chapter 2 : Literature Review ………………………………………………………………………………………….. 6 2.1 – Combustion Instability …………………………………………………………………………………………. 6 2.2 – Active Control Approach …………………………………………………………………………………….. 10 2.2.1 – Open Loop Control ……………………………………………………………………………………… 11 2.2.2 – Closed Loop Control ……………………………………………………………………………………. 11 2.3 – Passive Control Approach …………………………………………………………………………………… 13 Chapter 3 : Van Der Pol …………………………………………………………………………………………………. 14 3.1 – Non-linear Limit Cycle Oscillator………………………………………………………………………….. 14 3.2 – Lindstedt-Poincaré Method ……………………………………………………………………………….. 18 3.3 – Multiple Scale Method ………………………………………………………………………………………. 23 3.4 – Method of Averaging ………………………………………………………………………………………… 27 3.5 – Runge-Kutta Approximation ……………………………………………………………………………….. 31 Chapter 4 : Online Identification Algorithm……………………………………………………………………… 33 4.1 – Standard Least Means Square Algorithm ……………………………………………………………… 33 4.2 – X-Filtered Least Means Square Algorithm …………………………………………………………….. 36 Chapter 5 : Results & Discussion…………………………………………………………………………………….. 39 5.1 – Standard LMS ……………………………………………………………………………………………………. 39 5.2 – X-Filtered LMS …………………………………………………………………………………………………… 49 Chapter 6 : Conclusion ………………………………………………………………………………………………….. 52 6.1 – Further Work…………………………………………………………………………………………………….. 54 Chapter 7 : References………………………………………………………………………………………………….. 55 7.1 – Other useful links & papers…………………………………………………………………………………. 57 Chapter 8 : Appendix…………………………………………………………………………………………………….. 58
1|Page

Project MP4079/B427 | Calum Alexander Clarke
Table of Figures
Figure 2-1 : Thermo-Acoustic Feedback Loop………………………………………………………………………………….. 6 Figure 2-2 : Typical Combustion Instability Behavior [23] …………………………………………………………………. 8 Figure 2-3 : Open Loop Control System ………………………………………………………………………………………… 11 Figure 2-4 : Closed Loop Control System ………………………………………………………………………………………. 11 Figure 3-1 : Limit Cycle μ = 10 (Large)…………………………………………………………………………………………… 15 Figure 3-2 : Phase plot μ = 5 (Large) …………………………………………………………………………………………….. 15 Figure 3-3 : Phase Plot μ = 0.2 …………………………………………………………………………………………………….. 16 Figure 3-4 : Van der Pol (μ = 0.2) vs Real Combustion System …………………………………………………………. 16 Figure 3-5 : Lindstedt Poincare – A vs T (μ = 0.2) ……………………………………………………………………………. 20 Figure 3-6 : Lindstedt Poincare – A vs T (μ = 1) ………………………………………………………………………………. 21 Figure 3-7 : Lindstedt Poincare – A vs ω (μ = 0.2)……………………………………………………………………………. 21 Figure 3-8 : Lindstedt Poincare – A vs ω (μ = 1)………………………………………………………………………………. 22 Figure 3-9 : Multiple Scales, Initial 2.5, μ = 0.2 ………………………………………………………………………………. 25 Figure 3-10: Multiple Scales, Initial 2.5, μ = 1 ………………………………………………………………………………… 25 Figure 3-11 : Multiple Scale, Initial 0.5, μ = 1…………………………………………………………………………………. 26 Figure 3-12 : Method of averaging – a vs t (μ=0.2)…………………………………………………………………………. 29 Figure 3-13 : Method of averaging – a vs t (μ=1)……………………………………………………………………………. 30 Figure 3-14 : Sinusoidal Reconstruction ……………………………………………………………………………………….. 32 Figure 3-15 : Runge-Kutta Van der Pol Limit Cycle (μ = 0.2) …………………………………………………………….. 32 Figure 4-1 : LMS Adaptive Filter…………………………………………………………………………………………………… 33 Figure 4-2 : FIR Filter Approximation ……………………………………………………………………………………………. 34 Figure 4-3 : Black Box Controller Diagram …………………………………………………………………………………….. 35 Figure 4-4 : X-Filtered LMS Adaptive Filter ……………………………………………………………………………………. 36 Figure 4-5 : ANC Sound wave Illustration ……………………………………………………………………………………… 37 Figure 5-1 : Estimation Accuracy …………………………………………………………………………………………………. 39 Figure 5-2 : Instant Controller Activation ……………………………………………………………………………………… 40 Figure 5-3 : Controller Activation @ t = 2500 iterations, μ = 0.2 ………………………………………………………. 41 Figure 5-4 : Step size = 0.001 ………………………………………………………………………………………………………. 42 Figure 5-5 : Lower Frequency Limit Cycle (b = 62.5) ……………………………………………………………………….. 42 Figure 5-6 : Higher Frequency Limit Cycle (b=1000) ……………………………………………………………………….. 43 Figure 5-7 : Higher Frequency Limit Cycle Zoomed Graph ………………………………………………………………. 43 Figure 5-8 : Limit Cycle Control (μ = 1)………………………………………………………………………………………….. 44 Figure 5-9 : Limit Cycle Control (μ = 0.02)……………………………………………………………………………………… 45 Figure 5-10 : Limit cycle control with additional noise (μ=0.2)…………………………………………………………. 45 Figure 5-11 : Limit cycle control with additional noise (μ=0.2) Zoomed Graph…………………………………… 46 Figure 5-12 : Increased noise until loss of control ………………………………………………………………………….. 47 Figure 5-13 : 0.8*Additional Noise, 0.001 Stepsize ………………………………………………………………………… 47 Figure 5-14 : 10*Additional Noise, 0.001 Stepsize …………………………………………………………………………. 48 Figure 5-15 : Accuracy of Secondary Prop Path Estimation – X-Filtered Algorithm ……………………………… 49 Figure 5-16 : Results of Attenuation (X-Filtered) ……………………………………………………………………………. 50 Figure 5-17 : Results of Attenuation (X-Filtered) Zoomed Graph ……………………………………………………… 51
2|Page

Project MP4079/B427 | Calum Alexander Clarke
Abstract
Combustion instability is a serious problem and its control presents many challenges to the engineers who face it. It can have major negative effects on the systems it is included in, which is many, because combustion is arguably the most useful chemical reaction to humans, especially in industry. Fortunately these instabilities can be controlled. A common and practical control strategy is to introduce a secondary signal into the system which acts in a deconstructive manner to the instabilities. In this thesis an online identification algorithm is described, based on the Least Means Square (LMS) method, which is used to stabilise the unstable pressure oscillations within a simulated combustion chamber. Two separate versions of the LMS algorithm are used, the standard version and the X-filtered method. These are compared under the influence of additional noise being added to the input signal. The combustion chamber pressure oscillations are simulated by solving the van der pol equation to produce a limit cycle similar to that of an unstable combustion system. The algorithm makes an accurate estimate of the limit cycle and acts to attenuate it, performing the attenuation quickly and effectively as would be required in a real life system. The accurate estimation is required as in a real life system it is likely that the instability characteristics of the combustor would not be known previous to system activation. The procedure followed and the results obtained are all presented within this report.
3|Page

Project MP4079/B427 | Calum Alexander Clarke
Acknowledgements
The author would like to take this opportunity to express his gratitude to his home university, the University of Strathclyde, and also Nanyang Technological University for allowing this exchange to take place and with it this final year project.
The author would also like to express severe thanks to A/P Zhao Dan for his guidance and assistance with regards to the project.
Also, thanks go out to peers and colleagues who have supported the author in numerous respects throughout his time at NTU.
4|Page

Project MP4079/B427 | Calum Alexander Clarke
Chapter 1 : Introduction
This project sets the student the task of developing an online identification algorithm of which the main functionality is to stabilise an unstable combustion system. The main objectives are as follows: Simulate the inside of an unstable combustion chamber by representing the limit cycles within. Design a controller which will observe and estimate the future activity of the incoming signal from the combustion chamber with a prediction as accurate as is possible. Finally, use the controller to control an actuator to send a deconstructive signal which will counteract the result of instability to as low a value as achievable. In the case of this project the controller will also act as the actuator, sending a signal directly back into the system. The system will be tested at various frequencies and amplitudes of input, and also with the addition of external noise. The separate LMS algorithms will then be compared in this respect. If these objectives are a success the project can be deemed complete. For the coding and simulations involved MATLAB will be used. It is chosen for its compatibility with the problem and also because some previous experience of the software is possessed by the student.
After this initial introductory chapter the report will cover the literature review where background theory of relevant material to the project will be discussed. This section presents the results of the research undertaken at the beginning of and throughout the project. Chapter 3 will discuss how the combustion chamber is simulated, concentrating on the van der pol equation and various methods for accurately estimating limit cycle performance. Next, the online identification algorithm will be looked at in detail including its construction, implementation and the theory behind it. Following this in chapter 5, the results from the work undertaken will be presented and discussed. Subsequently, there will be a short conclusive chapter. An appendix is attached at the end of the report in the form of a CD containing all the matlab code used throughout the project, the figures in this report, and the report itself in pdf and word document format.
5|Page

Project MP4079/B427 | Calum Alexander Clarke
Chapter 2 : Literature Review
2.1 – Combustion Instability
Combustion systems have many uses in industry and other application areas, primarily aerospace and power generation. However many of these systems come at the cost of being vulnerable to combustion instabilities. Combustion instability within combustion systems can seriously reduce the working life of the engine involved. Besides the obvious implications of uncontrollable increases in pressure oscillations on the structure itself factors such as thermal output, thrust, emission levels and overall efficiency must also be considered. In certain cases it can also lead to extreme failure of the system. The first occurrences of combustion instability were discovered in rich burning devices such as ramjets, rockets and afterburners [1, 2] but it was not scientifically defined until Rayleigh in 1878 [3].
Combustion instability generally refers to large amplitude pressure or heat release fluctuations or oscillations within a combustor [4]. The cause is the excitement of one or more of the acoustic modes present in the system. In many systems, and specifically the one dealt with in this report, the instability can be simplified down to the positive feedback coupling of these two dominant processes; heat-release and pressure. This is sometimes known as thermo-acoustic instability. At this condition the excited acoustics effect the heat release which in turn influences the pressure causing a feedback loop between the two. This causes the oscillations which grow continuously until a saturation point is reached whereby the system continues to oscillate in the manner of a limit cycle. A simple figure is shown below to help illustrate this.
Figure 2-1 : Thermo-Acoustic Feedback Loop
6|Page

Project MP4079/B427 | Calum Alexander Clarke
As shown above in figure 2-1, the feedback causes two main branches which return to positively feedback each other. The heat release oscillations are denoted by 𝑄̇′ and the pressure oscillations by 𝑃′. Rayleigh’s index can be used to describe this situation.
𝑇1 ∮ 𝑃 ′ 𝑄 ̇ ′ 𝑑 𝑡 > 0 ( 1 . 0 )
Where the integral is taken over one time period, T. When the index is positive thermo-acoustic instability occurs by cause of acoustic modes and it can be known that the heat release and pressure forces are in phase. For negative values of the index thermo-acoustic damping is present and the pressure and heat release forces can be assumed to be out of phase. Essentially, the Rayleigh index can be used to find the total acoustic energy by integrating it over one full cycle. If the result from this integration is greater than the rate of acoustic loses then the system is unstable as there is excess acoustic power which can excite modes and cause instability.
Rayleigh’s criterion [5] is the basis for thermo-acoustic instability analysis and can be evaluated by using the above Rayleigh index over one cycle of instability. Rayleigh’s criterion states that if the heat-release rate and pressure are in-phase the system is considered unstable. However, if the heat-release rate and pressure can be made out of phase a cancelling out effect occurs reducing the oscillations and creating a stable combustion environment. To be able to accurately predict these conditions for instability the physical mechanisms which drive it must be understood, there has been considerable effort placed into this [6, 7, 8, 9, 10]. Unfortunately the information and understanding achieved thus far is not enough to design a system which is stable by itself. It is nearly impossible to predict if combustion instabilities will occur within a system or not. Hence the need for alternative control strategies.
A common example of combustion instability is shown in figure 2-2. In this case stable operation is present until a change in operation around 1.9 seconds. At this point the change in operation causes the acoustic mode to excite and hence starts the feedback loop described earlier and shown in figure 2-1. From here until around 3.8 seconds, the oscillations are growing exponentially until the saturation point is reached. The system then shows the resultant limit cycle continuing on indefinitely. These oscillations would have to be attenuated as swiftly and efficiently as possible in a real life system.
7|Page

Project MP4079/B427 | Calum Alexander Clarke
Figure 2-2 : Typical Combustion Instability Behavior [23]
When combustion instability occurs it can drastically increase the forces present inside the combustion chamber and potentially cause severe damage or breakage. This is a major concern as these systems cost a serious amount of money to produce and there is also the potential risk to life, especially when considering aero-engines, rockets and other systems which involve carrying human transports. As mentioned earlier, extreme failure is also a possibility of intense instability and this must be avoided at all costs for obvious reasons.
Although thermo-acoustic instability is specifically looked at earlier there are many factors which can cause instability, including specific operating conditions and combustor geometry. For example, with reference to a propulsion device such as a ramjet engine it is usually desirable to run this process at near stoichiometric conditions as this leads to a greater heat release and in turn, greater performance. In doing so however, the chances of instability are increased. The case is also the same for gas turbines where burning at lean ‘blow-out’ operating conditions with an equivalence ratio below 1 is preferable due to reduced exhaust emission [11], specifically NOx formation. But again this condition increases the chances of instability within the system.
8|Page

Project MP4079/B427 | Calum Alexander Clarke
To look at the cause of instability with an overall perspective, there are many different possible causes. These can include acoustics, heat-release, flame & chemical kinetics, convection, vaporisation, fluid dynamics and pressure. All of these play a part in combustion efficiency and stability, which in turn play a vital role in defining an engine’s operating capabilities and parameters including engine efficiency and emission levels. Ideally a situation where all of these separate parameters could be controlled would be implemented. Unfortunately this is not yet possible. In order to maximise combustion efficiency and system operation safety it is required that an active control strategy is adopted. Ideally a stable combustion system would be built but the technology for this has not yet been found, combustion instabilities cannot be predicted before a system is under operation. The system must be observed and when instabilities occur action must be taken to notice the oscillations, predict their future activity accurately and attenuate them effectively.
In terms of control, there are two main approaches. More recently there is active control, which is the control method looked at in this project and there is also passive control techniques which can be used. These are both discussed further in the following sections.
9|Page

Project MP4079/B427 | Calum Alexander Clarke
2.2 – Active Control Approach
Active combustion control is still a topic under active research. However, online identification algorithms are a popular solution for combustion instability due to their success in real-life systems. See [12] for a few examples of non-linear control models. They monitor the happenings of the combustor taking measurements of many of the combustion system dynamics and properties and use this information paired with a controller and an actuator to attenuate the unstable oscillations. This usually involves sending energy back in to the system in a controlled manner. The key components of an active control system are sensors, actuators and most importantly a control algorithm. Active control strategies have been implemented since as early as the 1950’s [13, 14].
Active control approaches force the combustion system to act in such a way that the system damps its own instabilities. This is usually achieved by using a phase difference technique (closed- loop) or by breaking the feedback cycle (open-loop) caused by excited acoustic modes. The phase difference technique is similar to that described earlier, Rayleigh [5], where equal amplitude waves to that of the instable oscillations are re-entered into the system 180o out of phase in order to cancel each other out. Active controllers can be broken into the two main categories of closed-loop or open-loop control, these are examined in more detail below.
An ideal control algorithm should be able to perform over a wide range of conditions without losing its level of control and while also maintaining robust stability. A good example of this is described in [15] where an observer based adaptive controller is used to stabilise multiple unstable modes without knowing the characteristics of the actuation or combustion system. However, as mentioned in [15] control & stability are not compatible with each other and a compromise must be made between achieving better control performance and minimizing the model uncertainties. As well as this there are other complications to consider. For example, the occurrence of instabilities is difficult to predict due to the unknown main causes, any change in stable operating condition could send the system into an unstable zone. Also, many combustors have many naturally occurring modes [16, 17] that can be easily excited under many operating conditions, causing difficult control. Sometimes quelling one unstable mode can cause another
10 | P a g e

Project MP4079/B427 | Calum Alexander Clarke
to be unstable so additional control is required. In certain cases control can be rendered completely ineffective [18].
2.2.1 – Open Loop Control
Open loop controllers do not use feedback, and are therefore more stable and simpler to build than closed loop controllers. The aim of an open loop active control system is to disrupt the positive feedback cycle which causes the unstable oscillations. It does this through many different ways such as those described in [19, 20].
Figure 2-3 : Open Loop Control System
Commonly, a trial & error approach is used to decide upon the ideal operating conditions of open loop controllers as the controller characteristics tend to be independent of the unstable parameter values.
2.2.2 – Closed Loop Control
Closed loop controllers make use of feedback to control the unwanted signal perturbations. These controllers add an out of phase signal to the system in the manner described earlier by the Rayleigh criterion. A typical closed loop control system is shown in figure 2-4. This is the type of active control strategy adopted in this project.
Figure 2-4 : Closed Loop Control System
11 | P a g e

Project MP4079/B427 | Calum Alexander Clarke
In the above diagram, y denotes the output from the combustor where y = A*sin(ωt + Φ). This is measured by the sensor and the information passed onto the algorithm which processes the data and filters it passing on the deconstructive wave, y’, to the actuator which produces the signal and feeds it back into the system. The input to the combustor is denoted by x’, which is the combination of y’ from the actuator and x, the original input to the system. In terms of sensors and actuators used, the actuator is most commonly a loudspeaker to counteract acoustics or a fuel injector to counteract the heat release perturbations. The sensor used could be a pressure transducer or microphone. See [11, 20, 21] for more details on specific actuators and sensors.
A major advantage of using an active control approach is that once perfected it can be adapted easily for use in different systems, saving cost and time. As well as this, preliminary use of an active control system to allow control in emergencies or during times of refurbishment to a system can be accomplished.
Another method which has been looked into is a model based algorithm approach where the algorithm attempts to predict and account for actuator & flame dynamics, combustor geometry and time delays. This method has found some success in the laboratory where the parameters it sets out to predict are already known [24] but in terms of a practical combustor this method is not yet feasible [25]. However, a fast multiple-mode observer was developed at Georgia Tech [26] in which the main problems of the above methods are bypassed. This observer is able to control one unstable mode without causing the destabilisation of another. More information on this particular observer can be found at [26].
12 | P a g e

Project MP4079/B427 | Calum Alexander Clarke
2.3 – Passive Control Approach
Before the use of these highly effective algorithms there were other alternative solutions, mainly in the form of physical adaptations to the systems geometry. These included methods such as modifying the fuel spray or injection positions, damping the acoustic energy using acoustic liners or attempting to stabilise the flame through methods such as bluff bodies, backward–facing steps or V-gutters. This kind of strategy is known as passive control. These methods were not as successful as online algorithms but they did have positive effects. The main aim of passive control is to prevent the acoustic excitation from happening or to increase the damping of the unstable modes as much as possible. Sometimes it is used to alter the natural frequency of certain modes in the system in an attempt to stop them occurring [27, 28].
For specific examples, in liquid jet engines the droplet size and distribution can be used to attenuate instabilities [29]. Another example of passive control could be that of the Helmholtz resonators which can be used as passive dampers, adapting the maximum damping frequency by changes in the resonator geometry [30]. Or, an alternative method to discovering and stopping instabilities passively could be one similar to that adopted with the Rocketdyne F1 engine [31, 32, 33] where small explosions are set off outside the chamber when it is off and the result of the explosion passing through the system is analysed in order to see where the ideal stabilising fuel injectors positions would be.
However, as for this project an online identification algorithm is required these options are not considered in detail and are only looked at marginally. Passive techniques are usually more difficult to implement than active ones, and active strategies have wider working ranges and don’t require expensive and timely modifications to the combustion systems in question. Although, every system is different and each requires its own specific strategy to ensure no instabilities occur, it would seem that active control is the cheaper and more compatible way to advance.
13 | P a g e

Project MP4079/B427 | Calum Alexander Clarke
Chapter 3 : Van Der Pol
3.1 – Non-linear Limit Cycle Oscillator
Balthasar van der Pol [34] was a Dutch physicist who is well known for his studies of radio. He is also well known for his discovery of what he called ‘relaxed oscillations’ or which is now more commonly known as a limit cycle. In fact the oscillator named after him is one of the most widely used models of non-linear self-oscillation today. It is used with regards to this project to simulate combustion instability within a combustion chamber. The desired limit cycle waveform is achieved by solving the following nonlinear second order differential equation;
𝑥′′ + μ (𝑥2 − 1)𝑥′ + 𝑥 = 0 (1.1)
This equation describes self-sustaining oscillations where energy is removed from larger oscillations and fed into small ones. The form shown above is the unforced version of the van der pol equation. For the forced equation the right hand side would contain an expression similar to F*cos(ωt), Where F is an external force with natural frequency, ω. This equation can be also derived from the Rayleigh differential equation [35].
When the constant, μ, in the equation is equal to zero the middle region of the equation is removed and it becomes the equation for simple harmonic motion. However, when μ is not equal to zero this middle section describes non-linear damping. This leads to the following occurrences; when ‘x’ is large the system acts as a damped harmonic oscillator but when ‘x’ is small an ‘anti- damping’ effect occurs. This ‘anti-damping’ effect is sometimes known as ‘pumping’ and basically means the signal will be amplified instead of attenuated. It is useful to note that for all positive values of μ a stable limit cycle will be produced. In this thesis different values of μ are used to obtain varying limit cycle inputs to test the controller’s ability to control different frequencies and amplitudes. Some graphs of the van der Pol equation at varying values of μ are shown below.
14 | P a g e

Project MP4079/B427 | Calum Alexander Clarke
Figure 3-1 : Limit Cycle μ = 10 (Large)
The above figures show the results of the van der pol equation at large values of μ. Figure 3.1 shows u(1) and u(2) which are the two components x & x’ of the equation solution respectively. Figure 3.2 shows the phase plot at μ = 5, which is x against x’. The graphs were obtained using the runge-kutta approximation which will be described later in chapter 3.5.
Figure 3-2 : Phase plot μ = 5 (Large)
15 | P a g e

Project MP4079/B427 | Calum Alexander Clarke
Figure 3-3 : Phase Plot μ = 0.2
Above is shown the phase plot for small values of μ. The van der Pol equation is used as the simulated combustion chamber as it results in a limit cycle similar to that of a real combustion system when certain values of μ are used. This is shown below in figure 3.4.
Figure 3-4 : Van der Pol (μ = 0.2) vs Real Combustion System
16 | P a g e

Project MP4079/B427 | Calum Alexander Clarke
The upper graph shows an example of the inside of a combustion chamber reaching instability and the lower figure shows a solved van der Pol equation plot with a value of μ = 0.2. The oscillator produced from the van der Pol equation doesn’t specifically follow the model of a combustion chamber but it does provide the basis to model an active control system. As can be seen from figure 3.4 (above) the limit cycles are very similar. So it can be deduced that if control can be obtained on the van der Pol limit cycle then it can also be achieved in the real life combustion chamber limit cycle. The following methods are going to be used to look more closely at the van der Pol limit cycles, the methods looked at include the Lindstedt-Poincare method, the Multiple Scales method and the method of Averaging.
These methods rely on perturbation theory which is defined as is the mathematical method of finding an approximate to one problem by starting with the exact known solution of a similar problem. It is only used if the problem is very difficult or impossible to solve exactly. The starting point is ‘perturbed’ by adding a ‘small’ term or ‘perturbation series’ which is used to quantify the difference the current solution has from the exactly solved solution. The technique breaks the problem into solvable and perturbation parts as it works through. This leads to an expression for the desired result in a power series format. The appropriate derivations are shown below as well as some figures to help illustrate the results.
17 | P a g e

Project MP4079/B427 | Calum Alexander Clarke
3.2 – Lindstedt-Poincaré Method
This method is a technique used in perturbation theory to approximate periodic solutions to ordinary differential equations. It is a single perturbation method and is usually applied when regular perturbation approaches fail. The Lindstedt-Poincare method is used here to analyse the limit cycle caused by the van der pol equation. For this method it is further assumed that the van der pol constant is positive and small:
0 < μ << 1 The first step is to introduce a new time scale in the form 𝜏 = 𝜔𝑡 which means the new period of the system is 2π. This further implies that: 𝑥′ = 𝜔−1𝑥̇ 𝑥′′ = 𝜔−2𝑥̈ Applying this to the van der Pol equation (1.1) it is found that; 𝜔2𝑥′′ + 𝜇𝜔(𝑥2 − 1)𝑥′ + 𝑥 = 0 A power series can then be formed; 𝑥𝜇(𝜏) = 𝑥0(𝜏) + 𝜇𝑥1(𝜏) + 𝜇2𝑥2(𝜏) + ⋯ 𝜔𝜇 = 1 + 𝜇𝜔1 + 𝜇2𝜔2 + ⋯ (2.1) (2.2) (2.3) (2.4) (2.5) Where it is found that 𝜔0 is equal to unity because the period, T, is 2π (𝜔 = 2𝜋⁄𝑇) when 𝜇 = 0. Each term added to the series is of a power less than the term before it. The key idea is to expand 𝜔 as a power series in 𝜇 in terms of unknown coefficients 𝜔1, 𝜔2 etc. These unknown coefficients are found by using conditions which cause the expansion to be uniform. So, subbing the expansion equations 2.4 & 2.5 into 2.3 it is found that: (1+𝜇𝜔1 +𝜇2𝜔2 +⋯)2(𝑥′′0(𝜏)+ 𝜇𝑥′′1(𝜏)+𝜇2𝑥′′2(𝜏)+⋯)+⋯ 𝜇(1+𝜇𝜔1 +𝜇2𝜔2 +⋯)[(𝑥0(𝜏)+ 𝜇𝑥1(𝜏)+𝜇2𝑥2(𝜏)+⋯)2 −1]+⋯ (𝑥0(𝜏)+ 𝜇𝑥1(𝜏)+𝜇2𝑥2(𝜏)+⋯)=0 18 | P a g e Collecting terms of the same order, in 𝜇, and tidying up to second order terms; 𝑥′′0 + 𝑥0 = 0 𝑥′′1 + 𝑥1 = 2𝜔1𝑥′′0 + (𝑥02 − 1)𝑥′0 𝑥′′2 + 𝑥2 = −(𝜔12 + 2𝜔2)𝑥′′0 − 2𝜔1𝑥′′1 − (𝑥02 − 1)(𝑥′1 + 𝜔1𝑥′0) − 2𝑥0𝑥1𝑥′0 (2.6) (2.7) (2.8) (2.9) (2.10) So initial conditions are; 𝑥(0) = 𝐴 𝑥̇(0) = 0 Assuming A to be on the first order it is possible to write the following initial conditions: 𝑥0(0)=0, 𝑥1(0)=0, 𝑥2(0)=0... 𝑥̇0(0)=0, 𝑥̇1(0)=0, 𝑥̇2(0)=0... This gives equation (2.11) by solving the solving the equation at the first order: 𝑥0(𝑡) = 𝐴𝑐𝑜𝑠(𝑡) Subbing this into the order of equations at order 𝜇; 𝑥′′1 + 𝑥1 = 12 𝐴2 − 1 − 2𝜔1𝐴𝑐𝑜𝑠(𝜏) + 12 𝐴2cos(2𝜏) (2.11) (2.12) Project MP4079/B427 | Calum Alexander Clarke The time origin can be taken from any point as the limit cycle oscillation is autonomous. Solving equation 2.12 for 𝑥1(𝜏) results in a term that grows secularly with time ‘t’. As this term grows, it will eventually break down when 𝜇𝑥1(𝜏) becomes as big as 𝑥0(𝜏). Therefore, to make the expansion valid the coefficient of cos(τ) is set to 0. This step is known as removal of secular terms and is a key step in perturbation methods which results in; 𝜔1 = 0 19 | P a g e And after removing the secular term it is found that; 𝑥1(𝜏)=1−12𝐴2 +(13𝐴2 −1)𝑐𝑜𝑠(𝜏)+16𝐴2cos(2𝜏) Substituting 2.13 and ?? into 2.8; 𝑥′′2 + 𝑥2 =(−2𝐴+2𝜔2𝐴+56𝐴3)𝑐𝑜𝑠(𝜏)+(𝐴−13𝐴3)(1+𝑐𝑜𝑠(𝜏))+16𝐴3cos(3𝜏) Again the secular term is to be removed, this results in the following expression for 𝜔2; 𝜔2 = 1 − 5 𝐴2 12 And from this the period and frequency can be deduced up to the second order: 𝜔 = 1 + 𝜇2 (1 − 5 𝐴2) 12 𝑇 = 1 − 𝜇2 (1 − 5 𝐴2) 12 (2.13) (2.14) (2.15) (2.16) (2.17) Project MP4079/B427 | Calum Alexander Clarke Given these two equations (2.16 & 2.17) matlab can be used to plot graphs and gain information about the van der pol limit cycle. The Lindstedt-Poincare method supplies information for frequency(ω) and time period(T) which is dependent on the initial condition for x(t) = A. Figure 3-5 : Lindstedt Poincare - A vs T (μ = 0.2) 20 | P a g e Project MP4079/B427 | Calum Alexander Clarke Figure 3-6 : Lindstedt Poincare - A vs T (μ = 1) From the previous two figures it can be seen that as the initial guess, A, increases so too does the resultant time period of the limit cycle, T. It can also be noticed, by comparison of the two figures, that as μ increases so too does the rising rate of T with A, and hence larger periods are found for larger values of the van der Pol constant μ. As well as this, for initial condition of A = 0, the time period is very close to zero. For higher values of μ it converges closer to zero but for low positive values of μ, as it approaches zero, the time period approaches 1. As is suggested from equations 2.16 & 2.17 the time period ‘T’ and frequency ‘ω’ of the limit cycle grow in different directions when plotted against ‘A’. The results of plotting frequency against initial condition ‘A’ is shown in the following two figures. Figure 3-7 : Lindstedt Poincare - A vs ω (μ = 0.2) 21 | P a g e Project MP4079/B427 | Calum Alexander Clarke Figure 3-8 : Lindstedt Poincare - A vs ω (μ = 1) So as the initial condition increases the frequency of the limit cycle decreases. The same occurs when μ is increased except as μ grows the faster ω reacts. It can also be stated that for higher values of μ and initial condition of zero, the frequency will also have a higher value. For huge values of μ the frequency tends towards μ2 for initial condition A = 0. 22 | P a g e 𝑥 ( 𝑡 ) = 𝑋 ( 𝑇0 , 𝑇1 , 𝑇2 ... ) Which can be further expanded to yield; 𝑋(𝑇0,𝑇1,𝑇2...)= 𝑋0(𝑇0,𝑇1,𝑇2...)+𝑋1(𝑇0,𝑇1,𝑇2...)+𝑋2(𝑇0,𝑇1,𝑇2...) Resulting in the derivatives with respect to t becoming; (3.1) (3.2) (3.3) (3.4) (3.5) (3.6) (3.7) Project MP4079/B427 | Calum Alexander Clarke 3.3 – Multiple Scale Method This method is also a single perturbation method used to approximate the solution to perturbation problems such as the one presented by the 2nd order non-linear van der Pol equation. The Lindstedt-poincare method cannot be used to obtain solutions that evolve aperiodically on a slow time-scale. This is where the multiple scale method comes in, in this method ‘slow’ variables are introduced with respect to each time scale of interest. For this method the assumption is made that the solution sought is dependent on various time scales: 𝑇0 =𝑡,𝑇1 =𝜇𝑡,𝑇2 =𝜇2𝑡... Where T0 represents the fast scale and the other variables are the slow ones. So, to write this as power series in ‘x’; Substituting; 𝜕(.) = 𝜕(.) + 𝜇 𝜕(.) + 𝜇2 𝜕(.) 𝜕𝑡 𝜕𝑇0 𝜕𝑇1 𝜕𝑇2 𝜕2(.)=𝜕2(.)+2𝜇𝜕2(.) +⋯ 𝜕 𝑡 2 𝜕 𝑇0 2 𝜕 𝑇1 𝜕 𝑇0 𝜕2𝑋0+𝑋 =0 𝜕𝑇02 0 𝜕2𝑋1+𝑋=−2𝜕2𝑋0 −(𝑋2−1)𝑑𝑋0 𝜕𝑇02 1 𝜕𝑇1𝜕𝑇0 0 𝑑𝑇0 Eq. ?? can be solved for X0 as follows; Subbing in to 3.6; 𝑋0 = 𝐴(𝑇1, ... )𝑐𝑜𝑠𝑇0 + 𝐵(𝑇1, ... )𝑠𝑖𝑛𝑇0 23 | P a g e 𝜕2𝑋1 +𝑋 = (2𝜕𝐴 −𝐴+𝐴(𝐴2 +𝐵2))𝑠𝑖𝑛𝑇 +(−2𝜕𝐵 +𝐵−𝐵(𝐴2 +𝐵2))𝑐𝑜𝑠𝑇 + 𝜕𝑇02 1 𝜕𝑇1 4 0 𝜕𝑇1 4 0 𝐴4 (𝐴2 − 3𝐵2)𝑠𝑖𝑛3𝑇0 + 𝐵4 (𝐵2 − 3𝐴2)𝑐𝑜𝑠3𝑇0 (3.8) Project MP4079/B427 | Calum Alexander Clarke To obtain the ‘slow flow’ the coefficients of the resonant terms are equated to zero, aka the secular terms are removed. This results in the following: 𝜕𝐴 =𝐴−𝐴(𝐴2+𝐵2) (3.9) 𝜕𝑇1 2 8 𝜕𝐵 =𝐵−𝐵(𝐴2+𝐵2) (3.10) 𝜕𝑇1 2 8 The next step involves transforming the equation to polar coordinates and using the following equations; Which results in; 𝐴 = 𝑅𝑐𝑜𝑠𝜑 𝐵 = 𝑅𝑠𝑖𝑛𝜑 𝐴̇ = 𝜇 𝜕𝐴 𝜕𝑇1 𝐵̇ = 𝜇 𝜕𝐵 𝜕𝑇1 𝑅̇ = 𝜇 (𝑅 − 𝑅3) 28 𝜑̇ = 0 (3.11) (3.12) (3.13) (3.14) (3.15) (3.16) Slow flow equilibria correspond to periodic solutions of the system, in the case of van der Pol these are R = 0 and R = 2. Solving equation 3.15 with the initial condition R(0) to obtain equation 3.17, this is shown below. This reveals that the limit cycle is stable because as t tends to infinity, R(t) tends to 2. √1+( 4 2−1)𝑒−𝜇𝑡 𝑅(0) 𝑅(𝑡) = 2 (3.17) 24 | P a g e Project MP4079/B427 | Calum Alexander Clarke The above equation (3.17) means the amplitude of the limit cycle should be equal to 2, and should remain stable for the period of its duration. Below figures are shown to support this theory. Figure 3-9 : Multiple Scales, Initial 2.5, μ = 0.2 Figure 3-10: Multiple Scales, Initial 2.5, μ = 1 Figures 3-9 & 3-10 show that for initial value above 2 the amplitude of the limit cycle will eventually converge to 2. The speed at which this occurs depends on the van der Pol constant μ, for higher values the amplitude will reach 2 faster. This is shown by the differences in the two above graphs. 25 | P a g e Project MP4079/B427 | Calum Alexander Clarke Figure 3-11 : Multiple Scale, Initial 0.5, μ = 1 Figure 3-11 above shows the results of the multiple scale method on the van der Pol equation with regards to the limit cycle amplitude with an initial value of 0.5. This value was chosen because it is near the other periodic solution value found earlier just previous to equation 3.17. It can be seen that the amplitude tends to 2 for initial values smaller than 2 as well. So the limit cycle is stable with amplitude 2. 26 | P a g e In this case; 𝑥̈+𝑥= 𝜇𝐹(𝑥,𝑥̇) 𝜇𝐹(𝑥,𝑥̇)= −(𝑥2−1)𝑥̇ (4.1) (4.2) (4.3) (4.4) It is required to solve Eq. 4.1 into the following form: 𝑥 = 𝑎(𝑡)cos(𝑡 + 𝜑(𝑡)) 𝑥̇ = −𝑎(𝑡)sin(𝑡 + 𝜑(𝑡)) When 𝜇 = 0, equation 4.1 has a solution of the form equation (4.3) with ‘a’ and ‘𝜑’ acting as constants. For small values of 𝜇 it is expected the same form will be achieved but with ‘a’ and ‘𝜑’ now slowly varying function of ‘t’. Differentiate (4.3) and (4.4) to yield: 𝑎̇ ∗cos(𝑡+ 𝜑(𝑡))−𝑎∗𝜑̇ ∗sin(𝑡+𝜑(𝑡))=0 −𝑎̇ ∗ sin(𝑡 + 𝜑(𝑡)) − 𝑎 ∗ 𝜑̇ ∗ cos(𝑡 + 𝜑(𝑡)) = 0 Sub (4.6) into (4.1) gives; −𝑎̇ ∗ sin(𝑡 + 𝜑) − 𝑎 ∗ 𝜑̇ ∗ cos(𝑡 + 𝜑) = 𝜇𝐹(𝑎(𝑡) cos(𝑡 + 𝜑) , −𝑎(𝑡) sin(𝑡 + 𝜑)) The next step is to solve (4.6) and (4.7) for 𝑎̇ & 𝜑̇ ; 𝑎̇ = −𝜇𝐹(𝑎(𝑡) cos(𝑡 + 𝜑) , −𝑎(𝑡) sin(𝑡 + 𝜑))sin(𝑡 + 𝜑) 𝜑̇ =−𝜇𝑎𝐹(a∗cos(𝑡+𝜑),−𝑎∗sin(𝑡+𝜑)cos(𝑡+𝜑) 𝐹(... ) = 𝑎[𝑎2 ∗ 𝑐𝑜𝑠2(𝑡 + 𝜑) − 1]sin(𝑡 + 𝜑) (4.5) (4.6) (4.7) (4.8) (4.9) (4.10) Project MP4079/B427 | Calum Alexander Clarke 3.4 – Method of Averaging The third and final method looked at in analysing the van der pol limit cycle is known as the method of averaging. This method begins by using a system of a more general form such as the following: Where: 27 | P a g e 𝑎̇=−𝜇1 ∫2𝜋𝑑𝜃∗𝑎∗(𝑎2∗𝑐𝑜𝑠2(𝜃)−1)sin2(𝜃) 2𝜋 0 𝜑̇ =−𝜇 1 ∫2𝜋𝑑𝜃∗𝑎∗(𝑎2 ∗𝑐𝑜𝑠2(𝜃)−1)sin(𝜃)cos(𝜃) 2𝜋 0 Where the averages are denoted by; 1 ∫2𝜋𝑑𝜃... 2𝜋 0 The right hand side of equation 4.14 is zero. The averaging for 4.13 can be done using the following expressions: 𝑐𝑜𝑠2(𝜃) = 12 (1 + cos(2𝜃)) 𝑠𝑖𝑛2(𝜃) = 12 (1 − cos(2𝜃)) 1 ∫2𝜋 𝑑𝜃 𝑐𝑜𝑠2(𝑛𝜃) = 1 ∫2𝜋 𝑑𝜃 𝑠𝑖𝑛2(𝑛𝜃) = 1 ,𝑛 = 1,2,3... 2𝜋 0 2𝜋 0 2 1 ∫2𝜋 𝑑𝜃 𝑐𝑜𝑠(𝑛𝜃) = 1 ∫2𝜋 𝑑𝜃 𝑠𝑖𝑛(𝑛𝜃) = 0,𝑛 = 1,2,3... 2𝜋 0 2𝜋 0 12𝜋 12𝜋1 2 2𝜋∫ 𝑑𝜃𝑐𝑜𝑠4(𝑛𝜃)=2𝜋∫ 𝑑𝜃 (2(1+cos(2𝜃))) =⋯ 00 =1 1 ∫2𝜋𝑑𝜃(1+2cos(2𝜃)+cos2(2𝜃))=1(1+1)=3 42𝜋0 428 (4.13) (4.14) (4.15) (4.16) (4.17) (4.18) (4.19) (4.20) Project MP4079/B427 | Calum Alexander Clarke 𝑎̇ = −𝜇𝑎[𝑎2 ∗ 𝑐𝑜𝑠2(𝑡 + 𝜑) − 1]sin2(𝑡 + 𝜑) (4.11) 𝜑̇ = −𝜇[𝑎2 ∗ 𝑐𝑜𝑠2(𝑡 + 𝜑) − 1]𝑠𝑖𝑛(𝑡 + 𝜑)cos(𝑡 + 𝜑) (4.12) So far, the procedure has been a variation of parameters in attempt to obtain a particular solution to a non-homogeneous linear differential equation. It can now be said that if 𝜇 is small then so too must 𝑎̇ & 𝜑̇ be small. So, over one cycle the right hand side of equations 4.11 & 4.12 are basically constant. This leads to the option of replacing them by their averages, so these equations then become: 28 | P a g e 1 2𝜋 2𝜋∫ 𝑑𝜃∗𝑎∗(𝑎2 ∗𝑐𝑜𝑠2(𝜃)−1)𝑠𝑖𝑛2(𝜃)=... 0 00 = 𝑎3 ∫2𝜋 𝑑𝜃 (𝑐𝑜𝑠2(𝜃) − 𝑐𝑜𝑠4(𝜃)) − 𝑎 = 𝑎3 (1 − 3) − 𝑎 = 1 𝑎(𝑎2 − 4) 2𝜋0 2 2828 𝑎̇ = 𝜇8 𝑎(4 − 𝑎2) Project MP4079/B427 | Calum Alexander Clarke 𝑎3 2𝜋 = 2𝜋 ∫ 𝑑𝜃 𝑐𝑜𝑠2(𝜃)(1 − 𝑐𝑜𝑠2(𝜃))𝜃 − 2𝜋 ∫ 𝑑𝜃 𝑠𝑖𝑛2(𝜃) = ⋯ 𝑎 2𝜋 (4.21) (4.22) The above equation 4.22 is solvable by separating the variables which results in a final solution shown in equation 4.23. For intermediate steps see [22]. 𝑎 ≈ 2 − 𝑒−𝜇𝑡 (4.23) This equation presents a formula for the amplitude of the van der Pol limit cycle. When plotted against a time, t, it should be seen that the graph rises from 1 until the value of 2 is reached and then proceed at this value indefinitely. This is because for t = 0 it can be seen that the formula will simplify down to a = 2-1 which equals 1, the smallest obtainable value. For large values of ‘t’ the equation will converge on 2 and the expression involving time and the van der Pol constant will be settled on very small values close to zero. The formula was presented in code form into matlab and the results attained are shown below. Figure 3-12 : Method of averaging – a vs t (μ=0.2) 29 | P a g e Project MP4079/B427 | Calum Alexander Clarke It can be seen from figure 3-12 that at around t=30 the graph has reached its final value of 2 as predicted. By increasing μ it is possible to increase the convergence time of the graph. A figure of this is shown below for μ = 1. The same works in the opposite manner whereby decreasing μ the graph takes much longer to reach its ultimate value of 2. Figure 3-13 : Method of averaging – a vs t (μ=1) From the three perturbation methods looked at the van der Pol limit cycle is now less unknown. It is now known that the equation causes a stable limit cycle of amplitude 2. The frequency and time period of the limit cycle depends on the initial condition of x(t) or what is also referred to as U(1). These methods are useful for finding out information about the oscillations, especially if computer software is unavailable then such methods can be used. Of course in this case matlab was used to produce graphs of information for much quicker results than the traditional method. The next section will discuss he runge-kutta approximation which is used to gain a very accurate approximation of the van der Pol limit cycle. 30 | P a g e Project MP4079/B427 | Calum Alexander Clarke 3.5 - Runge-Kutta Approximation The Runge-Kutta method is a numerical technique used to solve first order ordinary differential equations. The fourth order runge-kutta is one of the most accurate numerical methods available and hence chosen in this project. The runge-kutta method will be used to approximate the van der Pol limit cycle and act as the simulated combustor. As stated previously the van der pol equation is of second order and hence must be adapted into the form of a first order differential equation to be solved using this method. The runge-kutta is an iterative procedure of the form shown on the following page. Equation 5.1 represents the iterative equation which finds the next estimate of the root given the previous (or initial) guess and the coefficients k1, k2, k3 & k4, while equation 5.2 updates the value of time or steps taken. The following four equations are the iterative steps required to update equation 5.1. The step size is denoted by ‘h’. For the first run an initial value for yn and tn is required. For n = 0, 1, 2, 3... 𝑦 𝑛 + 1 = 𝑦 𝑛 + h6 ∗ ( 𝑘 1 + 2 𝑘 2 + 2 𝑘 3 + 𝑘 4 ) 𝑡𝑛+1 = 𝑡𝑛 + h 𝑘1 =𝑓(𝑡𝑛,𝑦𝑛) 𝑘 2 = 𝑓 ( 𝑡 𝑛 + h2 , 𝑦 𝑛 + h2 𝑘 1 ) 𝑘 3 = 𝑓 ( 𝑡 𝑛 + h2 , 𝑦 𝑛 + h2 𝑘 2 ) 𝑘4 =𝑓(𝑡𝑛 +h,𝑦𝑛 +h𝑘3) (5.1) (5.2) (5.3) (5.4) (5.5) (5.6) Other runge-kutta techniques exist such as the 1st, 2nd and 3rd order methods but the fourth order is used in this project as it is the most accurate. By online repeating the above procedures for a number of samples, n, of the input signal, its characteristics can be deciphered and then reconstructed accurately. A diagram is shown here to help visualize what is meant by this using a simple sinusoidal as example. 31 | P a g e Project MP4079/B427 | Calum Alexander Clarke Figure 3-14 : Sinusoidal Reconstruction So a number of samples, n, will be taken of the signal before t = t0, and used as a basis to predict the proceeding movements of the signal at t > t0. The samples taken are shown by the black dots on the sinusoidal wave previous to t = t0. The prediction is shown afterwards by the dotted line. The more samples taken the more accurate the reconstructed wave will be but with a larger number of samples also comes longer estimation time. So a trade-off between accuracy and estimation time is required. Of course a simple sinusoid is not used in this project so the sampling will be more complex than that shown in the previous figure. Below in figure 3-15 is shown the van der pol limit cycle used in this project as approximated by the runge-kutta method.
Figure 3-15 : Runge-Kutta Van der Pol Limit Cycle (μ = 0.2)
32 | P a g e

Project MP4079/B427 | Calum Alexander Clarke
Chapter 4 : Online Identification Algorithm 4.1 – Standard Least Means Square Algorithm
The Least Means Square algorithm has many forms, commonly there is the standard LMS, normalized LMS, leaky LMS and the sign LMS. The two explored in this project are the standard LMS and the Normalized (X-Filtered) LMS. The standard LMS algorithm deploys an adaptive filter whose coefficients are chosen to minimize the error that a key variable has from its desired value. In the case of this project the algorithm takes the output signal, ŷ, from the filter and compares it with the real (desired) output, y. The error is found from the following equation.
𝐸(𝑛) = 𝑌(𝑛) − 𝑌̂(𝑛) (6.1)
Where ‘n’ denotes the sample number at that current time iteration. The algorithm then acts to reduce this error with each sample batch. The controller works by using a sensor to measure the happenings inside the combustor and hence gaining the input for the algorithm, x(n). In this case, the input to the controller is the limit cycle caused by the van der pol equation. The algorithm then calculates an estimate of the signal. Although in this case the input is known, in a real life system it is not and hence an accurate estimation of the input value is required.
Figure 4-1 : LMS Adaptive Filter
33 | P a g e

Project MP4079/B427 | Calum Alexander Clarke
Figure 4-1 shows the basic layout of the LMS adaptive filter. The box on the left is the adaptive filter which produces the estimate result ŷ(n). The box on the right is an ordinary filter which the input data is sent through to produce Y which is the value compared with ŷ to produce the error output e(n). Each iteration of the input is sent through both filters and compared and with each iteration the result is fed back into the adaptive filter in order to produce a closer estimate than the previous iteration. This is continued even after convergence is reached to ensure future input activity is accurately predicted. The input to the algorithm is the runge-kutta samples obtained earlier. This convergence method is based on the use of a Finite Impulse Response (FIR) filter.
A FIR filter is a digital filter with a finite impulse response. They are also sometimes known as non-recursive filters as they do not involve any feedback, unlike the Infinite Impulse Response (IIR) filter which does adopt the use of feedback. When designing a FIR filter the ideal characteristics are sought but never achieved, as close as this to possible is desirable though of course. As the filter order increases so does the models likeness to that of an ideal filter. But as well as this the complexity and processing time of the filter is also increased. Figure 4-2 below shows an example of an ideal low-pass filter approximation.
Figure 4-2 : FIR Filter Approximation
Where the gold line shows the ideal result and the red and blue lines show ideal approximations.
34 | P a g e

Project MP4079/B427 | Calum Alexander Clarke
For a FIR filter the output is related to the input as follows;
𝑦̃𝑛 = ∑𝑀𝑘=0 𝑎𝑘 ∗ 𝑥𝑛−𝑘 (6.2)
Where xn denotes the sample inputs and ỹn denotes the output. The sampling time is denoted by ‘n’ and ak represents filter coefficients from k = 0 until k = M, where M = L – 1, and L is the filter length. Equation 6.2 is sometimes expressed in matrix form for ease of implementation, this is shown below.
𝑦̃𝑛 = 𝐴𝑡 ∗ 𝑋𝑛 (6.3)
Where At = (a0, a1… aM) the filter coefficients and Xn = (xn, x n-1 … xn-M)t the inputs. This form is the form used in the matlab code to successfully implement the algorithm. This can be observed in appendix B. Once it is known that the prediction is accurate enough the next step is to attenuate the instabilities by sending the relative information to the actuator so it can produce the required deconstructive signal. A basic diagram of the controller is shown and described below.
Figure 4-3 : Black Box Controller Diagram
Figure 4-3 shows a ‘black box’ controller figure where the input is the measured value inside the combustion chamber obtained from the sensor and the output the deconstructive wave required to cease unstable oscillations. This deconstructive wave will be the same amplitude of the estimate obtained from the algorithm but sent in 180o out of phase.
35 | P a g e

Project MP4079/B427 | Calum Alexander Clarke
4.2 – X-Filtered Least Means Square Algorithm
In the case now of the X-Filtered LMS algorithm, it differs slightly from the standard LMS in that the secondary propagation path is now considered in more detail and included in the algorithm. This describes the path taken from the actuator device back into the system. If the device has a delay or an offset this will have an effect on the results and so any discrepancies must be considered. To do this the actuator is usually tested with the system off and the resultant wave analysed for any errors, however slight. These are then incorporated into the algorithm to make them void. The algorithm still works in the same way as described previously for the standard method in that it acts to reduce the error. A diagram of an adaptive filter using the x-filtered algorithm is shown in figure 4-4. x(n) & y(n) denote the input and output respectively. Ys(n) is the signal from the secondary propagation path and fx(n) is the estimate. D(n) denotes external noise entering the system. E(n) is as described earlier, the resultant signal of the various inputs.
Figure 4-4 : X-Filtered LMS Adaptive Filter
Combustion instability can be compared to active noise control (ANC) problems in many ways. They are both problems caused by unstable oscillations within their relative systems which can cause harm and/or damage to machinery. Both are desired to be avoided at all costs and are commonly solved using identification algorithms. The X-filtered LMS algorithm is a very popular algorithm choice for solving ANC problems because of its robustness and consideration of the secondary propagation path. Hence why it was considered for use in this combustion instability problem.
36 | P a g e

Project MP4079/B427 | Calum Alexander Clarke
Active noise control is the method of reducing noise in a system by introducing ‘anti-noise’, whereby the anti-noise is simply another sound source acting in such a way that it cancels out the unwanted noise. The figure below shows an illustration of ANC. In theory the resultant wave produced from this process should be zero. As can be seen from the figure speaker 1 and speaker 2 are out of phase by 180 degrees, so the waveforms peak at the exact opposite times causing the overall result to be equal to zero. This is very similar to what happens in active control of combustion systems. This method is very useful when operating loud machinery that can be damaging to the ears or even the equipment and the unwanted noise is required to be lessened or removed.
Figure 4-5 : ANC Sound wave Illustration
The resultant output from this cancellation is zero, where the dotted lines in figure 4-5 indicate the waves are still there but can no longer be heard as they cancel out each other.
The mathworks example using the x-filtered algorithm to attenuate noise control is very informative; it was this example which was used as the basis for the x-filtered algorithm within this project [36].
37 | P a g e

Project MP4079/B427 | Calum Alexander Clarke
Step 1 using the X-Filtered LMS algorithm is to generate a synthetic secondary propagation path. When the system is on the secondary propagation path is not known until it is occurring so if there is a time delay or other offset error it has to be taken into consideration before the system is in operation. A signal must be generated from the actuator before the system is in use to discover these discrepancies.
Step 2 is estimating and adapting the secondary propagation path. As this signal is not known until it is happening an accurate estimate must be made. This step will ‘observe’ the synthetic signal created in step 1 and produce an estimate of its output. The estimate must be as accurate as possible to the real value of the wave so that when the system is in action it is known how the secondary propagation path will perform. Using this information the actuators signal is now altered to account for any results found. For example, an appropriate filter length will be chosen to accommodate for any time delay present. If the filter length is too long or short the delay could throw the results off initially causing some problems.
Step 3 is to introduce the primary propagation path which is simply the van der pol limit cycle described earlier. Following this the final step is to activate the control algorithm and run the program. It should be noted that the level of detail into this algorithm will be much less than that of the standard LMS. This algorithm is mainly used to compare the results of extra noise being added on top of the input signal. As it is popular in ANC applications it suggests it will deal with the additional noise better. The results in chapter 5 shall reveal if this is true.
38 | P a g e

Project MP4079/B427 | Calum Alexander Clarke
Chapter 5 : Results & Discussion
5.1 – Standard LMS
With the control algorithms complete and the combustion instability simulated the code was run and the following results were produced. Firstly, the standard LMS algorithm will be analyzed with the results from this section shown before the results achieved from the X-filtered algorithm are shown. Figure 5.1 shows the estimation accuracy of the standard LMS algorithm.
Figure 5-1 : Estimation Accuracy
The figures above show the estimation function of the controller in action. It can be seen that the algorithm converges very quickly on the input waveform and matches it more or less exactly by around 50 iterations. The upper figure shows the limit cycle as the output from the FIR filter (Y) in green, the Estimate (Ŷ) in blue and the difference between the two, the error (E) in red. Recalling equation 6.1 from section 4.1 which describes the error value in relation to the estimation and the output. It can be seen that while the algorithm converges at the start of the graph the error begins to increase before falling when the values are estimated accurately. The lower figure shows specific iterations and their estimates. It can be seen that the estimates are very close to the real values and only differ by a maximum of around 0.1. This estimation will increase as the time goes on and the error is reduced.
39 | P a g e

Project MP4079/B427 | Calum Alexander Clarke
Figure 5-2 : Instant Controller Activation
Figure 5-2 shows the controller in action from the instant the limit cycle begins. It can be seen that the control of the limit cycle is very efficient and only takes just over 50 iterations to fully converge and attenuate the signal. The estimate continues to measure the would-be value of the unstable oscillations throughout the whole of the above graph. During the first 50 iterations it can seen that the output is fighting against the algorithm to try and build into its full cycle. The waveform reaches its maximum value at around 30 iterations where the algorithm begins to take control and the actuation kicks in. From here the signal slowly depreciates to zero. For the following examples a much higher number of iterations are used to ensure the curves are measured as accurately as possible. A value of 5000 iterations is used from now on.
In the next graph, figure 5-3, the controller is activated at t = 2500 which is the midpoint of the simulation. This is to show that the limit cycle reaches its full potential if left undeterred and also how quickly the controller can attenuate the full limit cycle signal. As can be seen, the controller estimates the values very quickly. This is because once the limit cycle is in full form it is very repetitive and the algorithm can notice this and estimate it accurately within short time. The attenuation takes roughly 1000 to get completely to zero. However, safe levels are reached long before this, probably nearer to 100 iterations after controller activation.
40 | P a g e
Signal Amplitude

Project MP4079/B427 | Calum Alexander Clarke
Figure 5-3 : Controller Activation @ t = 2500 iterations, μ = 0.2
The x-axis shows the number of iterations passed and the y-axis shows the signal amplitude. As is the way of the LMS algorithm the output finally converges to the value of the error, E, which is close to zero. The exact values for Y, Ŷ and E can be viewed using the matlab code, a command has been included in the code to perform this if necessary. The attenuation presented here is very acceptable and would stabilize the system in enough time to avoid damage, of course assuming the controller is activated from t=0 and not 2500!
In order to get the controller to start at 2500 iterations a second vector had to be defined and the controller has to work through this. These details can all be found in the code attached in appendix B. Note: the function ‘f.m’ must be in the same path as the main code for it to work.
The filter length used in the controller was 10. If a higher value the filter looks back too far and can cause the estimation and attenuation to be inaccurate. The same applies for a small filter length where not looking back far enough causes the estimation to be not as accurate as it could be. A value of 10 was found to be the ideal length. In terms of step size, the ideal value was found to be 0.05. Too high a value causes a complete misrepresentation of the signal as the sampling points can be skipped completely. Too low results in slower response time as shown in figure 5- 4 below but also a better final convergence value.
41 | P a g e

Project MP4079/B427 | Calum Alexander Clarke
Figure 5-4 : Step size = 0.001
The above graph shows the attenuation with a step size value of 0.001 in the LMS adaptive filter. This causes the response time of the system to increase, now resulting in the algorithm taking less time to attenuate the unwanted limit cycle while also increasing the accuracy of the final convergence. But what about if the input oscillations are of a higher or lower amplitude or frequency? Will this affect the controller’s performance? The following figures show some results from putting these questions to the test. Figures 5-5 & 5-6 display the results from changing the frequency of the limit cycle. This was done by altering the upper edge limits (b) of the iterations to a smaller vaue for a lower frequency and vice versa. Figure 5-5 shows the lower frequency result.
Figure 5-5 : Lower Frequency Limit Cycle (b = 62.5)
42 | P a g e

Project MP4079/B427 | Calum Alexander Clarke
It is clear the limit cycle is taking a lot longer to reach its full amplitude at a lower frequency, however the controller converges its estimation in a very short amount of iterations, roughly 100. The longer estimation time is due to the lower frequency sampling the filter with more spaced out points making it harder to acurately estimate quickly. It can also be seen that the estimate continues to accurately predict the resultant limit cycle even after attenutation. This was the result of an upper edge four times less than previously used. Figure 5-6 below shows the higher frequency waveform (upper edge four times greater, b = 1000) and as expected the controller reacts quicker than before to this specific waveform. The higher frequency wave gives the filter more points closer together to measure and hence a convergence can occur in less time than that of lower frequencies. It is difficult to see just how quickly the controller reacts so a zoomed in graph is supplied in figure 5-7.
Figure 5-6 : Higher Frequency Limit Cycle (b=1000)
Figure 5-7 : Higher Frequency Limit Cycle Zoomed Graph
43 | P a g e

Project MP4079/B427 | Calum Alexander Clarke
The following figures 5-8 & 5-9 show the results from altering the amplitude of the ‘μ′ constant in the van der Pol equation, which changes the shape and amplitude of the limit cycle. Note the original value for μ is used is 0.2.
Figure 5-8 : Limit Cycle Control (μ = 1)
In figure 5-8 the adapted value of μ causes the limit cycle to enter its relaxed oscillatory state almost immediately and as of such is quite predictable early on leading to a very fast estimation. This also leads to a fast convergence with the signal being attenuated completely by around 150/200 iterations.
In figure 5-9 the value for μ is 0.02, which is a factor of ten smaller than the value used for the initial results shown in figures 5-2 & 5-3. The limit cycle now increasing constantly and at a lower rate than before. The attenuation is still performed very quickly given the change in signal amplitude. It should also be noted that the overall amplitude of the signal is much smaller than previously when the controller kicks in.
44 | P a g e

Project MP4079/B427 | Calum Alexander Clarke
Figure 5-9 : Limit Cycle Control (μ = 0.02)
Once the properties of the limit cycle had been adapted and the results of the controllers reactions analyzed an additional noise source was added to the input signal to act as an external disturbance. This will test the controller’s ability to estimate and attenuate a noisy limit cycle input. This is also the result which will be compared to the X-filtered algorithm. Figure 5-10 shows the resultant graph. The additional noise was added by using the matlab code ‘KK = 0.5*randn(size(X))’ which produced a matrix KK of random values of size ‘X’. This matrix was then added to the original input ‘X’ to produce a noisy version of the original limit cycle
Figure 5-10 : Limit cycle control with additional noise (μ=0.2)
45 | P a g e

Project MP4079/B427 | Calum Alexander Clarke
Figure 5-11 : Limit cycle control with additional noise (μ=0.2) Zoomed Graph
As can be seen from figure 5-10 the resultant signal is the image of a noisy limit cycle. The limit cycle can be seen as less stable than before and hence causes some problems for the controller. After controller activation there is two major peaks of activity which are of concern. Figure 5-11 shows that these occur at the peaks of the limit cycle values. It can be seen from figure 5-10 that the oscillations grow beyond previously reached values. The attenuation is complete successfully after 1000 iterations and seems to retain stability after this point. If the initial peaks could be reduced to lower values the control could be deemed acceptable. At this point the spikes present after control activation are a worry, however this is because the controller is activated when the limit cycle is already in full form. Activation from t=0 would avoid this issue. Also increasing the number of iterations or perhaps reducing step size could resolve this issue more. This figure shows the noise added multiplied by a half. Following this noise was added in greater quantities to the signal until the controller lost control completely. This occurred around the point where noise was added with a standard deviation greater 0.75 and a mean of 0. The resultant graph is shown in figure 5-12. It can be observed that the controller has lost control of the oscillations and cannot actuate it accurately. Any values of noise added above this are completely out with the controller capabilities and the result is intense oscillations instead of accurate attenuation.
46 | P a g e

Project MP4079/B427 | Calum Alexander Clarke
Figure 5-12 : Increased noise until loss of control
In the case of above figure, the oscillations occurring after controller activation can still be seen reaching values as high as 15, which is actually 7.5 times the amplitude of the original limit cycle. This is not safe for the system at all, this occurs after a long stretch of attenuation where the threat seems to be all but contained as well. This confirms the controller has lost control and the signal amplitudes are unpredictable from here onwards. The controller is no longer effective.
Figure 5-13 : 0.8*Additional Noise, 0.001 Stepsize
However, as can be seen from figure 5-13, by reducing the step size of the LMS filter one can increase the controller’s stability. Although the controller now takes longer to completely attenuate the instabilities it handles the additional noise much better. Noise was again added until controller failure with the new step size value of 0.001. It was found that the controller now successful upto values as high as 10 times of random noise added. This is shown below.
47 | P a g e

Project MP4079/B427 | Calum Alexander Clarke
Figure 5-14 : 10*Additional Noise, 0.001 Stepsize
From figure 5-14 it can be seen that despite the high levels of external noise attached to the input limit cycle signal the LMS controller with a step size of 0.001 is standing strong and performing very effectively to attenuate the signal in a short period of time. It takes the controller around 150 iterations to actuate the signal completely in this case.
Overall, it can be said that the standard LMS is an effective control algorithm strategy for combustion instability. With these results it proves itself efficient enough to ensure safe operation of the system at hand. Of course further slight modifications could be made to enhance the algorithm further, if time allowed. Especially with considerations towards increasing the algorithm’s effectiveness with additional noise in place. The controller performed poorly with noise added to the system and this should be a consideration of any future work involved with this algorithm. Although reducing the step size increased the capabilities of dealing with noise it affected the attenuation time in a negative way, a suitable trade-off is required for the particular system it is to function within. Of course the results found from figure 5-14 would indicate that the controller is sufficient enough for use in a real life system. It can also be seen from the results that the controller works best with high frequency limit cycles. However, care must be taken to ensure the correct filter length and step size are used. Once again, additional noise is a problem which must be addressed properly. However, the X-filtered algorithm is expected to work better under these conditions so this was put to the test and the results are shown and analyzed in the following pages.
48 | P a g e

Project MP4079/B427 | Calum Alexander Clarke
5.2 – X-Filtered LMS
The results from the X-filtered LMS algorithm will now be shown and discussed with comparison to those achieved using the standard LMS method. As mentioned previously this algorithm was used as it is a popular algorithm choice for ANC applications and it is similar to the standard LMS algorithm already implemented. This algorithm will not be looked at it in as much detail as the previous one. The algorithm will be tested with the addition of external noise and compared to the result obtained from the standard LMS algorithm with the addition of external noise. The following figures show the accuracy of estimation of the second algorithms secondary propagation path.
Figure 5-15 : Accuracy of Secondary Prop Path Estimation – X-Filtered Algorithm
49 | P a g e

Project MP4079/B427 | Calum Alexander Clarke
The signal does not start instantly as a delay was added intentionally to see if the estimation would accurately pick it up, which it did. As for the rest of the signal the estimation is very accurate, apart from at the tail of the signal which is acceptable as the value as too low to cause concern. Figure 5-15 shows two separate figures regarding the estimation of the secondary path and it can be seen from both figures that the estimation of the X-filtered algorithm is very accurate. In this case the sampling frequency (Fs) is 8kHz, and the number of iterations (N) taken is 500 which leads to a filter length of 0.0625 seconds. It can be seen from the graph that this is the value where the waveform ceases to continue. Such a small value for iterations is used as this section only concerns the secondary propagation path. The time scale of this algorithm covers the range from 1/Fs to N/Fs in steps of 1/Fs.
In designing an estimate of the secondary propagation path the quickest convergence and best final result was sought after. To ensure this a step size of 1 was used and a filter length of 250 was chosen. Higher values for filter length would mean a more accurate final result but a longer convergence time and a lower value would result in a quicker convergence but with a less accurate final result. A compromise had to be made and the value of filter length which seemed to provide the best of both results was 250.
Figure 5-16 : Results of Attenuation (X-Filtered)
Figure 5-16 shown above displays the results of the x-filtered controller being activated from t=0 with additional noise present in the system over 5000 iterations. In this case it seems as though the controller has complete control over the noise apart from the initial spurt of peak oscillations. Figure 5-16 shows a zoomed in graph over the first 1000 iterations.
50 | P a g e

Project MP4079/B427 | Calum Alexander Clarke
Figure 5-17 : Results of Attenuation (X-Filtered) Zoomed Graph
So it can now be seen clearly that the problem lies with the rate of estimation for this algorithm. The time taken to accurately converge and implement the algorithm is too great and as a result the limit cycles rein out of control for the first 500 iterations. The oscillations are eventually controlled and from figure 5-16 it is known that they are controlled for good, with no more unstable oscillations further along the line.
It is quite obvious that the x-filtered algorithm deals with the additional noise much more effectively than the standard LMS algorithm without its modified step size. Although its convergence rate is much less efficient the extra noise doesn’t cause too much of an issue, whereas for the standard algorithm the additional noise causes major problems if the correct step size is not used. A reference for future work involving these algorithms would be to tweak the x-filtered algorithm to increase convergence rate and actuation time while still keeping the stability gained here. Or to continue to tweak the standard LMS algorithm so as to obtain better noise control. Overall, the standard LMS algorithm is perfectly suitable with a small step size value. It manages to control additional noise added as well as attenuate the instabilities in good time. The x-filtered algorithm seems to cope with the noise better but the current convergence rate is not fast enough and further modifications are required. The standard LMS is preferred as is it simpler and easier to modify.
51 | P a g e

Project MP4079/B427 | Calum Alexander Clarke
Chapter 6 : Conclusion
This project required the student to write an algorithm, using matlab, which would maintain a combustion system at a stable working condition by controlling the heat-pressure oscillations present inside the simulated combustor. The combustion instabilities were simulated by solving the second order non-linear van der pol equation and using this to create a limit cycle similar to that which occurs during combustion instabilities. The Least Means Square algorithm was then introduced to measure the properties this input wave and estimate them as accurately as possible. The algorithm then sent information to an actuator allowing it to send out a deconstructive wave to quell the instabilities. The standard and the X-filtered LMS algorithms were both used for the case of a noisy input and compared.
Chapter 1 gave a brief introduction to the project and the work done before chapter 2 covered a literature review of relevant background information and theory. Chapter 3 described how the inside of an unstable combustion chamber was simulated by utilizing the van der Pol equation. This section also covers various numerical techniques which are used to analyze the limit cycles caused by the non-linear second order differential equation. Chapter 4 discussed the LMS algorithms and how they were implemented for use in this project. In chapter 5 the results from various tests and procedures were shown and discussed. It was found that although the standard LMS algorithm performed very well under normal conditions when additional noise was included some tweaks were required to maintain sufficient control. The X-filtered algorithm performed well with additional noise but lacked the quick attenuation time and high accuracy of the standard algorithm. Overall, the algorithms were implemented successfully and the system created was managed to be kept within stable working conditions, hence the project was a success.
On a personal note I enjoyed the challenge presented by this project assignment. I feel that by accepting this project at the beginning of the semester I was setting myself up for a challenge and although I do feel I have completed that challenge I also feel I could have done more, as is suggested in the further work section. There was many engineering topics covered that I had not ventured into before, mainly combustion instability itself. This is a topic which I have not the
52 | P a g e

Project MP4079/B427 | Calum Alexander Clarke
privilege of exploring before and I am glad that this project has given me such an opportunity. In terms of matlab as well, I feel much more competent with its use now than I did at the beginning of the project. On top of this topics such as thermodynamics, modelling and control, engineering mathematics and dynamics were covered throughout the course of the project and my horizons have definitely been widened in all these respects. It was especially interesting to read about specific real life events such as the Rocketdyne F-1 engine mentioned previously in the report (Chapter 2, Section 2.3).
The different technique areas covered by this project include literature review/lab investigation, software development & application and numerical analysis & mathematical formulation. As well as this the project relates strongly to aerodynamics and control theory.
53 | P a g e

Project MP4079/B427 | Calum Alexander Clarke
6.1 – Further Work
For initial ideas towards further work regarding this subject another algorithm such as the one described by Tongxun Yi and Ephraim J. Gutmark in Adaptive Control of Combustion Instability Based on Dominant Modes Reconstruction[15] could be attempted and then compared to the LMS methods. This algorithm is more complex than the LMS methods but also much more accurate. It concentrates on a controller which can identify and attenuate multiple unstable modes of unknown characteristics. Its procedures for frequency estimation are similar in a way to the discrete Fourier transform and controller stability is not an issue as feedback is avoided. This algorithm would be useful for study as it differs from the simpler LMS methods explored in this project.
Another point in terms of future work would be setting up a practical experiment to test the algorithms on a real-life system. Perhaps a simple set up involving a flame and a loudspeaker at first before advancing to more complex apparatus.
The future of combustion instability control seems to be remaining model based. The current methods of using online identification algorithms as active control are very successful and can be adapted for almost any system. This leads to many obvious advantages for engineers working with these systems.
54 | P a g e

Project MP4079/B427 | Calum Alexander Clarke
Chapter 7 : References
[1] – Oyediran, A. Darling, D. Radhakrishnan, K. ”Review of combustion-acoustic instabilities,” NASA Tech Memo, Aug 1995.
[2] – Culick, F.E.C., Yang, V., “Overview of combustion instabilities in liquid-propellant rocket engines,” Progress in Astronautics and Aeronautics.
[3] – Lord Rayleigh. Strutt, J.W. 1878 “On the instability of jets” proceedings of the London mathematical society Vol. 10.
[4] – “Combustion instability due to the nonlinear interaction between sound and flame” Xuesong Wu†, Meng Wang and Parviz Moin
[5] – J.W.S. Rayleigh, “The Theory of Sound”, vol 2. New York, Dover, 1945
[6] – Torres, H., Lieuwen, T.C., Johnson,C.E., Daniel, B.R., Zinn, B.T., “Experimental Investigation of Combustion Instabilities in a Gas Turbine Combustor Simulator,” AIAA Paper # 99-0712, 37th Aerospace Sciences Meeting, Reno, NV, January 1999
[7] – Harrje, D.J., Reardon, F.H.E., Liquid Propellant Rocket Combustion Instability, NASA report SP-194, 1972.
[8] – Peracchio, A.A., Proscia, W.M., “Nonlinear heat-release/acoustic model for thermoacoustic instability in lean premixed combustors,” ASME Paper 98-GT-269.
[9] – Cohen, J.M., Anderson, T.J., “Experimental Investigation of Near-blowout Instabilities in a Lean, Premixed Step Combustor,” AIAA Paper 96-0819, 34th Aerospace Sciences Meeting and Exhibit, Reno, NV, Jan. 1996.
[10] – Richards, G.A., Janus, M.C., “Characterization of Oscillations During Premix Gas Turbine Combustion,” ASME Paper 97-GT-244.
[11] – Active Control of Combustion Instability: Theory and Practice by Anuradha M. Annaswamy and Ahmed F. Ghonie MIT
[12] – A.M. Annaswamy, “Nonlinear modeling and control of combustion dynamics,” in Fluid Flow Control P. Koumoutsakos, I. Mezic, and M. Morari, eds, New York: Springer Verlag, 2002.
[13] – Marble, F.E., Cox, D.W., Jr., “Servo-Stabilization of Low-Frequency Oscillations in a Liquid BiPropellant Rocket Motor,” ARS Journal, 23, 63, 1953.
[14] – Lee, Y.C., Gore, M.R., Ross, C.C., “Stability and Control of Liquid Propellant Rocket Systems,” ARS Journal, 23, 75, 1953.
55 | P a g e

[15] – Adaptive Control of Combustion Instability Based on Dominant Acoustic Modes Reconstruction; Tongxun Yi and Ephraim J. Gutmark
[16] – Billoud, G., Galland, M.A., Huynh Huu, C., Candel, S., “Adaptive Active Control of Combustion Instabilities”, Combust. Sci. and Tech., 1992, Vol. 81 pp. 257-283
[17] – Neumeier, Y., Markopolous, N., Zinn, B.T., “A Procedure for Real-time Mode Decomposition, Observation, and Prediction for Active Control of Combustion Instabilities”, IEEE Paper 97318, Conference on Control Applications, Hartford, CT, Oct 5-7, 1997.
[18] – Fleifil, M., Annaswamy, A.M., Hathout, J.P., Ghoniem, A.F., “The Origin of Secondary Peaks with Active Control of Thermoacoustic Instability,” Combustion Science and Technology, v133, pp. 227-260, June 1998.
[19] – McManus, K.R., Vandsburger, U., Bowman, C.T., “Combustor Performance Enhancement Through Direct Shear Layer Excitation,” Combust. Flame, v82, pp. 75-92, 1990.
[20] – Richards, G.A., Yip, M.J., “Combustion Oscillation Control by Cyclic Fuel Injection”, J Eng Gas Turbines Power Trans ASME v119 n2 pp. 340-343, 1997
[21] – Seume, J.R. Vortmeyer, N. Krause, W. Hermann, J. Hantschk, C.-C. Zangl, P. Gleis, S. Vortmeyer, D. Orthmann, A. “Application of active combustion instability control to a heavy duty gas turbine,” J Eng Gas Turbines Power Trans ASME v 120 n 4 Oct 1998.
[22] – Van der Pol Oscillator Draft – https://docs.google.com/viewer?url=http%3A%2F%2Fwww.phys.uconn.edu%2F~rozman%2FC ourses%2FP2200_13F%2Fdownloads%2Fvanderpol%2Fvanderpol-oscillator-draft.pdf
[23] – Adaptive control of combustion instabilities using real-time modes observation, academic thesis, Clifford Edgar Johnson
[24] – Hathout, J.P., Annaswamy, A.M., Fleifil, M., Ghoniem, A.F., “Active Control of Thermoacoustic Instability Using a Model-Based Approach”, Proceedings of the 1997 ASME International Mechanical Engineering Congress and Exposition, November, 1997, Dallas, TX
[25] – Johnson,C.E., Neumeier, Y., Lieuwen, T.C., Zinn, B.T., “Experimental Determination of the Stability Margin of a Combustor Using Exhaust Flow and Fuel Injection Rate Modulations,” 28th International Symposium on Combustion, Edinburgh, 2000.
[26] – Neumeier, Y., Zinn, B.T., AIAA Paper 96-0758, 34th Aerospace Sciences Meeting and Exhibit, Jan 15-18, 1996, Reno, NV.
[27] – McManus, K.R. Poinsot, T. Candel, S.M., “Review of active control of combustion instabilities,” Prog Energy Combust Sci v 19 n 1 1993 p 1-29
Project MP4079/B427 | Calum Alexander Clarke
56 | P a g e

Project MP4079/B427 | Calum Alexander Clarke
[28] – Schadow, K.C. Gutmark, E., “Combustion instability related to vortex shedding in dump combustors and their passive control,” Energy Combust Sci v 18 n 2 1992 p 117132.
[29] – http://www.nvbv.org/index.php?option=com_edocman&task=document.viewdoc&id=318&Ite mid=271
[30] – Tuned passive control of combustion instabilities using multiple Helmholtz resonators; Dan Zhao; A.S. Morgans
[31] – http://www.rocketdynearchives.com/engines/f1.html
[32] – http://www.thespacerace.com/forum/index.php?topic=1054.0 [33] – http://www.astronautix.com/engines/f1.htm
[34] – http://www-history.mcs.st-and.ac.uk/Biographies/Van_der_Pol.html
[35] – http://mathworld.wolfram.com/RayleighDifferentialEquation.html
[36] – Mathworks – X-filtered LMS algorithm for Active Noise Control.
Weblink: http://www.mathworks.com/help/dsp/examples/active-noise-control-using-a- filtered-x-lms-fir-adaptive-filter.html
7.1 – Other useful links & papers
[37] – https://docs.google.com/viewer?url=http%3A%2F%2Fieeecss.org%2Fsites%2Fieeecss.org%2Ffil es%2Fdocuments%2FIoCT-Part4-11CombustionInstability-HR.pdf
[38] – Adaptive Active Control of Combustion Instabilities – G. Billoud and M.A.Galland
[39] – https://docs.google.com/viewer?url=http%3A%2F%2Fdspace.mit.edu%2Fbitstream%2Fhandle %2F1721.1%2F88841%2F45993716.pdf%3Fsequence%3D1
[40] – Minimizing transient energy growth of non-linear thermoacoustic oscillations – Dan Zhao, X.Y. Li
57 | P a g e

Chapter 8 : Appendix
The appendix for this FYP report is found on an attached CD.
Project MP4079/B427 | Calum Alexander Clarke
END OF APPENDIX
58 | P a g e