A N A N A L Y S I S O F T H R E E M O D E L - B A S E D E S T I M A T I O N M E T H O D S F O R D I E S E L E N G I N E C O N D I T I O N M O N I T O R I N G B y Raluca F. Constantinescu B. Eng. (Control Engineering) Bucharest Polytechnic Institute A THESIS S U B M I T T E D IN PART IAL F U L F I L L M E N T OF T H E R E Q U I R E M E N T S F O R T H E D E G R E E OF M A S T E R OF A P P L I E D S C I E N C E S in T H E F A C U L T Y OF G R A D U A T E STUDIES E L E C T R I C A L E N G I N E E R I N G We accept this thesis as conforming to the required standard T H E U N I V E R S I T Y OF BRIT ISH C O L U M B I A Apr i l 1995 © Raluca F. Constantinescu, 1995 In presenting this thesis in partial fulfilment of the requirements for an advanced degree at the University of British Columbia, I agree that the Library shall make it freely available for reference and study. I further agree that permission for extensive copying of this thesis for scholarly purposes may be granted by the head of my department or by his or her representatives. It is understood that copying or publication of this thesis for financial gain shall not be allowed without my written permission. Department of fyjiheeOrlQ The University of British Columbia. Vancouver, Canada Date ApHl^ mr DE-6 (2/88) Abst rac t In the context of engine fault detection and isolation, we focus on three main aspects: engine modeling and validation, engine parameter identification, and cylinder pressure waveform reconstruction. The problem of diesel engine modeling is solved using Euler's equation, under the assumption that the crankshaft is perfectly rigid. The model inputs are represented by the cylinder pressures and the output is the flywheel angular velocity fluctuation. Simulation results obtained with M A T L A B show a root mean square (RMS) error of 0.0891 rad/sec between the estimated and the actual crankshaft angular velocity fluctu-ation for the normal operating condition. The elimination of a strong sinusoidal trend for the faulty condition results in a RMS error range of 0.0973 rad/sec to 0.1836 rad/sec. The identification methods involved the off-line standard least-squares technique, the recursive gradient estimator and the on-line least-squares estimators with exponential for-getting. The parameters of interest are engine inertia, and torque fluctuation. The RMS velocity error for the normal operation has a value of 0.0559 rad/sec, which represents approximately 30% improvement over the initial result, before identification. The issue of cylinder pressure waveform reconstruction is addressed. The inverse dynamics are solved by redefining the system input as the torque due to gas pressure. The cylinder pressure waveform is approximated by an impulse-like periodic function. We considered the problem of fault detection and isolation. The procedure uses 6 pressure templates. The estimated pressure variations are obtained using a standard least-squares technique. An under-fueling fault in the i-th firing cylinder can be determined exactly by the minimum value of the estimated pressure variation. ii Using the pressure correction we are able to improve the estimation of the gas pressure torque. The RMS torque error for the normal operation reduces to 99.24 Nm. The case of an under-fueling fault is characterized by a reduced RMS torque error range of 85.7 Nm to 198.1 Nm. The pressure waveform reconstruction is characterized by a RMS pressure error range of 0.155 M P a to 0.277 MPa for the normal operating condition. For each of the six under-fueling faults, the pressure waveform corresponding to the faulty cylinder is reconstructed. The RMS error range is of 0.155 MPa to 0.386 MPa. iii Table of Contents Abs t rac t i i Lis t of Tables v i L is t of Figures v i i Acknowledgements v i i i 1 In t roduct ion 1 2 Background 5 2.1 Engine Model ing 6 2.2 Diesel Engine Characteristics 7 2.2.1 Data Acquisit ion 9 2.2.2 Engine Control 9 2.3 Fault Detection and Isolation in Control Systems 10 2.4 Trends in Engine Condit ion Monitor ing 13 2.5 Residual Generation and Decision Making 16 2.5.1 State Estimation-Based Techniques 17 2.5.2 Parameter Estimation-Based Techniques 18 3 Engine M o d e l i n g and Val ida t ion 21 3.1 Geometrical Consideration 22 3.2 Model Development 24 iv 3.3 Simulation Results 29 4 Parameter Identification 34 4.1 Identification Model 35 4.2 Performance Index 37 4.3 Simulation Results 39 5 Pressure Waveform Reconstruct ion 45 5.1 Gas Pressure Torque Estimation 46 5.2 Pressure Waveform Characteristics 49 5.3 Simulation Methods 51 5.3.1 Torque Estimation and Pressure Approximation 52 5.3.2 Fault Detection and Isolation Results 55 6 Summary, Discussion, and Conclusions 60 6.1 Modeling and Validation 61 6.2 Parameter Identification 61 6.3 Pressure Waveform Reconstruction 62 6.3.1 Individual Cylinder Pressure Reconstruction 63 6.4 Towards Condition Monitoring 64 6.5 Conclusions 64 Nomencla ture 66 Bib l iography 70 A Test C e l l Descr ip t ion 74 B P r o g r a m List ings 76 v Lis t of Tables 2.1 Specifications for the D D C 6V 92TA Detroit Diesel Engine 8 4.1 Estimates of the Engine Inertia 41 4.2 Estimates of the Torque Fluctuation 43 4.3 Residual to Signal Ratio 44 5.1 Pressure Variation Parameters 56 vi List of Figures 2.1 Ideal Diesel Cycle 8 2.2 Fault-Tolerant Control System 12 3.1 Crank Arm and Connecting Rod Mechanism 23 3.2 Piston Position, Speed and Acceleration 24 3.3 Connecting Rod Angular Acceleration 25 3.4 Torque Due to Gas Pressure 26 3.5 Torque Coefficients 27 3.6 Normal Operating Condition 30 3.7 Sine Wave Residual 32 3.8 Faulty Operation Due to Cylinder Under-Fueling 33 4.1 Output of the Identification Model 38 4.2 Velocity Fluctuation Estimation for the Normal Operation 42 4.3 Velocity Fluctuation Estimation for the Faulty Operation . . 42 5.1 Gas Pressure Torque Estimation for the Normal Operation 53 5.2 Gas Pressure Torque Estimation for the Faulty Operation 54 5.3 Typical Cylinder Pressure Waveform 55 5.4 Pressure Waveform Estimation for the Faulty Operation 57 5.5 New Torque Estimate for the Normal Operation 58 5.6 New Torque Estimate for the Faulty Operation 58 A . l Diesel Engine Test Cell 75 vii Acknowledgements The ^ opportunity to pursue a study of diesel engine rotational dynamics was provided by a research assistantship in the Department of Electrical Engineering, University of British Columbia, a short time after my arrival in Canada. This project started as a collaboration between the Robotics and Control Laboratory and the Institute for Ma-chinery Research, National Research Council of Canada. Given my initial background in control engineering, the transition to the new, for me, field of condition monitoring and fault diagnosis of internal combustion engines was challenging and sometimes full of dif-ficulties. Without the help and insight of my research supervisor, Dr. Peter D. Lawrence this accommodation might have been impossible to achieve. During our multiple discus-sions, concepts like fault detection and isolation, engine inverse dynamics, and parametric identification were linked to each other, thus providing me with the necessary tools to approach this project. In the spring of 1994, I had the opportunity to meet Dr. Phillip G. Hil l , who, through many helpful disscusions, put my work into the mechanical engineering perspective. He also offered me the possibility to develop the pressure measurement experiments using the test cell in the Engine Test Laboratory. In the same laboratory, I was welcomed and helped by K . Bruce Hodgins, and Peter Mtui. Dr. Terrence S. Brown and Devin Ostrom from the Institute for Machinery Research showed an active interest in this project. We discussed together many of the possible engine modeling approaches, and they helped me obtain the experimental data for model testing and validation. My final acknowledgement, to my husband, Matthew Palmer is of another nature. viii His patience in proofreading the manuscript, constant encouragement, and interest in my work, were invaluable and deeply appreciated. ix Chapter 1 Introduct ion dust and diesel rise like incense from the road— Bruce Cockburn, Dust and Diesel The reciprocating internal-combustion engine is one of the most common power pro-duction systems. The diesel engine plays an essential role in commercial automotive applications. Engine condition monitoring systems have the potential to improve reli-ability, lower maintenance costs, and improve safety by early detection of operational faults. Model-based condition monitoring permits the localization of faults within the engine and therefore aids in fault diagnosis. Our research is directed at developing and validating an effective diesel engine model for condition monitoring applications. In this project we analyze 6 different actuator faults and the normal operating condition. The actuator fault is defined as the phenomenon of cylinder under-fueling which has as a direct consequence, the reduction of cylinder peak pressure. We use both pressure and flywheel angular velocity measurements. Data correspond to a two-stroke diesel engine, D D C 6V 92TA, manufactured by the Detroit Diesel Corporation. The engine specifica-tions and much of the test data were furnished by the Institute for Machinery Research, National Research Council. Some measurements were done using the test cell in the Engine Test Laboratory in the Department of Mechanical Enginering at the University of British Columbia. The concept of analytical redundancy emphasizes the use of accurate dynamic and 1 Chapter 1. Introduction 2 static models for data processing and analysis. The major benefit is realized in the low cost and flexibility of a software versus a hardware implementation. A combination of analytical and physical redundancy is almost always necessary for a fault-tolerant system to maintain its function in the presence of certain failures. On the other hand, the paradigm of analytical redundancy has the advantage of fully exploiting the engine model and thus extracting information that otherwise might be difficult to obtain. This is the motivation for pursuing three directions of investigation: engine modeling and validation, parameter identification, and pressure waveform reconstruction. Thus, the second chapter of this thesis outlines the theory behind modeling tech-niques and estimation procedures to be applied in our study and relates them to the field of fault detection and isolation. The analytical model-based approach is compared to knowledge-based methods and to neural network-based techniques. Parameter identifi-cation techniques can be employed in the detection and isolation of incipient faults. The literature review provides the basis for our investigations by emphasizing the importance of preserving a physical relationship to the engine throughout the modeling and identifi-cation processes. In this context, we also present the basic terms and specifications of the D D C 6V 92TA diesel engine and the D D E C II electronic engine controller, manufactured by Detroit Diesel Corporation. The third chapter of this thesis aims to present the model for the dynamics of an TV-cylinder diesel engine, emphasizing its features and properties from a system theory point of view. The system inputs are given by cylinder combustion pressures, Pi(9), i = 1,... TV. Its output is represented by the flywheel angular velocity fluctuation, Sio(0). Both are expressed as functions of the crank angle, 0. The model is then verified using the M A T L A B simulation package. The model effectiveness is assessed by measuring the standard deviation of the residual (root mean square error) and the ratio between the standard deviation of the residual and the standard deviation of the measured signal Chapter 1. Introduction 3 (residual to signal ratio). In each of the 6 under-fueling faults, the velocity fluctuation residual is found to exhibit a strong sinusoidal trend which accounts for unmodeled torques. We attributed this harmonic torque fluctuation, A T , to the crankshaft vibration and the oscillatory behaviour of the dynamometer loading. The main goal of the fourth chapter is to describe a meaningful procedure for iden-tification of engine parameters. Identification methods attempt to determine values for the unknown or uncertain system parameters using measurements of input and output signals. This concept is transferred to the field of engine fault detection and isolation through the assumption that the occurrence of a fault has as a direct consequence, changes of process parameters, which modify the system output. Therefore, parameter estimation techniques can be employed to detect incipient process faults. If the residual generator is implemented as a parametric identification procedure, then the residuals are repre-sented by the output prediction errors. The decision making stage employs statistical comparisons between residuals and known fault signatures. Changes in parameter means and variances can be used as decision criteria. The identification problem is solved using the recursive gradient estimator and least-squares estimator with exponential forgetting. The results are presented in comparison with those obtained from the off-line standard least-squares estimator. The model parameters are functions of the crank-angle 9. The harmonic torque fluctuation, AT(6), and the total engine inertia, JE, are estimated. In the fifth chapter we analyze the model properties and the potential for solving the inverse problem of reconstructing the cylinder pressure waveform from noisy mea-surements of crankshaft/fly wheel angular velocity. Because our model implements a multi-input single output (MISO) system, it is not possible to explicitly decouple the inverse dynamics. We have therefore considered two approaches. First, we redefine the system input and we rewrite the engine model as a single input single output (SISO) system. If the input pressure waveforms did not overlap, then one could imagine that Chapter 1. Introduction 4 the system is SISO with sharp discontinuities at the end of each input period. The new input is then the torque due to gas pressure, Tp{0). Secondly, we approximate the pressure waveform with a periodic impulse-like continuous function. Using this latter template-based approach we are able to explicitly identify the condition of each cylinder and improve the torque estimation procedure. The sixth chapter summarizes the results and outlines the use of these three methods in a model-based engine condition monitoring system. Chapter 2 Background Again, nothing, and again the machine asked me politely: "Do you have the password?" Umberto Eco, Foucault's Pendulum Several methods of verifying the correct functioning of the diesel engine have been used by engineers from the early stages of the engine's development. The majority of those techniques relied on experience and close observation of some important engine operating factors such as: temperatures, pressures, flows, noise and vibration. Deviations from the normal operating characteristics were then recognized and classified. Advances in transducers technologies and, in particular, the availability of low-cost embedded computer hardware has rendered feasible, the application of much more so-phisticated condition monitoring schemes. In this chapter we introduce the definitions associated with fault diagnosis of systems and outline some methods applicable to diesel engines. The review of some recent engine condition monitoring applications puts the model-based approach into perspective by comparing it with knowledge-based methods and neural network-based techniques. In this context, we also present the basic terms and specifications of the D D C 6V 92TA diesel engine and the D D E C II electronic engine controller. We conclude the chapter by specifying the objective of our research in the general framework of model-based design for diesel engine fault diagnosis. 5 Chapter 2. Background 6 2.1 Engine M o d e l i n g The diesel engine is a complex device for which it is difficult to write a comprehensive mathematical model. For this reason, most models presented in the literature are ap-plication dependent. The laws of physics that govern a system represent an important modeling tool. Our modeling approach illustrates this aspect and is based on Euler's equations [1] applied to a generic internal combustion engine. These equations charac-terize the rotational dynamics of a diesel engine and can be subsequently used to derive the relationship between cylinder pressures and flywheel/crankshaft angular velocity fluc-tuation. The theory behind this approach is often applied to determine the equations of motions of any mechanical system [2]. Similar approaches dedicated to internal combus-tion engines can be found in references [3], [4] and [5]. The use of such a model in condition monitoring applications seems feasible as long as the correspondence between the model and the physical engine is transparent. We are interested in three main aspects. The first deals with validating the direct engine model, i.e. given the actual pressure inputs we are interested in measuring the error between the actual and the estimated crankshaft angular velocity. A second goal for our investigation aims at estimating some important engine parameters (inertia, and load torque) given input-output measurements. The third goal of this project is to present a method for pressure waveform reconstruction using the model and noisy measurements of the crankshaft angular velocity. This offers the possibility of replacing the direct measurement of cylinder pressures. The estimated pressures could then be used for a variety of calculations including heat-release analysis, net engine torque estimation, and cylinder power output calculation. Chapter 2. Background 7 2.2 Diesel Engine Characterist ics The design of a model-based fault diagnosis system takes into account the nature of faults to consider, the complexity of implementation, the robustness with respect to model inaccuracies, and the diagnosis performance index (number of false alarms). The model design is a critical step in the development of such a system. In this section we introduce some general diesel engine definitions and present some features of the D D C 6V 92TA diesel engine and D D E C II electronic controller. A diesel engine is a compression-ignition, reciprocating, internal-combustion engine. Four stroke diesel engines operate on a mechanical cycle [6, pp. 33-34] that has the idealized form shown in Figure 2.1. Each stroke is denned by the piston travel between top dead centre (TDC) and bottom dead centre (BDC). The induction stroke is characterized by constant pressure. The compression phase is followed by combustion at constant volume. The fuel is injected at the end of the compression stroke. Autoignition is made possible by a high compression ratio. After the expansion stroke the exhaust valve opens and the blow-down occurs at practically constant volume. The exhaust stroke takes place at constant pressure. Four-stroke diesel engines are used in many automotive applications. In a two-stroke engine the induction (entry of fresh air) and the exhaust (exit of burned gas) occur at the same time [6, pp. 274-323]. This phenomenon is called scavenging [7, pp. 213-215]. Some motorcycles, buses and locomotives are equipped with two-stroke diesel engines. The use of large two-stroke diesel engines is also common in marine applications. The advantages of two-stroke diesel engines are increased power output and high firing frequency. In our experiments we used the D D C 6V 92TA Detroit Diesel engine. Data were obtained from [3] and are illustrated in Table 2.1. This engine is turbocharged and uses in-cylinder fuel injection, therefore the scavenging is very efficient. Chapter 2. Background 8 induction/exhaust V r,Vd vol) luKe F i g u r e 2 . 1 : I l l u s t r a t e s a P - V d i a g r a m of a n i d e a l f our s t r o k e - D i e s e l c y c l e ( f r om [6, p p . 34]). T h e s t r okes are : i n d u c t i o n (0—1) , c o m p r e s s i o n ( 1 - 2 ) a n d c o m b u s t i o n (2 — 3) , e x p a n s i o n (3 — 4) a n d b l o w - d o w n (4 — 1), a n d e x h a u s t (1 — 0) . Vd is t h e v o l u m e at t o p d e a d cen t re , a n d rv is t h e c o m p r e s s i o n r a t i o . T a b l e 2.1: G e n e r a l features o f the t u r b o c h a r g e d D D C 6 V 9 2 T A D e t r o i t D i e s e l E n g i n e . Type: two stroke Number of cylinders: 6 ( V ) Bore: 123 mm Stroke: 127 mm Displacement: 9.05 liters Compression ratio: 17 : 1 Gross rated power output: 224 kW at 2100 rpm Friction power loss: 44.073 kW at 1800 rpm Overall mechanical efficiency: 84% Chapter 2. Background 9 2.2.1 D a t a Acqu i s i t ion The test cell consists of a fully instrumented two-stroke diesel engine, D D C 6V 92TA, manufactured by Detroit Diesel Corporation (see Table 2.1). The engine is electronically controlled by two Detroit Diesel electronic controllers (DDEC II). Data used in our simu-lation were furnished by the Institute for Machinery Research (IMR), National Research Council (NRC) and consist of cylinder pressures and flywheel angular velocity waveforms recorded at 1,200 rpm. The dynamometer torque was 990 Nm. There are 7 data sets, each corresponding to a certain engine condition. One set was obtained under normal operating conditions (baseline) and the other 6 were from operation with one cylinder at a time under-fueled by 10%. We also performed pressure measurements in the Engine Test Laboratory in the De-partment of Mechanical Engineering at U B C using a dedicated, completely instrumented test cell. This system was not equipped to measure the crankshaft speed fluctuations. The sensing system is very important for the accuracy of the condition monitoring sys-tem. Water-cooled piezoelectric pressure transducers were employed for measuring com-bustion pressure development [8]. For model validation purposes, both cylinder pressures and instantaneous crankshaft speed must be obtained simultaneously during the test. We measured the mean engine speed and the crank angle (CA) in the Engine Test Labo-ratory. A n update of the test cell is proposed in [9] and its description is detailed in Appendix A. 2.2.2 Engine Con t ro l A comprehensive description of the Detroit Diesel Electronic Control (DDEC) system was reported in [10]. The fuel volume flow and timing are controlled via a solenoid actuator. Each cylinder receives approximately equal quantities of fuel. The main control Chapter 2. Background 10 outputs are the beginning of injection (BOI) and the pulse-width (PW), i.e. duration of injection. Both depend on engine speed and torque. The top dead centre (TDC) corresponds to the minimum piston displacement, and the bottom dead centre (BDC) corresponds to the maximum piston displacement. The crank angle reference assumes 0 degrees at T D C and 180 degrees at B D C . The variables to be controlled are the engine speed and throttle position. Magnetic induction sensors are used to pick-up the time reference signal (TRSs), at 73.5 degrees B D C for each cylinder. The control action is not synchronized with changes in engine speed, and this can generate injection delays. The gain-scheduling (open-loop) solution is chosen for mid-range engine speeds. The engine controller functions in closed-loop for idle and rated speeds (when maximum output power is produced). The control time is 13 msec. The compressor boost pressure value is used to limit the allowable PW, and thus to compensate for smoke emissions caused by an inefficient scavenging process. A fast injection reduces the smoke production, while a torque increase has the opposite effect. Smoke control considerations require the acquisition of the following engine data: oil pressure and temperature, water flow and temperature, and air inlet temperature. Noise is caused mainly by the combustion process, and is reduced by delaying the BOI, which results in moving the combustion towards the end of the compression stroke. 2.3 Faul t Detect ion and Isolation in Con t ro l Systems Assume that the diesel engine can be described by a multi-input single output (MISO) model in the crank angle domain, 0. In this case, N defines the total number of cylinders, the model input is given by the N cylinder pressures, and the output is the crankshaft angular velocity fluctuation. We also assume that this dynamic system, V is characterized Chapter 2. Background 11 by the equations: x'(9) dx{9) A(p, 9)x(9) + B(P, 9)u{9), x(0) = x0, (2.1) y(9) = c(p,x) (2.2) where x,x' G Rn, are the state and its derivative, respectively, p £ Rq describes the process parameters, u G RN represents the input, and y £ R the output. The tracking, or servo problem, can be stated as follows: design a control law, uc(9) such that V XQ G X C Rn, a measure of the tracking error, e(9) — yd{9) — y{9) is minimized, while x is maintained stable (bounded) [11, pp. 192-197]. The initial state is XQ, and ya{9) is the desired output trajectory. The stabilization or regulation problem is equivalent to the tracking problem for con-stant yd{9). The control system refers, in general, to the following components: controller (C), process to be controlled (V), actuator (.4.) and transducer (T). A fault-tolerant control system is designed to take into consideration modeling inac-curacies, measurement noise and possible sensor or actuator faults [12]. It incorporates a fault diagnosis (FD) module which has the ability to detect, isolate and identify com-ponent malfunctions from observed symptoms. From a systems theory point of view, an abrupt fault is defined as a sudden jump in system response or parameters that do not necessarily correspond to a physical failure [13]. Biases or drifts are included in the category of incipient faults [12]. The process of observing (supervising) certain variables for diagnosis is called condition monitoring (CM). Figure 2.2 illustrates a block-diagram of a fault-tolerant control system. It is hierarchically structured, and the top level per-forms the tasks of condition monitoring, fault diagnosis and fault accommodation. A fault-tolerant control system has to respond rapidly when a failure occurs. This design requirement has as a consequence, an increased sensitivity to measurement noise. Fault detection and isolation (FDI) are important steps in a fault diagnosis procedure. The Chapter 2. Background 12 Figure 2.2: Illustrates a block-diagram of a supervised control system. The condition monitor-ing and fault diagnosis module allows the detection and isolation of component malfunctions. The controller is designed to exhibit robustness with respect to modeling inaccuracies and faulty sensors. first implies the acknowledgement .that a system component is defective (alarm) and the second locates the fault, i.e. differentiates between certain possible faults and determines the source of the failure [13]. The extent of the failure is evaluated as tolerable, condi-tionally tolerable or intolerable in the identification stage. System reorganization (fault accommodation) employs the substitution of the faulty component with a healthy one, if the fault is intolerable. Fault diagnosis and control could be carried out using a minimum set of measurements. In order to satisfy the requirements of high reliability and safety it is necessary to provide a degree of redundancy. The concept of physical (hardware) re-dundancy refers to the use of arrays (duplex, triplex, etc) of sensors for the measurement of the same variable. The concept of analytical redundancy refers to the use of process models and redundant functional relationships between plant variables of interest. The analytical (software) redundancy has the advantage of low cost and weight, flexibility Chapter 2. Background 13 and portability. Its reliability is strongly dependent on model accuracy. A combination of analytical and physical redundancy is almost always necessary for a fault-tolerant con-trol system to maintain its function in the presence of certain system failures. There are model-based, knowledge-based and neural network-based fault detection and isolation schemes. Model-based fault detection and isolation algorithms make use of measurable process inputs and outputs, nonmeasurable state variables, nonmeasurable process pa-rameters and nonmeasurable characteristic quantities [14]. The detection and isolation tasks are accomplished with the help of estimation algorithms developed from a priori knowledge of the system to be controlled. Thus, a central issue is the correct definition of the normal process. A l l these estimation techniques focus on the model accuracy. There are two types of model uncertainties: structured or parametric and unstruc-tured or unmodeled dynamics. An adaptive or a sliding-mode control scheme can com-pensate for the parametric uncertainties [11, pp. 276]. H°° optimization takes into ac-count a general perturbation model [15]. The principle of adaptive control (AC) is the continuous modification of the plant model and control law in response to parameter vari-ation. The parameters can be identified explicitly (self-tuning AC) or implicitly (model reference AC) . 2.4 Trends in Engine Cond i t ion M o n i t o r i n g The spread of electronic engine control modules (ECMs) in the 1980s, led mainly by the United States Clean Air Act, facilitated the incorporation of computers in cars and made possible the notion of up-integration (powertrain and vehicle control modules) and fault-diagnosis with the benefit of higher reliability [16]. In the U.S., the California Air Resource Board On-Board Diagnostic-II ( C A R B OBD-II) legislation had a powerful effect' on new developments for automotive digital electronics. Thus, "zero emission" vehicles, Chapter 2. Background 14 with sophisticated diagnosis tools and better human interfaces are expected by the end of this decade. The airline industry is a pioneer in the field and many ideas have now been transferred to portable commercial power machinery applications in fields such as naval, agriculture, mining, oil-well, road-building. The new trends are also encouraged by a more demanding consumer in the automotive market and by the increased degree of automation in manufacturing and process control. Most recent engine fault detection and isolation applications use: model, knowledge or neural network-based techniques [17]. This classification isn't exhaustive but repre-sents the methods with which we are most familiar. The model-based approach employs techniques for state and/or parameter identification. The state can be estimated using a Luenberger observer (the deterministic case) or a Kalman filter (the stochastic case). Parameter identification methods are suitable for early detection of incipient faults, while state estimation techniques are commonly used for detection of abrupt faults. The con-cepts of state and parameter estimation are interchangeable to a certain extent. In the field of fault detection and isolation the estimation errors define the process of residual generation. The phase of decision making is accomplished by statistical threshold tests. A review of the state of the art of model-based fault diagnosis methods for jet engine systems can be found in [18, 19]. The authors emphasize engine sensor failures as the most critical problem to deal with. In the aerospace industry, there is a lot of interest in replacing the hydro-mechanical controllers with fault-tolerant electronic ones for aerojet engine systems. A n example is the so-called full authority digital electronic controller (FADEC) for gas turbines. It controls the fuel flow rate and the exhaust nozzle area. The measured variables include: pressures, flows, temperatures and compressor shaft speed. A model-based fault diagnosis scheme for this complex non-linear system has been tested at NASA. It uses a Kalman filter for residual generation and a generalized likelihood ratio for decision making. It is therefore capable of detecting abrupt faults, but not dealing Chapter 2. Background 15 with the issue of incipient faults. Research in the field of engine fault detection and isola-tion is currently connected under the auspices of an international Technical Cooperation Panel with the participation of the Commonwealth countries and the United States [18]. In [20], the authors focus on a robust solution to the fault detection and isolation problem using the "disturbance decoupling" principle. The residual generator is a Lu-enberger state observer and the weighted output estimation error represents the fault indicator. This approach doesn't take into consideration the presence of measurement noise, and assumes a measurable system state. A combination of model and knowledge-based methods is used in [21] for the devel-opment of a condition monitoring system for a Spey SMI A marine gas turbine engine. There are three components that define the overall structure: the signal processing sub-system (SPS) which collects and processes the raw sensor data, the signal identification subsystem (SIS) which makes use of linear engine models and SPS information to calcu-late engine parameters and performance, and the monitoring and diagnosis system (MDS) which gives a symbolic interpretation of the SIS output. This approach doesn't exploit the possibility of using dynamic models and their ability to predict system behaviour. It focuses mainly on the issue of explicit fault classification by means of rule-based tech-niques. The subject of engine sensor faults or failures in automotive applications is treated in [22]. The authors introduce a model-based fault detection and isolation method for internal combustion engines. Residuals are generated using the eigenstructure assignment method, which is equivalent to the unknown input observer technique. These detection niters are used to localize failures in the throttle position and in the manifold absolute pressure sensors. A n attempt at diagnosis of the supercharger of a diesel engine is reported in [23]. The expert-system is tested on board a French Navy ship. There are four variables to be Chapter 2. Background 16 monitored: atmospheric pressure, temperature, accelerator position and engine speed. A "qualitative state" is associated with each of these parameters. The observed and the expected qualitative states are compared and then a decision is made using forward chain inference rules. A neural network-based method is used in [24] to develop the classification of me-chanical faults in newly manufactured engines. Data acquisition includes the following measurements: intake and exhaust manifold pressures, the crankcase air pressure and the oil pressure. These five waveforms are monitored by a neural network-based super-visor. The 29 classes consist of 28 different faults and the normal operating condition. The training is performed by the error back-propagation algorithm [25]. The informa-tion concerning the different faults is embedded in the network weights, but there is no explicit way of translating these values into physically meaningful parameters. Our approach to the problem of engine fault diagnosis is based on the concept of parameter estimation. We are interested in early detection of faults that manifest as parameter biases or drifts. The engine model is developed as a continuous-time dynamic system, with parameters that have a physical relevance. 2.5 Res idua l Generat ion and Decision M a k i n g A model-based fault detection and isolation procedure consists of two phases: residual generation (RG) and decision making (DM) [12]. First, the effect of the fault is am-plified in order to make it recognizable (determine a fault signature). This is obtained by processing the sensor signals. The processed measurements are called residuals or fault indicators [26]. If the residual generator is designed as a parametric identifica-tion procedure, then the residuals are represented by the output prediction errors. The Chapter 2. Background 17 decision can be accomplished by a simple threshold test or by statistical methods. Con-sider the system described by (2.1)—(2.2) and assume that u(9) and y(9) are measurable, zero-mean signals. We also assume that a fault has as a direct consequence, changes of process parameters, p and/or state, x(9), which modify the system output y(9) [17]. Therefore, either state or parameter estimation techniques respectively can be used to detect incipient or abrupt process faults. 2.5.1 State Estimation-Based Techniques State estimation-based residual generators include: the parity space method, dedicated observer approach, and fault detection filter approach [12]. The problem of discerning between different faults or between faults and other disturbances can be solved using the disturbance decoupling principle. If we take into account the effect of perturbations, (2.1)-(2.2) can be rewritten: x'{9) = A(p, 9)x(8) + B(p, 9)u{9) + d{9), x(0) = x0j (2.3) y(9) = c(p,x)+n(9), (2.4) where d G Rn and n G R describe the process and measurement noise, respectively. Both are considered unknown. In the stochastic case, d, and n are assumed zero-mean, indepen-dent white Gaussian processes. The estimator of the system normal state (see (2.3)-(2.4)) has the following form: x'{9) = A(p, 9)x{9) + B(p, 9)u{9) + H(9)(y - y), x(0) = x0, (2.5) y(9) = c(p,x), (2.6) where H(9) G Rn is the feedback gain. The state and the output estimation errors are given by: x' = A(p,9)x + d(9)- H(9)y, Chapter 2. Background 18 y = y -y = c(p, x) - c(p, x) + n{9). The condition monitoring system supervises the evolution of the residual (innovation) y. The fault detection can be the result of certain statistical tests performed on the residual: chi-squared test, or generalized likelihood ratio test [13]. Another approach consists of designing the gain H{6) with the objective of amplifying the effect of certain failures in y rather than providing good state tracking. Such observers are called Beard-Jones failure sensitive filters. Assume d{9) = dfa with 6 ; the i-th column of b, i.e. the fault corresponds to a bias in the i-th actuator. The expression for the residual becomes y' = (A(p, 9) — H(9))y + dibi, if y = x. For H(6) = A(p, 9) + a0In, y has the orientation of bi and a magnitude proportional to dj. The gain design procedure is similar to the disturbance decoupling problem [13]. 2.5.2 Parameter Est imat ion-Based Techniques Identification methods attempt to find suitable approximations to the parameters asso-ciated with real systems. Let p be an estimate of the process parameter vector p. The output prediction is y = c(p,x) and the prediction error is y = y — y — y — c(p,x). The parameter estimation problem can be formulated as the minimization of a generalized prediction error function, E [27]. The necessary zero gradient condition for the minimum is: VpE\p=p = 0. (2.7) From (2.7).the parameter estimates, p, have to be determined. An analytical solution is tractable if c is linear in p, i.e y(9) = c(p,x) = $(9)p with $ e Rlxq [27]. In this case y = <&(9)(p — p) = $(9)p, where p is the parameter estimation error. The gradient estimators and the least-squares estimators with exponential forgetting are characterized Chapter 2. Background 19 by good robustness with respect to noise and parameter variation [11, pp. 367-369, 380-381]. The simplest parameter estimator is the gradient estimator. In this case, the parameter adaptation law is given by: ^ = - * 0 V p £ = - * o $ T y , (2.8) where E = | || y | | 2 is the squared output prediction error, and \T/o € Rsxs > 0 is the gain matrix. The signal matrix, $ has to be persistently exciting in order to achieve exponential parameter convergence, i.e. 3ao, © > 0 such that V 0 > 0 / > Q 0 / „ (2.9) J 0 with Iq E Rqxq the idendity matrix [11, pp. 366]. The cost function used in least-squares estimation with exponential forgetting is: E = \ [8e-I>^dr || y(t) - $(t)p ||2 dt, (2.10) 2 Jo where A represents the forgetting factor. The parameters are updated using: fQ = -*mTy, (2.11) where tf"1^) = tf-^Oje-Zo A « d t + [* e-St8x{-r)dr$T(t)$(t)dt. Jo The gain matrix, \&(0) is calculated recursively using the formula: d\Sj — = - A ( 0 ) * + tf$T(0)$(0)* . dO The decision making stage employs statistical comparisons between actual residuals and known fault signatures. Changes in parameter means and variances can be used as decision criteria. Chapter 2. Background 20 Our research is directed at developing and validating an effective diesel engine model for fault diagnosis applications. The model is used in a parametric identification proce-dure to detect incipient faults in the fueling of engine cylinders [14]. The early detection and isolation of engine faults is accomplished using the gradient estimator and the least-squares estimator with exponential forgetting. The choice of these two identification methods is based on their robustness with respect to measurement noise and parameter variation [11, pp. 367,380]. The same model is used to obtain the cylinder pressure wave-form, reconstructed from noisy measurements of flywheel/crankshaft angular velocity fluctuations. Chapter 3 Engine Modeling and Validation No more friction, no more slowing. The length of the Earth day will then be more than fifty times as long as the present day; and the more distant Moon will turn in its orbit in twice the period it now turns. Isaac Asimov, Asimov on Astronomy The first step in designing a comprehensive diesel engine fault detection and isolation system is to develop and validate a mathematical model. The model we develop here is configured to correspond to the D D C 6V 92TA Detroit Diesel engine. The analytical methodology follows that of C. F. Taylor [28, pp. 240-305], and the testing was performed using data supplied by the Institute for Machinery Research, National Research Council. Our modeling approach is based on Euler's equation [1] applied to a generic internal combustion engine. These equations of motion characterize the rotational dynamics of a diesel engine and can be subsequently used to derive the relationship between cylinder pressures and flywheel/crankshaft angular velocity fluctuation. We define the total engine torque, TE(9), as the sum of all torques acting on the crankshaft taking into consideration the phase shift determined by the cylinder arrangements: TE(6) = Tp(9) + Tt{9) + Tr{9) + Tf(9) + T,(0), (3.1) where 9 is the crank angle. Tp(9) is the indicated torque due to the gas pressure forces, 21 Chapter 3. Engine Modeling and Validation 22 Tt{0) is the torque due to the inertia of the reciprocating parts, and Tr(9) is the torque due to the rotation of the connecting rods. Ty(0) refers to the total friction torque (assumed constant), and Ti(9) is the dynamometer (brake) load torque. The model is developed under the assumption that the crankshaft is a purely rotating rigid body. Similar approaches dedicated to internal combustion engines can be found in references [4] and [5]. The novelty of our model relies on the fact that it incorporates the torque due to angular acceleration of the connecting rods. We compare the measured angular velocity with the calculated one. The testing set consists of 6 different actuator faults and the normal operating condition. The actuator fault is defined here as the phenomenon of cylinder under-fueling which has as a direct consequence the reduction of cylinder peak pressure. We use both pressure and flywheel angular velocity measurements. 3.1 Geomet r ica l Considerat ion The engine crankshaft is assumed perfectly rigid. The connecting rod is modeled by two masses: mt which is included in the piston assembly and mr which rotates with the crank pin. Consider the diagram for the crank and connecting rod mechanism illustrated in Figure 3.1. The piston instantaneous position with respect to the top dead centre (TDC) is: where r is the crank radius, I is the connecting rod length, 9 is the crank angle, and a is the angle between the connecting rod axis and the cylinder axis. The piston instantaneous velocity is calculated by taking the time derivative of (3.2): s = I + r — (r cos 9 + I cos a) (3-2) with (3-3) s — ruoci{9) (3-4) Chapter 3. Engine Modeling and Validation 23 Figure 3.1: Crank arm and connecting rod mechanism: r is the crank radius, I is the connecting rod length, 9 is the crank angle, and o is the angle between the connecting rod axis and the cylinder axis. The piston instantaneous position with respect to the top dead centre is s — I + r — (r cos 6 + I cos a) . where u = 9. The geometrical coeffcient c\(9) has the following expression: . . / r cos 9 \ Cl{9)= 1 + - sin0, \ I cos a / (3.5) Similarly, the piston acceleration is given by: a = r (ci(0)w + c2(9)u2) (3.6) The coefficient C2(9) is given by: dcx(9) d,9 1 + r cos 9 I cos a cos 9 + r s in r cos I I cos a \ I cos a (3.7) Figure 3.2 illustrates the variation of piston position (see (3.2)), speed (see (3.4)), and acceleration (see (3.6)) during one revolution, assuming a nominal rotational velocity of 1, 200 rpm. The angular acceleration of the connecting rod is calculated from Figure 3.1 \ Chapter 3. Engine Modeling and Validation 24 -0.06 Figure 3.2: Instantaneous piston position, s (left), speed, s (middle) and acceleration, s (right) during one revolution, assuming a nominal crankshaft angular velocity of 1, 200 rpm. as a function of us and us2: ix = c3(0)d) + c4(0)u;2 (3.8) The coefficients c3{9) and c^{9) have the following form: <*(*) = dc3(8) r cos 9 I cos a' r sin 9 d,9 I cos a v ' (3-9) (3.10) Figure 3.3 illustrates the angular acceleration of the connecting rod (see (3.8)) during one revolution, assuming a nominal rotational velocity of 1, 200 rpm. 3.2 M o d e l Development A modified engine model is developed in this thesis by finding explicit formulas for the torques involved in (3.1). Tp(9) is defined under the assumption that the work done on Chapter 3. Engine Modeling and Validation 25 T 1 1 r degrees Figure 3.3: Angular acceleration of the connecting rod, a, during one revolution, assuming a constant crankshaft angular velocity of 1, 200 rpm. the piston is equal to the work done on the crankshaft [28, pp. 268]. It depends directly on each cylinder pressure: N Tp(9) = AprY,Pi{8)ci(e - , (3.11) i=i where ./V is the total number of cylinders, and Ap is the piston crown area. Pi(0) refers to the pressure in the i-th cylinder and 0,; is the phase shift corresponding to the i-ih firing cylinder. Figure 3.4 illustrates the measured torque due to gas pressure for the D D C 6V 92TA operating at 1,200 rpm. The phase angles 0?; are the angles of the crank arm relative to the crank angle when cylinder 1 is at top dead centre. For a two-stroke, even-firing, in-line, six-cylinder engine, 0, = | ( i — 1), i = 1... 6 rad. In the case of the D D C 6V 92TA Detroit Diesel engine, the V-angle between cylinder banks is 63.5 degrees. The phase origin is the top dead centre of the first cylinder firing, and 4>\ — 0, 02 = 56.5, 03 = 120, 0 4 = 176.5, 05 = 240, 06 = 296.5 degrees. The work done by the crankshaft on the connecting rod is equal to the change in Chapter 3. Engine Modeling and Validation 26 1 8 0 0 i • i i r 1 6 0 0 \ T p _ m = 8 5 6 . 5 N m / \ 1 4 0 0 -g- 1 2 0 0 § 1 ooo c r f 8 0 0 | 6 0 0 4 0 0 2 0 0 O - 2 0 0 O 5 0 1 0 0 1 5 0 2 0 0 d e g r e e s 2 5 0 3 0 0 3 5 0 Figure 3.4: Torque due to gas pressure for the case of the six-cylinder DDC 6V 92TA Detroit Diesel Engine operating at 1,200 rpm. kinetic energy of the piston [28, pp. 264]. Tt(0) is due to the inertia of the translating parts and depends on all piston accelerations (see (3.6) and Figure 3.2): N Ti{9) = -myY,ci(0-^)(ci(d-^ + c^9-^)u)2)^ (3-12) ?;=i where mp is the mass which is considered to reciprocate with the piston, i.e. piston, piston rings, piston pin and the upper end of the connecting rod. The torque acting on the crankshaft due to the angular acceleration of all connecting rods is (see (3.8) and Figure 3.3): N Tr{e) = J c r V ^ C 3 ^ - ^ ) ( C 3 ^ - ^ ) w + c 4 ( f 5 - ^ ) ^ 2 , (3-13) i=l where J C T is the connecting rod moment of inertia. We define the following variable parameters: N cn(0) = £ c ? ( 0 - , (3-14) ?;=i Chapter 3. Engine Modeling and Validation 27 3 . 1 1 0 O 2 0 0 d e g r e e s 1 0 O 2 0 0 d e g r e e s 1 0 0 2 0 0 d e g r e e s F i g u r e 3.5: Var ia t ion of the torque coefficients: c\i(6), 012(0), £ 3 3 ( 0 ) , and 034(6) w i t h j = 0.247, dur ing one crankshaft revolution. TV N C3 3(0) = £ C 3 ( 0 - & ) i=l (3.15) (3.16) (3-17) i=l These parameters are illustrated in the graphs of Figure 3.5. We can rewrite (3.12) and (3.13) as: Tt(9) = - m y 2 (cn(^)w + c12(e)u2) , (3.18) T r(0) - Jcr (c33(0)u + cu(6)u2) . (3.19) Let JE be the engine moment of inertia. The Euler equation for the diesel engine, viewed as a rotating rigid body is: d(JEio) dt TE(0), (3.20) Chapter 3. Engine Modeling and Validation 28 or taking into consideration (3.1), (3.11), (3.18), (3.19): (jE ~ JcrC33(0) + m p r 2 c n ( 0 ) ) u = (j c rc 3 4(0) - mpr2c12{e)) u2 + TF{9) + TT(9) + N V E W - & ) c i (3-21) i=l Let x = UJ2, then, dw dw dd f 1 da; 1 , u — — : = to LO — — —x . (3.22) dt dO dt 2 d£ 2 v ' Therefore, (3.21) becomes: \ (JE ~ JcrC33(e) + m p r 2 C l l ( 0 ) ) rr' = ( j c r c 3 4 (0) - mpr2c12(9J) x + TF{6) + TT(9) + N AprY.PiifyiV-fc). (3.23) i=l Let = describe the engine operating point, with UQ the nominal engine speed. The engine speed fluctuation is Su = u — UQ, and this corresponds to a fluctuation of the variable x of the form: Sx = Su(2uQ + Su) « 2UQSU . (3.24) At the same time, the speed fluctuation can be determined from the variation Sx using: Su = yjsx + u2 — UQ . (3.25) With these considerations, the state equation, given by (3.23), is rewritten in a form similar to (2.1): Sx' = A(e)Sx + bT(e)u(e) + d(e), (3.26) and (3.25) becomes the output equation similar to (2.2): Su = c(Sx), (3.27) where A (o\ - J c r c 3 4 ( g ) ~ mpr2c12(9) Ay&) - TT~r r 77VT 2 7oY\ ' (3.28) 2 W£ - ^ 3 3 ( 0 ) + rripr^cu^e)) Chapter 3. Engine Modeling and Validation 29 b{9) Apr \{JE~ Jcrcs3(0)+ m,pr2cu(8)) (J c r c 3 4 (0) - m,pr2c12(8)) x0 + Tf(9) + Ti(0) ci{e-i) ' Pi(e) ' , u(9) = _ c\{9 - 4>N) _ _ PN(8) _ d(9) \ (JE ~ JcrC33{9) + mpr2cn(9)) and c(6x) = \JSx + UQ — U>Q (3.29) (3.30) (3.31) 3.3 S imula t ion Results The model described by (3.26) is validated using the M A T L A B simulation package, and the code is listed in Appendix B. The most important M-files are: • C-l.m, C-2.m, cS.m, C-4-m which implement the coefficients described by (3.5), (3.7), (3.9), (3.10), respectively, • system.m which implements the differential equation (3.26), • demol.m which calculates the residual and the estimation for a certain cylinder condition (faulty or not), and • rmse.m which calculates the root mean square error (see (3.33)). The model parameters, from [3] are: N = 6, r = 0.0635 m, I = 0.2571 m, J C R = 0.0745 kgm 2 , mp = 6.03 kg, Ap = 0.01188 m 2 . The engine effective moment of inertia includes the crankshaft and the flywheel moments of inertia and has the lumped value JE = 4.044 kgm 2 . The phase origin is the top dead centre of the first cylinder to fire, and i — 0, 02 = 56.5, 03 = 120, 04 = 176.5, 0 5 = 240, 0 6 = 296.5 degrees. The firing order is: IL, 3R, 3L, 2R, 2L, 1R, where the letters "L" and "R" refer to the left and the right bank of cylinders, respectively. Data acquisition is performed for a nominal speed UJQ = 125.7 Chapter 3. Engine Modeling and Validation 30 degrees Figure 3.6: Presents the actual (solid line), 6u>, and the estimated (dashed line), <5u), flywheel angular velocity fluctuation for the normal operation of the diesel engine. rad/sec (1, 200 rpm) and a dynamometer torque, T/ = 990 Nm. This load represents 75% of the engine peak torque. Tf = 156 Nm is the value of the mechanical friction torque corresponding to a crankshaft angular velocity of coo = 125.7 rad/sec and an overall damping coefficient of 1.24 Nmsec. The testing set consists of 6 different actuator faults and the normal operating con-dition. The actuator fault is represented by a 10% drop in cylinder fueling. In order to measure the modeling error we calculate the output residual: 8u = Su — Su , (3.32) where 6u> and SUJ are the actual and the estimated flywheel angular velocity fluctuation, respectively. Figure 3.6 illustrates the measured and estimated flywheel angular velocity fluctuation for the normal operating condition (baseline). The root mean square (RMS) error is defined as the standard deviation of the zero-mean residual: Chapter 3. Engine Modeling and Validation 31 RMS velocity error = \ ^ ^ 0 ) ) dB . (3.33) V Z7T The ratio between the RMS error and the standard deviation of the measured fluctuation is the residual to signal ratio (RSR) and is used here to express the effectiveness of the model (3.26)-(3.27): R g R = RMS error i o Q % ) a(ooj) where a(8u>) is the standard deviation of 6u>. The first step of the validation procedure is concerned with the analysis of the model behaviour in the case of engine normal operation. In the case of this experiment (see Fig-ure 3.6), the RMS error is 0.0891 rad/sec. The standard deviation of the measured signal is 0.1956 rad/sec. This result yields an RSR value of 46%. The model is then tested for the case of an under-fueling fault cylinder 1L. The difference between the actual and the estimated velocity waveforms for this situation is plotted in Figure 3.7. In each of the 6 under-fueling faults (10% down condition for one cylinder at a time), this residual is found to exhibit a strong sinusoidal trend. A number of factors can explain this harmonic trend. The first is crankshaft vibration [29, pp. 64-74] which is not taken into consideration in our model since it assumes a perfectly rigid crankshaft. The ampli-tude of the vibration increases in the presence of cylinder under-fueling. Secondly, the dynamometer and the dynamometer controller form a feedback system which attempts to maintain a constant engine speed. This can induce an oscillatory behaviour in the dynamometer loading. In order to eliminate this trend, the following technique is employed: 6u>(9) = 6ou(B) - 6u>(9) = Qhsin(0 + 7) + n, (3.35) Chapter 3. Engine Modeling and Validation 32 2.5 O 50 100 150 200 250 300 3 5 0 degrees Figure 3.7: The difference between the estimated and measured velocity fluctuation, 6ui, in the case of 10% under-fueling in cylinder IL. with 6u> the estimation error, fi/j the unknown amplitude, 7 the unknown phase, and n the measurement noise. We added the data to a version of itself shifted by 180 degrees. Since the estimation and the actual SLO are periodic with period 360 degrees, then the residual, Su), has the same property, and the shift becomes rotation: 6u>(9) + 6u(0 + TT) » 0. (3.36) The resulting estimate is shown in Figure 3.8 The RMS error reduces to 0.1368 rad/sec. This corresponds to a 38% RSR. The RMS error for the six under-fueled conditions lies between 0.0973 and 0.1836 rad/sec. This corresponds to a RSR of 26% to 54%. The best results are obtained for an under-fueling fault in cylinder 1R, and the worst case corresponds to an under-fueling fault in cylinder 2L. Chapter 3. Engine Modeling and Validation 33 degrees Figure 3.8: The actual (solid line) and the estimated (dashed line) flywheel angular velocity fluctuation for the faulty operation of the first cylinder on the left bank of the engine. These results are posterior to the elimination of a vibrational sine trend as per (3.36). The results reported here can be improved via a parameter identification procedure. The uncertain system parameters are engine inertia, JE and the unknown torque fluctu-ation, A T = Th sin(# + 7), where T/,, is the torque amplitude, and 7 the phase. This is the main objective of the following chapter. C h a p t e r 4 P a r a m e t e r I d e n t i f i c a t i o n The one that has the most will be the greatest? Now I understand, Sir, you are equating quality with quantity. Eugene Ionesco, The Lesson Identification methods attempt to determine values for the unknown or uncertain system parameters using measurements of input and output signals. This concept is transferred to the field of engine fault detection and isolation by assuming that the occurrence of a fault has as a direct consequence, changes of process parameters, which modify the system output. Therefore, parameter estimation techniques can be employed to detect incipient process faults. If the residual generator is implemented as a parametric identifi-cation procedure, then the residuals are represented by the output prediction errors. The decision making stage employs statistical comparisons between residuals and known fault signatures. Changes in parameter means and variances can be used as decision criteria. The parameter identification procedure can be accomplished on-line or off-line. In this chapter we are focusing on two on-line identification methods: the gradient estimator and the least-squares estimator with exponential forgetting. The results are compared with those obtained from the off-line standard least-squares technique. The parameters of interest are the engine inertia, JE, and the torque fluctuation, AT(9), which we attributed to unmodeled crankshaft vibration and oscillatory loading of the dynamometer. The reason for choosing these two estimators is found in their robustness with respect to measurement noise, and parameter variation [11, pp. 367, 380]. The analysis of the 34 Chapter 4. Parameter Identification 35 estimated parameters allows the classification of an engine cylinder as faulty or not. 4.1 Identification M o d e l There are a number of papers that attempt internal combustion engine parametric iden-tification using least squares approaches. Many of them have as a final goal, the adaptive control of systems that have a diesel engine as the principal component [30, 31, 32]. Three such applications—adaptive speed control of a stationary diesel engine for power generation, self-adaptive idle speed control of an automotive diesel engine and adaptive performance optimization—are discussed in reference [30]. In that work the author uses the recursive least squares (RLS) method for parametric identification with engine models in the autoregressive (AR) form: where A(q~x) and B(q~1) are polynomials in the shift operator q~l, u is the fuel rack position, and u> is the engine speed. The RLS identification of a diesel prime-mover with unknown dead-time is treated in references [31, 32]. Here the authors accelerate the algorithm convergence by imposing limits on model parameters. The model used in [33] for gain scheduling control is a fifth order A R type (see (4.1)) and describes the dynamics relating the fuel rack to engine speed. The model parameters (the coefficients of polynomials A(q~x) and B(q"1)) are nonlinear functions of engine speed and power output. A l l models mentioned above are suitable for adaptive control applications, but they all have the disadvantage that they don't preserve any physical relationship to the engine throughout the identification process and therefore they are not ideally suited to condition monitoring applications. (4.1) Chapter 4. Parameter Identification 36 A parametric identification approach to engine fault detection and isolation is reported in [4]. The authors assume that the influences of cylinder faults on the flywheel angular velocity can be decoupled from each other. This is a simplification that may lead to misclassification as faulty, a normal cylinder adjacent to the faulty cylinder. This problem is solved in [34], but the solution involves pattern recognition techniques. A combination of parameter estimation methods and expert-system interpretation techniques is used in [17]. This approach is applied to 4 case study experiments in [35]. The authors recommend the use of an analytical knowledge paradigm for well understood processes, followed by a heuristic analysis for aspects that cannot be fully explained by mathematical means. Our identification model is obtained from the engine model developed in the previous chapter (see (3.26)-(3.27)): 8x = A{B)8x + bT(B)u(8) + d(9), 8u> = c(8x), where A(9), b(9), d(B), c(8x) are the system coefficients defined in (3.28)-(3.31), the variable 8x is defined in (3.24), and 6u> is the angular velocity fluctuation. The sinu-soidal trend of the residual, depicted in Figure 3.7 which we attributed to the crankshaft vibration and oscillatory behaviour of the dynamometer loading is: AT(B) = Thsm(8 + 1), (4.2) where 7),, is the unknown torque amplitude, and 7 the unknown phase. In order to obtain a linear parametrization model we rewrite (4.2) as: AT (9) = Ts sine? + TC cos B, (4.3) where Ts = T/,, cos 7, and Tc = T/,sin7 are parameters to be identified, and B is the crank angle. Taking into consideration (4.3), we can represent the system described by Chapter 4. Parameter Identification 37 (3.26)-(3.27) in the linear parametrization form: Y(9) = $(9)p, (4.4) where $(#) is the signal matrix: -^6x'(9) sinB cos6 p is the vector of unknown parameters: (4.5) p = [JE TS Tcf , (4.6) and Y(9) the modified torque output for the identification model: Y(9) = -mpr2 \\Cll(0)6x'(9) + c12(9) (6x + x0) + Jc \czz(9)6x'(9) + cM{B) (Sx + xQ)] + Tp{9) + Tf{6) + T,(0), (4.7) with mp the mass of the reciprocating parts, r the crankshaft radius, Jcr the connecting rod mass moment of inertia, Tp(0) the torque due to gas pressure forces (see (3.11)), Tj(/9) the nominal load torque, and Tf(8) the friction torque. The torque coefficients, c n(#), ci2(8), £33(0), 034(8) were defined in (3.14)-(3.17), and XQ = u;2, with UJQ the nominal engine speed. Y(6) is a zero-mean variable and is plotted in Figure 4.1. Assume that the number of available measurements is equal to NM > 3. Then, the objective of the identification procedure is to find the solution of the algebraic system of NM equations of the form (4.4) with 3 unknowns: JE, TS, and Tc. 4.2 Performance Index Assume that p is a solution of the following algebraic system: Y(9k) = $(8k)p, k = l,.:.NN (4.8) Chapter 4. Parameter Identification 38 1 0 0 0 degrees Figure 4.1: The variation of the corrected engine torque, Y (6), during one crankshaft revolution for the normal operating condition. where 0 k defines the fc-th measurement (sampling instant). The output prediction is Y{0) = $(9)p and the output prediction error is Y{9) = $(0)p, with p = p — p. The parameter estimate, p is found by minimizing a measure of the output prediction error, Y{9). For standard least-squares estimation the cost function is given by: E = \f\\ Y(f<) ~ * (* )P II' d t (4-9) I Jo Geometrically, this is equivalent to determining the projection of the unknown parameter vector p on the hyperplane defined by the columns of $(9). The off-line solution to the standard least-squares estimation problem is: p = $+(O)Y(0), (4.10) where $+(6>) is the pseudoinverse of <&(9). The pseudoinverse exists if the matrix QT(9)$(9) Chapter 4. Parameter Identification 39 is nonsingular: $+(f?) = ($ T (0 )$ (0 ) ) - 1 $ T (0 ) . (4.11) The cost function used for the least-squares estimator with exponential forgetting is (see (2.10)): E = - /V^ A ( r ) d r || Y(t) - $(t)p ||2 dt, (4.12) 2 Jo where A represents the forgetting factor. The on-line implementation of the least-squares estimator with exponential forgetting is described by the following equations: Y9 = -m)*T(8i)Y(0t), (4.13) W) = A 0 ( l - " 1 1 ) , (4.14) ^ = -X(9i)^(9i) + ^(9^(9^(9^(8,), (4.15) d,9 where ^(9,) represents the estimator gain at the i-th iteration, || \P(0j) || is the norm denned as the maximum singular value of the matrix ^ (#,;), Ao is the maximum allowable value for the forgetting factor, and k0 > 0 [11, pp. 374-377]. The performance index used by the gradient estimator is the squared output predic-tion error: E = ^ || Y | | 2 . (4.16) The parameters are estimated recursively following the law: d™ „ ~ = -Vo&MYiOi), (4.17) d,9 where v&o > 0 is the descending step (estimator gain) [11, pp. 364-367] 4.3 S imula t ion Results We analyzed the performance of the on-line gradient estimator and least-squares with exponential forgetting in comparison with the off-line standard least-squares estimator. Chapter 4. Parameter Identification 40 The simulation is performed using the M A T L A B package and the code is listed in Ap-pendix B. The main functions are: • sls.m which implements the standard least-squares estimator, • itgr.m which implements the on-line gradient estimator, • itls.m which implements the recursive least-squares estimator with exponential for-getting, and • est.m which calculates the velocity fluctuation using the parameter estimates. The angular velocity is sampled at a period of ir/59 rad. This sampling period corre-sponds to the total number of flywheel teeth, 118. The cylinder pressure data is char-acterized by a sampling period of 7r/360 rad. The model parameters are given the same numerical values as in Section 3.3. These are: TV = 6, r = 0.0635 m, I = 0.2571 m, Jcr = 0.0745 kgm 2 , mp = 6.03 kg, Ap = 0.01188 m 2 . The phase origin is the top dead center of the first cylinder to fire, and 0 i = 0, 4>2 = 56.5, 3 = 120, 4>4 = 176.5, Q = 125.7 rad/sec and an overall damping coefficient of 1.24 Nmsec. The parameters to be identified are the engine effective moment of inertia, JE, which includes the crankshaft and the flywheel moments of inertia, and the torque fluctuation, A T . The testing set consists of 6 different actuator faults and the normal operating condition. The actuator fault consists of a 10% drop in cylinder fueling. The estimation is performed for Nm = 118 measurements. The effect of the parameter identification procedure is measured using the same criteria as in Section 3.3: the root mean square Chapter 4. Parameter Identification 41 Table 4.1: Engine inertia estimates, JE, depending on the number of measurements, Nm = 118, the engine operating condition (normal or faulty), the identification method: off-line standard least-squares (SLS), on-line gradient estimator (GE), and least-squares with exponential for-getting (LSEF). Faulty Cylinder JE (kgm2) SLS GE LSEF None 3.86 3.90 3.85 1 4.89 4.90 4.86 2 5.06 5.06 4.97 3 5.17 5.13 5.08 4 5.01 5.04 4.97 5 5.01 5.01 4.97 6 4.96 4.96 4.91 (RMS) error, defined as the standard deviation of the velocity residual and the residual to signal ratio (RSR), defined by the ratio between the standard deviation of the residual and the standard deviation of the actual angular velocity fluctuation. The different engine inertia estimates, J E , are illustrated in Table 4.1. Estimates for the other two parameters, Ts and Tc, characterize the harmonic torque fluctuation (see (4.3)). Table 4.2 presents the values Th and 7 (see (4.2)) for different engine conditions. Figures 4.2 and 4.3 present the estimation results for flywheel angular velocity fluc-tuation corresponding to two test cases—normal operating condition, and under-fueling fault in cylinder 1R. In the case of normal operation, we used the estimated parameters: JE — 3.90 kgm 2 and AT = 60sin(# + 210) Nm, obtained with the recursive gradient estimator. The RMS error has a value of 0.0559 rad/sec, which represents approximately 30% improvement over the result reported in Section 3.3. This identification procedure results in a RSR of 27%. This result is put in perspective by the fact that the occurrence of a 10% under-fueling fault in any of the six cylinders results in a RSR value between 55% and 159%. This RSR range is characteristic of a faulty condition. In addition, Chapter 4. Parameter Identification 42 0.5i r degrees Figure 4.2: Actual (solid line), and the estimated (dashed line) velocity fluctuation during one crankshaft revolution for the normal operation condition. The results correspond to the gradient estimator. degrees Figure 4.3: Actual (solid line), and estimated (dashed line) velocity fluctuation during one crankshaft revolution corresponding to an under-fueling fault in the first cylinder of the right bank (cylinder 6). The estimation is carried out by the on-line least-squares estimator with exponential forgetting. Chapter 4. Parameter Identification 43 Table 4.2: The parameter estimates, Th and 7 for Nm — 118. The estimates are obtained using the on-line gradient estimator. Faulty Cylinder Th (Nm) 7 (degrees) None 60 210 1 600 160 2 540 129 3 520 29 4 700 303 5 800 258 6 700 223 the detection of the faulty condition can be accomplished by observing the change in a parameter's mean and standard deviation. Our experiments show that the under-fueling of cylinder 1L causes a 6% change in the mean, and a 620% increase in the standard deviation of the estimated inertia. This result is obtained using the gradient estimator. For the same case, the least-squares estimator with exponential forgetting indicates a 16% change in the mean and a 327% change in the standard deviation of JE-Table 4.3 illustrates the RSRs for all 6 actuator faults. The RMS error for the case presented in Figure 4.3 is 0.0829 rad/sec. This calculation is carried out using the estimated parameters: JE = 4.91 kgm 2 and A T = 700sin(# + 223) Nm, obtained from the on-line estimator with exponential forgetting. The prediction of this under-fueling fault is therefore characterized by a RSR of 14% (see Table 4.3). The overall RMS error range for the faulty condition is between 0.0804 and 0.1315 rad/sec. The best results are obtained for faulty conditions in cylinders 3L and 1R. We considered the problem of fault isolation which deals with the correct classification of an engine cylinder as faulty or not. We considered the case when an under-fueling fault in cylinder 2 (3R) is misclassified as a fault in cylinder 1 (1L). This resulted in a RSR value of 36%, using the gradient estimator. Likewise, similar results were obtained Chapter 4. Parameter Identification 44 Table 4.3: Residual to signal ratio (RSR) for the 6 different under-fueling faults. Faulty Cylinder RSR (%) 1 20 2 24 3 20 4 13 5 16 6 14 for other misclassified cylinders. Alternatively, large changes in estimated parameter standard deviation can be used for the elimination of this misclassification. The gradient estimator shows a 227% change in the standard deviation of JE- The least-squares estimator with exponential forgetting records a 302% change for the same case. On-line parameter identification methods are powerful tools for increasing the predic-tion ability of the engine model and, at the same time, can be employed in fault detection and isolation by monitoring parameter estimate changes. The analytical redundancy provided by the engine model is enhanced further by the study of the inverse dynamics which allows the pressure waveform reconstruction from angular velocity fluctuation. This subject is treated in the next chapter. o Chapter 5 Pressure Waveform Reconst ruct ion It's here they got the range and the machinery for change and it's here they got the spiritual thirst. Leonard Cohen, Democracy In this chapter we investigate the procedure for pressure waveform reconstruction using the engine model and angular velocity measurements. This problem has multiple practical applications in engine control and diagnosis. It provides an indirect method for cylinder pressure measurement, and permits the implementation of a technique for estimating the torque due to gas pressure [36, 37, 38, 39] or the power contribution of each cylinder [40], and therefore allows the detection and identification of a variety of cylinder faults [41, 34]. The dynamics of a given system are represented by a collection of differential equations that permit the process output determination using the input history. The term "inverse dynamics" is used to characterize the calculation of the system input given the output history [11, pp. 263]. In our case, the diesel engine model is crank angle based and has N inputs, Ui(0) = Pj{6),i = 1,... iV, represented by the cylinder pressures. The system output is given by the angular velocity fluctuation, Sto(d). Because our model implements a multi-input single output (MISO) system, it is not possible to explicitly decouple the inverse dynamics. We have therefore considered two approaches. First, we redefine the system input and we rewrite the engine model as a single input single output (SISO) 45 Chapter 5. Pressure Waveform Reconstruction 46 system. The new input is then the torque due to gas pressure, TP(0), defined in (3.11) as: N Tp(e) = AprY/Pi(8)ci(8-i), i=\ where TV is the total number of cylinders, A p is the piston crown area, Pi(9) refers to the pressure in the i-th cylinder, fa is the phase shift corresponding to the i-th firing cylinder, and c\(9) is defined in (3.5) as: / x ( r c o s 8 \ ci(0) = 1 + - sinfl. y I cos a J Secondly, we approximate the pressure waveform with a periodic impulse-like continuous function. Using this latter template-based approach we are able to explicitly identify the condition of each cylinder. 5.1 Gas Pressure Torque Es t imat ion The estimation of cylinder power from engine speed fluctuations is considered in [40]. The authors use a linear second order model to describe the relationship between engine speed and crankshaft torque. An inverse filtering technique is employed for pressure torque calculation. The power contribution is determined as the area under this torque curve, corresponding to the power stroke of each cylinder. Similarly, the approach taken in [37] involves a time-based 4 degree-of-freedom vibration model for the system composed by the vibration damper, engine, flywheel, and dynamometer. The authors assume a purely harmonic solution for the crank angle: 9 = GejtJ\ where 0 is the amplitude, and to is the engine speed. The model is rewritten in the frequency domain as a matrix equation which is then solved for different harmonics of Chapter 5. Pressure Waveform Reconstruction 47 the engine fundamental frequency. The cylinder pressure is determined locally, in the vicinity of the pressure peak. This approach relies on the assumption that each cylinder contributes to the gas pressure torque only during its power stroke. In both these two papers, the engine model is developed in the time domain, and crankshaft/flywheel speed is sampled in the crank angle domain. The authors point, out that this inadvertence represents a source of error. The model used in [39] is written in the crank angle domain. The authors employ a deconvolution technique for the estimation of the pressure torque, Tp(9) from noisy measurements of the crankshaft angular acceleration. Their analysis is based on the assumption that the torque fluctuation is concentrated at a frequency equal to the firing frequency. This results in the neglect of higher harmonics. The solution adopted in [38] employs the design of an estimator in the crank angle domain. The authors use the following relationship to relate the crank angle and time bases: t = where 9 is the crank angle, t is the time, and LOQ is the nominal engine speed. In terms of the Laplace transform (£) this relationship becomes: sg = where st, and SQ describe the Laplace transform variable in the time domain and the crank angle domain, respectively. This is correct under the assumption of constant speed. A comprehensive analysis of the crank angle domain versus the time domain is re-ported in [42]. Here the authors model the engine rotational dynamics in the crank angle domain by a first order system characterized by the following transfer function: u>(sg) 1 H{sf)) = —-—- = -Tp(sg) JE^OSQ + b where u)(sg) — £{u>(8)}, Tp(se) = C{TP(9)}, JE is the engine mass moment of inertia, and b is the viscous friction coefficient. The authors point out that the assumption of constant speed introduces a fractional error of 6U/LOQ due to this approximation for the pole of the transfer function. Chapter 5. Pressure Waveform Reconstruction 48 A pattern recognition technique used for detection and isolation of faults due to cylinder under-fueling is reported in [34]. The authors constructed classes based on the fault signatures as observed in the flywheel angular velocity fluctuation data. A drawback of this approach is the assumption of reproducibility and linear scaling of the fault signatures with respect to the stored patterns. The method allows the estimation of the pressure peak, but doesn't allow the reconstruction of the actual pressure waveform. Our diesel engine model was defined in (3.26)-(3.27). The equations are repeated here: Sx' = A{9)Sx + bT{9)u{9) + d(9), Sto = c(6x), where the system coefficients A(9), b(B), d(9), and c(8x) were denned in (3.28)-(3.31). The system input is the vector of cylinder pressures: u{9) = PN{9) By redefining the system input as u(9) = Tp(9) we obtain the equivalent SISO system: Sx' = A(9)6x + b'(9)u(9) + die), where b'{9) is: b'(9) = -Apr (5.1) (5.2) 2 (JE - JcrC33{6) + m,pr2cu(9))' where Ap is the piston head area, r is the crank radius, mp is the mass of the reciprocating parts, and JE and J C R are the mass moments of inertia of the engine and the connecting rod, respectively. The torque coefficient cn(9) was defined in (3.14). The system state, 6x was defined in (3.24) as: Sx = 5LO(2LUQ + 6u) f a 2LUQ6U>. Chapter 5. Pressure Waveform Reconstruction 49 Using the approximation we obtain: fiW = f , w = ^ ^ W + ^ ) f a W ] - < » ( » ) , ( 5 , 3 ) which allows the estimation of the gas pressure torque using the measurement of angular velocity fluctuation. 5.2 Pressure Waveform Characterist ics Cylinder pressure is uniquely connected to the combustion process. Cycle-by-cycle and cylinder-by-cylinder pressure variations are correlated with the injection timing, the quantity of fuel injected, and the rate of mixing between the injected fuel and the air. The most important pressure-related parameters are: maximum cylinder pressure, the crank angle at which this maximum pressure occurs, the maximum rate of pressure rise, and the corresponding crank angle, and the indicated mean effective pressure (IMEP) [43, pp. 415]. A measure of the pressure cyclic variability is the coefficient of variability (COV) defined in [43, pp. 417]: GOV = x 100, imep where o"i m e p is the standard deviation in indicated mean effective pressure. The author asserts that vehicle driveability problems are characterized by a COV> 10%. The cycle-by-cycle variation in cylinder pressure was modeled by a Gaussian proba-bility distribution in [44]. A similar approach was taken in [5]. Here the authors modeled the cyclic pressure variation as a raised-cosine window amplitude-modulated by a white Bernoulli-Gaussian random sequence. The same authors reported their experimental re-sults in [45]. The proposed stochastic model describes the gross pressure waveforms, but doesn't accurately characterize the instantaneous shape of the actual pressure waveforms. Chapter 5. Pressure Waveform Reconstruction 50 Rather than trying to determine a stochastic description of the cycle-by-cycle and cylinder-by-cylinder pressure variation, our approach is to find a functional approxima-tion for the cylinder pressure waveform, and subsequently relate the pressure variation to a minimum number of parameters. This approach is based on the observation that the cylinder pressure is an impulse-like periodic function. We consider a general description for a family of functions, related to the <5-function, given in [46, pp. 487-492]: \ / W d r i . ] . . 8(w,z) = = — tan (wz) , (5.4) wLz* + 1 az 1 J where z € R, and w is a variable parameter. When z = sin with 0 the crank angle, and 4>i the phase of the i-th cylinder, the proposed impulse functional approximation for the pressure waveform is: w,sm—^)-P0 = . 2V, -Po, (5.5) l ) wz s in —f1- + 1 where p P a x = kjW represents the maximum pressure value, fc7; > 0 is a factor related to the pulse height, w > 0 is inversely proportional to the width of the pressure pulse, and Po is a pressure offset. Taking into consideration (5.3), and (5.5), we can write the estimated gas pressure torque in a form similar to (3.11): N fp(9) = Apr 6PiPi(0)Cl(0 - fa), (5.6) ?:=i where 8Pi,i = 1,. . . N are parameters that describe the pressure variation with respect to the template given by Pi{8). The estimation of these parameters allows us to determine the condition of each cylinder (faulty or not). The standard least-squares estimator is: 8P = ( $ T $ ) - 1 $ T f p , (5.7) fp = $8P, (5.8) Chapter 5. Pressure Waveform Reconstruction 51 where 6P1 SP = SPN and $ = Apr P^cxiO - fa) PN(6)Cl{9 - 0/\r) We analyzed the use of the engine model in solving the inverse dynamics, and thus pro-viding an estimate for the torque due to gas pressure. We investigated a meaningful procedure for pressure waveform reconstruction based on the functional approximation introduced in (5.5). We are able to detect and isolate cylinder under-fueling faults us-ing a standard least-squares identification method. The steps of our procedure can be summarized as follows: 1. Calculate the derivative Sui', using the measured angular velocity fluctuation, Sui. 2. Calculate the ini t ia l estimated pressure torque, Tp, using (5.3). 3. Obta in the least-squares estimates of the pressure variations, SP,, using (5.5) and 4. Determine the new estimate of the gas pressure torque, Tp, using (5.8). 5. Reconstruct the cylinder pressure waveforms using (5.5) and (5.7). 5.3 S imula t ion Methods The simulation is performed with the M A T L A B package and the programs are listed in Appendix B . The main functions are: • diffc.m which implements a smooth differentiator, (5.7). Chapter 5. Pressure Waveform Reconstruction 52 • tpe.m which implements (5.3), • dlt.m which implements (5.5), and • pwr.m which calculates the pressure variation vector, SP. The angular velocity is sampled at a period of 7r /59 rad. This sampling period corre-sponds to the total number of flywheel teeth, 118. The cylinder pressure data is char-acterized by a sampling period of 7r/360 rad. The model parameters are given the same numerical values as in Sections 3.3 and 4.3. These are: N = 6, r = 0.0635 m, I = 0.2571 m, Jcr = 0.0745 kgm 2 , m,p = 6.03 kg, Ap = 0.01188 m 2 . The phase origin is the top dead center of the first cylinder to fire, and 0 i = 0, 02 = 56.5, 03 = 120, 04 = 176.5, 05 = 240, 06 = 296.5 degrees. The firing order is again: IL, 3R, 3L, 2R, 2L, 1R. Data acquisition is performed for a nominal speed LUQ = 125.7 rad/sec (1,200 rpm) and a dynamometer torque, T/ = 990 Nm. This load represents 75% of the engine peak torque. Tf = 156 Nm is the value of the mechanical friction torque corresponding to a crankshaft angular velocity of UQ = 125.7 rad/sec and an overall damping coefficient of 1.24 Nmsec. 5.3.1 T o r q u e E s t i m a t i o n a n d P r e s s u r e A p p r o x i m a t i o n The testing set consists of 6 different actuator faults and the normal operating condition. The actuator fault consists of a 10% drop in cylinder fueling. The root mean square (RMS) error is defined as the standard deviation of the torque residual, Tp = Tp — Tp: The gas pressure torque is estimated from (5.3). The residual to signal ratio (RSR) is: (5.9) RSR = RMS error 100% (5.10) °(TP) Chapter 5. Pressure Waveform Reconstruction 53 degrees Figure 5.1: Actual (solid line), and the estimated (dashed line) gas pressure torque fluctuation during one crankshaft revolution for the normal operation condition. where cr(Tp) is the standard deviation of the actual gas pressure torque. The same performance indicators are used to assess the magnitude of the pressure residual, P = P - P . Figures 5.1 and 5.2 present the estimated gas pressure torque, Tp(6), corresponding to two test cases—normal operating condition, and an under-fueling fault in cylinder 1R. The results are obtained after the elimination of the harmonic torque fluctuation, A T = r^sin(i9 + 7 ) , where Th is the amplitude of the torque fluctuation, and 7 is the phase (see Chapter 4). In the case of normal operation (see Figure 5.1), the engine moment of inertia is JE — 3.85 kgm 2 , and the mean value for the gas pressure torque is Tp = 856 Nm. The resulting RMS error has a value of 196.2 Nm, which corresponds to a RSR of 35%. In the case of an under-fueling fault in cylinder 1R (see Figure 5.2), the engine moment of inertia is JE = 5 kgm 2, and the mean value for the gas pressure torque is Tp = 1215 Nm. The resulting RMS error has a value of 271.4 Nm, which corresponds Chapter 5. Pressure Waveform Reconstruction 54 d e g r e e s Figure 5.2: Actual (solid line), and estimated (dashed line) gas pressure torque fluctuation during one crankshaft revolution corresponding to an under-fueling fault in the first cylinder of the right bank (cylinder 6). to a RSR of 38%. The gas pressure torque for the faulty operation is estimated by a RMS torque error range of 271.4 Nm to 308.7 Nm, which corresponds to a RSR range of 38% and 43%. The worst results are associated with a fault in cylinder 2L. The best results are obtained for an under-fueling fault in cylinder 1R. Figure 5.3 illustrates the actual and the estimated pressure waveform for the normal operation of cylinder 3L. We used (5.5) and the following values: P 3 m a x = 11.8 MPa, w = 5.5, and P0 = 0.5 MPa for the approximation of P3(9). The RMS pressure error is 0.149 MPa. This corresponds to a RSR of 5%. For the normal operating condition, this functional approximation is characterized by a RMS pressure error range of 0.145 M P a to 0.177 MPa, which corresponds to a RSR range of 5% to 6%. The smallest error is obtained for cylinder IL. Chapter 5. Pressure Waveform Reconstruction 55 1 2 , r I i I I I I 1 1 O 5 0 1 0 0 1 5 0 2 0 0 2 5 0 3 0 0 3 5 0 crank angle (degrees) F i g u r e 5.3: Illustrates a typical cylinder pressure waveform (solid line) and its approximat ion (circles) as per (5.5). The phase corresponds to the th i rd firing cylinder. 5.3.2 Fault Detect ion and Isolation Results We considered the problem of fault detection and isolation. The procedure uses 6 pressure templates determined by taking into consideration each cylinder firing phase, fa. A l l the 6 pressure waveforms are characterized by (5.5) with the same values p P a x = p m a x = n MPa, and uii = w = 5.5, i = 1,... 6. The estimated parameters, <5P7; are obtained using a standard least-squares technique, and their values are presented in Table 5.1. A n under-fueling fault in the i-th firing cylinder can be determined exactly by the value of the estimated pressure variation, <5P, using the decision rule: 6Pi = min 6Pk. (5.11) fc=l,...6 This detection problem can be formulated also from the perspective of the Bayes classifier [47, pp. 221-233]. The two classes are determined from the estimates presented in Table 5.1: faulty cylinder, which corresponds to a mean pressure variation, Spf — 0.83 MPa, and healthy Chapter 5. Pressure Waveform Reconstruction 56 Table 5.1: Estimated pressure variation parameters, SP, (MPa) for the normal condition, and the 6 different under-fueling faults using the impulse like template described by (5.5). Faulty Cylinder SPi 6P2 SP3 8P4 SP5 6P6 None 1.01 1.09 1.09 1.10 1.15 1.06 1 0.78 1.05 1.08 1.11 1.10 0.93 2 0.90 0.79 1.03 1.09 1.14 0.95 3 0.96 0.89 0.80 0.96 1.16 1.07 4 1.05 1.04 0.95 0.87 1.10 1.02 5 1.01 1.10 1.02 0.99 0.93 0.99 6 0.97 1.10 1.07 1.05 1.07 0.78 cylinder, which corresponds to a mean pressure variation, Sph = 1-04 MPa. The noise is assumed Gaussian with zero mean and a calculated standard deviation of crn — 0.07. If the data is assumed Gaussian, 95% of the estimated pressure variation samples, <5P,; lie within ±2er n of the mean estimate, Sph-0.9 < SPi < 1.18. Large changes in the value of the estimated parameter can be used for the detection and isolation of the faulty engine cylinder. For example, a fault in the first firing cylinder (IL) results in a 23% decrease in 6Pi, and a fault in the second firing cylinder (3R) results in a 28% decrease in SP2. Under-fueling faults are thus uniquely determined by the values of the estimated pressure variation, Sp, i — 1,... 6. The same parameter allows the reconstruction of the cylinder pressure waveform by correcting the value of pP3* initially assigned in (5.5): The pressure waveform reconstruction is characterized by a RMS pressure error range of 0.155 MPa to 0.277 MPa for the normal operating condition. This corresponds to a RSR range of 5% to 10%. Chapter 5. Pressure Waveform Reconstruction 57 12 i i 10 r \ 8 -«* 6 Q_ / " 4 2 O o 50 100 150 200 d e g r e e s 250 300 350 Figure 5.4: Actual (solid line), and estimated (dashed line) gas pressure waveform during one crankshaft revolution corresponding to an under-fueling fault in cylinder 1R. For each of the six under-fueling faults, the pressure waveform corresponding to the faulty cylinder is reconstructed. The RMS error range is of 0.155 MPa to 0.386 MPa. The corresponding RSR range is 7% and 18%. The best results (RSR from 5% to 9%) are obtained in the case of an under-fueling fault in the six firing cylinder (1R). This situation is illustrated in Figure 5.4- Using this correction (see (5.12)) we are also able to improve the estimation of the gas pressure torque, Tp{6). The new estimate is plotted in Figure 5.5, for the normal operation. The RMS torque error is reduced to 99.24 Nm, and this corresponds to a RSR of 14%. The case of an under-fueling fault in cylinder 1R is illustrated in Figure 5.6. The RMS torque error reduces to 85.7 Nm, and the RSR is 12%. These results are compatible with those reported in [34] for the calculation of only the peak cylinder pressure. Our procedure has the advantage that the overall pressure waveform can be reconstructed using (5.5) and the correction (5.12). This approach is Chapter 5. Pressure Waveform Reconstruction 58 -iooo\-I I I I I I I I O 50 100 150 200 250 300 3 5 0 degrees Figure 5.5: Actual (solid line), and the estimated (dashed line) gas pressure torque during one crankshaft revolution for the normal operating condition, using the correction described in (5.12). T 1 1 1 1 1 r degrees Figure 5.6: Actual (solid line), and estimated (dashed line) gas pressure torque during one crankshaft revolution corresponding to an under-fueling fault in cylinder 1R, taking into con-sideration (5.12) Chapter 5. Pressure Waveform Reconstruction 59 simpler and less computationally demanding than the stochastic techniques reported in [5, 44]. Chapter 6 Summary, Discussion, and Conclusions What we call the beginning is often the end And to make an end is to make a beginning The end is where we start from. T. S. Eliot, Little Gidding The concept of analytical redundancy emphasizes the use of accurate dynamic and static models for data processing and analysis. The major benefit is realized in the low cost and flexibility of a software implementation versus a hardware implementation. A combina-tion of analytical and physical redundancy is almost always necessary for a fault-tolerant system to maintain its function in the presence of certain failures. On the other hand, the paradigm of analytical redundancy has the advantage of fully exploiting the engine model and thus extracting information that otherwise might be difficult to obtain. This was the motivation for pursuing three directions of investigation: engine modeling and validation, parameter identification, and pressure waveform reconstruction. Chapter 2 established a global framework for our analysis. The field of model-based fault detection and isolation was linked to the problem of system state and parameter identification. The latter option, parametric identification, was pursued in our study because of its ability to detect early system faults. 60 Chapter 5. Pressure Waveform Reconstruction 51 where SPi 8P = 5PN and $ = v PMciie - fa), PN{e)Cl{e - N) We analyzed the use of the engine model in solving the inverse dynamics, and thus pro-viding an estimate for the torque due to gas pressure. We investigated a meaningful procedure for pressure waveform reconstruction based on the functional approximation introduced in (5.5). We are able to detect and isolate cylinder under-fueling faults us-ing a standard least-squares identification method. The steps of our procedure can be summarized as follows: 1. Calculate the derivative 6u>', using the measured angular velocity fluctuation, 8to. 2. Calculate the initial estimated pressure torque, Tp, using (5.3). 3. Obtain the least-squares estimates of the pressure variations, 8Pi, using (5.5) and 4. Determine the new estimate of the gas pressure torque, Tp, using (5.8). 5. Reconstruct the cylinder pressure waveforms using (5.5) and (5.7). 5.3 S imula t ion Methods The simulation is performed with the M A T L A B package and the programs are listed in Appendix B. The main functions are: • diff cm which implements a smooth differentiator, (5.7). Chapter 5. Pressure Waveform Reconstruction 52 • tpe.m which implements (5.3), • dlt.m which implements (5.5), and • pwr.m which calculates the pressure variation vector, SP. The angular velocity is sampled at a period of -zr/59 rad. This sampling period corre-acterized by a sampling period of 7r/360 rad. The model parameters are given the same numerical values as in Sections 3.3 and 4.3. These are: TV = 6, r = 0.0635 m, I = 0.2571 m, Jcr = 0.0745 kgm 2 , mp = 6.03 kg, Ap = 0.01188 m 2 . The phase origin is the top dead center of the first cylinder to fire, and fa = 0, fa — 56.5, 03 = 120, 04 = 176.5, fa = 240, 06 = 296.5 degrees. The firing order is again: IL, 3R, 3L, 2R, 2L, 1R. Data acquisition is performed for a nominal speed U>Q — 125.7 rad/sec (1,200 rpm) and a dynamometer torque, T/ = 990 Nm. This load represents 75% of the engine peak torque. Tf = 156 Nm is the value of the mechanical friction torque corresponding to a crankshaft angular velocity of UJQ = 125.7 rad/sec and an overall damping coefficient of 1.24 Nmsec. 5.3.1 T o r q u e E s t i m a t i o n a n d P r e s s u r e A p p r o x i m a t i o n The testing set consists of 6 different actuator faults and the normal operating condition. The actuator fault consists of a 10% drop in cylinder fueling. The root mean square (RMS) error is defined as the standard deviation of the torque residual, Tp = Tp — Tp: The gas pressure torque is estimated from (5.3). The residual to signal ratio (RSR) is: sponds to the total number of flywheel teeth, 118. The cylinder pressure data is char-RMS torque error = (5.9) RSR = RMS error 100% (5.10) Chapter 5. Pressure Waveform Reconstruction 53 T 1 1 1 1 r d e g r e e s Figure 5.1: Actual (solid line), and the estimated (dashed line) gas pressure torque fluctuation during one crankshaft revolution for the normal operation condition. where cr(Tp) is the standard deviation of the actual gas pressure torque. The same performance indicators are used to assess the magnitude of the pressure residual, P = P-P. Figures 5.1 and 5.2 present the estimated gas pressure torque, Tp(9), corresponding to two test cases—normal operating condition, and an under-fueling fault in cylinder 1R. The results are obtained after the elimination of the harmonic torque fluctuation, A T = ThSm(0 + 7), where Th is the amplitude of the torque fluctuation, and 7 is the phase (see Chapter 4). In the case of normal operation (see Figure 5.1), the engine moment of inertia is JE = 3.85 kgm 2, and the mean value for the gas pressure torque is T p = 856 Nm. The resulting RMS error has a value of 196.2 Nm, which corresponds to a RSR of 35%. In the case of an under-fueling fault in cylinder 1R (see Figure 5.2), the engine moment of inertia is JE = 5 kgm 2, and the mean value for the gas pressure torque is Tp = 1215 Nm. The resulting RMS error has a value of 271.4 Nm, which corresponds Chapter 5. Pressure Waveform Reconstruction 54 degrees Figure 5.2: Actual (solid line), and estimated (dashed line) gas pressure torque fluctuation during one crankshaft revolution corresponding to an under-fueling fault in the first cylinder of the right bank (cylinder 6). to a RSR of 38%. The gas pressure torque for the faulty operation is estimated by a RMS torque error range of 271.4 Nm to 308.7 Nm, which corresponds to a RSR range of 38% and 43%. The worst results are associated with a fault in cylinder 2L. The best results are obtained for an under-fueling fault in cylinder 1R. Figure 5.3 illustrates the actual and the estimated pressure waveform for the normal operation of cylinder 3L. We used (5.5) and the following values: p™** = H.8 MPa, w = 5.5, and Po = 0.5 MPa for the approximation of Pz(9). The RMS pressure error is 0.149 MPa. This corresponds to a RSR of 5%. For the normal operating condition, this functional approximation is characterized by a RMS pressure error range of 0.145 MPa to 0.177 MPa, which corresponds to a RSR range of 5% to 6%. The smallest error is obtained for cylinder 1L. Chapter 5. Pressure Waveform Reconstruction 55 crank angle (degrees) Figure 5.3: Illustrates a typical cylinder pressure waveform (solid line) and its approximation (circles) as per (5.5). The phase corresponds to the third firing cylinder. 5.3.2 Fault Detect ion and Isolation Results We considered the problem of fault detection and isolation. The procedure uses 6 pressure templates determined by taking into consideration each cylinder firing phase, fa. A l l the 6 pressure waveforms are characterized by (5.5) with the same values p™3* = p m a x — \ \ MPa, and = w = 5.5, i = 1,.. . 6. The estimated parameters, SPi are obtained using a standard least-squares technique, and their values are presented in Table 5.1. A n under-fueling fault in the i-th firing cylinder can be determined exactly by the value of the estimated pressure variation, 6Pi using the decision rule: SPi = min SPk. (5.11) fc=l,...6 This detection problem can be formulated also from the perspective of the Bayes classifier [47, pp. 221-233]. The two classes are determined from the estimates presented in Table 5.1: faulty cylinder, which corresponds to a mean pressure variation, Spf = 0.83 MPa, and healthy Chapter 5. Pressure Waveform Reconstruction 56 Tab l e 5.1: Es t imated pressure variat ion parameters, SPi (MPa ) for the normal condit ion, and the 6 different under-fueling faults using the impulse like template described by (5.5). F a u l t y C y l i n d e r SPi SP2 SP3 8P4 SP5 SP6 None 1.01 1.09 1.09 1.10 1.15 1.06 1 0.78 1.05 1.08 1.11 1.10 0.93 2 0.90 0.79 1.03 1.09 1.14 0.95 3 0.96 0.89 0.80 0.96 1.16 1.07 4 1.05 1.04 0.95 0.87 1.10 1.02 5 1.01 1.10 1.02 0.99 0.93 0.99 6 0.97 1.10 1.07 1.05 1.07 0.78 cylinder, w h i c h corresponds to a mean pressure va r i a t i on , Sph = 1-04 M P a . T h e noise is assumed G a u s s i a n w i t h zero m e a n and a ca l cu la t ed s t a n d a r d dev i a t i on of ern = 0.07. If the d a t a is assumed G a u s s i a n , 9 5 % of the es t imated pressure va r i a t i on samples , SPi l ie w i t h i n ±2cr n of the m e a n est imate , Sph'-0.9 < SPi < 1.18. La r ge changes i n the va lue of the es t imated parameter can be used for the de t ec t i on a n d i so l a t i on of the fau l ty engine cy l inder . Fo r example , a fault i n the first f i r ing c y l i nde r ( I L ) resul ts i n a 2 3 % decrease i n SP\, and a fault i n the second firing cy l inde r (3R) resul ts i n ' a 2 8 % decrease i n SP2. Under - fue l i ng faul ts are thus un ique l y de te rmined by the values of the e s t ima ted pressure v a r i a t i o n , SPi, i — 1 , . . . 6. T h e same parameter a l lows the r e cons t ruc t i on of the cy l inde r pressure wave form by cor rec t ing the value of _ p . m a x i n i t i a l l y ass igned i n (5.5): p ™ a x = 8PiPma*. (5.12) T h e pressure wave form recons t ruc t i on is character i zed by a R M S pressure-error range of 0.155 M P a to 0.277 M P a for the n o r m a l ope ra t ing cond i t i on . T h i s corresponds to a R S R range of 5% to 10%. / Chapter 5. Pressure Waveform Reconstruction 57 12 10 8 - 'A ™ e 4 -2 - J O 0 so 100 150 200 degrees 250 300 350 Figure 5.4: Actual (solid line), and estimated (dashed line) gas pressure waveform during one crankshaft revolution corresponding to an under-fueling fault in cylinder 1R. For each of the six under-fueling faults, the pressure waveform corresponding to the faulty cylinder is reconstructed. The RMS error range is of 0.155 MPa to 0,386 MPa. The corresponding RSR range is 7% and 18%. The best results (RSR from 5% to 9%) are obtained in the case of an under-fueling fault in the six firing cylinder (1R). This situation is illustrated in Figure 5.4. Using this correction (see (5.12)) we are also able to improve the estimation of the gas pressure torque, Tp(9). The new estimate is plotted in Figure 5.5, for the normal operation. The RMS torque error is reduced to 99.24 Nm, and this corresponds to a RSR of 14%. The case of an under-fueling fault in cylinder 1R is illustrated in Figure 5.6. The RMS torque error reduces to 85.7 Nm, and the RSR is 12%. These results are compatible with those reported in [34] for the calculation of only the peak cylinder pressure. Our procedure has the advantage that the overall pressure waveform can be reconstructed using (5.5) and the correction (5.12). This approach is Chapter 5. Pressure Waveform Reconstruction 58 T 1 1 i 1 • r d e g r e e s F i g u r e 5.5: A c t u a l ( so l id l ine ) , a n d th e e s t i m a t e d ( dashed l ine ) gas p r essure t o r q u e d u r i n g one c r a n k s h a f t r e v o l u t i o n for the n o r m a l o p e r a t i n g c o n d i t i o n , u s i n g t h e c o r r e c t i o n d e s c r i b e d i n (5.12) . d e g r e e s F i g u r e 5.6: A c t u a l ( so l id l ine ) , a n d e s t i m a t e d ( dashed l ine ) gas p r essure t o r q u e d u r i n g one c r a n k s h a f t r e v o l u t i o n c o r r e s p o n d i n g to a n u n d e r - f u e l i n g f au l t i n c y l i n d e r 1 R , t a k i n g i n t o c o n -s i d e r a t i o n (5.12) Chapter 5. Pressure Waveform Reconstruction 59 simpler and less computationally demanding than the stochastic techniques reported in [5, 44]. Chapter 6 Summary, Discussion, and Conclusions What we call the beginning is often the end And to make an end is to make a beginning The end is where we start from. T. S. Eliot, Little Gidding The concept of analytical redundancy emphasizes the use of accurate dynamic and static models for data processing and analysis. The major benefit is realized in the low cost and flexibility of a software implementation versus a hardware implementation. A combina-tion of analytical and physical redundancy is almost always necessary for a fault-tolerant system to maintain its function in the presence of certain failures. On the other hand, the paradigm of analytical redundancy has the advantage of fully exploiting the engine model and thus extracting information that otherwise might be difficult to obtain. This was the motivation for pursuing three directions of investigation: engine modeling and validation, parameter identification, and pressure waveform reconstruction. Chapter 2 established a global framework for our analysis. The field of model-based fault detection and isolation was linked to the problem of system state and parameter identification. The latter option, parametric identification, was pursued in our study because of its ability to detect early system faults. 60 Chapter 6. Summary, Discussion, and Conclusions 61 6.1 M o d e l i n g and Val ida t ion The model we developed in Chapter 3 was configured to correspond to the D D C 6V 92TA Detroit Diesel engine. The analytical methodology followed that of C. F. Taylor [28, pp. 240-305], and the testing was performed using data supplied by the Institute for Machinery Research, National Research Council. The system inputs were given by cylinder combustion pressures, and its output is represented by the flywheel angular velocity fluctuation. Both are expressed as functions of the crank angle. A drawback of our modeling technique derives from the assumptions that the crankshaft is perfectly rigid and the engine operates at steady state. This was compensated by the accuracy of the model which, after a change of state variable, is a first order linear ^-variant multi-input single output system. The diesel engine model was validated using the M A T L A B simulation package, and its effectiveness was measured by the root mean square (RMS) error, and the residual to signal ratio (RSR). The use of an unified performance index allowed us to assess each method using a similar scale, and thus provided the first step towards the integration of these three techniques in a condition monitoring system. 6.2 Parameter Identification The prediction capability of the model was improved via a parametric identification procedure in Chapter 4. The estimated parameters were the engine mass moment of inertia, and the harmonic torque fluctuation, which we attributed to the crankshaft vibration and oscillatory behaviour of the dynamometer loading. The estimation of this unmodeled torque reduced the error created by those modeling assumptions. We used three identification methods: off-line standard least-squares, on-line gradient estimator, and recursive least-squares with exponential forgetting. The results obtained Chapter 6. Summary, Discussion, and Conclusions 62 were similar with these three methods. For the normal operation we obtained an esti-mated engine inertias between 3.85 kgm 2 and 3.90 kgm 2. For the six fault conditions, the estimated engine inertia had values between 4.88 kgm 2 and 5.17 kgm 2 . In the case of cylinder under-fueling faults, the amplitude of the estimated torque fluctuation is approximately 10 times larger than in the case of normal operation. The correct prediction of the six under-fueling faults was characterized by a RSR range of 13% to 24%. In addition, the detection of a cylinder under-fueling fault condition was accomplished by observing the change in a parameter's mean and standard deviation. For example, the under-fueling of cylinder IL caused a 6% change in the mean, and a 620% increase in the standard deviation of the inertia estimated with the gradient estimator. For the same case, the least-squares estimator with exponential forgetting indicated a 16% change in the mean and a 327% change in the standard deviation of the estimated engine inertia. Alternatively, large changes in estimated parameter standard deviation can be used for the elimination of fault misclassification. We considered the case when an under-fueling fault in cylinder 3R is misclassified as a fault in cylinder IL. The gradient estimator showed a 227% change in the standard deviation of estimated inertia. The least-squares estimator with exponential forgetting recorded a 302% change for the same variable. 6.3 Pressure Waveform Reconstruct ion In Chapter 5 we investigated the procedure for pressure waveform reconstruction using the engine model and angular velocity measurements. Because our model implements a multi-input single output system, it was not possible to explicitly decouple the inverse dynamics. We therefore considered two approaches. First, we redefined the system input and we rewrote the engine model as a single input single output system. The new input Chapter 6. Summary, Discussion, and Conclusions 63 was then the torque due to gas pressure. Secondly, we approximated the cylinder pressure waveforms with a periodic continu-ous impulse-like function. Using this template-based approach we were able to exactly identify the condition of each cylinder, and to improve the estimation procedure for the gas pressure torque. Initial simulation of normal operating condition resulted in a RMS error of 196.2 Nm, while the template-based approach resulted in a RMS torque error of only 99.24 Nm. In the case of an under-fueling fault, the gas pressure torque was initially estimated with a RMS error between 271.4 Nm and 308.7 Nm, while the pressure-template approach reduce the final RMS torque error to the range of 85.7 Nm to 198.1 Nm. 6.3.1 Ind iv idua l Cy l inder Pressure Reconstruct ion Cylinder pressure waveforms were approximated by an impulse-like periodic function. For the normal operating condition, this functional approximation was characterized by a RMS pressure error between 0.145 MPa and 0.177 MPa, which correspond to a RSR range of only 5% to 6%. The torque due to gas pressure waveform reconstruction was obtained from the super-position of the six individual pressure waveform templates. A l l 6 pressure waveforms are initially characterized by the same maximum value and width. The estimated pressure variations were obtained using a standard least-squares technique. A n under-fueling fault in the i-th cylinder could be determined exactly by the minimum value of the estimated pressure variation. The RMS pressure error was between 0.155 MPa and 0.277 M P a for the normal operating condition corresponding to a RSR between 5% and 10%. For each of the six under-fueling faults, the pressure waveform corresponding to the faulty cylinder was reconstructed. The RMS error was between 0.155 M P a and 0.386 MPa. The RSR was between 7% and 18%. Chapter 6. Summary, Discussion, and Conclusions 64 6.4 Towards Cond i t ion M o n i t o r i n g In the foregoing work, we have laid the groundwork for a diesel engine condition moni-toring system. In short, we have developed and validated a model, and shown its utility in detecting faults using two separate strategies—parameter estimation, and cylinder pressure reconstruction. In the latter strategy, we have introduced a new template-based scheme for estimating individual cylinder pressures. These results suggest a condition monitoring scheme which is highly sensitive to input variation due to fueling faults. Fur-thermore, reconstructed torque due to gas pressure by superposition of the templates results in very accurate estimates—RSR of about 5%—which are as good as the best results reported in the literature. The integration of these three methods (modeling and validation, parametric iden-tification, and pressure waveform reconstruction) exploits the use of the diesel engine model from the perspective of the analytical redundancy paradigm. Our simulation re-sults suggest that parametric identification methods represent important tools for solving engine fault detection and isolation problems. Parameter means and standard deviations were used to differentiate between healthy and faulty conditions, and to eliminate mis-classifications of faulty cylinders. Using a combination of functional approximation and parametric identification techniques we were able to determine the condition of each cylinder and to reconstruct close approximations to the cylinder pressure waveforms. j 6.5 Conclusions A new method for computing individual cylinder pressures from crankshaft velocity vari-ation has been developed. The pressure template-based method provides estimates of pressure to an accuracy of 5% residual to signal ratio. The single cylinder model of C. F. Taylor [28, pp. 240-305] has been transformed Chapter 6. Summary, Discussion, and Conclusions 65 to a model based on engine speed fluctuation as a function of pressure input. The engine model was further modified to multiple cylinders, and was linearized in state and parameters. The resulting model was validated using recorded data to produce a RMS estimation error of 0.0559 rad/sec for the normal operating condition. The above model, improved by parametric identification allowed the correction for harmonic torque fluctuation, and the identification of total engine inertia. Nomenclature We have used the following notation: m JE = engine mass moment of inertia (kgm2); • Jcr = connecting rod mass moment of inertia (kgm2); • mt + mr = total connecting rod mass (kg); • rn,p = mass of the reciprocating parts (kg); • Ap — piston crown area (m2); • r — crankshaft radius (m); • I = connecting rod length (m); • TE = total engine torque (Nm); • Tp = torque due to gas pressure (Nm) and thus, Tp = estimated gas pressure torque (Nm), Tp = torque residual, and Tp — mean pressure torque (Nm); • Tt = torque due to the inertia of the reciprocating masses (Nm); • Tr = torque due to the angular acceleration of the connecting rod (Nm); • Tf = friction torque (Nm); • Ti = load torque (Nm); • A T = unmodeled harmonic torque fluctuation (Nm) and thus, A T = estimated torque fluctuation (Nm); 66 Nomenclature 67 • Th = torque fluctuation amplitude (Nm); • 7 = torque phase (degrees); • N = number of cylinders, or number of system inputs; • Nm = number of measurements (data points); • 9 = crank angle (degrees or rad), and thus, 9 = LO angular velocity (rad/sec), and 9 = co angular acceleration (rad/sec2); • LOQ = nominal engine speed (rpm or rad/sec); • 8u> = measured angular velocity fluctuation, and thus 8to = estimated crankshaft angular velocity fluctuation, and SCo = angular velocity residual (rad/sec); • Qh = velocity residual amplitude (rad/sec); • a = angle between the connecting rod axis and the cylinder axis (degrees or rad), and thus, a = connecting rod angular velocity (rad/sec), and a = connecting rod angular acceleration (rad/sec2); • s = piston position (m), and thus, i = piston velocity (m/sec), and s — piston acceleration (m/sec2); • Pi = pressure in the i-th cylinder (MPa); • 8Pi = pressure variation and thus, 8Pi = estimated pressure variation; • 4>i = phase shift corresponding to the i-th cylinder (degrees); • x G Rn system state with XQ = initial condition, and thus, x = estimated state, and x = state estimation error. In the case of the engine model, x = LO2, and the state is described by the variable 8x G R; Nomenclature 68 • u E RN system input, u estimated system input, and uc = control law. For the engine model the input is represented by N = 6 cylinder pressures; • y & R system output, and thus, y — predicted output, and y = prediction error. For the engine model, y = Scu; • Y = corrected engine torque, output of the identification model; • yd = desired trajectory, and e = tracking error; • p G Rq vector of system parameters and thus, p = parameter estimates, and p = parameter estimation error; • A, c, and B, b or b' = engine model coefficients; • i] — measurement noise; • d = process disturbance; • H = gain of the Beard-Jones filter; • c i ) c2> c3, c4 = geometrical coefficients; • cn, C12, C33, C34 = torque coefficients; • E = performance index; • $ = signal matrix; • XJ/Q = gain matrix for the gradient estimator; • \T/ = variable gain matrix for the least-squares estimator with exponential forget-ting; • A = forgetting factor; Nomenclature • a — standard deviation. Bibl iography [1] S. E. Salcudean, Introduction to Robotics (course notes). Dept. of E E — U B C , 1993. [2] P. E. Wellstead, Physical System Modelling. Academic Press, 1979. [3] Diesel Dynamic Model Mathematical Analysis, Technical Report # GTL-19-3^-TR.l. GasTOPS Ltd., March 1991. [4] A . K . Sood, A . A . Fahs, and N . A Henein, "Engine fault analysis: Part Il-Parameter estimation approach," IEEE Trans. Automat. Contr., vol. 32, pp. 301-307, Nov. 1985. [5] F. T. Connolly and A. F. Yagle, "Modeling and identification of the combustion pressure process in internal combustion engines using engine speed fluctuations," in Transportation Systems-92: Winter Annual Meeting of the American Society of Mechanical Engineers, pp. 191-206, Nov. 1992. [6] R. Stone, Introduction to Internal Combustion Engines. Macmillan Press, 2nd ed., 1992. [7] C. F. Taylor, The Internal-Combustion Engine in Theory and Practice, Volume I: Thermodynamics, Fluid Flow, Performance. The Technology Press of the MIT and John Willey & Sons, 1960. [8] K . B. Hodgins and P. Mtui, Evaluation of PCB Pressure Transducers, Internal Report. Dept. of M E — U B C , June 1994. [9] P. G. Hil l and P. D. Lawrence, Diesel Engine Advanced Condition Monitoring System (Research Proposal—DRAFT). U B C , 1994. [10] M . J. Jennings, P. N . Blumberg, and R. W. Amann, "A dynamic simulation of the Detroit Diesel electronic control system in heavy duty powertrains," SAE Trans., J. Engines, vol. 95, pp. 942-966, 1986. Section 3, Paper No. 861959. [11] J. J. Slotine and W. L i , Applied Nonlinear Control. Prentice Hall, 1991. [12] P. M . Frank, "Fault diagnosis in dynamic systems using analytical and knowledge-based redundancy-A survey and some new results," Automatica, vol. 26, pp. 459-474, Sep. 1990. 70 Bibliography 71 [13] A . S. Willsky, "A survey of design methods for failure detection in dynamic sytems," Automatica, vol. 12, pp. 601-611, Dec. 1976. [14] R. Isermann, "Process fault detection based on modeling and estimation methods-A survey," Automatica, vol. 20, pp. 387-404, Dec. 1984. [15] H. Kwakernaak, "Robust control and #°°-optimization-Tutorial paper," Automat-ica, vol. 29, pp. 255-273, Apr. 1993. [16] J. Auzins and R. V . Wilhelm, "Automotive electronics: Steering toward the next frontiers," IEEE Potentials, vol. 13, pp. 32-36, Oct. 1994. [17] R. Isermann and B. Freyermuth, "Process fault diagnosis based on process model knowledge-Part I: principles for fault diagnosis with parameter estimation," Trans. ASME, J. Dynamic Systems, Measurement and Control, vol. 113, pp. 620-626, Dec. 1991. [18] R. J. Patton and J. Chen, "Detection of faulty sensors in aero jet engine systems us-ing robust model-based methods," in IEE Colloquim No. 156: Condition Monitoring for Fault Diagnosis, pp. 2/1-2/22, IEE, 1991. [19] R. J. Patton, "Fault detection and diagnosis in aerospace systems using analytical redundancy," Computing and Control Engineering Journal, vol. 2, pp. 127-135, May 1991. [20] R. J. Patton, J. Chen, and H. Y . Zhang, "Modeling methods for improving robust-ness in fault diagnosis of jet engine system," in Proc. of the 31st Annual Conference on Decision and Control, pp. 2330-2335, IEEE, 1992. [21] M . N . Brown, R. W. Stewart, T. S. Durrani, and T. W. Buggy, "DSP subsystem for knowledge based health monitoring of gas turbine engines," in ICASSP'92: Acous-tics, Speech & Signal Processing Conference, pp. 69-72, IEEE, 1992. ( [22] G. Rizzoni and P. S. Min, "Detection of sensor failures in automotive engines," IEEE Trans. Veh. Technol, vol. 40, pp. 487-500, May 1991. [23] L. Dubost and J. N . Heude, "On-line diagnosis of a superchargher in a noisy en-vironment," in Proc. of the IEE Conference on Intelligent Systems Engineering, pp. 65-70, 1992. [24] K . A . Marko, B. Bryant, and N . Soderborg, "Neural network application to compre-hensive engine diagnostics," in Proc. of the International Conference on Systems, Man, and Cybernatics, pp. 1016-1022, IEEE, 1992. Bibliography 72 [25] K . S. Narendra and K . Parthasarathy, "Identification and control of dynamical sys-tems using neural networks," IEEE Trans. Neural Networks, vol. 1, pp. 4-27, March 1990. [26] R. J. Patton and J. Chen, "Advances in fault diagnosis using analytical redundancy," in IEE Colloquim No. 019: Plant Optimisation for Profit, pp. 6/1-6/12, IEE, 1993. [27] P. Eykhoff, System Identification: Parameter and State Estimation. John Wiley & Sons, 1974. [28] C. F. Taylor, The Internal-Combustion Engine in Theory and Practice, Volume IT. Combustion, Fuels, Materials, Design. The Technology Press of the MIT and John Willey & Sons, 1968. [29] J. S. Rao, Advanced Theory of Vibration. John Wiley Sz Sons, 1st ed., 1992. [30] P. E. Wellstead, "Application of adaptive techniques to internal combustion en-gine control," in The benefits of electronic control systems for internal combustion engines—Seminar Proceedings, pp. 11-23, Inst. Mech. Engineers, Jan. 1989. [31] S. Roy, O. P. Malik, and G. S. Hope, "Adaptive control of speed and equivalence ratio dynamics of a diesel driven power plant," IEEE Trans. Energy Conversion, vol. 8, pp. 13-19, March 1993. [32] S. Roy, O. P. Malik, and G. S. Hope, "Real time test results with adaptive speed con-trollers for a diesel prime mover," IEEE Trans. Energy Conversion, vol. 8, pp. 499-505, Sep. 1993. [33] J. Jiang, "Design and implementation of an optimal gain scheduling controller for a diesel engine," in Proc. of the Second IEEE Conference on Control Applications, pp. 667-672, 1993. [34] T. Brown and W. Neill, "Determination of engine cylinder pressures from crankshaft speed fluctuations," SAE Trans., J. Engines, vol. 101, pp. 10-19, 1992. Section 3, Paper No. 920463. [35] R. Isermann and B. Freyermuth, "Process fault diagnosis based on process model knowledge-Part II: case study experiments," Trans. ASME, J. Dynamic Systems, Measurement and Control, vol. 113, pp. 627-633, Dec. 1991. [36] T. H. B. Jewitt and B. Lawton, "The use of speed sensing for monitoring the condi-tion of military vehicles engines," in International Conference on Vehicle Condition Monitoring and Fault Diagnosis, pp. 67-72, Inst. Mech. Engineers, 1985. Bibliography 73 [37] S. J. Citron, J. E. O'Higgins, and L. Y Chen, "Cylinder by cylinder engine pres-sure and pressure torque waveform determination utilizing speed fluctuations," SAE Trans., J. Engines, vol. 98, pp. 933-947, 1989. Section 3, Paper No. 890486. [38] M . T. K Srinivasan, G Rizzoni and G. C. Luh, "On-line estimation of net engine torque from crankshaft angular velocity measurement using repetitive estimators," in Proc. of 1992 American Control Conference, pp. 516-520, 1992. [39] W. B. Ribbens and G. Rizzoni, "Applications of precise crankshaft position measure-ments for engine testing, control and diagnosis," SAE Trans., J. Engines, vol. 98, pp. 1582-1596, 1989. Section 3, Paper No. 890885. [40] J. W. Freestone and E. G. Jenkins, "The diagnosis of cylinder power faults in diesel engines by flywheel speed measurement," in International Conference on Vehicle Condition Monitoring and Fault Diagnosis, pp. 15-24, Inst. Mech. Engineers, 1985. [41] G. Rizzoni, "Diagnosis of individual cylinder misfires by signature analysis of crankshaft speed fluctuations," SAE Trans., J. Engines, vol. 98, pp. 1572-1581, 1989. Section 3, Paper No. 890884. [42] Y . K . Chin and F. E. Coats, "Engine dynamics: Time-based versus crank-angle based," SAE Trans., J. Engines, vol. 95, pp. 937-956, 1986. Section 3, Paper No. 860412, [43] J. B. Heywood, Internal Combustion Engine Fundamentals. McGraw Hil l , 1st ed., 1988. [44] G. Rizzoni, "A stochastic model for the indicated pressure process and the dynamics of the internal combustion engine," IEEE Trans. Veh. Technol., vol. 38, pp. 180-192, Aug. 1989. [45] F. T. Connolly and A. E. Yagle, "Modeling and identification of the combustion process in internal combustion engines: II—Experimental results," Trans. ASME, J. of Engineering for Gas Turbines and Power, vol. 115, pp. 801-809, Oct. 1993. [46] Y . Z. Tsypkin, Relay Control systems. Cambridge University Press, 1984. [47] M . Schwartz and L. Shaw, Signal Processing: Discrete Spectral Analysis, Detection, and Estimation. McGraw-Hill, 1st ed., 1975. A p p e n d i x A Test C e l l Descr ipt ion A n I B M - P C computer equipped with analog to digital (A/D) conversion board and a pulse-timer board is used for high speed data acquisition. In order to reduce equipment cost, only two cylinder pressure signals are acquired (one "healthy" cylinder and one under-fueled cylinder). Data from 2 channels (2 cylinder pressures) are sampled simul-taneously by the sample and hold preamplifiers; a multiplexer directs each sample to the A / D converter. The A / D converter is triggered by an index signal, bottom dead centre (BDC) of the first left (1L) cylinder, from an optical shaft encoder connected to the crankshaft. The optical shaft encoder is used for measuring the crank angle. The resolu-tion is 360 pulses per revolution. The flywheel proximity sensor is a sensing-transduction coil which uses the magnetic pick-up principle. The signal conditioning system includes charge amplifiers with very high input impedance suitable for use with piezoelectric cylinder pressure gauges. The signal from the crankshaft/flywheel proximity sensor is passed through a zero-crossing circuit to produce digital pulses which are then fed to the pulse-timer input of the data acquisition system. The variations of angular velocity are measured indirectly by recording the number of clock cycles (absolute time) between two successive crank angle sampling intervals. Each measurement is repeated for different brake (dynamometer load) torque values and with several levels of fuel starvation. The proposed data acquisition system is very similar to the one used by the Institute for Machinery Research (IMR). 74 Appendix A. Test Cell Description 75 DDC6V92TA proximity signal sensor conditioning TJ C + LH flywheel S T 1 charge amps. * • IP Pulse/Timer • T PC AT A/D • „ CLX * T Figure A . l : Illustrates a data acquisition system for measuring cylinder pressures and crankshaft/flywheel speed variations. Appendix B Program Listings cl.m function y = c_l(theta,r,l) M '/.FUNCTION '/.ARGUMENTS : 7. 7. 7.RETURNS : x=r/l; fll=sin(theta)."2; fl2=sqrt(ones(size(theta))-(x~2)*f11); y=sin(theta)+x*sin(theta).*cos(theta)./fl2; C_l r 1 theta y = crankshaft radius (m); = connecting rod length (m); = crank angle (rad); = the ratio v/(r.w). c2.m function y = c_2(theta,r,l) 7.FUNCTI0N C_2.M 7.ARGUMENTS: r 7. 1 7. theta 7.RETURNS: y x=r/l; calpha=sqrt(ones(size(theta))-(x~2)*sin(theta).~2); xx=x*cos(theta)./calpha; y=(ones(size(theta))+xx).*(cos(theta)+... x.*(xx-ones(size(theta))).*(sin(theta)."2)./calpha); = crankshaft radius (m); = connecting rod length (m); = crank angle (rad); = acceleration coefficient. c3.m 3(theta,r,l) C_3.M r = crankshaft radius; 1 = connecting rod length; theta = crank angle; y = c3 acceleration coefficient. function y = c 7.FUNCTI0N 7.ARGUMENTS: 7. 7. 7.RETURNS: x=r/l; calpha=sqrt(ones(size(theta)) (x~2)*sin(theta)."2); y=x*cos(theta)./calpha; 76 Appendix B. Program Listings 77 c4.m function y = c_4(theta,r,1) '/.FUNCTION C_4.M '/.ARGUMENTS: r '/. 1 '/„ theta "/.RETURNS: y x=r/l; calpha=sqrt(ones(size(theta))-(x~2)*sin(theta).~2); xx=x*cos(theta)./calpha; y=x*(sin(theta)).*(xx.~2- . . . ones(size(theta)))./(calpha); crankshaft radius; connecting rod length; crank angle; c4 acceleration coefficient. cmpal l .m '/.CMPALL.M '/.Calls: demol.m, rmse.m, rot.m dispC NORMAL OPERATING CONDITION') wOhat=demol(0,wO) ; dispCROOT MEAN SQUARE velocity error (rad/sec):') eO=rmse(wO-0.8325*w0hat) dispC RESIDUAL TO SIGNAL RATIO ('/„):') RSRO=eO/rmse(wO)*100 figure(1) plot(ca0*180/pi,w0,ca0*180/pi,.8325*w0hat,':') xlabel('degrees') ylabel('rad/sec') title('Estimated (..) and Actual (—) Velocity Fluctuation') dispCPress any key to continue...') pause dispC TESTING 6 ACTUATOR FAULTS') disp('The actuator fault is defined as') dispClO 1/. under-fueling of each cylinder at a time.') wlhat=demol(1,wll); w2hat=demol(2,w21); w3hat=demol(3,wlr) ; w4hat=demol(4,wll); w5hat=demol(5,w3r); w6hat=demol(6,w31); wlx = wlhat + rot(wlhat,59); w2x = w2hat + rot(w2hat,59); w3x = w3hat + rot(w3hat,59); w4x = w4hat + rot(w4hat,59); w5x = w5hat + rot(w5hat,59); w6x = w6hat + rot(w6hat,59) ; wllx=wll+rot(wll,59) ; w3rx=w3r+rot(w3r,59) ; w31x=w31+rot(w31,59); w2rx=w2r+rot(w2r,59); w21x=w21+rot(w21,59); Appendix B. Program Listings w l r x = w l r + r o t ( w l r , 5 9 ) ; w l x=w lx -mean (w lx ) ; w2x=w2x-mean(w2x); w3x=w3x-mean(w3x); w4x=w4x-mean(w4x); w5x=w5x-mean(w5x); w6x=w6x-mean(w6x); w l l x = w l l x - m e a n ( w l l x ) ; w21x=w21x-mean(w21x); w31x=w31x-mean(w31x); w l r x = w l r x - m e a n ( w l r x ) ; w2rx=w2rx-mean(w2rx ) ; w3rx=w3rx-mean(w3rx ) ; d i s p ( ' E l i m i n a t i o n o f t h e S i n u s o i d a l T r e n d ' ) d i s p C T h e RMS v e l o c i t y e r r o r s ( r ad/ sec ) a n d ' ) d i s p C t h e c o r r e s p o n d i n g RSRs (%)') d i s p C a r e d i s p l a y e d i n t h e c y l i n d e r o r d e r : ' ) e l = r m s e ( w l l x - . 6 6 4 7 * w l x ) R S R l = e l / r m s e ( w l l x ) * 1 0 0 e2=rmse (w3rx - .6154*w2x ) RSR2=e2/rmse (w3rx ) *100 e3=rmse (w31x- .5599*w3x) RSR3=e3/rmse(w31x)*100 e4=rmse (w2rx - .5897*w4x ) RSR4=e4/rmse (w2rx ) *100 e5=rmse (w21x- .5820*w5x) RSR5=e5/rmse(w21x)*100 e6= rmse (w l r x - . 7686*w6x ) RSR6=e6/rmse (w l rx ) *100 d i s p C T h e a ve r a ge RMS e r r o r ( r ad/ sec ) i s : ' ) e_m=(el+e2+e3+e4+e5+e6)/6 d i s p C T h e a ve r a ge RSR (*/.) i s : ' ) RSR_m=(RSR1+RSR2+RSR3+RSR4+RSR5+RSR6)/6 f i g u r e ( 2 ) s u b p l o t ( 3 , 2 , 1 ) p l o t ( c a 0 * 1 8 0 / p i , w l l x , c a 0 * 1 8 0 / p i , . 6 * w l x , ' : ' ) x l a b e l ( ' d e g r e e s ' ) y l a b e l ( ' r a d / s e c ' ) s u b p l o t ( 3 , 2 , 2 ) p l o t ( c a 0 * 1 8 0 / p i , w 3 r x , c a 0 * 1 8 0 / p i , . 6 * w 2 x , ' : ' ) x l a b e l ( ' d e g r e e s ' ) y l a b e l ( ' r a d / s e c ' ) s u b p l o t ( 3 , 2 , 3 ) p l o t ( c a 0 * 1 8 0 / p i , w 3 1 x , c a 0 * 1 8 0 / p i , . 6 * w 3 x , ' : ' ) x l a b e l ( ' d e g r e e s ' ) y l a b e l ( ' r a d / s e c ' ) s u b p l o t ( 3 , 2 , 4 ) p l o t ( c a 0 * 1 8 0 / p i , w 2 r x , c a 0 * 1 8 0 / p i , . 6 * w 4 x , ' : ' ) x l a b e l ( ' d e g r e e s ' ) y l a b e l ( ' r a d / s e c ' ) s u b p l o t ( 3 , 2 , 5 ) p l o t ( c a 0 * 1 8 0 / p i , w 2 1 x , c a 0 * 1 8 0 / p i , . 6 * w 5 x , ' : ' ) Appendix B. Program Listings xlabel('degrees') yl a b e l ( ' r a d / s e c ' ) subplot(3,2,6) plot(ca0*180/pi,wlrx,ca0*180/pi,.6*w6x,':') xlabel('degrees') yl a b e l ( ' r a d / s e c ' ) disp('Press any key to continue... ') pause demol.m f u n c t i o n dwhat = demol(n,dw) "/.FUNCTION DEM01.M "/.ARGUMENTS: n 7. 7. dw 7.RETURNS: dwhat 7.Calls: geom.m, ldown.m gl o b a l ca caO xO w_0 [pi p2 p3 p4 p5 p6] = ldown(n); [u,q,m] = geom(ca,pl,p2,p3,p4,p5,p6); uu = (x0*q + u); d = -mean(uu); dT = c a ( 2 ) - c a ( l ) ; dxhat = dw(l)*(2*w_0+dw(D) +. . . i g r ( ( u u + d)./m,dT); dwhat = sqrt(xO + dxhat) - w_0; dwhat=dwhat-mean(dwhat); dwhat=rot(interpl(ca,dwhat,caO),1); integer corresponding to the down c y l i n d e r ; a c t u a l angular v e l o c i t y ; estimated angular v e l o c i t y . demo2.m f u n c t i o n [JEhatOl,JEhat02,JEhat03,TLhat01,... TLhat02,TLhat03,dwhat01,dwhat02, . . . dwhat03,rms,rsr]=demo2(n,dw,PARO) 7.FUNCTI0N DEM02.M 7.ARGUMENTS: n = c y l i n d e r number; %' dw = angular v e l o c i t y waveform; 7. PARO = parameter i n i t i a l i z a t i o n ; 7.RETURNS: JEhatOl = estimated engine i n e r t i a using 7. standard least-squares; 7. JEhat02 = estimated engine i n e r t i a using 7. the re c u r s i v e gradient estimator; 7, JEhat03 = estimated engine i n e r t i a using 7. least-squares with exponential f o r g e t t i n g 7. TLhatOl = estimated torque f l u c t u a t i o n using SLS; 7. TLhat02 = estimated torque f l u c t u a t i o n using GE; 7. TLhat03 = estimated torque f l u c t u a t i o n using LSEF; 7. dwhatOl = estimated speed f l u c t u a t i o n using SLS; 7, dwhat02 = estimated speed f l u c t u a t i o n using GE; 7. dwhat03 = estimated speed f l u c t u a t i o n using LSEF; Appendix B. Program Listings 80 % rms = vector of root mean square velocity errors; 7. rsr = vector of residual to signal ratios. '/.Calls: ldown.m, diffc.m, sls.m, itgr.m, est.m, itls.m, rmse.m. p0=0.3; P0=diag([l 1 1 ] ) ; lambda=15; %crank angle for pressure sampling ca = [0:719]' /720 * 2 * p i ; '/.crank angle for velocity sampling caO = [0:117]'/118 * 2 * p i ; [pl,p2,p3,p4,p5,p6]=ldown(n); pl=interpl(ca,pl,caO); p2=interpl(ca,p2,caO); p3=interpl(ca,p3,caO) ; p4=interpl(ca,p4,caO) ; p5=interpl(ca,p5,caO) ; p6=interpl(ca,p6,caO) ; dca0=ca0(2)-ca0(l) ; dwp=diffc(dw,dcaO); [Y01,PHI01,PAR01,ERR01,JEhatOl,TLhatOl,epOl,dwhatOl]=... sIs(caO,p1,p2,p3,p4,p5,p6,dw,dwp); [ERR02,PAR02]=itgr(Y01,PHI01,pO,PARO); [JEhat02,TLhat02,ep02,dwhat02]=... est(PAR02,caO,pi,p2,p3,p4,p5,p6,dw,dwp); [ERR03,PAR03]=itls(Y01,PHI01,lambda,PO,PARO); [JEhat03,TLhat03,ep03,dwhat03]=... est(PAR03,caO,pl,p2,p3,p4,p5,p6,dw,dwp); rms=[rmse(epOl) rmse(ep02) rmse(ep03)]; rsr=[rmse(epOl) rmse(ep02) rmse(ep03)]/rmse(dw)*100; diffc.m function y = diffc(x,dT) '/.FUNCTION DIFFC.M '/.ARGUMENTS: x 7. dT '/.RETURNS: y ll=length(x); for i=3:ll-2 y(i)=(x(i+2)/4 + x(i+l)/2-... x(i-l)/2-x(i-2)/4)/2; end y(l)=(x(3)/4+x(2)/2-x(ll)/2 x ( l l - l ) / 4 ) / 2 ; y(2)=(x(4)/4+x(3)/2-x(l)/2-... x(ll)/4)/2; y(ll)=(x(2)/4+x(l)/2-x(ll-l)/2 -... x(ll-2)/4)/2; y(ll-l)=(x(l)/4+x(ll)/2-x(ll-2)/2 -x(ll-3)/4)/2; y=yVdT; = vector of samples; = crank angle interval; = smooth derivative. Appendix B. Program Listings d l t .m f u n c t i o n y=dlt(alpha,beta,gamma,x) '/.FUNCTION 7.ARGUMRNTS: 7. 7. 7. 7.RETURNS: DLT.M alpha beta gamma x y y=beta./(beta~2*sin(x/2).~2+l); y=y*alpha-gamma; = impulse gain; = impulse width; = o f f s e t ; = crank angle vector; = p e r i o d i c impulse. est.m f u n c t i o n [JEhat,TLhat,ep,dwhat]=... est(PAR,caO,pl,p2,p3,p4,p5,p6,dw,dwp) 7.FUNCTI0N EST.M 7.ARGUMENTS: PAR = vector of estimated parameters 7. caO = crank angle vector; 7. p i = c y l i n d e r pressure, i = l , . . . , 6; 7. dw = angular v e l o c i t y waveform; 7, dwp = angular a c c e l e r a t i o n waveform; 7oRETURNS: JEhat = vector of estimated i n e r t i a s ; 7. TLhat = estimated torque f l u c t u a t i o n ; 7. ep = v e l o c i t y r e s i d u a l s ; 7. dwhat = estiamted angular v e l o c i t y ; 3.m, c_4.m, igr.m. dwhat 7.Calls: c_l.m, c_2.m, 7.Engine data correspond to DDC 6V 92TA: w_0=125.5; 7.engine speed (rad/sec) b_c=1.24; 7.damping (Nm/(rad/sec)) m_p=6.03; 7.piston mass (kg) r=0.0635; 7.crankshaft radius (m) 1=0.2571; 7oconnecting rod length (m) Ap=0.01188; 7.piston area (m"2) J_C=0.1544; 7.crankshaft i n e r t i a (kg m"2) J_CR=0.0745; 7oconnecting rod i n e r t i a (kg m~2) J_F=3.85; 7.flywheel i n e r t i a (kg nT2) x0=w_0~2; 7.Pressure inputs: p l p2 p3 p4 p5 p6 7.Firing order ( r i g h t hand ro t a t i o n ) : 7.1L 3R 3L 2R 2L 1R 7.1 2 3 4 5 6 dca0=ca0(2)-ca0(l); cal=ca0; c _ l l = c _ l ( c a l , r , l ) ; c_21=c_2(cal,r,l); c_31=c_3(cal,r,l); c_41=c_4(cal,r,l); ca2=cal-(60-3.5)*pi/180; c_12=c_l(ca2,r,l); c_22=c_2(ca2,r,l); Appendix B. Program Listings c _ 3 2 = c _ 3 ( c a 2 , r , l ) ; c _ 4 2 = c _ 4 ( c a 2 , r , l ) ; c a 3 = c a 2 - ( 6 0 + 3 . 5 ) * p i / 1 8 0 ; c _ 1 3 = c _ l ( c a 3 , r ) l ) ; c _ 2 3 = c _ 2 ( c a 3 , r , l ) ; c _ 3 3 = c _ 3 ( c a 3 , r , l ) ; c _ 4 3 = c _ 4 ( c a 3 , r , l ) ; c a 4 = c a 3 - ( 6 0 - 3 . 5 ) * p i / 1 8 0 ; c _ 1 4 = c _ l ( c a 4 , r , l ) ; c _ 2 4 = c _ 2 ( c a 4 , r , l ) ; c _ 3 4 = c _ 3 ( c a 4 , r , l ) ; c _ 4 4 = c _ 4 ( c a 4 , r , l ) ; c a 5 = c a 4 - ( 6 0 + 3 . 5 ) * p i / 1 8 0 ; c _ 1 5 = c _ l ( c a 5 , r , l ) ; c _ 2 5 = c _ 2 ( c a 5 , r , l ) ; c _ 3 5 = c _ 3 ( c a 5 , r , l ) ; c _ 4 5 = c _ 4 ( c a 5 , r , l ) ; c a 6 = c a 5 - ( 6 0 - 3 . 5 ) * p i / 1 8 0 ; c _ 1 6 = c _ l ( c a 6 , r , l ) ; c _ 2 6 = c _ 2 ( c a 6 , r , l ) ; c _ 3 6 = c _ 3 ( c a 6 , r , l ) ; c _ 4 6 = c _ 4 ( c a 6 , r , l ) ; s l l = c _ l l . * c _ l l + c _ 1 2 . * c _ 1 2 + c _ 1 3 . * c _ 1 3 + . . . c _ 1 4 . * c _ 1 4 + c _ 1 5 . * c _ 1 5 + c _ 1 6 . * c _ 1 6 ; s l 2 = c _ l l . * c _ 2 l + c _ 1 2 . * c _ 2 2 + c _ 1 3 . * c _ 2 3 + . . . c _ 1 4 . * c _ 2 4 + c _ 1 5 . * c _ 2 5 + c _ 1 6 . * c _ 2 6 ; s 3 3 = c _ 3 1 . * c _ 3 1 + c _ 3 2 . * c _ 3 2 + c _ 3 3 . * c _ 3 3 + . . . c _ 3 4 . * c _ 3 4 + c _ 3 5 . * c _ 3 5 + c _ 3 6 . * c _ 3 6 ; s 3 4 = c _ 3 1 . * c _ 4 1 + c _ 3 2 . * c _ 4 2 + c _ 3 3 . * c _ 4 3 + . . . c _ 3 4 . * c _ 4 4 + c _ 3 5 . * c _ 4 5 + c _ 3 6 . * c _ 4 6 ; s p = ( l e + 6 ) * ( c _ l l . * p l + c _ 1 2 . * p 2 + c _ 1 3 . * p 3 + . . . c _ 1 4 . * p 4 + c _ 1 5 . * p 5 + c _ 1 6 . * p 6 ) ; dx=2*w_0*dw; dxp=2*w_0*dwp; p h i l = 0 . 5 * d x p ; f a c t = m a x ( p h i l ) ; l h a t = l e n g t h ( P A R ) ; t h = ( 0 : l h a t - l ) / l h a t * 2 * p i ; J E h a t = P A R ( l , : ) . / f a c t ; T L h a t = P A R ( 2 , : ) . * s i n ( t h ) + P A R ( 3 , : ) . * c o s ( t h ) ; J E = P A R ( 1 , l h a t ) . / f a c t ; T L = P A R ( 2 , l h a t ) . * s i n ( c a 0 ) + P A R ( 3 , l h a t ) . * c o s ( c a 0 ) y 3 = A p * r * s p ; m= ( J E + m _ p * r ' " 2 * s l l - J _ C R * s 3 3 ) / 2 ; i n p = y 3 - T L ; °/,TLhat ( 1 : l h a t - 1 ) ' ; i n p = i n p - m e a n ( i n p ) ; d x h a t = d w ( l ) * ( 2 * w _ 0 + d w ( l ) ) + i g r ( i n p . / m . d c a O ) dwhat = s q r t ( x O +• d x h a t ) - w_0 ; d w h a t = d w h a t - m e a n ( d w h a t ) ; °/.dwhat=0. 9 * r o t ( dwha t , 1 ) ; ep=dw-dwhat ; Appendix B. Program Listings 83 geom.m f u n c t i o n [ u , q , m ] = g e o m ( c a . p l , p 2 , p 3 , p 4 , p 5 , p6) "/.FUNCTION: GEOM.M "/.ARGUMENTS: c a = c r a n k a n g l e ( r a d ) ; °/„ p i = p r e s s u r e i n t h e i - t h c y l i n d e r (MPa) ; "/.RETURNS: [u ,q ,m] = mode l c o e f f i c i e n t s . " / . C a l l s : c _ l . m , c _ 2 . m , c _ 3 . m , c _ 4 . m "/.Engine d a t a c o r r e s p o n d t o DDC 6V 92TA: °/.w_0=125.5; "/.nominal e n g i n e s p e e d ( r a d / s e c ) b _ c = 1 . 2 4 ; "/.damping (Nm/ ( rad/sec ) ) ' m_p=6 .03 ; "/ .p is ton mass (kg ) r = 0 . 0 6 3 5 ; "/ . c ranksha f t r a d i u s (m) 1 = 0 . 2 5 7 1 ; "/ .connec t ing r o d l e n g t h (m) A p = 0 . 0 1 1 8 8 ; "/ .p iston a r e a (m~2) J _ C = 0 . 1 5 4 4 ; "/ . c ranksha f t i n e r t i a (kg m~2) J _ C R = 0 . 0 7 4 5 ; "/ .connect ing r o d i n e r t i a ( k g m~2) J _ F = 3 . 8 5 ; "/ . f lywheel i n e r t i a ( kg m~2) "/.Pressure i n p u t s : p l p2 p3 p4 p5 p6 " / . F i r i n g o r d e r ( r i g h t hand r o t a t i o n ) °/.lL 3R 3L 2R 2L 1R 7.1 2 3 4 5 6 c a l = c a ; c _ l l = c _ l ^ c a l , r , l ) c _21=c_2 ^ c a l , r , l ) c _31=c_3 ^ c a l , r , l ) c _41=c_4 ' c a l , r , l ) c a 2 = c a l - : 6 0 - 3 . 5 ) * p i / 1 8 0 ; c _ 1 2 = c _ l ; c a 2 , r , l ) c _22=c_2 ; c a 2 , r ) l ) c _32=c_3 ; c a 2 , r , l ) c _42=c_4 ^ c a 2 , r , 1 ) c a 3 = c a 2 - ; 6 0 + 3 . 5 ) * p i / 1 8 0 ; c _ 1 3 = c _ l ^ c a 3 , r , l ) c _23=c_2 ^ c a 3 , r , 1 ) c_33=c_3 ' c a 3 , r , 1 ) c_43=c_4 ' c a 3 , r , 1 ) c a 4 = c a 3 - : 6 0 - 3 . 5 ) * p i / 1 8 0 ; c 14=c 1 ' c a 4 , r , l ) c _24=c_2 ' c a 4 , r , l ) c _34=c_3 ^ c a 4 , r , l ) c _44=c_4 ' c a 4 , r , l ) c a 5 = c a 4 - : 6 0 + 3 . 5 ) * p i / 1 8 0 ; c _ 1 5 = c _ l ' c a 5 , r , l ) c _25=c_2 ' c a 5 , r , D c_35=c_3 ( c a 5 , r , l ) c _45=c_4 ^ c a 5 , r , l ) c a 6 = c a 5 - : 6 0 - 3 . 5 ) * p i / 1 8 0 ; c _ 1 6 = c _ l [ c a 6 , r , D c_26=c_2 ( c a 6 , r , 1 ) c_36=c_3 ( c a 6 , r , l ) Appendix B. Program Listings 84 c_46=c_4(ca6,r,l); s l l = c _ l l . * c _ l l + c_12.*c_12 + c_13.*c_13 +. c_14.*c_14 + c_15.*c_15 + c_16.*c_16; sl2=c _ l l . * c _ 2 1 + c_12.*c_22 + c_13.*c_23 +. c_14.*c_24 + c_15.*c_25 + c_16.*c_26; s33=c_31.*c_31 + c_32.*c_32 + c_33.*c_33 +. c_34.*c_34 + c_35.*c_35 + c_36.*c_36; s34=c_31.*c_41 + c_32.*c_42 + c_33.*c_43 +. c_34.*c_44 + c_35.*c_45 + c_36.*c_46; s p = ( l e + 6 ) * ( c _ l l . * p l + c_12.*p2 + c_13.*p3 + c_14.*p4 + c_15.*p5 + c_16.*p6); q = J_CR*s34-m_p*r"2*sl2; u = Ap*r*sp; m = 0.5*(J_C+J_F-J_CR*s33+m_p*r~2*sll); igr .m f u n c t i o n II=igr(y,dT) "/.FUNCTION IGR.M "/.ARGUMENTS: y = vector of samples; "/. dT = crank angle i n t e r v a l ; "/.RETURNS: II = numerical i n t e g r a t i o n . I K D = 0 ; cs=cumsum(y)-y(1); f o r i=2:length(y) I L ( i ) = y ( D / 2 + c s ( i - l ) + y ( i ) / 2 ; end II=dT*II'; i tgr .m f u n c t i o n [err,PAR]=itgr(y,phi,pO,PARO) "/.FUNCTION "/.ARGUMENTS: ITGR.M y = output vector; 7. phi = s i g n a l matrix; 7. pO = gain f o r the gradient estimator 7. PARO = i n i t i a l i z a t i o n ; "/.RETURNS: err = vector of r e s i d u a l s ; 7. PAR = estiamted parameters. Nit=length(y); h=2*pi/118; PARC:,1)=PAR0; f o r i = l : N i t e r r ( i ) = p h i ( i , : ) * P A R ( : , i ) - y ( i ) ; E ( i ) = 0 . 5 * e r r ( i ) ' * e r r ( i ) ; PARC:,i+l)=PAR(:,i)-h*p0*phi(i,:)'*err(i) ; end err=err'; Appendix B. Program Listings i t l s .m f u n c t i o n [err,PAR]=itIs(y,phi,lambda,PO,PARO) '/.FUNCTION ITLS.M '/.ARGUMENTS: y = output vector; 7. phi = s i g n a l matrix; 7. lambda = f o r g e t t i n g f a c t o r ; 7. PO = gain i n i t : L a l i z a t i o n ; 7. PARO = parameter i n i t i a l i z a t i o n '/.RETURNS: err = vector of re s i d u a l s ; 7. PAR = estimated parameters. Nit=length(y); h=2*pi/118; P=P0; PARC:,1)=PAR0; f o r i = l : N i t e r r ( i ) = p h i ( i , : ) * P A R ( : , i ) - y ( i ) ; PAR(:,i+l)=PARC:,i)-h*P*phi(i,:)'*err(i); lambda=lambda*Cl-norm(P)/70); p=p-h*P*(lambda+phi(i,:)'*phi Ci , : ) * P ) ; end err=err'; ldown.m f u n c t i o n [ y l , y2, y3, y4, y5, y6] = ldown(n) g l o b a l p l c o l p5col presO '/.FUNCTION: LDOWN.M '/.ARGUMENTS: n = integer; '/.RETURNS: [Yl, Y2, Y3, Y4, Y5, Y6] = s i x c y l i n d e r 7, pressure waveforms, '/.Load 5 pressures from healthy c y l i n d e r s '/.originating from PRES.DAT, f i r s t 6*720 '/.values and a s i n g l e 'down' c y l i n d e r corresponding °/,to number n. If n i s not i n the range 0 through 6 '/.then a l l healthy c y l s are loaded. If n i s zero then '/.the presO waveforms are loaded. '/.Note v a r i a b l e p l c o l should be loaded (healthy c y l s ) , '/.and p5col (down c y l s ) , and presO (baseline) . '/.Calls: separate.m [yl y2 y3 y4 y5 y6] = s e p a r a t e ( p l c o l ) ; [dl d2 d3 d4 d5 d6] = separate(p5col); i f n==0 [yl y2 y3 y4 y5 y6] = separate(presO); e l s e i f n==l y l = d l ; e l s e i f n==2 y2 = d2; e l s e i f n==3 y3 = d3; e l s e i f n==4 y4 = d4; e l s e i f n==5 Appendix B. Program Listings 86 y5 = d5; e l s e i f n==6 y6 = d6; end pia l l .m 7.PIALL. M 7 0Calls: demo2.m, rmse.m input('Enter option number: ' ) ; n=ans; 7 0 i n i t i a l i z a t i o n f o r the normal condition PARO [970 -50 -30] ' ; 7 0 i n i t i a l i z a t i o n f o r 1L down PAR011=[1510 -580 204] '; 7 0 i n i t i a l i z a t i o n f o r 3R down PAR03r=[1380 -350 400]'; 7 0 i n i t i a l i z a t i o n f o r 3L down PAR031=[1650 450 250]'; 7 0 i n i t i a l i z a t i o n f o r 2R down PAR02r=[1780 400 -615]'; 7 0 i n i t i a l i z a t i o n f o r 2L down PAR021=[1700 -170 -790]'; 7 o i n i t i a l i z a t i o n f o r 1R down PAR01r=[1730 -510 -480]'; i f n==0 PARO = PARO; dw=1.05*rot(w0,117); disp('NORMAL OPERATING CONDITION') e l s e i f n==l PARO = PAR011; dw=wll; disp('FAULT 1L') e l s e i f n==2 PARO = PAR03r; dw=w3r; d i s p C FAULT 3R') e l s e i f n==3 PARO = PAR031; dw=w31; d i s p C FAULT 3L') e l s e i f n==4 PARO = PAR02r; dw=w2r; d i s p C FAULT 2R') e l s e i f n==5 PARO = PAR021; dw=w21; d i s p C FAULT 2L') e l s e i f n==6 PARO = PAROlr; dw=wlr; disp('FAULT 1R') e l s e i f n==10 PARO = PARO; dw=wO; n=l; Appendix B. Program Listings disp('FAULT IL vs. e l s e i f n==20 PARO = PARO; dw=wO; n=2; disp('FAULT 3R vs. e l s e i f n==30 PARO = PARO; dw=wO; n=3; disp('FAULT 3L vs. e l s e i f n==40 PARO = PARO; dw=wO; n=4; disp('FAULT 2R vs. e l s e i f n==50 PARO = PARO; dw=wO; n=5; disp('FAULT 2L vs. e l s e i f n==60 PARO = PARO; dw=wO; n=6; disp('FAULT 1R vs. e l s e i f n==21 PARO = PAROll; dw=wll; n=2; disp('FAULT 3R vs. e l s e i f n==31 PARO = PAROll; dw=wll; n=3; disp('FAULT 3L vs. e l s e i f n==41 PARO = PAROll; dw=wll; n=4; disp('FAULT 2R vs. e l s e i f n==51 PARO = PAROll; dw=wll; n=5; disp('FAULT 2L vs. e l s e i f n==61 PARO = PAROll; dw=wll; n=6; disp('FAULT 1R vs. end "/gradient step p0=0.3; "/ogain matrix P0=diag([l 1 1 ] ) ; "/oforgetting f a c t o r lambda=15; NORMAL OPERATING CONDITION') NORMAL OPERATING CONDITION') NORMAL OPERATING CONDITION') NORMAL OPERATING CONDITION') NORMAL OPERATING CONDITION') NORMAL OPERATING CONDITION') FAULT IL') FAULT IL') FAULT IL') FAULT IL') FAULT IL') Appendix B. Program Listings °/0crank angle f o r pressure sampling ca = [0:719]' /720 * 2 * p i ; °/0crank angle f o r v e l o c i t y sampling caO = [0:117] '/118 * 2 * p i ; [pl,p2,p3,p4,p5,p6]=ldown(n); p l = i n t e r p l ( c a . p l , c a O ) ; p2=interpl(ca,p2,ca0); p3=interpl(ca,p3,ca0); p4=interpl(ca,p4,ca0); p5=interpl(ca,p5,caO); p6=interpl(ca,p6,caO); dca0=ca0(2)-ca0(l); dwp=diffc(dw.dcaO); [Y01,PHI01,PAR01,ERR01,JEhatOl,TLhatOl,epOl,dwhatOl] sIs(ca0,pl,p2,p3,p4,p5,p6,dw )dwp); [ERR02,PAR02]=itgr(Y01,PHIO1,pO,PARO); [JEhat02,TLhat02,ep02,dwhat02]=... est(PAR02,caO,pi,p2,p3,p4,p5,p6,dw,dwp); [ERR03,PAR03]=itls(Y01,PHI01,lambda,PO,PARO); [JEhat03,TLhat03,ep03,dwhat03]=... est(PAR03,ca0.pl,p2,p3,p4,p5,p6,dw,dwp); rms=[rmse(epOl) rmse(ep02) rmse(ep03)]; rsr=[rmse(epOl) rmse(ep02) rmse(ep03)]/rmse(dw)*100; subplot(3,1,1) plot(ca0*180/pi,dw,ca0*180/pi,dwhat01,':',... ca0*180/pi,dwhat02,'-.',ca0*180/pi,dwhat03,'--') t i t l e ( ' A c t u a l and Estimated V e l o c i t y Fluctuation') y l a b e l ( ' r a d / s e c ' ) subplot(3,1,2) sn=(l:length(JEhat02))'/length(JEhat02) * 2 * p i ; plot(sn*180/pi,JEhat02,'-.',sn*180/pi,JEhat03,'—') t i t l e ( ' E s t i m a t e d Engine I n e r t i a ' ) ylabel('kgm"2') subplot(3,1,3) plot(sn*180/pi,TLhat02,'-.',sn*180/pi,TLhat03,'--') t i t l e ( ' E s t i m a t e d Torque Fluctuation') ylabel('Nm') xlabel('degrees') disp('INERTIA VALUES') J_sls=JEhat01 J_ge=JEhat02(length(JEhat02)) J_lsef=JEhat03(length(JEhat03)) disp('PARAMETERS MEAN AND STANDARD DEVIATION') disp('GRADIENT ESTIMATOR') Jm=me an(JEhat 02) Js=rmse(JEhat02) Tm=mean(TLhat02) Ts=rmse(TLhat02) disp('LEAST-SQUARES WITH EXPONENTIAL FORGETTING') Jm=mean(JEhat03) Appendix B. Program Listings 89 Js=rmse(JEhat03) Tm=mean(TLhat03) Ts=rmse(TLhat03) dispCROOT MEAN SQUARE ERROR (rad/sec)') rms d i s p C RESIDUAL TO SIGNAL RATIO (°/„)') r s r prsal l .m '/.PRSALL.M "/.Calls: tpe.m, pwr.m, rmse.m. input('Enter c y l i n d e r number (0-7): ' ) ; n=ans; i f n==0 dw=wO; d i s p C NORMAL OPERATING CONDITION 1') e l s e i f n==l dw=wll; d i s p C FAULT 1L') e l s e i f n==2 dw=w3r; d i s p C FAULT 3R') e l s e i f n==3 dw=w31; disp('FAULT 3L') e l s e i f n==4 dw=w2r; disp('FAULT 2R') e l s e i f n==5 dw=w21; disp('FAULT 2L') e l s e i f n==6 dw=wlr; disp('FAULT 1R') e l s e i f n==7 • dw=wO; disp('NORMAL OPERATING CONDITION 2') end [pl p2 p3 p4 p5 p6]=ldown(n); disp('TORQUE DUE TO GAS PRESSURE') [Tp )Tpm,Tphat,deltaT,err]=... tpe(n )ca,dw Jpl,p2 )p3,p4,p5,p6); disp('RMS torque er r o r (Nm):') e_l=rmse(Tp-Tphat-deltaT) dispCRSR (%) : ') r_l=100*e_l/rmse(Tp) plot(ca*180/pi,Tp,ca*180/pi,Tphat+deltaT,':') a x i s ( [ 0 360 -1500 1500]) xlabel('degrees') ylabel('Nm') t i t l e ( ' E s t i m a t e d (..) and Actual (—) Pressure -Torque') disp('Press any key to continue...') pause [par,Tpx,plx,p2x,p3x,p4x,p5x,p6x]=... Appendix B. Program Listings pwr(n,ca,Tphat+deltaT+Tpm); disp('RMS torque e r r o r (Nm):') e_2=rmse(Tp-Tpx) disp ('RSR (7.) : ' ) r_2=100*e_2/rmse(Tp) plot(ca*180/pi,Tp,ca*180/pi,Tpx-mean(Tpx),':') a x i s ( [ 0 360 -1500 1500]) xlabel('degrees') ylabel('Nm') t i t l e ( ' E s t i m a t e d (..) and Actual (—) Pressure Torque') disp('Press any key to continue...') pause disp('PRESSURE WAVEFORM RECONSTRUCTION') disp('6 RMS pressure (MPa) er r o r s , and RSRs (7 . ) : ' ) err=[rmse(pl-plx) rmse(p2-p2x) rmse(p3-p3x)... rmse(p4-p4x) rmse(p5-p5x) rmse(p6-p6x)] rsr=100*[rmse(pl-plx)/rmse(pi) rmse(p2-p2x)/rmse(p2)... rmse(p3-p3x)/rmse(p3) rmse(p4-p4x)/rmse(p4)... rmse(p5-p5x)/rmse(p5) rmse(p6-p6x)/rmse(p6)] plot(ca*180/pi,pl,ca*180/pi,plx,':',ca*180/pi, ... p2,ca*180/pi,p2x,':',ca*180/pi,p3,ca*180/pi,p3x,':', ... ca*180/pi,p4,ca*180/pi,p4x,':',ca*180/pi,p5,ca*180/pi, .. p5x,':',ca*180/pi,p6,ca*180/pi,p6x,':') a x i s ( [ 0 360 -1 12]) xlabel('degrees') ylabel('MPa') t i t l e ( ' E s t i m a t e d (..) and Actual (—) Pressure Waveform') disp('Press any key to continue...') pause pwr .m f u n c t i o n [par,Tpx,plx,p2x,p3x,p4x,p5x,p6x]=. pwr(n,theta,Tphat) 7.FUNCTI0N PWR.M 7.ARGUMENTS: n 7, theta 7. Tphat 7.RETURNS: par 7. Tpx 7. pix 7oCalls: c_l.m, dlt.m 7.Engine data correspond to DDC 6V 92TA: r=0.0635; 7oCrankshaft radius (m) 1=0.2571; 7oconnecting rod length (m) Ap=0.01188; 7.piston area (m"2) cal=theta; c _ l l = c _ l ( c a l , r , 1 ) ; ca2=cal-(60-3.5)*pi/180; c_ 1 2=c_l ( c a 2,r,l); ca3=ca2-(60+3.5)*pi/180; c y l i n d e r number; crank angle; estimated torque; estimated pressure v a r i a t i o n ; new estimate of the pressure torque; estimated pressure i n the i - t h c y l i n d e r . Appendix B. Program Listings 91 c_13=c_l(ca3,r,l); ca4=ca3-(60-3.5)*pi/180; c_14=c_l(ca4,r,l); ca5=ca4-(60+3.5)*pi/180; c_15=c_l(ca5,r,l); ca6=ca5-(60-3.5)*pi/180; c_16=c_l(ca6,r,1); IX=[13 125 253 366 493 606]; i f n==0 IX=[17 130 256 370 496 610]; e l s e i f n==l IX(1)=17; e l s e i f n==2 IX(2)=131; e l s e i f n==3 IX(3)=253; e l s e i f n==4 IX(4)=368; e l s e i f n==5 IX(5)=496; e l s e i f n==6 IX(6)=611; end plx=dlt(2,5.5,0.45,theta-theta(IX(1))); p2x=dlt(2,5.5,0.45,theta-theta(IX(2))); p3x=dlt(2,5.5,0.45,theta-theta(IX(3))); p4x=dlt(2,5.5,0.45,theta-theta(IX(4))); p5x=dlt(2,5.5,0.45,theta-theta(IX(5))); p6x=dlt(2,5.5,0.45,theta-theta(IX(6))); phi=(le+6)*Ap*r*... [ c _ l l . * p l x c_12.*p2x c_13.*p3x c_14.*p4x c_15.*p5x c_16.*p6x]; par=inv(phi'*phi)*phi'*Tphat; par=par-0.15; i f n==0 par=par+0.15; end Tpx=phi*par; plx=plx*par(1); p2x=p2x*par(2); p3x=p3x*par(3); p4x=p4x*par(4); p5x=p5x*par(5); p6x=p6x*par(6); rmse.m f u n c t i o n y=rmse(x) "/.FUNCTION RMSE.M "/.ARGUMENTS: x = vector of r e s i d u a l s ; "/.RETURNS: y = root mean square e r r o r . y=norm(x-mean(x),2)/sqrt(length(x)); Appendix B. Program Listings ro t .m f u n c t i o n y=rot(x,n) "/.FUNCTION ROT.M "/.ARGUMENTS: x = vector of samples; 7, n = integer; "/.RETURNS: y = x s h i f t e d by n; 11 = length(x); y = x; f o r i = l : l l - n y ( i ) = x(i+n); end f o r i = l l - n + l : l l y ( i ) = x ( i - l l + n ) ; end separate.m f u n c t i o n [ p i , p2, p3, p4, p5, p6]=. separate(pres) "/.FUNCTION SEPARATE. M "/.ARGUMENTS: pres 7. 7. "/.RETURNS: [pi,p2,p3,p4,p5,p6] 7. p i = pres(l+3*720:4*720); p2 = pres(1+2*720:3*720) ; p3 = pres(1+5*720:6*720) ; p4 = pres(1+1*720:2*720); p5 = pres(1+4*720:5*720) ; p6 = pres(1+0*720:1*720); sls.m f u n c t i o n [Y,PHI,PAR,ERR,JEhat,TLhat,ep,dwhat]=... sls(caO,pl,p2,p3,p4,p5,p6,dw,dwp) "/.FUNCTION "/.ARGUMENTS: SLS.M caO = crank angle vector; 7. P i = c y l i n d e r pressure, i = l , . . . , 6; 7. dw = v e l o c i t y waveform; 7. dwp = v e l o c i t y d e r i v a t i v e ; '/.RETURNS: Y = vector of outputs; 7. PHI = s i g n a l matrix; • 7. PAR = estimated parameters; 7. JEhat = estimated engine i n e r t i a ; 7. TLhat = estimated torque f l u c t u a t i o n ; 7. ep = vector of r e s i d u a l s ; 7. dwhat = estimated v e l o c i t y f l u c t u a t i o n "/.Calls: c_l.m, c_2.m, c_3 m, c_4.m, igr.m. "/.Engine data correspond to DDC 6V 92TA: = a column vector with 6x720 elements which represent pressure data; = pressure data corresponding to each c y l i n d e r . Appendix B. Program Listings 93 w_0=125.5; b_c=1.24; m_p=6.03; r=0.0635; 1=0.2571; Ap=0.01188; J_C=0.1544; J_CR=0.0745; J_F=3.85; °/0engine speed (rad/sec) °/odamping (Nm/ (rad/sec)) 7 0pistori mass (kg) °/ 0crankshait radius (m) 7 0connecting rod length (m) °/ 0piston area (m~2) °/ 0crankshaft i n e r t i a (kg m"2) 7 0connecting rod i n e r t i a (kg m~2) 7oflywheel i n e r t i a (kg m"2) x0=w_0"2; 7«Pressure inputs: p l p2 p3 p4 p5 p6 dca0=ca0(2)-ca0(l); cal=caO; c _ l l = c _ l ( c a l , r , l ) c_21=c_2(cal , r , l ) c_31=c_3(cal , r , l ) c_41=c_4(cal , r , l ) ca2=cal-(60-3.5)*pi/180; c_12=c_l(ca2,r , l ) c_22=c_2(ca2,r,l) c_32=c_3(ca2,r,l) c_42=c_4(ca2,r,l) ca3=ca2-(60+3.5)*pi/180; c_13=c_l(ca3,r , l ) c_23=c_2(ca3,r,l) c_33=c_3(ca3,r,l) c_43=c_4(ca3,r,l) ca4=ca3-(60-3.5)*pi/180; c_14=c_l(ca4,r , l ) c_24=c_2(ca4,r,l) c_34=c_3(ca4,r,l) c_44=c_4(ca4,r,1) ca5=ca4-(60+3.5)*pi/180; c_15=c_l(ca5,r ,1); c_25=c_2(ca5,r , l) ; c_35=c_3(ca5,r , l) ; c_45=c_4(ca5,r , l) ; ca6=ca5-(60-3.5)*pi/180; c_16=c_l(ca6,r , l ) c_26=c_2(ca6,r,l) c_36=c_3(ca6,r,l) c_46=c_4(ca6,r,l) s l l = c _ l l . * c _ l l + c_12.*c_12 +. c_13.*c_13 + c_14.*c_14 +. . . c_15.*c_15 + c_16.*c_16; s l2=c_l l .*c_21 + c_12.*c_22 +. c_13.*c_23 + c_14.*c_24 +. . . c_15.*c_25 + c_16.*c_26; s33=c_31.*c_31 + c_32.*c_32 +. c_33.*c_33 + c_34.*c_34 +. . . c_35.*c_35 + c_36.*c_36; s34=c_31.*c_41 + c_32.*c_42 +. c_33.*c_43 + c_34.*c_44 +. . . c_35.*c_45 + c_36.*c_46; Appendix B. Program Listings 94 s p = ( l e + 6 ) * ( c _ l l . * p l + c_12.*p2 +... c_13.*p3 + c_14.*p4 +... c_15.*p5 + c_16.*p6); dx=2*w_0*dw; dxp=2*w_0*dwp; phil=0.5*dxp; fact=max(phil); p h i l = p h i l / f a c t ; phi2=sin(ca0); phi3=cos(caO); PHI=[phil,phi2,phi3]; yl=m_p*r~2*(-0.5*sll.*dxp-sl2.*(dx+xO)); y2=J_CR*(0.5*s33.*dxp+s34.*(dx+xO)); y3=Ap*r*sp; y=yl+y2+y3; my=mean(yl+y2+y3); Y=y-my; PAR=inv(PHI'*PHI)*PHI'*Y; ERR=PHI*PAR-Y; JEhat=PAR(l)/fact; TLhat=PAR(2)*sin(ca0)+PAR(3)*cos(ca0) ; m=(JEhat+m_p*r~2*sll-J_CR*s33)/2; inp=y3-TLhat; inp=inp-mean(inp); dxhat = dw(l)*(2*w_0+dw(l)) + igr(inp./m.dcaO); dwhat = sqrt(xO + dxhat) - w_0; dwhat=dwhat-mean(dwhat); °/.dwhat=0 . 9*rot (dwhat, 1) ; ep=dw-dwhat; system.m f u n c t i o n xp=system(t,x) °/„FUNCTION: SYSTEM, m '/.ARGUMENTS: t = crank angle; 7. x = engine state; "/.RETURNS: xp = state d e r i v a t i v e , g l o b a l p i p2 p3 p4 p5 p6 ca 7.Engine data correspond to DDC 6V 92TA: w_0=125.5; '/oengine speed (rad/sec) b_c=1.24; "/.damping (Nm/ (rad/sec)) m_p=6.03; "/.piston mass (kg) r=0.0635; "/.crankshaft radius (m) 1=0.2571; "/.connecting rod length (m) Ap=0.01188; "/.piston area (m~2) J_C=0.1544; "/.crankshaft i n e r t i a (kg m'2) J_F=3.85; "/.flywheel i n e r t i a (kg m~2) J_CR=0.0745; "/.connecting rod i n e r t i a (kg m~2) "/.Pressure inputs: p l = i n t e r p l ( c a , p l , t ) ; p2=interpl(ca,p2,t); Appendix B. Program Listings p3=in te rp l (ca ,p3 , t ) ; p4=in te rp l (ca ,p4 , t ) ; p5=in te rp l (ca ,p5 , t ) ; p6=in te rp l (ca ,p6 , t ) ; cal=t; c _ l l = c _ l ( c a l , r , 1 ) ; c_21=c_2(cal ,r ,1); c_31=c_3(ca l , r , l ) ; c_41=c_4(ca l , r , l ) ; ca2=cal-(60-3.5)*pi/180; c_12=c_l (ca2 , r , l ) ; c_22=c_2(ca2,r , l) ; c_32=c_3(ca2,r,1); c_42=c_4(ca2,r , l) ; ca3=ca2-(60+3.5)*pi/180; c _ i 3 = c _ l ( c a 3 , r , l ) ; c_23=c_2(ca3,r , l) ; c_33=c_3(ca3,r , l) ; c_43=c_4(ca3,r , l) ; ca4=ca3-(60-3.5)*pi/180; c _ i 4 = c _ l ( c a 4 , r ) l ) ; c_24=c_2(ca4,r , l) ; c_34=c_3(ca4,r , l) ; c_44=c_4(ca4,r , l) ; ca5=ca4-(60+3.5)*pi/180; c_15=c_l(ca5,r ,1); c_25=c_2(ca5,r,1); c_35=c_3(ca5,r , l) ; c_45=c_4(ca5,r , l) ; ca6=ca5-(60-3.5)*pi/180; c_16=c_l (ca6 , r ; i ) ; c_26=c_2(ca6,r , l) ; c_36=c_3(ca6,r , l) ; c_46=c_4(ca6,r , l) ; s l l = c _ l l * c _ l l + c_12*c_12 +. . . c_13*c_13 + c_14*c_14 +. . . c_15*c_15 + c_16*c_16; s l2=c_ l i*c_2 l + c_12*c_22 +. . . c_13*c_23 + c_14*c_24 +. . . c_15*c_25 + c_16*c_26; s33=c_31*c_31 + c_32*c_32 +. . . c_33*c_33 + c_34*c_34 +. . . c_35*c_35 + c_36*c_36; s34=c_31*c_41 + c_32*c_42 +. . . c_33*c_43 + c_34*c_44 +. . . c_35*c_45 + c_36*c_46; sp=(le+6)*(c_ll*pl + c_12*p2 +. . . c_13*p3 + c_14*p4 +. . . c_15*p5 + c_16*p6); m=0.5*(J_F+J_C-J_CR*s33+m_p*r~2*sll); A = J_CR*s34-m_p*r"2*sl2; fact = Ap*r; bTu=fact*sp; Tf=-b_c*w_0; Appendix B. Program Listings Tl=-958.7535; d=Tf+Tl; x0=w_0~2; xp = (A*(x+x0) + bTu + d)/m; tpe .m f u n c t i o n [Tp,Tpm,Tphat,deltaT,err]=... tpe(n,theta,dw,p1,p2,p3,p4,p5,p6) "/.FUNCTION TPE.M •/.ARGUMENTS: n = c y l i n d e r number; °/„ theta = crank angle vector; °/„ dw = angular v e l o c i t y ; °/. p i = pressure of the i - t h c y l i n d e r ; •/.RETURNS: Tp = actu a l pressure torque; •/„ Tpm = mean pressure torque; °/„ Tphat = estimated pressure torque; 7. deltaT = torque f l u c t u a t i o n ; 7, err = vector of r e s i d u a l s . "/.Calls: c_l.m, c_2.m, c_3.m, c_4.m, diffc.m. •/.Engine data correspond to DDC 6V 92TA: w_0=125.5; x0=w_0~2; b_c=l.24; m_p=6.03; r=0.0635; 1=0.2571; Ap=0.01188; J_C=0.1544; J_CR=0.0745; J_F=3.85; J_E=5; end cal=theta; c _ l l = c _ l ( c a l , r , 1 ) ; c_21=c_2(cal,r,l); c_31=c_3(cal,r,l); c_41=c_4(cal,r,l); ca2=cal-(60-3.5)*pi/180; c_12=c_l(ca2,r,l); c_22=c_2(ca2,r,l); c_32=c_3(ca2,r,l); c_42=c_4(ca2,r,l); ca3=ca2-(60+3.5)*pi/180; c_13=c_l(ca3,r,l); c_23=c_2(ca3,r,l); c_33=c_3(ca3,r,1); c_43=c_4(ca3,r,l); ca4=ca3-(60-3.5)*pi/180; c_14=c_l(ca4,r,l); c_24=c_2(ca4,r )l); Appendix B. Program Listings c_34=c_3(ca4,r,1); c_44=c_4(ca4,r,l); ca5=ca4-(60+3.5)*pi/180; c_15=c_l(ca5,r,l); c_25=c_2(ca5,r,l); c_35=c_3(ca5,r,l); c_45=c_4(ca5,r,l); ca6=ca5-(60-3.5)*pi/180; c_16=c_l(ca6,r,l); c_26=c_2(ca6,r,l); c_36=c_3(ca6,r,l); c_46=c_4(ca6,r,l); s l l = c _ l l . * c _ l l + c_12.*c 12 +... c_13.*c_13 + c_14.*c_14 +... c_15.*c_15 + c_16.*c_16; sl2=c _ l l . * c _ 2 1 + c_12.*c 22 +... c_13.*c_23 + c_14.*c_24 +... c_15.*c_25 + c_16.*c_26; s33=c_31.*c_31 + c_32.*c_32 +... c_33.*c_33 + c_34.*c_34 +... c_35.*c_35 + c_36.*c_36; s34=c_31.*c_41 + c_32.*c 42 +... c_33.*c_43 + c_34.*c_44 +... c_35.*c_45 + c_36.*c_46; s p = c _ l l . * p l + c_12.*p2 +... c_13.*p3 + c_14.*p4 +... c_15.*p5 + c_16.*p6; Tp=(le+6)*Ap*r*sp; Tpm=mean(Tp); Tp=Tp-Tpm; q = J_CR*s34-m_p*r~2*sl2; m = 0.5*(J_E-J_CR*s33+m_p*r"2*sll); ca0=[0:117] '/H8 * 2 * p i ; h=ca0(2)-ca0(l); dwp=diffc(dw,h); dw=interpl(caO,dw,theta,'spline'); dwp=interpl(caO,dwp,theta,'spline') dx=2*dw*w_0; dxp=2*dwp*w_0; dl=x0*q; Y=m.*dxp-q.*dx-dl; Tphat=Y-mean(Y); err=Tp-Tphat; phi= [sin(theta) c o s ( t h e t a ) ] ; par=inv(phi'*phi)*phi'*err; deltaT=phi*par; Index B Beard-Jones filter, 18 bottom dead centre, 7, 10, 23, 73 C connecting rod acceleration, 23, 24 angle, 22 axis, 23 inertia, 26, 37, 48 length, 22, 23 mass, 22, 26 mechanism, 22, 23 crankshaft radius, 22, 23, 37, 48 cylinder pressure template, i i , 4, 46, 50, 55, 56, 63 variation, i i , 49-52, 55, 56, 63 D decision rule, 3, 14, 16, 17, 19, 34, 55 E estimator gradient, i i , 3, 18, 19, 34, 39-44, 61, 62 least-squares, i i , 3, 18, 19, 34, 38-40, 42-44, 55, 61-63 F forgetting factor, 19, 39 K Kalman filter, 14 L Luenberger observer, 14, 15, 17 M misclassification, 36, 44, 62 P parameters mean, 3, 19, 34, 43, 62, 64 standard deviation, 34, 43, 44, 62, 64 piston acceleration, 23, 24, 26 crown area, 25, 46, 48 position, 22-24 velocity, 22-24 probability distribution, 49 R residual generation, 3, 14-18, 34 to signal ratio, 2, 31, 32, 41, 43, 44, 52-54, 56, 57, 61, 63 root mean square error, i i , i i i , 30-32, 41, 43, 52-54, 56, 57, 61, 63 T top dead center, 40, 52 98 Index top dead centre, 7, 8, 10, 22, 23, 25, 29 torque coefficients, 27, 37, 48 connecting rod, 22, 26 fluctuation, i i , 3, 31, 33, 34, 36, 40, 41, 53, 61, 62 friction, 22, 30, 37, 40, 52 gas pressure, i i , i i i , 4, 21, 37, 45-47, 49, 51, 53, 54, 57, 58, 63 load, 6, 9, 22, 30, 31, 37, 40, 52, 73 reciprocating parts, 22 total engine, 10, 21, 30, 40, 52