- Email: [email protected]
On the Application of the Empirical Mode Decomposition to Dynamic Positioning Systems
8th IFAC Conference on Control Applications in Marine Systems Rostock-Warnemünde, Germany September 15-17, 2010
On the Application of the Empirical Mode Decomposition to Dynamic Positioning Systems Leonardo Kubota*, Helio. M. Morishita**, Solenn Greuell, *** Lázaro Moratelli Jr.*** *Department of Naval Architecture and Ocean Engineering University of São Paulo Brazil (Tel: 55-11-30915345; e-mail: [email protected]). ** Department of Naval Architecture and Ocean Engineering University of São Paulo (e-mail: [email protected]) *** Department of Naval Architecture and Ocean Engineering University of São Paulo Abstract: In this work, a new filtering algorithm suitable for applications in Dynamic Positioning Systems (DPS) is presented. This new filtering method relies on the Empirical Mode Decomposition (EMD) and the reasons that make the EMD Filter really attractive are as follows: it can handle both non stationary and non linear signals, it is an adaptive filtering algorithm and it requires no previous knowledge on the system dynamics (model independent). In what follows, we provide details on the algorithm we designed, and then present the results of experimental tests. Keywords: Empirical Mode Decomposition, Dynamic Positioning System, Adaptive Filter, EMD Filter smoothed out, before they are fed into the control algorithm which calculates propeller forces, hence the need for a filtering algorithm. The purpose of the wave filter is to separate the high-frequency wave-induced motion from the motion caused by slowly varying disturbances. Feedback control action must be implemented using the filtered lowfrequency vessel motion, enabling thruster modulation and all related problems to be avoided. Commercially available DPSs feature the well-known Kalman Filter algorithm. The Kalman Filter relies on a previous knowledge of the (linear) system dynamics (in the problem at hand, the low frequency range of the vessel dynamics) and an assumed model for high-frequency motions of the vessel to perform the filtering of wave-induced motion. Actually, as the Kalman Filter is based on a linear mathematical model of the system, all its parameters need to be tuned for a range of operational conditions, to that end it is common practice to use an automatic gain-scheduling procedure.
1. INTRODUCTION Dynamic Positioning System (DPS) can be defined as a computer controlled system designed to keep a floating vessel on a specific position or pre-defined path by using her own propellers. Position and heading sensors, wind sensors, a set of algorithms and propellers are the members of a typical DPS. Those systems have found a wide array of applications in the offshore oil exploration setting, as evidenced by the deployment of such technology in operations such as drilling, pipe-laying, offloading as well as diving support. When in the ocean, floating vessels experience environmental forces that are anything but simple. Their complexity lies not only in the fact they can be seen as the sum of an endless amount of oscillatory components of varying amplitudes, frequencies, phases and directions, their mathematical description is also quite challenging, for each has its own intrinsic characteristics. Nonetheless, their action on floating vessels gives rise to rigid body motions made up of both high and low frequency motions. Wind-generated waves are at the root of large, high-frequency, oscillatory forces and moments that co-exist with their low frequency counterparts originated by wind and current loads, as well as the so-called wave-drift forces, associated with second order effects present in waveinduced forces.
In this work, the authors would like to introduce a novel approach to filtering which is not only model independent, but it is also appropriate for handling signals possessing a high degree of nonlinearity and non-stationarity, conditions which would be at the root of a poor performance of the Kalman Filter. We also intend to demonstrate how this filtering algorithm would perform in a tanker scale model fitted with DPS, where filtering will be performed in real time using our novel approach to filtering: the Empirical Mode Decomposition Filter [Huang et al, 1998], the EMD Filter.
For the design of DPS, it is generally those low-frequency and steady components that must be accounted for, because they are the ones whose physical nature would cause a ship to go astray and drift off its set-point or original path. The highfrequency motions, on the other hand, should not be actuated upon by the controller, since that would require an enormous amount of energy which could not be provided by the enginepropeller system. As a result, we arrive at the conclusion that high-frequency motion components must be suppressed or 978-3-902661-88-3/10/$20.00 © 2010 IFAC
In Section 2, we describe the concepts and ideas underlying the EMD Filter, and present some examples where the EMD Filter is applied to real signals obtained from GPS readings of a marine vessel position. Section 3 provides an overview of the mathematical model of a floating vessel under the 138
10.3182/20100915-3-DE-3008.00071
CAMS 2010 Rostock-Warnemünde, Germany, Sept 15-17, 2010
8.
Repeat iterations until the function become monotonic and there is no more IMFs to be extracted from it. A stopping criterion determines the number of sifting steps to produce an IMF. Huang et al (1998) suggest that the sifting should be stopped as soon the sum of the difference (SD):
external action of current, wind and wave forces, so the reader can have a glimpse of the multitude of parameters which must be estimated in most filtering algorithms. Then, the same reader is informed that the EMD Filter requires no previous knowledge on the system dynamics. Section 4 puts our EMD Filter to the test in an experimental setup where a scale model of an oil tanker fitted with DPS is run in a academic towing tank.
T
∑h
k −1
SD =
The main advantage of this novel approach to filtering, which is still in its early stages, is the fact it requires no previous knowledge regarding the system dynamics (such as inertia, damping, stiffness matrices), nor is it bounded by assumptions of linearity and/or stationarity. The main drawback of the EMD filter is the lack of a solid mathematical grounding; after all it is a heuristic algorithm, which makes it difficult to perform a stability analysis.
(t ) − hk (t )
2
0 T
∑h
2 k −1
(t )
0
2.2 A novel approach to filtering: EMD Filter The EMD claims to properly identify the time scales featuring the physical characteristics of any given system by analyzing a time series representative of the system dynamics. The original signal is then broken down into a set of functions, each displaying an oscillatory mode intrinsic to the system, the so-called Intrinsic Mode Function (IMF). Lower order IMFs retain the most oscillatory temporal modes.
2 THE EMPIRICAL MODE DECOMPOSTION FILTER (EMD FILTER) 2.1 The Empirical Mode Decomposition (EMD) algorithms The EMD algorithm decomposes any given data into a set of intrinsic mode functions (IMFs). An IMF is defined as a function that satisfies two conditions:
Mathematically, the EMD algorithm decomposes the function X (t ) into: N
X (t ) = ∑ hi (t ) + r (t )
1. In the whole data set, the number of extrema and the number of zero-crossings must either be equal or differ at most by one. 2. At any point, the mean value of the envelope defined by the local maxima and the envelope defined by the local minima is zero.
Where r(t) is the residual. Equation (1) contains the central idea behind the EMD Filter. It simply states that the raw signal can be mathematically written as the sum of N IMFs, note that IMFs isolate time scales present in the signal. From this equation, it can be seen that by retaining only the low-frequency terms in the reconstruction above, we obtain a filtered version of X (t ) , for each IMF is a representation of a temporal mode within the raw signal. That way, the first IMF tends to be a highly oscillatory mode, and as we move towards the higher order IMFs, they become less and less oscillatory, in other words, they are low frequency modes. That is the central idea behind the EMD Filter: use the EMD to break the raw signal down into IMFs, and then add up higher order (low frequency) IMFs. What comes out of this procedure is a filtered version of the raw signal, a version that contains only the low frequency temporal modes.
Huang et al (1998) suggests that an IMF represents a single oscillatory mode (not a structural mode, but a temporal mode) which generalizes a harmonic function in the way that its amplitude and frequency are time-varying. The process of extracting IMFs off a given signal is termed sifting, which is defined by the following algorithm: 1. 2. 3. 4. 5.
Between each successive pair of zero crossings, identify a local extremum in the test data. Connect all the local maxima by a cubic spline line as the upper envelope. Repeat the procedure for the local minima to produce the lower envelope. Calculate the envelope mean m1 Calculate the first estimate of IMF 1:
We again point out that the decomposition procedure requires no previous knowledge regarding the system dynamics. The EMD Filter corresponds to an automatic and adaptive (signaldependent) time-variant filter capable of handling both nonlinear and non stationary signals.
h1 = X (t ) − m1 6. Take h1 as the new signal and repeat iterations from 7.
(1)
i =1
1 to 4, until the stopping criteria is met (see below). We then have IMF 1. The first IMF h1 , contains the high frequency components of the original signal. From the relation r1 = X (t ) − h1 , we obtain the residue which is taken as the new signal to be sifted through.
2.3 The EMD Filter: offline filtering algorithm Now in the following example, we shall consider signals obtained from GPS measurements of a full scale oil tanker. The signal analyzed here corresponds to the vessel’s center of gravity calculated from positions captured from 2 GPS. The signals provided were fed into the algorithm and filtered offline, that is, in one single batch, as opposed to filtering it on139
CAMS 2010 Rostock-Warnemünde, Germany, Sept 15-17, 2010
out of the measured signal, thus producing a filtered point at any given time step, as opposed to producing a filtered vector of points (offline problem).
line using a moving window. In Section 2.4, we provide a more detailed explanation on how to interpret the online filtering problem as well as the offline filtering problem on the EMD Filter setting.
Our study has shown that filtering a noisy signal using an online and an offline approach do not lead to the same results. The reason for this can be grasped when we notice that during an offline analysis, at any given time instant, we have full knowledge of past and future values in the time series, in other words we benefit from a “look-ahead” bias, which is not present during an online filtering, where at any given moment in time we have no information whatsoever on future values. That will be the reason why online filtering performance will not be as good as its offline counterpart.
The original and filtered signals for the sway motion are shown in Fig. 1. Note that for an off-line analysis, the EMD filter shows a remarkably good performance without introducing any lags into the filtered signal. What kind of previous information regarding the system dynamics was required for the algorithm to produce such results? None. The method is signal dependent and model independent. Some shortcomings originating from this will be discussed in Section 5. In order to filter the raw signal, we added up IMFs from order N/2 (or the nearest integer to N/2) on.
During an online filtering analysis, at step k, one might wonder what the filtered value to be fed into the control law should be. Our tests have suggested that taking the last point in the filtered vector f fails to produce good results. The authors believe that end effects might be responsible for such poor filtering performance. As an attempt to circumvent such an issue, we suggest the following strategy for online filtering:
As it will be shown later on, in DPS applications the proper approach to EMD filtering should require the use of a moving window, as the filtering is to be performed in real time - the online filtering problem, which we should detail in the next section. -1.2575
x 10
6
Raw signal vs EMD Filtered Signal-SWAY yg filtered yg
1.
Update the noisy input vector s of size n by eliminating the oldest value and including the measured value at instant k; 2. Filter the input signal via EMD Filter. We have vector f of size n (it is important to see that both the filtered and noisy signal have the same number of elements), in our study we adopted n=100. 3. Take the n-j is the ilter delay time; 1 ) value in the f vector where the last observations correspond to the last obtained values, which allows us to move away from the filtered points contaminated by the end effect. 4. Take the quantity obtained in step 2 as the filtered value for j time steps. In order to avoid excessive phase lag j should be around 1. In fact there is a trade-off involved in the choice of parameter j , which is why we suggest the reader to adopt values around 1 . As we mentioned in the 3rd item above, the parameter j allows us to move away from the end-effects, but at the same time, values greater than 1 cause the filtered signal to lag behind the raw signal.
-1.2575
-1.2575
-1.2575
-1.2575
-1.2575
-1.2575
-1.2575
-1.2575
0
100
200
300
400
500
600
700
800
900
1000
Figure 1: Filtered sway motion in red.
2.4 The EMD Filter: online filtering algorithm Unfortunately, notwithstanding the great performance of the off-line filter, it is not feasible for application in DPS because its control loop requires a real time filtering of the signals of the position and heading that are captured at sampling instants k. The natural solution is to adopt a moving window with a pre-determined length in which at every sampling instant the oldest value is eliminated and the new value is inserted in the data vector.
To realize the influence of the parameter j, the spectrums of the filtered signals for j = 1, 2, and 3 are shown in Figure 2. The surge motion of the full scale ship captured by GPS, as previously explained, is taken as the raw signal. It is clear that values greater than 1 for the parameter j, which is associated with the phase lag, lead to better filtering, as evidenced by the spectrum, since we move away from the zone affected by the end effects.
Here it is important to distinguish the online filtering problem from the offline filtering problem, presented on Section 2.3. The offline filtering problem presents us with a noisy data, and demands us to find a filtered (noise-free) version of it. As we have shown in the Section 2.3, the EMD Filter is quite a suitable tool for tackling such problem. For one thing, it produces a smoothed version of a noisy signal without introducing any lag. Moreover, it requires no previous knowledge on the system dynamics.
3 MATHEMATICAL MODEL The modeling of dynamics of a vessel in the horizontal plane is formulated considering an earth-fixed frame OXEYEZE and a vessel fixed reference frame GXYZ, as indicated in Fig. 3. The origin of the vessel fixed frame is located in the centerline, in a distance from the center of gravity. The axes of each body-fixed co-ordinate system coincide with the
The online filtering problem, as seen in the DPS setting, consists of, at time step k, taking high frequency components 140
CAMS 2010 Rostock-Warnemünde, Germany, Sept 15-17, 2010
0 0 -) ) " .) ) /# - ) - ) " -) ) ! - ) - ) ! , (0 0 0 0 .) ) /! Where ) is the mass of the vessel, ) , 1, 1,2,6 are the added mass coefficients in surge, sway and yaw, respectively; ! and " are the surge and sway velocities, respectively;! and " are current speeds related to the body fixed frames, respectively; # is the yaw rate; + is the moment of inertia about vertical axis; %, & and ' represent the total external forces and moment in surge, sway and yaw directions, respectively; dot means time derivative of the variable. 4 5 is the control vector and 4 5 is a matrix defined by the configuration of the actuators. The position and heading of the vessel related to the earth-fixed coordinate system are obtained from the following equation
principal axes of inertia of the vessel. To represent the intricate motion of a vessel at sea the mathematical model is usually split in two terms namely the low frequency motion caused by wind, current and second-order wave forces, and the high frequency motions caused by the first order wave forces. The former results by applying the Newton law and the latter are time series generated from a spectrum of the motion for each degree of freedom.
Spectrum for j=1 Spectrum for j=2
0.15
Spectrum
Spectrum for j=3 0.1 0.05
6 76
0 3
6
0.3
2
0.2
j
;<=9 76 : =19 0
0.1 1
Frequency
0
Figure 2: Spectrum of the surge time series shown in Figure 2, where the parameter j associated with the phase lag is increased.
Figure 3 Earth fixed and vessel fixed references frames
% & '$
) ) 0 ( 0
0 ) ) ) )
0 ) ) , + )
-=19 ;<=9 0
0 0> 1
(6) (7)
?@ |BCD@, E| F@
The equations for low frequency motions, related to the body fixed reference frame, and considering added mass forces, are given by: ! " #$
9 $
where , 8 and 9 are the coordinates and heading of the vessel in the earth fixed frame respectively. The forces %, & and ' are due to the action of current, wind, and yaw hydrodynamic damping. Forces due to current and wind can be determined using drag equations (Leite et al., 1998; OCIMF, 1994). The second order wave forces are caused by the so-called slow and mean drift forces and they can be estimated as presented in Aranha (1994), and Aranha and Fernandes (1995). The high frequency motions (denoted as 6 in what follows) cannot be easily predicted by using differential equations because the calculations of the first order wave forces in real time are complex numerical processes. In light of this a time series generated from the high frequency motion spectrum is used to perform numerical analysis of the dynamics of the vessel. The spectrum for each degree of freedom is given by:
.
8
(5)
(2)
where RAO is the ship Response Amplitude Operator, β is the wave incidence angle related to the ship and S (ω ) is the sea spectrum.
(3)
3.1 Measurement System
(4)
For conventional ships, only position and heading measurements are available. Position is normally obtained from a DGPS while heading is usually measured by using a gyro compass. Hence, the position and heading measurement equation can be written as: G 6 6 H
(8)
Where d ∈ ℜ3 is a zero mean Gaussian white measurements noise.
141
CAMS 2010 Rostock-Warnemünde, Germany, Sept 15-17, 2010
The performance of the EMD filter was checked considering step variation in the reference values for the coordinates of the axes OXE, and OYE, and the heading 9, for different values of the parameter j. The test was performed adopting sampling period of 0.1 s. The model was disturbed by regular waves characterized by frequency of 1,0Hz and 1cm high. However, in this paper, only the influence of the parameter j in the performance of the EMD is pointed out. For that, results for surge motion for j=1, 2 and 3 are shown in Figure 4. Note that the EMD Filter in its online version outlined in Section 2.4 was capable of performing a real time filtering of the raw signal being measured. Besides, it also managed to properly capture the variation of the motion if the setpoint is altered, which higlights the adaptive nature of the EMD Filter. Again, we stress the fact we did not provide the filter with any prior information on the test model dynamics nor the wave characteristics, and yet the EMD Filter performed considerably well in all 3 conditions.
3.2 Conventional Filtering of the High Frequency Motions The high frequency motions of the vessel should be filtered because their averages are zero and an attempt to counteract them would require a huge power, leading to extra fuel consumption and propulsion degradation. The first generation of the SPDs were implemented by using notch filter to reduce the amplitude of the high frequency motions. The second generation of the SPDs have adopted Extended Kalman Filter that assumes an internal linear model based on equations (2) and (5) for estimating low frequency positions and velocities. For filtering of the high frequency motions a state-space realization of a transfer function that satisfies the following approximation is included: ?@ I |J@|
The transfer function for each direction is given by:
G( s) =
η w ( s) w(s)
=
σs s 2 + 2ξω0 s + ω02
Granted, in order for the EMD Filter to start, it does require an initial window of acquired data. In our tests, we collected 100 points during 10 seconds, and once we had our initial data window, the EMD Filter could be set off. As time goes by, the 100-point-window is continously updated, as it moves forward in time discarding older readings and incorporating new ones, as explained in Section 2.4. Figure 4 also brings Figure 2 to the time-domain, so we can see how larger values of parameter j lead to both better filtering (as evidenced by Figure 2) and – what is not evident from Figure 2 – to the introduction of a slight phase lag, the trade-off brought up in previous sections.
As a result, the Kalman Filter works with state vector of order 12 to filter the high frequency motions and estimate the low frequency motions and velocities. Therefore it is clear, from what we explained so far, that the application of the EMD Filter to solve the filtering problem in DPS would be a very interesting alternative indeed. For one, there would be no need to estimate dynamical properties of the vessel such as mass, added mass, damping coefficients and others which are not always easily obtained. Besides, should the sea state changes, the EMD Filter would automatically adapt, thus eliminating the need for a cumbersome gain adjustment procedure, which troubles the Kalman Filter.
As for the performance of the DPS system as a whole, we see in Figure 4 that the introduction of the EMD Filter does produce good results, as the test model manages to track the reference signal. A better adjustment procedure in the PID controller gains would have led to even better results. However, during the tests it was realized that if the input raw signal is not noisy enough (if it flatlines or oscillates at very low frequencies), the EMD procedure breaks down or iterates for too long until the stopping criterion is met. So in DPS applications, the more the floating vessel moves about the better the filtering will be. On the other hand, in case of a violent storm, our EMD Filter would quickly adapt to the external conditions and the filtering algorithm would carry on unflappably
4 APPLICATION OF THE EMD FILTER TO DPS In order to check the viability of the application of the EMD Filter in DPS using the online filtering procedure proposed in Section 2.4, a test with scale model was performed in the tank of the Department of Naval Architecture and Ocean Engineering of the University of São Paulo. The model corresponds to a typical offloading vessel and the actuators are the bow thruster, the stern thruster and the main propeller. The scale model main particulars in full load conditions and details of the laboratory facilities are described in Morishita et al (2009). A conventional PID (Proportional – Integral – Derivative) controller for each degree of freedom was used.
5 CONCLUSIONS
To perform tests with EMD Filter three vectors were considered, namely K 4 5 , L 4 5 and M 4 5 where K, L and M are vectors composed of measured data obtained from G (eq. 8) by using the procedure suggested above. At every time N those vectors are updated and the EMD filter is applied for each vector. Then the PID control law is calculated for each degree of freedom using the filtered signals. The parameters of the controller were tuned by means of genetic algorithm considering ! " 0 in equation (2). The computer code to perform the control law is based on Matlab but EMD algorithm is a C compiled program to save computer time.
This work presents a new filtering algorithm for application in DPS. We named it the EMD Filter, for the filtering strategy we created relies on a prior step where the raw signal is decomposed into Intrinsic Mode Functions by a procedure that came to be known as EMD proposed by Huang et al (1998). The filtering algorithm we created is capable of filtering out high frequencies of a raw signal by retaining only the low frequency Intrinsic Mode Functions (IMFs) [Huang et al.,1998] . We came forward with an offline filtering algorithm, which helps us solve the problem where given a raw signal, we are asked to produced a smoother, filtered version of it. And, we also devised an online filtering
4.1 Model test results 142
CAMS 2010 Rostock-Warnemünde, Germany, Sept 15-17, 2010
Workshop on Nonlinear Signal and Image Processing NSIP,v.3,pp.8-11.
algorithm, where the extraction of high frequency components happens in a real-time fashion. So given a moving window of fixed length containing noisy observations, we can provide at any given time step, a filtered point (as opposed to a filtered vector of points), that will be used in the control law. Both the offline filtering and the online filtering problem and its related algorithms are presented in detail in Section 2, where we also discuss the inevitability of introducing a phase lag in online filtering in order to escape the end effects, and the trade-off arising out of it.
(a) j = 1 -2.75
-2.8
surge [m]
-2.85
Perhaps, the most important claim we stake throughout this work is that the EMD Filter algorithm is model-independent, as it requires no information on the system dynamics. Unlike the Kalman Filter, where a multitude of gains and parameters need to be tuned (and readjusted once the external conditions, such as sea state, changes), the EMD Filter is also adaptive, as we have shown in Section 4, in a scale model run in a water tank.Some shortcomings found by the authors during the investigations were: 1) the input signal need be noisy to the EMD filter to work properly. 2) it is not yet possible to introduce a cutoff frequency in the EMD Filter design.
-2.9
-2.95
-3 70
75
80
85
90
95
100
100
110
120
75
80
85
time [s]
(b) j = 2 -2.7
The EMD Filter we propose in this work is still in its early stages. There certainly is room for further improvements and investigation on the algorithm. And the search for new fields of application of this inovative technique looks quite promissing.
-2.75
surge [m]
-2.8
ACKNOWLEDGMENTS The authors thank Petrobras for financial support in implementing the project. The second author thanks the CNPq, under the processes 558161/2009-0.
-2.85
-2.9
REFERENCES
-2.95 60
70
80
90
time [s]
Aranha, J.A.P. (1994) A Formula for Wave Damping in the Drift of a Floating Body”, Journal of Fluid Mechanics, vol. 272, pp.147-155. Aranha, J.A.P., Fernandes, A.C., (1995), Second-order slow drift force spectrum, Applied Ocean Research, v.17,n.5,p.311-313. Huang, N. E, Shen, Z., Long, S.R., Wu, M.C., Shih H.H., Zheng, Q., Yen, N., Tung, C.C., Liu, H.H. (1998) The empirical model decomposition and the Hilbert spectrum for nonlinear and non-stationary time series analysis, Proc.R. Soc. Lond. A454, 903-995. Leite, A.J.P., Aranha, J.A.P., Umeda, C. and de Conti, M.B.(1998), Current Forces in Tankers and Bifurcation of Equilibrium of Turret System: Hydrodynamic Model and Experiments, Applied Ocean Research, Vol. 20(3), pp. 145-156. Morishita, H.M., Tannuri, E.A., Saad, A.C., Sphaier, S.H., Lago, G.A., Moratelli Jr, L., (2009) Laboratory Facilities for Dynamic Positioning System, Proc. International Conference on Manoeuvring and Control of Marine Craft, Guarujá. OCIMF. Predictions of wind and current loads on VLCCs, (1994) Oil Companies International Marine Forum. Rilling,G.,Flandrin, P.,Gonçalvès,P., (2003), On Empirical Mode Decomposition and its Algorithms,IEEEURASIP
(c) j = 3 -2.8 -2.82 -2.84
surge [m]
-2.86 -2.88 -2.9 -2.92 -2.94 -2.96 -2.98 55
60
65
70
time [s]
Figure 4: Results of test on scale model DPS - surge motion (a) j=1; (b) j = 2; (c) j = 3 (blue:raw signal, black: reference signal, red: filtered signal)
143
On the Application of the Empirical Mode Decomposition to Dynamic Positioning Systems Leonardo Kubota*, Helio. M. Morishita**, Solenn Greuell, *** Lázaro Moratelli Jr.*** *Department of Naval Architecture and Ocean Engineering University of São Paulo Brazil (Tel: 55-11-30915345; e-mail: [email protected]). ** Department of Naval Architecture and Ocean Engineering University of São Paulo (e-mail: [email protected]) *** Department of Naval Architecture and Ocean Engineering University of São Paulo Abstract: In this work, a new filtering algorithm suitable for applications in Dynamic Positioning Systems (DPS) is presented. This new filtering method relies on the Empirical Mode Decomposition (EMD) and the reasons that make the EMD Filter really attractive are as follows: it can handle both non stationary and non linear signals, it is an adaptive filtering algorithm and it requires no previous knowledge on the system dynamics (model independent). In what follows, we provide details on the algorithm we designed, and then present the results of experimental tests. Keywords: Empirical Mode Decomposition, Dynamic Positioning System, Adaptive Filter, EMD Filter smoothed out, before they are fed into the control algorithm which calculates propeller forces, hence the need for a filtering algorithm. The purpose of the wave filter is to separate the high-frequency wave-induced motion from the motion caused by slowly varying disturbances. Feedback control action must be implemented using the filtered lowfrequency vessel motion, enabling thruster modulation and all related problems to be avoided. Commercially available DPSs feature the well-known Kalman Filter algorithm. The Kalman Filter relies on a previous knowledge of the (linear) system dynamics (in the problem at hand, the low frequency range of the vessel dynamics) and an assumed model for high-frequency motions of the vessel to perform the filtering of wave-induced motion. Actually, as the Kalman Filter is based on a linear mathematical model of the system, all its parameters need to be tuned for a range of operational conditions, to that end it is common practice to use an automatic gain-scheduling procedure.
1. INTRODUCTION Dynamic Positioning System (DPS) can be defined as a computer controlled system designed to keep a floating vessel on a specific position or pre-defined path by using her own propellers. Position and heading sensors, wind sensors, a set of algorithms and propellers are the members of a typical DPS. Those systems have found a wide array of applications in the offshore oil exploration setting, as evidenced by the deployment of such technology in operations such as drilling, pipe-laying, offloading as well as diving support. When in the ocean, floating vessels experience environmental forces that are anything but simple. Their complexity lies not only in the fact they can be seen as the sum of an endless amount of oscillatory components of varying amplitudes, frequencies, phases and directions, their mathematical description is also quite challenging, for each has its own intrinsic characteristics. Nonetheless, their action on floating vessels gives rise to rigid body motions made up of both high and low frequency motions. Wind-generated waves are at the root of large, high-frequency, oscillatory forces and moments that co-exist with their low frequency counterparts originated by wind and current loads, as well as the so-called wave-drift forces, associated with second order effects present in waveinduced forces.
In this work, the authors would like to introduce a novel approach to filtering which is not only model independent, but it is also appropriate for handling signals possessing a high degree of nonlinearity and non-stationarity, conditions which would be at the root of a poor performance of the Kalman Filter. We also intend to demonstrate how this filtering algorithm would perform in a tanker scale model fitted with DPS, where filtering will be performed in real time using our novel approach to filtering: the Empirical Mode Decomposition Filter [Huang et al, 1998], the EMD Filter.
For the design of DPS, it is generally those low-frequency and steady components that must be accounted for, because they are the ones whose physical nature would cause a ship to go astray and drift off its set-point or original path. The highfrequency motions, on the other hand, should not be actuated upon by the controller, since that would require an enormous amount of energy which could not be provided by the enginepropeller system. As a result, we arrive at the conclusion that high-frequency motion components must be suppressed or 978-3-902661-88-3/10/$20.00 © 2010 IFAC
In Section 2, we describe the concepts and ideas underlying the EMD Filter, and present some examples where the EMD Filter is applied to real signals obtained from GPS readings of a marine vessel position. Section 3 provides an overview of the mathematical model of a floating vessel under the 138
10.3182/20100915-3-DE-3008.00071
CAMS 2010 Rostock-Warnemünde, Germany, Sept 15-17, 2010
8.
Repeat iterations until the function become monotonic and there is no more IMFs to be extracted from it. A stopping criterion determines the number of sifting steps to produce an IMF. Huang et al (1998) suggest that the sifting should be stopped as soon the sum of the difference (SD):
external action of current, wind and wave forces, so the reader can have a glimpse of the multitude of parameters which must be estimated in most filtering algorithms. Then, the same reader is informed that the EMD Filter requires no previous knowledge on the system dynamics. Section 4 puts our EMD Filter to the test in an experimental setup where a scale model of an oil tanker fitted with DPS is run in a academic towing tank.
T
∑h
k −1
SD =
The main advantage of this novel approach to filtering, which is still in its early stages, is the fact it requires no previous knowledge regarding the system dynamics (such as inertia, damping, stiffness matrices), nor is it bounded by assumptions of linearity and/or stationarity. The main drawback of the EMD filter is the lack of a solid mathematical grounding; after all it is a heuristic algorithm, which makes it difficult to perform a stability analysis.
(t ) − hk (t )
2
0 T
∑h
2 k −1
(t )
0
2.2 A novel approach to filtering: EMD Filter The EMD claims to properly identify the time scales featuring the physical characteristics of any given system by analyzing a time series representative of the system dynamics. The original signal is then broken down into a set of functions, each displaying an oscillatory mode intrinsic to the system, the so-called Intrinsic Mode Function (IMF). Lower order IMFs retain the most oscillatory temporal modes.
2 THE EMPIRICAL MODE DECOMPOSTION FILTER (EMD FILTER) 2.1 The Empirical Mode Decomposition (EMD) algorithms The EMD algorithm decomposes any given data into a set of intrinsic mode functions (IMFs). An IMF is defined as a function that satisfies two conditions:
Mathematically, the EMD algorithm decomposes the function X (t ) into: N
X (t ) = ∑ hi (t ) + r (t )
1. In the whole data set, the number of extrema and the number of zero-crossings must either be equal or differ at most by one. 2. At any point, the mean value of the envelope defined by the local maxima and the envelope defined by the local minima is zero.
Where r(t) is the residual. Equation (1) contains the central idea behind the EMD Filter. It simply states that the raw signal can be mathematically written as the sum of N IMFs, note that IMFs isolate time scales present in the signal. From this equation, it can be seen that by retaining only the low-frequency terms in the reconstruction above, we obtain a filtered version of X (t ) , for each IMF is a representation of a temporal mode within the raw signal. That way, the first IMF tends to be a highly oscillatory mode, and as we move towards the higher order IMFs, they become less and less oscillatory, in other words, they are low frequency modes. That is the central idea behind the EMD Filter: use the EMD to break the raw signal down into IMFs, and then add up higher order (low frequency) IMFs. What comes out of this procedure is a filtered version of the raw signal, a version that contains only the low frequency temporal modes.
Huang et al (1998) suggests that an IMF represents a single oscillatory mode (not a structural mode, but a temporal mode) which generalizes a harmonic function in the way that its amplitude and frequency are time-varying. The process of extracting IMFs off a given signal is termed sifting, which is defined by the following algorithm: 1. 2. 3. 4. 5.
Between each successive pair of zero crossings, identify a local extremum in the test data. Connect all the local maxima by a cubic spline line as the upper envelope. Repeat the procedure for the local minima to produce the lower envelope. Calculate the envelope mean m1 Calculate the first estimate of IMF 1:
We again point out that the decomposition procedure requires no previous knowledge regarding the system dynamics. The EMD Filter corresponds to an automatic and adaptive (signaldependent) time-variant filter capable of handling both nonlinear and non stationary signals.
h1 = X (t ) − m1 6. Take h1 as the new signal and repeat iterations from 7.
(1)
i =1
1 to 4, until the stopping criteria is met (see below). We then have IMF 1. The first IMF h1 , contains the high frequency components of the original signal. From the relation r1 = X (t ) − h1 , we obtain the residue which is taken as the new signal to be sifted through.
2.3 The EMD Filter: offline filtering algorithm Now in the following example, we shall consider signals obtained from GPS measurements of a full scale oil tanker. The signal analyzed here corresponds to the vessel’s center of gravity calculated from positions captured from 2 GPS. The signals provided were fed into the algorithm and filtered offline, that is, in one single batch, as opposed to filtering it on139
CAMS 2010 Rostock-Warnemünde, Germany, Sept 15-17, 2010
out of the measured signal, thus producing a filtered point at any given time step, as opposed to producing a filtered vector of points (offline problem).
line using a moving window. In Section 2.4, we provide a more detailed explanation on how to interpret the online filtering problem as well as the offline filtering problem on the EMD Filter setting.
Our study has shown that filtering a noisy signal using an online and an offline approach do not lead to the same results. The reason for this can be grasped when we notice that during an offline analysis, at any given time instant, we have full knowledge of past and future values in the time series, in other words we benefit from a “look-ahead” bias, which is not present during an online filtering, where at any given moment in time we have no information whatsoever on future values. That will be the reason why online filtering performance will not be as good as its offline counterpart.
The original and filtered signals for the sway motion are shown in Fig. 1. Note that for an off-line analysis, the EMD filter shows a remarkably good performance without introducing any lags into the filtered signal. What kind of previous information regarding the system dynamics was required for the algorithm to produce such results? None. The method is signal dependent and model independent. Some shortcomings originating from this will be discussed in Section 5. In order to filter the raw signal, we added up IMFs from order N/2 (or the nearest integer to N/2) on.
During an online filtering analysis, at step k, one might wonder what the filtered value to be fed into the control law should be. Our tests have suggested that taking the last point in the filtered vector f fails to produce good results. The authors believe that end effects might be responsible for such poor filtering performance. As an attempt to circumvent such an issue, we suggest the following strategy for online filtering:
As it will be shown later on, in DPS applications the proper approach to EMD filtering should require the use of a moving window, as the filtering is to be performed in real time - the online filtering problem, which we should detail in the next section. -1.2575
x 10
6
Raw signal vs EMD Filtered Signal-SWAY yg filtered yg
1.
Update the noisy input vector s of size n by eliminating the oldest value and including the measured value at instant k; 2. Filter the input signal via EMD Filter. We have vector f of size n (it is important to see that both the filtered and noisy signal have the same number of elements), in our study we adopted n=100. 3. Take the n-j is the ilter delay time; 1 ) value in the f vector where the last observations correspond to the last obtained values, which allows us to move away from the filtered points contaminated by the end effect. 4. Take the quantity obtained in step 2 as the filtered value for j time steps. In order to avoid excessive phase lag j should be around 1. In fact there is a trade-off involved in the choice of parameter j , which is why we suggest the reader to adopt values around 1 . As we mentioned in the 3rd item above, the parameter j allows us to move away from the end-effects, but at the same time, values greater than 1 cause the filtered signal to lag behind the raw signal.
-1.2575
-1.2575
-1.2575
-1.2575
-1.2575
-1.2575
-1.2575
-1.2575
0
100
200
300
400
500
600
700
800
900
1000
Figure 1: Filtered sway motion in red.
2.4 The EMD Filter: online filtering algorithm Unfortunately, notwithstanding the great performance of the off-line filter, it is not feasible for application in DPS because its control loop requires a real time filtering of the signals of the position and heading that are captured at sampling instants k. The natural solution is to adopt a moving window with a pre-determined length in which at every sampling instant the oldest value is eliminated and the new value is inserted in the data vector.
To realize the influence of the parameter j, the spectrums of the filtered signals for j = 1, 2, and 3 are shown in Figure 2. The surge motion of the full scale ship captured by GPS, as previously explained, is taken as the raw signal. It is clear that values greater than 1 for the parameter j, which is associated with the phase lag, lead to better filtering, as evidenced by the spectrum, since we move away from the zone affected by the end effects.
Here it is important to distinguish the online filtering problem from the offline filtering problem, presented on Section 2.3. The offline filtering problem presents us with a noisy data, and demands us to find a filtered (noise-free) version of it. As we have shown in the Section 2.3, the EMD Filter is quite a suitable tool for tackling such problem. For one thing, it produces a smoothed version of a noisy signal without introducing any lag. Moreover, it requires no previous knowledge on the system dynamics.
3 MATHEMATICAL MODEL The modeling of dynamics of a vessel in the horizontal plane is formulated considering an earth-fixed frame OXEYEZE and a vessel fixed reference frame GXYZ, as indicated in Fig. 3. The origin of the vessel fixed frame is located in the centerline, in a distance from the center of gravity. The axes of each body-fixed co-ordinate system coincide with the
The online filtering problem, as seen in the DPS setting, consists of, at time step k, taking high frequency components 140
CAMS 2010 Rostock-Warnemünde, Germany, Sept 15-17, 2010
0 0 -) ) " .) ) /# - ) - ) " -) ) ! - ) - ) ! , (0 0 0 0 .) ) /! Where ) is the mass of the vessel, ) , 1, 1,2,6 are the added mass coefficients in surge, sway and yaw, respectively; ! and " are the surge and sway velocities, respectively;! and " are current speeds related to the body fixed frames, respectively; # is the yaw rate; + is the moment of inertia about vertical axis; %, & and ' represent the total external forces and moment in surge, sway and yaw directions, respectively; dot means time derivative of the variable. 4 5 is the control vector and 4 5 is a matrix defined by the configuration of the actuators. The position and heading of the vessel related to the earth-fixed coordinate system are obtained from the following equation
principal axes of inertia of the vessel. To represent the intricate motion of a vessel at sea the mathematical model is usually split in two terms namely the low frequency motion caused by wind, current and second-order wave forces, and the high frequency motions caused by the first order wave forces. The former results by applying the Newton law and the latter are time series generated from a spectrum of the motion for each degree of freedom.
Spectrum for j=1 Spectrum for j=2
0.15
Spectrum
Spectrum for j=3 0.1 0.05
6 76
0 3
6
0.3
2
0.2
j
;<=9 76 : =19 0
0.1 1
Frequency
0
Figure 2: Spectrum of the surge time series shown in Figure 2, where the parameter j associated with the phase lag is increased.
Figure 3 Earth fixed and vessel fixed references frames
% & '$
) ) 0 ( 0
0 ) ) ) )
0 ) ) , + )
-=19 ;<=9 0
0 0> 1
(6) (7)
?@ |BCD@, E| F@
The equations for low frequency motions, related to the body fixed reference frame, and considering added mass forces, are given by: ! " #$
9 $
where , 8 and 9 are the coordinates and heading of the vessel in the earth fixed frame respectively. The forces %, & and ' are due to the action of current, wind, and yaw hydrodynamic damping. Forces due to current and wind can be determined using drag equations (Leite et al., 1998; OCIMF, 1994). The second order wave forces are caused by the so-called slow and mean drift forces and they can be estimated as presented in Aranha (1994), and Aranha and Fernandes (1995). The high frequency motions (denoted as 6 in what follows) cannot be easily predicted by using differential equations because the calculations of the first order wave forces in real time are complex numerical processes. In light of this a time series generated from the high frequency motion spectrum is used to perform numerical analysis of the dynamics of the vessel. The spectrum for each degree of freedom is given by:
.
8
(5)
(2)
where RAO is the ship Response Amplitude Operator, β is the wave incidence angle related to the ship and S (ω ) is the sea spectrum.
(3)
3.1 Measurement System
(4)
For conventional ships, only position and heading measurements are available. Position is normally obtained from a DGPS while heading is usually measured by using a gyro compass. Hence, the position and heading measurement equation can be written as: G 6 6 H
(8)
Where d ∈ ℜ3 is a zero mean Gaussian white measurements noise.
141
CAMS 2010 Rostock-Warnemünde, Germany, Sept 15-17, 2010
The performance of the EMD filter was checked considering step variation in the reference values for the coordinates of the axes OXE, and OYE, and the heading 9, for different values of the parameter j. The test was performed adopting sampling period of 0.1 s. The model was disturbed by regular waves characterized by frequency of 1,0Hz and 1cm high. However, in this paper, only the influence of the parameter j in the performance of the EMD is pointed out. For that, results for surge motion for j=1, 2 and 3 are shown in Figure 4. Note that the EMD Filter in its online version outlined in Section 2.4 was capable of performing a real time filtering of the raw signal being measured. Besides, it also managed to properly capture the variation of the motion if the setpoint is altered, which higlights the adaptive nature of the EMD Filter. Again, we stress the fact we did not provide the filter with any prior information on the test model dynamics nor the wave characteristics, and yet the EMD Filter performed considerably well in all 3 conditions.
3.2 Conventional Filtering of the High Frequency Motions The high frequency motions of the vessel should be filtered because their averages are zero and an attempt to counteract them would require a huge power, leading to extra fuel consumption and propulsion degradation. The first generation of the SPDs were implemented by using notch filter to reduce the amplitude of the high frequency motions. The second generation of the SPDs have adopted Extended Kalman Filter that assumes an internal linear model based on equations (2) and (5) for estimating low frequency positions and velocities. For filtering of the high frequency motions a state-space realization of a transfer function that satisfies the following approximation is included: ?@ I |J@|
The transfer function for each direction is given by:
G( s) =
η w ( s) w(s)
=
σs s 2 + 2ξω0 s + ω02
Granted, in order for the EMD Filter to start, it does require an initial window of acquired data. In our tests, we collected 100 points during 10 seconds, and once we had our initial data window, the EMD Filter could be set off. As time goes by, the 100-point-window is continously updated, as it moves forward in time discarding older readings and incorporating new ones, as explained in Section 2.4. Figure 4 also brings Figure 2 to the time-domain, so we can see how larger values of parameter j lead to both better filtering (as evidenced by Figure 2) and – what is not evident from Figure 2 – to the introduction of a slight phase lag, the trade-off brought up in previous sections.
As a result, the Kalman Filter works with state vector of order 12 to filter the high frequency motions and estimate the low frequency motions and velocities. Therefore it is clear, from what we explained so far, that the application of the EMD Filter to solve the filtering problem in DPS would be a very interesting alternative indeed. For one, there would be no need to estimate dynamical properties of the vessel such as mass, added mass, damping coefficients and others which are not always easily obtained. Besides, should the sea state changes, the EMD Filter would automatically adapt, thus eliminating the need for a cumbersome gain adjustment procedure, which troubles the Kalman Filter.
As for the performance of the DPS system as a whole, we see in Figure 4 that the introduction of the EMD Filter does produce good results, as the test model manages to track the reference signal. A better adjustment procedure in the PID controller gains would have led to even better results. However, during the tests it was realized that if the input raw signal is not noisy enough (if it flatlines or oscillates at very low frequencies), the EMD procedure breaks down or iterates for too long until the stopping criterion is met. So in DPS applications, the more the floating vessel moves about the better the filtering will be. On the other hand, in case of a violent storm, our EMD Filter would quickly adapt to the external conditions and the filtering algorithm would carry on unflappably
4 APPLICATION OF THE EMD FILTER TO DPS In order to check the viability of the application of the EMD Filter in DPS using the online filtering procedure proposed in Section 2.4, a test with scale model was performed in the tank of the Department of Naval Architecture and Ocean Engineering of the University of São Paulo. The model corresponds to a typical offloading vessel and the actuators are the bow thruster, the stern thruster and the main propeller. The scale model main particulars in full load conditions and details of the laboratory facilities are described in Morishita et al (2009). A conventional PID (Proportional – Integral – Derivative) controller for each degree of freedom was used.
5 CONCLUSIONS
To perform tests with EMD Filter three vectors were considered, namely K 4 5 , L 4 5 and M 4 5 where K, L and M are vectors composed of measured data obtained from G (eq. 8) by using the procedure suggested above. At every time N those vectors are updated and the EMD filter is applied for each vector. Then the PID control law is calculated for each degree of freedom using the filtered signals. The parameters of the controller were tuned by means of genetic algorithm considering ! " 0 in equation (2). The computer code to perform the control law is based on Matlab but EMD algorithm is a C compiled program to save computer time.
This work presents a new filtering algorithm for application in DPS. We named it the EMD Filter, for the filtering strategy we created relies on a prior step where the raw signal is decomposed into Intrinsic Mode Functions by a procedure that came to be known as EMD proposed by Huang et al (1998). The filtering algorithm we created is capable of filtering out high frequencies of a raw signal by retaining only the low frequency Intrinsic Mode Functions (IMFs) [Huang et al.,1998] . We came forward with an offline filtering algorithm, which helps us solve the problem where given a raw signal, we are asked to produced a smoother, filtered version of it. And, we also devised an online filtering
4.1 Model test results 142
CAMS 2010 Rostock-Warnemünde, Germany, Sept 15-17, 2010
Workshop on Nonlinear Signal and Image Processing NSIP,v.3,pp.8-11.
algorithm, where the extraction of high frequency components happens in a real-time fashion. So given a moving window of fixed length containing noisy observations, we can provide at any given time step, a filtered point (as opposed to a filtered vector of points), that will be used in the control law. Both the offline filtering and the online filtering problem and its related algorithms are presented in detail in Section 2, where we also discuss the inevitability of introducing a phase lag in online filtering in order to escape the end effects, and the trade-off arising out of it.
(a) j = 1 -2.75
-2.8
surge [m]
-2.85
Perhaps, the most important claim we stake throughout this work is that the EMD Filter algorithm is model-independent, as it requires no information on the system dynamics. Unlike the Kalman Filter, where a multitude of gains and parameters need to be tuned (and readjusted once the external conditions, such as sea state, changes), the EMD Filter is also adaptive, as we have shown in Section 4, in a scale model run in a water tank.Some shortcomings found by the authors during the investigations were: 1) the input signal need be noisy to the EMD filter to work properly. 2) it is not yet possible to introduce a cutoff frequency in the EMD Filter design.
-2.9
-2.95
-3 70
75
80
85
90
95
100
100
110
120
75
80
85
time [s]
(b) j = 2 -2.7
The EMD Filter we propose in this work is still in its early stages. There certainly is room for further improvements and investigation on the algorithm. And the search for new fields of application of this inovative technique looks quite promissing.
-2.75
surge [m]
-2.8
ACKNOWLEDGMENTS The authors thank Petrobras for financial support in implementing the project. The second author thanks the CNPq, under the processes 558161/2009-0.
-2.85
-2.9
REFERENCES
-2.95 60
70
80
90
time [s]
Aranha, J.A.P. (1994) A Formula for Wave Damping in the Drift of a Floating Body”, Journal of Fluid Mechanics, vol. 272, pp.147-155. Aranha, J.A.P., Fernandes, A.C., (1995), Second-order slow drift force spectrum, Applied Ocean Research, v.17,n.5,p.311-313. Huang, N. E, Shen, Z., Long, S.R., Wu, M.C., Shih H.H., Zheng, Q., Yen, N., Tung, C.C., Liu, H.H. (1998) The empirical model decomposition and the Hilbert spectrum for nonlinear and non-stationary time series analysis, Proc.R. Soc. Lond. A454, 903-995. Leite, A.J.P., Aranha, J.A.P., Umeda, C. and de Conti, M.B.(1998), Current Forces in Tankers and Bifurcation of Equilibrium of Turret System: Hydrodynamic Model and Experiments, Applied Ocean Research, Vol. 20(3), pp. 145-156. Morishita, H.M., Tannuri, E.A., Saad, A.C., Sphaier, S.H., Lago, G.A., Moratelli Jr, L., (2009) Laboratory Facilities for Dynamic Positioning System, Proc. International Conference on Manoeuvring and Control of Marine Craft, Guarujá. OCIMF. Predictions of wind and current loads on VLCCs, (1994) Oil Companies International Marine Forum. Rilling,G.,Flandrin, P.,Gonçalvès,P., (2003), On Empirical Mode Decomposition and its Algorithms,IEEEURASIP
(c) j = 3 -2.8 -2.82 -2.84
surge [m]
-2.86 -2.88 -2.9 -2.92 -2.94 -2.96 -2.98 55
60
65
70
time [s]
Figure 4: Results of test on scale model DPS - surge motion (a) j=1; (b) j = 2; (c) j = 3 (blue:raw signal, black: reference signal, red: filtered signal)
143