Thorsten Johl
email: thorsten.johl@helmholtz-hzi.de

Manfred Nimtz

Lothar Jänsch

Frank Klawonn

Detecting Glycosylations in Complex Samples

Keywords: proteomics, glycosylation, mass spectrometry 1

Glycoproteins are the highly diverse key element in the process of cell cell recognition and hostpathogen interaction. It is this diversity that makes it a challenge to identify the glyco-peptides together with their modification from trypsin-digested complex samples in mass spectrometry studies. The biological approach is to isolate the peptides and separate them from their glycosylation to analyse both separately. Here we present an in-silico approach that analyses the combined spectra by using highly accurate data and turns previously established knowledge into algorithms to refine the identification process. It complements the established method, needs no separation, and works on the most readily available clinical sample of them all: Urine.

Introduction

Glycosylations are complex post-translational modifications of proteins and the resulting glycoproteins are often found as part of the cell membrane where the glycosylations extend from the plasma membrane into the intracellular space and form the glycocalyx. Due to a high variability in molecular structure, they play a vital role in cellcell interactions and immune responses [START_REF] Marino | A systematic approach to protein glycosylation analysis: a path through the maze[END_REF] comparable to a key that fits only a specific lock, and several pathogens exploit this system to gain entry into specific cell types. Other roles include helping in protein folding, as transport molecule, lubricant or providing frost resistance. Current estimates state that about half of the eukaryotic proteome is glycosylated, making it the most abundant of modifications and one of the most versatile. Glycosylations most commonly attach either to the asparagine (N) amino acid via nitrogen, or to serine (S) or threonine (T) via oxygen and are hence called N-or O-linked glycosylations. These two types are made up of the same basic mono-saccharides, which can be broken down into four groups of identical masses (Fig. 1 A). These basic components are arranged into link type-dependent core structures of highly variable overall complexity (Fig. 1 B). This complexity hampers mass spectrometric (MS) analyses of glycosylated samples, because MS relies on a known total peptide mass for identification before it compares the most intense signals of a fragmentation pattern (MS²) of this peptide to a number of theoretical signals produced by known peptides of similar mass from a reference database [START_REF] Perkins | Probability-based protein identification by searching sequence databases using mass spectrometry data[END_REF]. The best fitting spectrum is then assumed to be the correct interpretation. But these theoretical spectra do not contain any modification unless each is specified, and to predict all possible glycosylation patterns would soon reach incalculable complexity.

Fig. 1. Overview of saccharide makeup und structure

The general accepted solution to this problem is to separate the oligo-saccharide side chain from the peptide by using endoglycosidases specific to N-linked, and strong bases for both side chains [3] [4]. The resulting components are then analyzed separately, with a single mono-saccharide remaining attached to the peptide, which identifies the modified amino acid in MS experiments. This method works best on samples of low complexity where it is clear which oligo-saccharide originated from which peptide, but it is unsuited for large-scale proteomic studies. Yet glycosylations are ideally suited to be analyzed by MS workflows. First of all, they are readily ionized and thus result in good signal to noise ratios in the mass spectra. Secondly, the mono-saccharides produce distinct signals in the low mass region of such a spectrum. These identify spectra of glycosylated peptides with absolute certainty. And thirdly, the signals stemming from amino acid and saccharide residues in conjunction have a characteristic atomic mass that sets them apart from pure amino acid residues of similar mass. The reason for this is the greater amount of oxygen contained in saccharides which weighs in at just under sixteen Dalton, while every other atom prevalent in amino acids has a decimal value just above a full Dalton (C,H,N,S) [START_REF] Dodds | Enhanced peptide mass fingerprinting through high mass accuracy: Exclusion of non-peptide signals based on residual mass[END_REF]. It will be these criteria that are exploited by the program introduced herein in order to identify the peptides and their glycosylations from the same spectrum. This paper will demonstrate the feasibility of this new approach that does not rely on a reference database for glycosylated peptides like UniProt [14] or Glycosuite [START_REF] Cooper | GlycoSuiteDB: a curated relational database of glycoprotein glycan structures and their biological sources. 2003 update[END_REF].

Methods

Sample preparation and acquisition

3x 2.5ml urine were desalted in turn with a PD-10 column, united, and then concentrated to 400µl with a speedvac and prepared according to the "glycoprotein isolation kit, WGA" from Thermo Scientific with minor adjustments. Incubation time was increased to 30min, the elution time to 15min. Subsequent digestion was performed over night with 20µg trypsin. Two MS measurements were performed in turn, one injecting 1µl onto the loading column of a Dionex Ultimate 3000 HPLC, the other 3µl. The samples were separated over 60 min on a reverse phase LC at 350µl/min flow rate with 80% Acetonitrile as mobile phase. Ions were introduced into the Thermo Fisher MS Orbitrap Velos device via an electrospray ion source. The FTMS Orbitrap was used as the detector for all scans. Survey scans were performed at a resolution of 60000 in the range of 700 to 2000 and the 5 most intense signals were chosen with an isolation width of 4 Da for MS² fragmentation in the HCD collision cell at normalized collision energy of 40 for 40ms. The minimum signal threshold was set at 10000 and the default dynamic exclusion list was enabled, removing signals for 90 sec from the selection process after having been measured twice. Concluding MS² scans were performed at a resolution of 7500 and a range of 100 to 2000. All preprocessing steps like baseline estimation and noise removal were performed by Xcalibur 2.7 at the time of acquisition with the device inherent parameters.

2.2

Program setup

The program was created using Java programming language from Oracle and Eclipse as programming environment. Readw.exe [6] from the trans-proteomic pipeline [START_REF] Pedrioli | Trans-proteomic pipeline: a pipeline for proteomic analysis[END_REF] was used in conjunction with files from ThermoFinnigan for centroidization and to convert Thermo Fisher RAW files into mzxml. Mzxml in turn was converted to mgf using JRAP [START_REF]JRAP Extended[END_REF]. Alternative mzML parsing is provided by jmzML [START_REF] Jmzml | [END_REF]. Settings are stored as xml file using the Apache commons parser [START_REF]The Apache XML Project[END_REF]. The program requires the JRE and has been tested on Windows 32 and 64bit as stand-alone application. Note that the readw.exe only works on 32bit windows operating systems, and thus 64bit machines can only process mzxml, mzml and mgf files.

Spectra preparation and identification of glycosylated spectra

All signals of a deconvoluted spectrum were checked for isotopic C13 satellites allowing for 20ppm mass tolerance, and the appropriate charge state was deduced from the distance between the isotopes. All following steps were then performed on the mono-charged mass of the most intense isotopic signal. All spectra were searched in turn for any two of the following marker oxonium-ion masses: dehydro-N-acetyl hexose (204.079 Da), a dehydro di-saccharide of N-acetyl hexose and a hexose (366.132 Da) or di-dehydro sialic acid (274.085 Da). The search was performed with a mass tolerance of 20ppm and a minimum marker signal intensity of 10% of the most intense signal in the spectrum. Every signal that was found during identification and, if present, an additional mono-dehydro sialic acid (291.095 Da) signal were marked as signals of saccharide origin. All spectra that met the identification criteria were then considered glycosylated and passed on to the following processing steps.

Signal makeup

Any mass from top to bottom, beginning with the total mass of the modified peptide, was then considered as the starting point of a tree of saccharide masses. In order to qualify as such a seed, all mono-saccharide masses from Fig. 1 A where in turn checked to find a smaller signal in the estimated mass range of starting signal minus saccharide mass. The steps are repeated from these seeds onward until no further signals are found that can be explained by the loss of a saccharide. There are additional criteria that influence the signal selection process:

Firstly, the peaks are only considered if the mass delta between the actual and expected decimal fractions of a peptide (Mi) [START_REF] Dodds | Enhanced peptide mass fingerprinting through high mass accuracy: Exclusion of non-peptide signals based on residual mass[END_REF] increases from higher peak (Pi) mass to lower peak mass (Pi+1), denoting the loss of an oxygen-rich monosaccharide from the remaining molecule. The logic allows for a small mass-dependent inaccuracy (Δ) of 15ppm [START_REF] Marino | A systematic approach to protein glycosylation analysis: a path through the maze[END_REF].

( ) { (1)
Secondly, the program allows for a shift in the most abundant isotope form, searching for a second peak at +/-1.00335 Da distance with 15ppm tolerance if no primary peak can be found. Such shifts are marked, and block a further shift into the same direction, which is highly unlikely given the standard distribution of C13 isotopes. It may, however, shift back to the original isotope pattern once to accommodate for inaccuracies in the measured intensities.

Thirdly, signals without an isotope pattern that denote the charge state are considered to be of any charge state up to the charge state of the mother mass. Such signals are called jokers, and may be assigned only one charge during the enlargement of the saccharide trees. Signals that pass all criteria and thus carry a glycosylation are marked and the saccharide tree is stored for further analyses.

Typecasting spectra

The kind of glycosylation found in any given set of spectra can be determined depending on the completeness of the saccharide tree. Four different core structures as introduced in Fig. 1 B and three further N-linked sub types are known to be expressed in human tissue. At least four saccharides from the start of the tree are necessary to discern between N and O-linked glycosylations. The more complex N-linked core structure is checked first. This is done by a pattern recognition that performs a depthfirst traversal of the saccharide tree and checks first if there are one or two leading HexNAc saccharides, and then if there are three consecutive hexoses following that.

The traversal for any sub-tree is aborted if the child does not match expectations, and the next child is considered until the pattern is found, or no children are left. There is some allowance for a single extra saccharide attached to the leading HexNAc, and either the last hexose or the first and its connected HexNAc may be substituted by a Hex-HexNAc di-saccharide.

If the algorithm fails to find the key elements for an N-linked core structure, it checks the less complex O-linked core structures in the following order: 2 HexNAc and a hexose, 3 HexNAc, 2 HexNAc, HexNAc and a hexose. The leading HexNAc may be missing in all cases.

Determining true peptide mass

All signals that have lost a saccharide are removed from the spectrum. The remaining signals are then considered from top to bottom. If they have a relative intensity of 5% of the most intense signal ever present in the spectrum, the algorithm tries to fit any combination of saccharide masses into the mass delta of this signal and the total mother ion mass. This is achieved by beginning with just one saccharide, trying each in turn before adding a second saccharide and checking all possible mass combinations of these. The method stops adding more saccharide elements once it has found one combination that fits, or multiple combinations of the smallest mass have become larger than the mass delta it was trying to fit. Signals that can be explained by the loss of a combination of saccharides are stored. All spectra are then grouped by their mother mass and charge state, and the stored signal intensities are added up. The four most intense combined signals are then chosen in turn and proposed as the true peptide mass and new mother mass of the spectrum retaining the charge the spectrum originally carried, resulting in up to four new spectra for each that was identified as containing saccharides. The spectra are then stored separately and can be submitted to the Mascot server for identification.

Mascot search

The Mascot search (V 2.3.02, Matrix Science) was submitted via Mascot Daemon against a UniProt database (Release 2012_02) restricted to human proteins with a peptide mass tolerance of 20ppm, a fragment mass tolerance of 0.4 Da and no missed cleavages. Mother ion charge was limited to up to 4 charges. Allowed modifications were Carbamidomethyl (fixed) and Oxidation (M) (variable).

Results

8615 and 9577 spectra, respectively, were processed from the 1µl and 3µl samples. About a third of these were identified as glycosylated and resulted in more than 11000 proposed spectra per sample (Table 1). Roughly 34% of all spectra were found to carry a glycosylation. Only 5% of these could be identified as carrying an N-linked glycosylation tree. The majority is estimated to carry the less complex O-linked glycosylation. Less than 20% resulted in an unknown type of glycosylation. The final number of spectra proposed to Mascot by the program is nearly equal to four times the amount of glycosylated spectra, meaning that almost every spectrum had enough signals that qualified it as one of four new possible true peptide masses. Only 13 spectra in total could not be assigned a new mass and thus resulted in no new spectra (not shown). The results of the identification by Mascot are summarized in Table 2. The identity score cut offs for the 1µl and 3µl samples were 19 and 20, respectively, and the significance (homology) threshold was a p-value of 0.05. The false discovery rate (FDR) of Mascot was 5.13 (7.5) above identity (above homology) for the 1µl sample, and 1.83 (5.17) in the 3µl sample. The table shows all proteins that were found in both samples by Mascot with their respective score. Shown next are the number of matching spectra that were associated with this protein in total and the number of unique sequences. The affirmed glycosylations are the number of unique sequences already known and found in this analysis. The number in parenthesis for all three groups is the number of hits above homology. The numbers of predicted and known glycosylations are taken from the UniProt entry and show the total number of listed glycosylations followed by the number of confirmed ones in parenthesis. Asterisks mark entries where more than one glycosylation occurs on the same tryptic peptide and the location is not unambiguous. Values were taken from the protein overview tab of Mascot and UniProt Release 2012_02.

Sample Original Spectra

Glycosylated

Nlinked

Olinked

Unknown

EGF 140 125 9 (1) 4(3) 2(1) 3(2) 4(3)
) 0 0 3*(1) YIPF3 102 89 6 (6) 3(3) 2(1) 1(1
PGCB 107 21

CF072 26 48 2 (2) 2(2) 1(1) 1(1) 1(1) 1(1) 2(1)
VASN 18 52 1

IC1 46 22 5 (1) 2(2) 1(1) 1(1) 0 0 5(3)
P3IP1 25 29 4 (2) 4(1) 1(1) 1(1) 2*(2) 2*(2) 15(0) GOLM1 30 30 4(2) 4(3) 2(1) 1(1) 0 0 3(1)
OSTP 27 23 5

Table 2. Comparison of identified and known glycosylations

The 3µl sample proved to be the more reliable and was consecutively used for further analyses of the saccharide trees and glycosylation sites. 7 of the identified 19 proteins were detected with at least one peptide with a known glycosylation site above homology. All of these sites, except the one from CF072, were related to Oglycosylations. All of them were identified by the program as O-glycosylated at least once per sequence. Three spectra are shown in Fig. 2. They exemplify different O-glycosylations that were correctly identified by the program, and demonstrate the potential to reliably identify the glycosylation site from glycosylated spectra. The notable exception is CF072_HUMAN, which is predicted to have an N-linked glycosylation at the identified site 191, but was identified as O-type. The identified spectra in question only contain unlinked Hex-HexNAc and Hex, and thus do not permit the direct identification of an N-linked glycosylation, although the peptide could be identified with an adjusted mono-charged mother mass of 1461.77 Da (down from 1206.06, z=2). Fig. 2. Saccharide tree and peptide y-ion series (A) Is taken from P3IP1_HUMAN. A consecutive saccharide tree links the total mass of the mother ion over three steps with the true peptide mass of 1170.58 Da. Three further saccharide trees were found that stem from different dissociation events. One of them is linked directly to S39, the reported glycosylation site. Spectrum (B) shows a glyco-peptide that was associated with site T183 of Protein X3CL1_HUMAN. The lower half of the spectrum shows several saccharide ions from the b-series, in one case extending the mass 274 used for identification by NeuNAc and HexNAc. The true peptide mass 1346.67 Da is again linked to the total mother mass. (C) Shows two saccharide trees in the lower mass spectrum. They originate from the same fragment, which dissociated from the peptide and were found complementary forward and backward. Note the double charged di-saccharide that was found and extends the mono charged saccharide that is attached to 1252.67 Da. Neither the total mother mass nor the actual site is directly attached to a saccharide. It was nonetheless identified as IGF2_HUMAN, which could be modified either at T96 or at T99. The HexNac-Hex tree connected to proline (P) suggests that it is site T96. Spectra were taken from Thermo Fisher Xcalibur -QualBrowser [13].

Discussion

This paper describes a new approach to analyze peptides together with their glycosylation. As the current standard is to analyze both separately, there is just one other project that uses a similar approach [START_REF] Neue | Straightforward N-glycopeptide analysis combining fast ion trap data acquisition with new ProteinScape functionalities[END_REF]. That approach uses the ProteinScape 3.0 and GlycoQuest programs in conjunction. It compares the unadjusted glycosylated spectra with a list of theoretical spectra of dissociated glyco-peptides, very similar to the database-based approach performed by Mascot. Although this method is less likely to provide false positives, it will not be able to identify new glycosylations. Indeed the high number of predicted glycosylations apparent from Table 2, especially with regards to N-linked glycosylations, demonstrates there is need for an approach similar to de-novo sequencing for peptides. The example of an unidentified N-linked glycosylation described above demonstrates the problems when identifying large saccharide trees in combination with large peptides in mass spectrometric studies. The variety of possible dissociation events and multiple fragmentations within the tree means that no specific fragment reaches a detectable quantity in the highly accurate, but not as sensitive Orbitrap detector [START_REF] Scigelova | Orbitrap Mass Analyzer -Overview and Applications in Proteomics[END_REF]. This results in detectable amino acids and mono-saccharides, but no larger consecutive saccharide elements. Peptide identification itself is based largely on correct total mass and some preferred dissociation locations like proline. The situation may be improved by allowing for di-and tri-saccharides during peak analysis, however, this could decrease the specificity by increasing the chance of false positives. This problem may, in turn, be addressed by using collision-induced dissociation (CID), which generates only two fragments per ion and thus can generate a more gradual distribution of signals. The results may be even better when coupled to the more sensitive but less accurate ion trap detector. On the other hand, this would produce even less pure amino acid ions, resulting in less ion intensity. As is, many of the 3000+ original spectra did not lead to peptide identification, even when allowing four possible interpretations of the remaining non-glyco signals. The reason for this may be found in the more readily ionized saccharides. The remaining amino acid part is less likely to be ionized, provides weaker signals when fragmented, and eventually vanishes into the background noise. A possible solution to this is to perform a third dissociation and detection cycle (MS³) of the most likely true peptide mass. The evident workflow for this would be as follows:

After the first data acquisition cycle, the peptides would have to be analyzed by the program and the mother mass, estimated true peptide mass and retention time recorded. A second data acquisition cycle would follow the first, using the recorded data as an inclusion list for MS³ experiments. The results would then be sent to Mascot for identification. It is, however, presently unclear whether such a general approach is possible and would result in spectra of the required intensity.

As demonstrated above, the basic principle works. The program uses a shotgun approach to reliably identify O-linked glycosylations and their sitebut further development is necessary to extend this to N-linked glycosylations, which is already underway at the time of this writing.

The program is available from the author upon request.

Table 1 .

 1 Overview of glycosylated Spectra

	Spectra


Jacek Kabziński
email: jacek.kabzinski@p.lodz.pl

Fuzzy friction modeling for adaptive control of mechatronic systems

Keywords: fuzzy modeling, adaptive control, motion control, friction compensation

We discuss several fuzzy models to approximate friction and other disturbances in mechatronic systems, especially linear and rotarional electrical drives. Some methods of experimental identification of disturbance forces are presented. We consider several fuzzy models to compromise between model accuracy and complexity. Fuzzy model is used in an adaptive control loop. Several adaptive control algorithms are discussed and the influence of fuzzy model accuracy on the system performance is investigated.

Introduction

It is well recognized that the presence of friction often destructs a performance of a precision motion control systems, especially servo drives realizing tracking tasks. The friction phenomenon is rather complicated and not yet completely understood, so existing friction models are also far from universality and accuracy. In this contribution we propose to connect fuzzy modeling with adaptive control. This approach allows to connect the simplicity of static friction models with the accuracy offered by adaptation to changing conditions, such as a lubricant temperature for example. As the experimental information about friction is usually corrupted and inaccurate we believe that using a flexible fuzzy model connected with adaptation of its parameters is an effective approach.

We consider the motion dynamics given by

v dt dx  ext friction e F F F dt dv m    . ( 1
)
where m is a forcer mass, F e is a thrust force, F ext is external force, load (usually constant or slow-varying) and F friction represents all kinds of friction forces. The mover speed is v and position x. Although equation ( 1) is written according to linear motion convention, it may be also used for rotational movement description if we read m as a moment of inertia and consider torques instead of forces.

In this paper we shortly present the basic friction models and discuss the problem of experimental acquisition of the friction force data. We propose a TSK fuzzy model-based friction estimation structure that can be used for real-time nonlinear friction identification. We introduce a procedure to automatically decide the fuzzy model rules and starting parameters according to the desired modeling accuracy. Finally we apply friction fuzzy model in adaptive backstepping control assuring the position tracking stability without exact knowledge of all plant parameters, including the control gain coefficient. The presented contribution may be placed among many other concerning fuzzy adaptive control in presence of friction [START_REF] Wang | Extraction and Adaptation of Fuzzy Rules for Friction Modeling and Control Compensation[END_REF][START_REF] Wang | Modeling and control compensation of nonlinear friction using adaptive fuzzy systems[END_REF][START_REF] Lin | Stable Adaptive Fuzzy Control with TSK Fuzzy Friction Estimation for Linear Drive Systems[END_REF], but it develops a new and simpler (than for example in [START_REF] Wang | Extraction and Adaptation of Fuzzy Rules for Friction Modeling and Control Compensation[END_REF]) fuzzy model construction procedure and investigates new adaptive control approach.

Friction Models

Several models where proposed for friction forces. An excellent review is provided in [START_REF] Wojewoda | Hysteretic effects of dry friction: modelling and experimental studies[END_REF]. As we claim that an approximated model should be connected with adaptive control approach, we mention only basic ideas here. Usually it is assumed that friction forces are speed dependent and are roughly described by the formula

                 s v v c s c friction e v g Bv v sign v g f f f F ) ( , ) ( ) ( ) ( . ( 2
)
where f s is the level of static friction , f c is the minimum level of Coulomb friction v s is the lubricant parameter (so called Stribeck velocity), B -viscous friction parameter and δ is an even constant. The function g(v) describing a characteristic of the Stribeck curve is only one of possibilities -several other are reported [START_REF] Wojewoda | Hysteretic effects of dry friction: modelling and experimental studies[END_REF]. All parameters of this model are unknown and should be determined by empirical experiments, and still the model accuracy is doubtful. The simplified version of (2) takes into account only Coulomb and viscous friction:

Bv v sign f F c friction   ) ( . (3)
So called LuGre [START_REF] Wang | Extraction and Adaptation of Fuzzy Rules for Friction Modeling and Control Compensation[END_REF] dynamic friction model is supposed to capture most of the real friction behaviour, like Stribeck effect, hysteresis, spring-like characteristics, varying brake-away force. It is based on 'elastic bristles' model of contact surfaces. The average deflection z of the bristles is given by

) ( ) ( v g z v v t z    . ( 4
)
where g(v) is a positive function. To describe Stribeck effect g(v) is usually chosen as

                        s v v c s c e f f f v g ) ( 1 ) ( . ( 5
)
Friction force is given by

Bv z z F friction       . ( 6
)
where σ is the equivalent stiffness coefficient and τ is the equivalent damping coefficient of bristles. Several another models (more complicated, with bigger number of parameters and more difficult to identify) of friction forces are reported in literature [START_REF] Wang | Modeling and control compensation of nonlinear friction using adaptive fuzzy systems[END_REF][START_REF] Lin | Stable Adaptive Fuzzy Control with TSK Fuzzy Friction Estimation for Linear Drive Systems[END_REF]. As it follows from the above discussion friction and ripple forces are of very complicated nature, difficult to analyse and to model. In this paper we suggest modelling the sum of ripple and friction forces by a fuzzy inference system.

Acquisition of the Data for Friction Modeling

It is necessary to conduct some experiments to collect the data for the fuzzy model training.

One of possibilities is so called constant speed test. If we are able to produce a constant speed movement, it means that all the forces are balanced. If we are can measure or estimate the external force, calculate the thrust force (from measurement of motor currents for example), we are able to estimate the friction force.

Sporadically it is possible to apply a constant external force (from an another drive, or from a gravitational load), while the thrust force is zero. In this case we may try to tune the friction model parameters by curve fitting comparing measured position history with numerical solution of equation [START_REF] Wang | Extraction and Adaptation of Fuzzy Rules for Friction Modeling and Control Compensation[END_REF].

Both above methods are theoretically straightforward but difficult to implement in practice. Another possibility is to use a simple observer described by:

) ( ) ( 0 0 v v F dt d v v K F F v dt d m est est fric est est fric est        . ( 7
)
where K and Γ are design parameters,

m m m    0 and F F F F ext e   
                                              dt dF v K dt dv m F t e t e K t e t e dt d m m friction F v F v ) ( ) ( 0 1 ) ( ) ( 1 0 0 . ( 8
)
As we see error dynamics is described by a linear system with disturbances. The eigenvalues s 1 , s 2 of this system are connected with design parameters:

) ( , 2 1 2 1 s s m m K s s m m          . ( 9
)
and so we may choose values of s 1 , s 2 to obtain desired observer dynamics. We may tune observer parameters v  on the estimation error and to plan measurements properly. Special care must be taken to minimize v  , as it is multiplied by K in (7) and its influence cannot be decreased by increasing K.

As we conclude from the above discussion the obtained triples (position -velocity -estimated friction), denoted by

    m k f x v x x k k k k k ,..., 1 , , 2 , 1     . ( 10
)
will be corrupted by estimation method error and subject to estimation/measurement noise and outliers. We will develop special procedure to extract fuzzy rules to construct Takagi-Sugeno-Kang fuzzy model. Figure 1 presents about 300 triples of the data collected from an exemplary linear permanent magnet motor.

Fuzzy Model Construction

The proposed method of fuzzy friction modelling is based on One-dimensional Linear Local Prototypes (1dLLP) approach proposed in [START_REF] Kabziński | One-Dimensional Linear Local Prototypes for Effective Selection of Neuro-Fuzzy Sugeno Model Initial Structure[END_REF]. First we have to recognize if position is an irrelevant input or not. We consider two single-input fuzzy models (SIFM) described below:

 input -x i (i=1 -position, i=2 -velocity), output -c i ,  input linguistic categories:   m k x IS x k i i ,..., 1 ,  ,  membership functions:  i k i k b x x x a , , ( )          1 1 2 , ( 11
)
 rules:  . ,..., 1 0 , 0 0 0 , , , , , , , , , , , , m
k x if f q p x if q x f p q x p c c THEN x IS x IF k i k k i k i k i k i k i k k i k i i k i k i i k i i           . ( 12
)
The action curve given by the output c i of this system for the input data -x i,k generalises information coded by

k k i f x  , i=1
,2 and the degree of this generalisation depends on membership function parameter a . Recommendations for the choice of a and b are given in [START_REF] Kabziński | One-Dimensional Linear Local Prototypes for Effective Selection of Neuro-Fuzzy Sugeno Model Initial Structure[END_REF]. The shape of each action curve is robust to outliers in the measured data and to the measurement noise. If the i-th input is inessential the curve generated by corresponding SIFM will be flat, if it is meaningful the curve will cover significant part of the range of

    m k f k ,..., 1 , 
. Fig. 2 depicts action curves for position and speed for the data presented in fig. 1. Its visible that in this case position was the irrelevant input for friction modelling.

Selection of membership functions for each input is based on piece-wise linear approximation of action curves derived above. Uniform or mean-square approach are both applicable. As the result of piece-wise linear approximation for the i-th significant input we obtain m i linear local prototypes (LLP) defined on intervals

    i j i j i j i m j x x I , , 2 , 1 , , max , min ,    ( 13
)
by linear polynomials

    i j i j i i j i i j i i j i m j x x x p x p x P , , 2 , 1 , , ) ( , max , min , 0 , 1 ,
     (14)
for each interval. The design parameter  , which defines the approximation accuracy, governs the number of linear local prototypes. For each j we construct a bell-shaped functions

) ( , x j i  spanned over j i I
, and centred at the middle point of j i I

, . The choice of the third parameter b i,j is arbitrary -usually 1.5< b i,j <2 gives good results. The rules for the proposed TSK fuzzy model will be:

 , ,..., 2 , 1 , , 0 , 1 , , i j i i j i j i i j i i m j p x p c c THEN IS x IF      (15)
where starting values of parameters are taken from piecewise linear approximation results (14). The model will be trained by any suitable algorithm, we may choose a neural representation of the fuzzy system -ANFIS [START_REF] Jang | ANFIS: Adaptive-network-based fuzzy inference system[END_REF] and the appropriate training algorithm. Fig. 3 presents linear local prototypes obtained by linear mean-square approximation of the curve from fig. 2 with  =0.005 -in this case two "sticks" were enough. Smaller  will impose bigger number of LLPs and so bigger number of rules. In this paper we concentrate on the case when the friction is velocity-dependent and it was possible to eliminate the position as the inessential input (although the procedure of final rule selection and model construction if both inputs are important is possible as presented in [START_REF] Kabziński | One-Dimensional Linear Local Prototypes for Effective Selection of Neuro-Fuzzy Sugeno Model Initial Structure[END_REF]). So starting from this point we model friction by m velocity-dependant rules  ,

,..., 2 , 1 , , 0 , 1 m j p v p F THEN IS v IF j j j j     (16)
with generalised bell-shaped functions

       m j j T m j j j mdl v v v v v F v v F 1 1 ) ( : ) ( ), ( ) ( ) ( ) ( ) (       , (17) where . ) ( ) ( ) ( ) ( ) ( 1 ) ( , ] , , , , [ 1 1 , 0 , 1 1 , 0 1 , 1                   v v v v v v v v p p p p m m m m T          (18)
Several adaptive control techniques may be proposed for motion control of systems with friction, depending on details of the model description and the chosen method to prove the stability. In this contribution we describe adaptive backstepping position tracking. We also assume that the thrust force is proportional to the control variable (the motor current)

) ( ) ( t i t F e    , (19)
and coefficient  is not known exactly. As the external disturbance ext F may be compensated exactly the same way as the friction, we assume that 0  . The velocity will be 'virtual control' for position tracking. If we choose the desired velocity

d v according to 1 1 e k x v d d     , (20)
where

1 1 0  k
is a design parameter, we will be able to describe the tracking error dynamics as

2 1 1 2 1 1 1 e e k e e k x x e d d             , ( 21
) v e k x v v e d d       1 1 2  (22)   m m F i v m v m v m e m o friction d o o d o o           , 1 2     , ( 23
)   2 1 1 1 e e k k x v d d          , (24)
The control variable i has to compensate function

friction d o F v m D  1     , (25)
and to assure fast tracking. We will use a model D ˆ for D, incorporating the fuzzy friction model (17,18). The general structure of D ˆ will be given by ξ

A T D  , ( 26
)
where T A ˆ is a vector of adaptive parameters and ξ is known. We have several possibilities to choose the number of adaptive parameters, for example:

1 adaptive parameter:

a T  A ,          mdl N d oN F v m  1  ξ (27)
  k m o T ,  A          mdl N d F v  1  ξ , ( 28
)
2m+1 adaptive parameters:

  T o T m  ,  A            N d v 1  ξ , (29)
In ( 27

        A  . ( 31
)
To investigate the tracking stability we propose Lyapunov function

  A Γ A 2 1 1 2 2 2 1      T o e m e V . ( 32
)
with positive definite symmetric matrix Γ . Taking any of the adaptation laws

   2 ˆe A  , or ΓA A      2 ˆe  , or ΓA A 2 2 2 1 2 ˆe e e        . (33a,b,c)
with small positive  , we are able to prove that the system derivative of (32) is negative outside a certain, bounded set, and so

1 ,e e are uniformly ultimately bounded. For example with adaptation performed according to (33a) we get

2 2 2 2 2 1 1 2 2 2 2 2 1 1 2 1 2 1                      e k e k e e k e k V  (34)
and is negative outside

2 max 2 2 2 1 1  k e e   ,         2 1 , min 2 1 k k k (35)
Several experiments with various adaptive controllers were conducted with the linear motor investigated in the previous sections, with fuzzy model presented in fig. 4 and adaptive laws (33a,b,c). All controllers perform correctly. More accurate fuzzy model results in smaller control signals. In fig. 567we illustrate the performance of adaptive backstepping controller (28) with two adaptive parameters and adaptation law (33c).

the system was to track sinusoidal position trajectory

Conclusions

The proposed approach, connecting fuzzy modeling and adaptive control, may be used for any motion control problem with friction or other disturbances. The fuzzy model is robust against measurement noise and isolated outliers in the modeling data.

On-line gain adaptation by the adaptive laws allows to use a static friction model and to obtain good control performance. The proposed procedure based on linear local prototypes allows to build a fuzzy model as simple as possible. The same approach may be used if friction is function of position and speed. Some other concepts of fuzzy models -for example TSK models with nonlinear consequences were also applied by the author with promising results.

  v. Equation (7) allows also to estimate the influence of dt dF friction and

Fig. 1 . 2 .

 12 Fig. 1. Friction modelling data. Fig. 2. SIFM action curves for the data from fig. 1: + position, ٠speed.

Fig. 3 .

 3 Fig. 3. SIFM action curve for speed (solid) and piecewise linear approximation (dotted) by two LLP.

Fig. 4 .

 4 Fig. 4. Final fuzzy friction model (solid) and the data (++).

  . Let us denote the desired smooth position trajectory by d x , actual position by x and the tracking error by

  ) a ˆ is responsible for general model correction, in (28) o m ˆ is supposed to adapt the changing inertia and

  ' parameter separately. Without loss of generality we may assume existence of "the best" adaptive parameters

Fig. 5 .Fig. 6 .Fig. 7 .Fig. 8 .

 5678 Fig. 5. The tracking error history. Fig. 6. The tracking error history -an enlarged part of fig.5.


Adnene Arbi
email: adnen.arbi@enseignant.edunet.tn

Chaouki Aouiti
email: chaouki.aouiti@fsb.rnu.tn

Abderrahmane Touati
email: abder.touati@fsb.rnu.tn.

Uniform asymptotic stability and global asymptotic stability for time-delay Hopfield neural networks

In this paper, we consider the uniform asymptotic stability and global asymptotic stability of the equilibrium point for time-delays Hopfield neural networks. Some new criteria of the system are derived by using the Lyapunov functional method and the linear matrix inequality approach for estimating the upper bound of the derivative of Lyapunov functional. Finally, we illustrate a numerical example showing the effectiveness of our theoretical results.

Introduction

Hopfield ([1-2]) has proposed Hopfield neural networks (HNN) which have found applications in a broad range of discipline where the targeted problems can reduce to optimization problems. It has been extensively studied and developed in recent years, and it has attracted much attention in the literature on Hopfield neural networks with time delays, . They are now recognized as candidates for information processing systems and have been successfully applied to associative memory, pattern recognition, automatic control, model identification, optimization problems, etc. (we refer to reader [START_REF] Singh | On global robust stability of interval Hopfield neural networks with delays[END_REF][START_REF] Huang | On global asymptotic stability of recurrent neural networks with time-varying delays[END_REF][START_REF] Zhang | Delay-dependent global stability results for delayed Hopfield neural networks[END_REF][START_REF] Liu | Almost periodic solutions for Hopfield neural networks with continuously distributed delays[END_REF][START_REF] Zhou | Synchronization in complex delayed dynamical networks with impulsive effects[END_REF][START_REF] Zhang | Stability of impulsive neural networks with time delays[END_REF][START_REF] Chen | Global stability analysis of impulsive Cohen-Grossberg neural networks with delay[END_REF][START_REF] Xiang | Existence and global exponential stability of periodic solution for delayed high-order Hopfield-type neural networks[END_REF]). Therefore, the study of stability of HNN has caught many researchers' attention. HNN with time delays has been extensively investigated over the years, and various sufficient conditions for the stability of the equilibrium point of such neural networks have been presented via different approaches. In [START_REF] Zhang | Delay-dependent global stability results for delayed Hopfield neural networks[END_REF], [START_REF] Li | Stability properties for Hopfield Neural Networks with delays and impulsive perturbations[END_REF], some sufficient conditions of stability by utilizing the Lyapunov functional method, and linear matrix inequality approach for delayed continuous HNN are derived. In [START_REF] Zong | New Delay-dependent Global Asymptotic Stability Condition for Hopfield Neural Networks with Time-varying Delays[END_REF], G.Zong and J.Liu established a novel delay-dependent condition to guarantee the existence of HNN and its global asymptotic stability by resorting to the integral inequality and constructing a Lyapunov-Krasovskii functional. In [START_REF] Long | Delay-dependent stability analysis for impulsive neural networks with time varying delays[END_REF], S.Long and D.Xu got the sufficient conditions for global exponential stability and global asymptotic stability by using Lyapunov-Krasovskii-type functionally of negative definite matrix and Cauchy criterion. This paper is organized as follows: In section 2, a model of time-delay Hopfield neural network is described. In addition, we present some basic definitions and lemmas. New stability criteria for Hopfield neural network are derived in section 3. An example is given in section 4, to illustrate the advantage of the results obtained. Finally, some conclusions are drawn in section 5.

Preliminaries

Let R denote the set of real numbers, Z + denote the positive integers and R n denote the n-dimensional real space equipped with the Euclidean norm . . The identity matrix, with appropriate dimensions, is denoted by Id and diag(...) denotes the block diagonal matrix. Consider the following delayed HNN model with impulses

           ẋi (t) = -c i x i (t) + n j=1 a ij f j (x j (t)) + n j=1 b ij g j (x j (t -τ (t))) + I i if t = t k x i \ t=t k = d (i) k (x i (t - k ) -xi ) i = 1, ..., n, n, k ∈ Z + , (1)
where n ≥ 2 corresponds to the number of units in a neural network; the impulsive times t k satisfy:

0 ≤ t 0 < t 1 < ... < t k < ..., lim k-→+∞ t k = +∞;
x i corresponds to the state of the unit i at time t; c i is positive constant; f j , g j , denote respectively, the measures of response or activation to their incoming potentials of the unit j at time t and t -τ (t); constant a ij denotes the synaptic connection weight of the unit j on the unit i at time t; constant b ij denotes the synaptic connection weight of the unit j on the unit i at time t -τ (t); I i is the input of the unit i; τ (t) is the transmission delay such as 0 < τ (t) ≤ τ and τ (t) ≤ ρ < 1; t ≥ t 0 ; τ, ρ, d

k are constants; x is the equilibrium point of the first equation in the system [START_REF] Hopfield | Neural networks and physical systems with emergent collective computational abilities[END_REF]. The evolution of the neuron state i at time t k is described by the equation:

x i \ t=t k = x i (t k ) -x i (t - k ) = d (i) k (x i (t - k ) -xi ), i = 1, 2
, ..., n, k = 1, 2, .... The initial conditions associated with system (1) are of the form:

x(s) = φ(s), s ∈ [t 0 -τ, t 0 ] (2)
where

x(s) = (x 1 (s), x 2 (s), ..., x n (s)) T , φ(s) = (φ 1 (s), φ 2 (s), ..., φ n (s)) T ∈ P C([-τ, 0], R n ) P C([-τ, 0], R n ) = {ψ : [-τ, 0] -→ R n ,
) = ψ(t k )}. For ψ ∈ P C([-τ, 0], R n ), the norm of ψ is defined by ψ τ = sup -τ ≤θ≤0 ψ(θ) .
For any t 0 ≥ 0, let:

P C δ (t 0 ) = {ψ ∈ P C([-τ, 0], R n ) : ψ τ < δ}.
In this paper, we assume that some conditions are satisfied so that the equilibrium point of system (1) does exist, see ( [START_REF] Zhang | Delay-dependent global stability results for delayed Hopfield neural networks[END_REF], [START_REF] Zhang | Delay-dependent global stability condition for delayed Hopfield neural networks[END_REF]). Impulsive operator is viewed as perturbation of the equilibrium point x of such system without impulsive effects. Since x is an equilibrium point of system (1), one can derive from system (1) that the transformation y i = x i -xi , i = 1, 2, ..n transforms such system into the following system:

   ẏi (t) = -c i y i (t) + n j=1 a ij F j (y j (t)) + n j=1 b ij G j (y j (t -τ (t))) si t = t k y i (t k ) = (1 + d (i) k )y i (t - k ) i = 1, ..., n, n, k ∈ Z + , (3)
where

F j (y j (t)) = f j (x j + y j (t)) -f j (x j ) and G j (y j (t -τ (t))) = g j (x j + y j (t -τ (t))) -g j (x j ).
To prove the stability of x of system (1), it is sufficient to prove the stability of the zero solution of system (3).

In this paper, we assume that there exist constants L i , M i ≥ 0 such as

| F i (y) |≤ L i | y |, | G i (y) |≤ M i | y |, i ∈ Λ = {1, 2, ...n},
and we set

L max = max i∈Λ L i , M max = max i∈Λ M i , c max = max i∈Λ c i , c min = min i∈Λ c i , D k = diag(1 + d (1)
k , 1 + d (2) k , ..., 1 + d (n) k ).
So, the system (3) can be written as follows:

       ẏ(t) = -Cy(t) + A.F (y(t)) + B.G(y(t -τ (t))) si t = t k y(t k ) = D k y(t - k ) i = 1, ..., n, n, k ∈ Z + , y(t 0 + θ) = ϕ(θ), (4)
where

y(t) = (y 1 (t), ..., y n (t)) T ; y(t -τ (t)) = (y 1 (t -τ (t)), ..., y n (t -τ (t))) T ; C = diag[c 1 , ..., c n ]; A = (a ij ) n×n ; B = (b ij ) n×n ; F (y) = (F 1 (y 1 ), F 2 (y 2 ), ..., F n (y n )) T ; G(y) = (G 1 (y 1 ), G 2 (y 2 ), ..., G n (y n )) T .
Some definitions and lemma of stability for system (1) at its equilibrium point are introduced as follows:

Definitions 21 Assume y(t) = y(t 0 , ϕ)(t) be the solution of (3) through (t 0 , ϕ), then the zero solution of ( 3) is said to be [START_REF] Fu | Introduction of Impulsive Differential Systems[END_REF] P1 stable, if for any > 0 and t 0 ≥ 0, there exists some δ( , t 0 ) > 0 such as

ϕ ∈ P C δ (t 0 ) implies y(t 0 , ϕ)(t) < , t ≥ t 0 . P2 uniformly stable, if the δ in (P1) is independent of t 0 .
P3 uniformly attractive, if there exists some δ > 0 such as for any > 0, there exists some T = T ( , δ) > 0 such as t 0 ≥ 0 and ϕ ∈ P C δ (t 0 ) implies y(t 0 , ϕ)(t) < , t ≥ t 0 + T . P4 uniformly asymptotically stable, if (P2) and (P3) are held. P5 globally asymptotically stable, if (P1) holds and for any given initial value y 0 = ϕ, y(t 0 , ϕ)(t) -→ 0 as t -→ +∞. Now, we need the following basic lemmas used in our work.

Lemma 21 [START_REF] Liao | LMI approach for global periodicity of neural networks with time-varying delays[END_REF] For any a, b ∈ R n , the inequality

±2a T b ≤ a T Xa + b T X -1 b
holds, where X is any n × n matrix with X > 0.

Robust stability criteria

In this section, we shall establish some theorems which provide sufficient conditions for uniform asymptotic stability and global asymptotic stability of system [START_REF] Hopfield | Neural networks and physical systems with emergent collective computational abilities[END_REF].

Theorem 31 The system (1) is uniformly stable if there is * ∈ [0, 1], σ > 0 and positive n × n definite matrix Q such as:

(i) σ • max i,j ci•qij min i,j |bij | ≤ τ 2 +4-τ √ τ 2 +4 τ 2 +4+τ √ τ 2 +4 (ii) max 1≤i≤n { 1 ci n j=1 a ij } + max 1≤j≤n {L 2 j n i=1 aij ci } + max 1≤j≤n {M 2 j n i=1 |bij | ci } + * cmin + 1 σ λ max (C -1 BQ -1 B T C -1 ) < 2 (iii) t 0 <t k <t max{cmaxξ k ,1} 1+ * (t-t0) 2
< ∞, where ξ k is the largest eigenvalue of D k C -1 D k , and k ∈ Z + . In addition, if we have:

(iv) t 0 <t k <t max{cmaxξ k ,1} 1+ * (t-t0) 2
-→ 0, if t -→ +∞, then the system (1) is uniformly asymptotically stable and globally asymptotically stable.

Proof. First, we prove the equilibrium point of system (1) is uniformly stable. We consider this Lyapunov function:

V (y)(t) = [1+ * (t-t 0 ) 2 ] n i=1 1 c i y 2 i (t)+ n i=1 n j=1 | b ij | c i t t-τ (t) (1+(s-t 0 ) 2 )G 2 j (y j (s))ds.
(5) For condition (iii), there is a constant M * > 0, such as:

t0<t k <t max{c max ξ k , 1} 1 + * (t -t 0 ) 2 < M * , t ≥ t 0 .
For any t 0 ≥ 0, let y(t 0 , ϕ)(t) be a solution of system (1). So, ∀ > 0, we choose δ of the following manner:

δ = 1 c max [ 1 cmin + max 1≤i≤n { n j=1 M 2 j |bij | ci }(τ + τ 3 3 )]M *
.

By simple calculation:

∂V (y)(t) ∂t = [1 + * (t -t 0 ) 2 ]{-2 n i=1 y 2 i (t) + 2 n i=1 n j=1 a ij c i y i (t)F j (y j (t)) + 2 n i=1 n j=1 b ij c i y i (t)G j (y j (t -τ (t)))} + 2 * (t -t 0 ) n i=1 1 c i y 2 i (t) + n i=1 n j=1 | b ij | c i (1+(t-t 0 ) 2 )G 2 j (y j (t))- n i=1 n j=1 | b ij | c i (1+(t-τ (t)-t 0 ) 2 )G 2 j (y j (t-τ (t))).
Therefore,

∂V (y)(t) ∂t = -2(1 + * (t -t 0 ) 2 ) n i=1 y 2 i (t) + (1 + * (t -t 0 ) 2 ){ n i=1 n j=1 a ij c i y 2 i (t) + n i=1 n j=1 a ij c i F 2 j (y j (t)) + 1 σ λ max (C -1 BQ -1 B T C -1 ) n i=1 y 2 i (t)} + 2 * (t -t 0 ) n i=1 1 c i y 2 i (t) + σ(1 + * (t -t 0 ) 2 ) n i=1 n j=1 q ij G 2 j (y j (t -τ (t))) -(1+(t-τ (t)-t 0 ) 2 ) n i=1 n j=1 | b ij | c i G 2 j (y j (t-τ (t)))+(1+(t-t 0 ) 2 ) n i=1 n j=1 | b ij | c i G 2 j (y j (t)).
It follows that,

∂V (y)(t) ∂t ≤ (1 + (t -t 0 ) 2 ){-2 n i=1 y 2 i (t) + n i=1 n j=1 a ij c i y 2 i (t) + n i=1 n j=1 a ij c i F 2 j (y j (t))+ n i=1 n j=1 | b ij | c i G 2 j (y j (t))+ 1 σ λ max (C -1 BQ -1 B T C -1 ) n i=1 y 2 i (t) + 2 * (t -t 0 ) 1 + (t -t 0 ) 2 n i=1 1 c i y 2 i (t)} + n i=1 n j=1 [(1 + (t -t 0 ) 2 )σ • q ij -(1 + (t -τ (t) -t 0 ) 2 ) | b ij | c i ]G 2 j (y j (t -τ (t))). ( 6
)
We have:

2 n i=1 n j=1 1 c i b ij y i (t)G j (y j (t -τ (t))) = 2y T (t)C -1 BG(y(t -τ (t))) = 2G T (y(t -τ (t)))B T C -1 y(t) = 2[G(y(t -τ (t))) √ σ] T (B T C -1 y(t) 1 √ σ ) ≤ σG T (y(t -τ (t)))QG(y(t -τ (t))) + 1 σ y T (t)C -1 BQ -1 B T C -1 y(t) ≤ σ n i=1 n j=1 q ij G 2 j (y j (t -τ (t))) + 1 σ λ max (C -1 BQ -1 B T C -1 ) n i=1 y 2 i (t). ( 7
)
So, from ( 6) and ( 7):

∂V (y)(t) ∂t ≤ (1 + (t -t 0 ) 2 )[-2 + max 1≤i≤n { 1 c i n j=1 a ij } + max 1≤j≤n {L 2 j n i=1 a ij c i }+ max 1≤j≤n {M 2 j n i=1 | b ij | c i }+ 1 σ λ max (C -1 BQ -1 B T C -1 )+ * c min ] y(t) 2 + n i=1 n j=1 [(1 + (t -t 0 ) 2 )σ • q ij -(1 + (t -τ (t) -t 0 ) 2 ) | b ij | c i ]G 2 j (y j (t -τ (t))).
Then we obtain,

∂V (y)(t) ∂t < 0, (8)
if

(1 + (t -t 0 ) 2 )σ • q ij ≤ (1 + (t -τ (t) -t 0 ) 2 ) | b ij | c i .
Therefore, it is sufficient that:

σ•qij ci |bij | ≤ (1+(t-τ (t)-t0) 2 ) (1+(t-t0) 2 )
.

Let u(t) = 1+(t-τ (t)) 2 1+t 2
, next we show for t ≥ 0

u(t) ≥ τ 2 + 4 -τ √ τ 2 + 4 τ 2 + 4 + τ √ τ 2 + 4 . (9)
First, for t ∈ [τ, +∞[, we have:

u(t) ≥ 1+(t-τ ) 2 1+t 2 = v(t), it is easy to compute that for t ≥ 0 v min = v( τ + √ τ 2 + 4 2 ) = τ 2 + 4 -τ √ τ 2 + 4 τ 2 + 4 + τ √ τ 2 + 4 , also we obtain v(τ ) > v min , that is 1 1 + τ 2 ≥ τ 2 + 4 -τ √ τ 2 + 4 τ 2 + 4 + τ √ τ 2 + 4 . ( 10
)
Second, for t ∈ [0, τ [, we have:

u(t) = 1 + (t -τ (t)) 2 1 + t 2 ≥ 1 1 + τ 2 . ( 11
)
In view of ( 10) and ( 11), we obtain that (9) also holds for t ∈ [0, τ [. Then we obtain that we have proved ( 9) holds for all t ∈ [0, +∞[. We have ∀k ≥ 1:

V (y)(t k ) = [1+(t k -t 0 ) 2 ] n i=1 1 c i y 2 i (t k )+ n i=1 n j=1 | b ij | c i t k t k -τ (t k ) (1+(s-t 0 ) 2 )G 2 j (y j (s))ds = [1+(t k -t 0 ) 2 ]y T (t - k )D k C -1 D k y(t - k )+ n i=1 n j=1 | b ij | c i t - k t - k -τ (t - k ) (1+(s-t 0 ) 2 )G 2 j (y j (s))ds ≤ [1+(t k -t 0 ) 2 ]ξ k y T (t - k )y(t - k )+ n i=1 n j=1 | b ij | c i t - k t - k -τ (t - k ) (1+(s-t 0 ) 2 )G 2 j (y j (s))ds ≤ [1+(t k -t 0 ) 2 ] ξ k λ min (C -1 ) y T (t - k )C -1 y(t - k )+ n i=1 n j=1 | b ij | c i t - k t - k -τ (t - k ) (1+(s-t 0 ) 2 )G 2 j (y j (s))ds ≤ max{ ξ k λ min (C -1 ) , 1}V (t - k ). Therefore, V (y)
(t k ) ≤ max{ξ k c max , 1}V (t - k ). (12)
Then, we have:

1 c max (1 + * (t -t 0 ) 2 ) y(t) 2 ≤ V (t) ≤ V (t 0 ) × t0<t k ≤t max{ξ k c max , 1}. (13)
From ( 5), we have:

1 c max (1 + * (t -t 0 ) 2 ) y(t) 2 ≤ V (t) ≤ 1 c min (1 + * (t -t 0 ) 2 ) y(t) 2 + max 1≤i≤n { n j=1 M 2 j | b ij | c i }(τ (t) + (t -t 0 ) 3 -(t -t 0 -τ (t)) 3 3 ) y(t) 2 .
Therefore,

V (t) ≤ [ 1 c min (1 + * (t -t 0 ) 2 ) + max 1≤i≤n { n j=1 M 2 j | b ij | c i }(τ (t) + (t -t 0 ) 3 -(t -t 0 -τ (t)) 3 3 )] y(t) 2 . ( 14
)
For t = t 0 , we have:

V (t 0 ) ≤ [ 1 c min + max 1≤i≤n { n j=1 M 2 j | b ij | c i }(τ + τ 3 3 )] ϕ 2 . (15)
Then,

1 c max (1+ * (t-t 0 ) 2 ) y(t) 2 ≤ [ 1 c min + max 1≤i≤n { n j=1 M 2 j | b ij | c i }(τ + τ 3 3 )] ϕ 2 × t0<t k <t {ξ k c max , 1}.
This implies that:

y(t) 2 ≤ [ 1 c min + max 1≤i≤n { n j=1 M 2 j | b ij | c i }(τ + τ 3 3 )]c max ϕ 2 × t0<t k <t {ξ k c max , 1} (1 + * (t -t 0 ) 2 ) ≤ 2 .
Hence, the zero solution of system ( 1) is uniformly stable.

In view of condition (iv), it is obvious that : lim sup t-→+∞ y(t) 2 = 0, so the equilibrium point of system ( 1) is also uniformly asymptotically stable and globally asymptotically stable.

Which completes the proof.

If

t0<t k ≤t max{c max • ξ k , 1}
< ∞, then we can get the following criterion for stability with * = 0.

Corollary 32 Assume that there is a constant σ > 0 and n × n positive definite matrix Q such as:

(i) σ • max i,j ci•qij min i,j |bij | ≤ τ 2 +4-τ √ τ 2 +4 τ 2 +4+τ √ τ 2 +4 (ii) max 1≤i≤n { 1 ci n j=1 a ij } + max 1≤j≤n {L 2 j n i=1 aij ci } + max 1≤j≤n {M 2 j n i=1 |bij | ci } + 1 σ λ max (C -1 BQ -1 B T C -1 ) < 2.
Then, the equilibrium point of system (1) is uniformly asymptotically stable and globally asymptotically stable.

If σ = τ 2 +4-τ √ τ 2 +4 τ 2 +4+τ √ τ 2 +4 and Q = 1 60
Id in Corollary 32, then we can get the following criterion for stability.

Corollary 33 Assume that the following conditions are satisfied:

(i) max i ci 60 min i,j |bij | ≤ 1 (ii) max 1≤i≤n { 1 ci n j=1 a ij } + max 1≤j≤n {L 2 j n i=1 aij ci } + max 1≤j≤n {M 2 j n i=1 |bij | ci } + 60 σ λ max (C -1 BB T C -1 ) < 2.
Then, the equilibrium point of system (1) is uniformly asymptotically stable and globally asymptotically stable.

Application

In this section, we present a numerical example to illustrate that our conditions are more feasible than those given in some earlier references [START_REF] Zhang | Delay-dependent global stability results for delayed Hopfield neural networks[END_REF].

Example 41 Consider the two-neuron delayed neural network with impulses as follows:

                   ẋ1 (t) = -2.5x 1 (t) -0.5f 1 (x 1 (t)) + 0.1f 2 (x 2 (t)) -0.1g 1 (x 1 (t -τ )) + 0.2g 2 (x 2 (t -τ )) -1 ẋ2 (t) = -2x 2 (t) + 0.2f 1 (x 1 (t)) -0.1f 2 (x 2 (t)) +0.2g 1 (x 1 (t -τ )) + 0.1g 2 (x 2 (t -τ )) + 4 x i | t=t k = x i (t k ) -x i (t - k ) = d (i) k (x i (t - k ) -xi ) k ∈ Z + , (16)
where τ = 0.87, the activation functions are the following:

f 1 (x) = f 2 (x) = g 1 (x) = g 2 (x) = 0.5(|x + 1| -|x -1|) and d (1) k = 1 + 1 5k 2 -1, d (2)
k = 1 + 1 6k 2 -1, t k = k, k ∈ Z + . So, the matrix A, B and C are: A = -0.5 0.1 0.2 0.1 , B = -0.1 0.2 0.2 0.1 , C = 2.5 0 0 2 .
By Matlab, we note

C -1 BB T C -1 = 0.008 0 0 0.0125 so, λ max (C -1 BB T C -1 ) = 0.0125.
Considering the activation functions f 1 , f 2 , g 1 and g 2 , we can choose

L i = 1, M i = 1, i = 1, 2.
On the other hand, by Mathematica software, we note

∞ s=1 max i=1,2 (1 + d (i) s ) 2 = ∞ s=1 max i=1,2 (1 + 1 5s 2 ) 2 < 1.4
We show the equilibrium point of ( 16) is uniformly stable and globally asymptotically stable. From Corollary 33: σ = τ

| b ij | c i }+ 60 σ λ max (C -1 BB T C -1 ) < 2.
Hence, by Corollary 33, the equilibrium point (-0.2258, 1.9548) T of system ( 16) is uniformly asymptotically stable and globally asymptotically stable.

Remark 41 In [START_REF] Zhang | Delay-dependent global stability results for delayed Hopfield neural networks[END_REF], the authors proved that system ( 16) is globally asymptotically stable with upperbound of delays τ = 0.17. For this example, we additionally get that the equilibrium point of system ( 16) is uniformly asymptotically stable and globally asymptotically stable with upperbound of delays τ = 0.87 > 0.17 . However, the criteria given in [START_REF] Zhang | Delay-dependent global stability results for delayed Hopfield neural networks[END_REF] are invalid for (τ ≥ 0.87). Therefore, our results are less conservative and more efficient those that given in [START_REF] Zhang | Delay-dependent global stability results for delayed Hopfield neural networks[END_REF].

Conclusion

In this paper, a class of HNN with delays is considered. We obtain some new criteria ensuring the global and uniform asymptotic stability of the equilibrium point for such system by using the Lyapunov method and linear matrix inequality. Our results show the effects of delay on the stability of HNN. The results here are compared to earlier results. Our criteria are more simpler to verify. An example is given to illustrate the efficiency of the results.


Ruth Janning
email: janning@ismll.uni-hildesheim.de

Tomáš Horváth
email: horvath@ismll.uni-hildesheim.de

Andre Busche
email: busche@ismll.uni-hildesheim.de

Lars Schmidt-Thieme
email: schmidt-thieme@ismll.uni-hildesheim.de

GamRec: a Clustering Method Using Geometrical Background Knowledge for GPR Data Preprocessing

Keywords: Ground penetrating radar (GPR), Object detection, Clustering, Sweep line algorithm, Preprocessing

GPR is a nondestructive method to scan the subsurface. On the resulting radargrams, originally interpreted manually in a time consuming process, one can see hyperbolas corresponding to buried objects. For accelerating the interpretation a machine shall be enabled to recognize hyperbolas on radargrams autonomously. One possibility is the combination of clustering with an expectation maximization algorithm. However, there is no suitable clustering algorithm for hyperbola recognition. Hence, we propose a clustering method specialized for this problem. Our approach is a directed shape based clustering combined with a sweep line algorithm. In contrast to other approaches our algorithm finds hyperbola shaped clusters and is (1) able to recognize intersecting hyperbolas, (2) noise robust and (3) does not require to know the number of clusters in the beginning but it finds this number. This is an important step towards the goal to fully automatize the buried object detection.

Introduction

Ground penetrating radar (GPR) is a method to scan the shallow subsurface without destroying the road surface. It can be used to detect buried objects like pipes, cables, ducts and sewers. To get a 2D GPR radargram image (fig. 1 (b), 2 (a)) of a cross-sectional area of the subsurface, electromagnetic waves are transmitted into the ground and the reflected signals are caught by an antenna. The vector of reflections (intensities on the radargram image) measured at one certain point for different answer times is called an A-Scan. The radargram image, the B-Scan, is a sequence of consecutive A-Scans. If the antenna is moved in a line perpendicular to a pipe, then the signals reflected from the pipe have the shape of a hyperbola branch in the radargram image, because (as fig. 1 (a) shows) the more the antenna nears the pipe the shorter is the time in which the corresponding signal is reflected. Originally, GPR radargram images were interpreted manually by human experts, but the pipe localization takes much time in this way (cp. [START_REF] Simi | Hough transform based automatic pipe detection for array GPR: algorithm development and on-site tests[END_REF]). Hence, an automatization of the radargram interpretation which can accelerate the whole process is required. However, in the beginning it is not known how much hyperbolas there are and where the hyperbola branches are located. Additionally, there might be noise in the radargram image and the different hyperbola branches might intersect (see fig. 2,3).

General Problem

For A searched upper hyperbola branch can be defined as f a,b,j0,i0 : R → R, i →

n ∈ N let [n] := {x ∈ N | 0 ≤ x ≤ (n -1)}. Given is a dense 2D image B ⊆ [c] × [r] × [ι min , ι max ], c, r ∈ N, ι min , ι max ∈ R , (1) a B
(1 + (j-j0) 2 b 2
) • a 2 + i 0 where (j 0 , i 0 ) is the center of the hyperbola, a is the distance between the center and the apexes of the branches and b influences the curvature of the hyperbola. Let the distance between a point and a hyperbola branch be D((i, j), f ) := d((j, i), (j, f (j))) for any distance measure between points, e.g. the Euclidean distance d((j, i), (j , i )) := (j -j ) 2 + (i -i ) 2 . The hyperbola recognition problem can then be formalized as follows:

Definition 1 (Hyperbola recognition problem). Given B T and D, find a set of hyperbola branches

F := {f 1 , ..., f K } and a decomposition B T = ( f k ∈F h f k ) ∪ B noise , ∀f k ∈ F : B noise ∩ h f k = ∅, of B T into clusters h f k and a noise cluster B noise with |B noise | ≈ ν |F | • |B T | for some ν ∈ [0, 1
) such that the following sum of distances of all points ∈ B T \ B noise to their assigned hyperbola branch is minimal:

err := K k=1 (j,i)∈h f k D((j, i), f k ) (2)
3 Related Work

In the last few years different approaches were applied to the hyperbola recognition problem on GPR radargrams. The Hough Transform (HT) (e.g. used in [START_REF] Simi | Hough transform based automatic pipe detection for array GPR: algorithm development and on-site tests[END_REF], [START_REF] Windsor | The Estimation of Buried Pipe Diameters by Generalized Hough Transform of Radar Data[END_REF]) is an often used method. However, the HT is a brute force method and is hence computationally intensive (see e.g. [START_REF] Chen | Probabilistic Robust Hyperbola Mixture Model for Interpreting Ground Penetrating Radar Data[END_REF]). A more efficient approach is presented in [START_REF] Chen | Probabilistic Robust Hyperbola Mixture Model for Interpreting Ground Penetrating Radar Data[END_REF]. This method uses clustering combined with a classification expectation maximization algorithm. In an iterative process the algorithm fits hyperbola branches to clusters and then rearranges the clusters by assigning each point to one cluster according to the maximum posterior probability of the given point to be in this cluster. The hyperbola fitting and a rearrangement of the clusters is done alternately until convergence. However, for getting the initial clusters a K-means algorithm is used. The K-means clustering has some important drawbacks in relation to the hyperbola recognition problem:

1. It is not resistant against noise (even if a noise cluster is used which contains the q% most distant points of the other clusters, see fig. 2 (b),(d)). 2. It can recognize only convex clusters (see fig. 2 (c)). 3. The number K of clusters to detect has to be known beforehand.

In [START_REF] Chen | Probabilistic Robust Hyperbola Mixture Model for Interpreting Ground Penetrating Radar Data[END_REF] a Bayesian information criterion (BIC) is proposed to estimate K, but also for this procedure a range for the possible number of hyperbolas has to be estimated beforehand and for each of the values in this range the whole approach has to be applied. Afterwards the model with the highest overall likelihood is selected as the final result. A clustering method which avoids the above mentioned drawbacks of Kmeans is the shape based clustering method DBSCAN [START_REF] Ester | A density-based algorithm for discovering clusters in large spatial databases with noise[END_REF]. This method is able to treat noise, to find clusters of arbitrary shapes and the number of occurring clusters has not to be known beforehand. Nevertheless, the original DBSCAN 1. is not able to recognize intersecting hyperbola branches, which are caused by nearby buried objects and often occur in GPR radargrams (it finds one cluster for two or more intersecting hyperbolas (see fig. 3 (b))), and 2. does not differentiate between clusters with hyperbola shape and clusters with other shapes.

Our Contribution

The drawbacks of DBSCAN in relation to the hyperbola recognition problem can be treated by integrating geometrical background knowledge. For this purpose we propose GamRec (Geometric approach for multi-hyperbola Recognition), a kind of a directed DBSCAN combined with a sweep line algorithm. Sweep line algorithms (see e.g. [START_REF] Boissonnat | Algorithmic Geometry[END_REF]) originate from algorithmic geometry. The idea is to imagine that a line moves (sweeps) over the considered plane from one side to the other one and pauses if it touches a special point to apply some action to this point. Such points are visited just once, i.e. every point behind the sweep line is never visited again. The characteristic of sweep line algorithms is that they reduce the computational complexity by translating an n dimensional static problem into an n -1 dimensional dynamic problem. A clustering algorithm for spatial data which uses a sweep line algorithm is presented in [START_REF] Zalik | A sweep-line algorithm for spatial clustering[END_REF], but it is not able to separate overlapping clusters, i.e. it can also not recognize intersecting hyperbola branches. In GPR radargram images apexes are the highest points of hyperbola branches and if an apex is once found one has to search for the left and right hyperbola branch side below it in left-down and right-down direction. Hence, a sweep line can be moved top down over the image and collect hyperbola shaped clusters composed of an apex, a right and a left side. In doing so intersections are treated by considering that a point of an intersection must belong to the right side of one hyperbola branch and to the left side of another one. The usual sorting phase of a sweep line algorithm is not needed, as the considered points are already sorted by their occurrence in the input image. Our approach treats all the above mentioned problems of K-means and DBSCAN in relation to the hyperbola recognition problem, i.e. it (see also fig. The shapes of the resulting clusters are very close to the shapes of the original hyperbola branches (see fig. 3 (c)). That indicates that a hyperbola fitting has to be applied just once after running GamRec instead of applying it several times in an iterative process to rearranged clusters until the clusters are close enough to the hyperbola branches, as done in the approach of [START_REF] Chen | Probabilistic Robust Hyperbola Mixture Model for Interpreting Ground Penetrating Radar Data[END_REF]. To strengthen this assumption in our experiments in section 6 we compare the original initialization (K-means with K given) of the approach in [START_REF] Chen | Probabilistic Robust Hyperbola Mixture Model for Interpreting Ground Penetrating Radar Data[END_REF] with GamRec as an initialization and show that if we use GamRec instead of K-means with K given just one hyperbola fitting step is needed for achieving already very good results. But first we will introduce GamRec in the following section 5.

The GamRec Algorithm

Our approach (see also fig. 4) consists of three consecutive steps: the preparation, the sweep step and the noise removing step.

Preparation

Before applying our algorithm (and any other clustering method) the B-scan B has to be translated into B T . Hence, for every row i of the radargram image an intensity threshold τ i is computed for

B T := {(j, i) | (j, i, ι) ∈ B ∧ ι > τ i }.
Threshold τ i corresponds to the lower limit of the upper p% quantile of the intensities ι of the points (j, i, ι) ∈ B in row i. Parameter p depends on the contrast of the image but can be considered as a constant for the same kind of images (for our experiments we used p = 18).

The Sweep Step

GamRec searches for a set of clusters H

:= {h 1 , h 2 , ..., h K | h k ⊆ B T , k = 1, ..., K} with h k := {(j, i) ∈ B T | ∃(j , i ) ∈ h k = (j, i) : (|j -j | ≤ e) ∧ (|i -i | ≤ e)
} where e is the maximal distance between two neighboring points in a cluster. Such clusters h k shall possess the following properties:

-h k = h k L ∪ h k C ∪ h k R is a union of point sets.
• h k C contains the apex (j a , i a ).

• h k L := {(j, i) ∈ h k | j < j a , i ≥ i a } contains all points of the left side. • h k R := {(j, i) ∈ h k | j > j a , i ≥ i a }
contains all points of the right side.

-For every point (j, i)

∈ h k must hold |( (1 + (j-j0) 2 b 2
) • a 2 + i 0 ) -i| < δ with appropriate a, b, i 0 , j 0 and δ. That means the points in the cluster are positioned around and near to a hyperbola branch.

Sweep: H := {}; for i from 1 to r do estimate all middle points in row i;

for every found middle point poi = (j mid , i) do if poi is not used right and not used left then do

h k := directedNeighborhoodSearch(poi); if |h k L | > e 2 ∧ |h k R | > e 2 then do H := H ∪ {h k };
for all points poi L ∈ h k L do mark poi L as used left; for all points poi R ∈ h k R do mark poi R as used right; Our algorithm (fig. 5,6) applies a sweep line algorithm by moving top down over the thresholded image

T := {(j, i, ι) | (j, i, ι ) ∈ B ∧ ((j, i) ∈ B T ⇒ ι = 1) ∧ ((j, i) / ∈ B T ⇒ ι = 0)}.
With the sweep line all middle points (see def. 2) are visited.

Definition 2 (Middle point). Let s = ((j 1 , i), ..., (j n , i)) with j 1 < j 2 < ... < j n , |j w+1 -j w | ≤ e, w = 1, ..., n be a sequence of points ∈ B T in row i. For every such sequence s with length n in row i a middle point (j mid , i) is the point in s with the index mid = div( n 2 ). For every middle point a directed neighborhood search is started (fig. 6), if the point was not used left or used right before (a point is assigned as used left (used right), if it belongs to the left (right) side of a hyperbola cluster already found). For the points ∈ B T of the left side h k L of the appendant cluster is searched in the neighborhood in left and down direction and for the points ∈ B T of the right side h k R in right and down direction. A neighboring point ∈ B T is inserted into h k L , (1) if the point was not used left and not used right before or (2) if it was yet used right for another cluster and is now located left from the current apex. h k R is treated in the same way but with left and right transposed. In this way GamRec is able to treat intersections of different hyperbola branches. It allows that points may belong to more than one cluster if they are located in the intersection of two hyperbola shaped clusters. GamRec uses one hyper parameter e for the maximal distance between neighboring cluster points (equivalent to in DBSCAN). Parameter e depends on the closeness of directedNeighborhoodSearch((j mid , i)) the points belonging to a hyperbola cluster. For our experiments we used a constant value of e = 8, estimated by a grid search. The second hyper parameter of DBSCAN, the minimal number of neighbors min N , is not needed in GamRec, as noise is identified mainly by considering the hyperbola shape as described in the following section 5.3.

1. h k C := {(j mid , i)}; 2. h (0) k L := {poi L = (j mid -e1, i + e2) ∈ B T | e1 = 1, ..., e, e2 = 0, ..., e, poi L is not used left ∨ poi L is used right}; until h (t) k L = h (t-1) k L do h (t) k L := h (t-1) k L ∪ {poi L = (j, i) ∈ B T | ∃(jv, iv) ∈ h (t-1) k L : j = jv -e3, i = iv + e4, e3 = 0, ..., e, e4 = 0, ..., e, poi L is not used left ∨ poi L is used right}; h k L := h (t) k L ; 3. h (0) k R := {poi R = (j mid + e1, i + e2) ∈ B T | e1 = 1, ..., e, e2 = 0, ..., e, poi R is not used right ∨ poi R is used left}; until h (t) k R = h (t-1) k R do h (t) k R := h (t-1) k R ∪ {poi R = (j, i) ∈ B T | ∃(jv, iv) ∈ h (t-1) k R : j = jv + e3, i = iv + e4, e3 = 0, ..., e, e4 = 0, ..., e, poi R is not used right ∨ poi R is used left}; h k R := h (t) k R ; 4. return h k L ∪ h k C ∪ h k R ;

The Noise Removing Step

Finally, clusters which do not correspond to hyperbola branches are removed from H (see fig. 7). GamRec identifies such a non-hyperbola cluster by considering geometrical properties of hyperbolas. It investigates if the side length of one side of the rectangle around all of its points is less than 2e (to eliminate larger background noise), in this rectangle the points in the upper left corner, in the upper right corner or in the lower center belong to the cluster and the height of one side is less than half of the height of the other side.

In this cases the cluster cannot have the shape of a hyperbola branch.

Complexity

Let n be the number of considered points for a clustering. In our experiments we compare the original initialization (K-means with K given) used in [START_REF] Chen | Probabilistic Robust Hyperbola Mixture Model for Interpreting Ground Penetrating Radar Data[END_REF] with GamRec as an initialization. Hence, we consider the complexity of GamRec as well as of the K-means clustering. The complexity O(n • K • t) of K-means clustering results from the t iterations multiplied by the n investigations for every point of its distances to all K cluster centers. The complexity of GamRec corresponds to the sum of the complexities of the sweep and the noise removing

Noise removing: for all h k ∈ H do if (jmax -jmin) ≤ 2e ∨ (imax -imin) ≤ 2e then do H := H \ h k ; else do if ({(jmin + e3, imin + e4) | e3 = 0, ..., e, e4 = 0, ..., e} ⊂ h k L ) ∨ ({(jmax -e3, imin + e4) | e3 = 0, ..., e, e4 = 0, ..., e} ⊂ h k R ) ∨ ({(j mid ± e3, imax -e4) | e3 = 0, ..., e, e4 = 0, ..., e} ⊂ h k ) then do H := H \ h k ; a else do if ((imax L -imin L ) < 0.5(imax R -imin R )) ∨ ((imax R -imin R ) < 0.5(imax L -imin L )) then do H := H \ h k ;
a jmin, imin are the minimal j and i values of all points ∈ h k , jmax, imax the maximal j and i values (correspondingly imin R , imax R and imin L , imax L for the right and left side) and j mid is the j value of the apex of h k . Fig. 7. GamRec algorithm: noise removing step.

step. As GamRec is applied to points of an image, the sorting of points for the sweep can be skipped and an efficient neighborhood search is enabled. The sweep consists of two main steps. The first step is the search for the m middle points, which takes O(n) time, as every point in every row is visited. The second step includes a directed neighborhood search with a complexity equal to the size of the corresponding cluster for every found middle point times 2e 2 , as for every point 2e ), because then every point has either no neighbors (two middle points of the same row cannot have a distance less or equal than e, see def. 2) and size i = 1 for every i, or a middle point possesses neighbors below it which are middle points and belong to its cluster but these points are then marked as used left or used right. That means that from these points no neighborhood search will be started in the following. Additionally, these points will either be visited never again or only once again if it belongs to an intersection.

Experiments

We simulated a cross-sectional area of 2.50 meter length and 0.60 meter depth. The radargrams (see 2 (a)) of this cross-sectional area were produced by GprMax [START_REF] Giannopoulos | Modelling ground penetrating radar by GprMax[END_REF], an electromagnetic wave simulator for modeling GPR data based on the Finite-Difference Time-Domain numerical method. GprMax is an often used tool in this field, as real data are rare and expensive to achieve. Additionally, with this data we can easily compare the results of GamRec and of the original initialization of the approach in [START_REF] Chen | Probabilistic Robust Hyperbola Mixture Model for Interpreting Ground Penetrating Radar Data[END_REF], as we know exactly the ground truth, i.e. where the pipes are located. This means we can compare the positions found with the real positions of the buried objects. In each of 45 scenarios, 5 perfect conductors (e.g. wires) with a radius of 3 centimeters are buried in a 2.50×0.60 meter area of wet sand at different locations with a distance of 25 centimeters between neighboring pipes. We used a 600 MHz antenna frequency for the measurements. In this way we investigated 225 pipes at different horizontal positions and different depths.

For the experiments we applied GamRec (with e = 8) as well as the original initialization (K-means with K given, K = 5, t = 20) of the approach in [START_REF] Chen | Probabilistic Robust Hyperbola Mixture Model for Interpreting Ground Penetrating Radar Data[END_REF] to the thresholded radargram images (with p = 18, see fig. 3 (a)). To the clusters found a hyperbola fitting algorithm was applied (from the K-means clusters the 10% most distant points were removed into the noise cluster). Subsequently, the positions found (fitted apex for K-means, depth of apex of hyperbola cluster and horizontal position of fitted apex for GamRec) were compared to the upper positions of the corresponding buried objects by computing the Euclidean distance between them (see table 1). Table 1 shows that after one run of GamRec the real positions of the buried objects can be estimated already with an error of just a few millimeters and at most about one centimeter. K-means with K given in contrast delivers position estimations with an error from a few centimeters up to more than one meter. These results indicate that GamRec delivers not just the better initialization but also that there is no need for a rearrangement of the clusters found by GamRec. One hyperbola fitting step after the application of GamRec suffices to reach already very good position estimation results.

Table 1. Average distance errors in meter of the 5 pipes in each of the 45 scenarios ordered from minimal to maximal Euclidean distance. The overall average Euclidean distance error of all 225 pipes is presented in the last column.

Conclusions

We presented GamRec, a clustering method specialized for an improved GPR data preprocessing. Just one run of GamRec and one hyperbola fitting step applied to the hyperbola shaped clusters found is needed to reach already very good pipe position estimation results. Furthermore, GamRec finds the number of occurring hyperbolas by itself, it is able to recognize intersecting hyperbola branches and it is resistant against background noise. This work is an initial study in the field of pipe detection by GPR. Our final goal is a full automatization of the whole pipe detection process, including pipe position estimation and pipe course detection.

Fig. 1 .

 1 Fig. 1. The antenna (rectangle) is moved perpendicular to the buried object (circle). The downwards directed arrows represent transmitted signals, the dotted lines indicate the position of the received signal values in the corresponding A-Scans. On the radargram (right) one can see a hyperbola branch corresponding to the buried object.

-

  Scan composed of c columns (A-scans) with intensities in a range from ι min to ι max . The r rows of B correspond to the signal travel times. In this B-scans upper branches (the origin is the upper left corner of the image) of North-South opening hyperbolas with extreme intensities are searched.B T := {(j, i) | (j, i, ι) ∈ B ∧ ι > τ } is a sparse representation of B containing only points with extreme intensities larger than a threshold τ .

Fig. 2 .

 2 Fig. 2. (a) Radargram image of 5 buried objects. (b) K-means does not recognize the noise points as noise and assigns them to a cluster (black, fitted hyperbola: white). (c) K-means can only find convex clusters (black, fitted hyperbola: white). (d) final result of K-means with K = 5, t = 20 and q = 10 for the noise cluster.

  3 (c), (d))1. ignores background noise, 2. finds hyperbola branch shaped clusters composed of apex, right and left side, 3. has not to know the number of occurring hyperbolas in the beginning, 4. recognizes intersecting hyperbola branches as different hyperbola branches, 5. distinguishes between hyperbola shaped clusters and non-hyperbola clusters.

Fig. 3 .

 3 Fig. 3. (a) Thresholded radargram image with p = 18. (b) One Single DBSCAN cluster (white) found for actually 5 hyperbolas. (c) Single GamRec cluster (white). (d) final resulting clusters (black ) and fitted hyperbola branches (white) of GamRec applied to the thresholded radargram image of the radargram in fig. 2 with e = 8.

Fig. 4 .

 4 Fig. 4. Architecture of GamRec (nbh-search = neighborhood search).

Fig. 5 .

 5 Fig. 5. GamRec algorithm: sweep step.

Fig. 6 .

 6 Fig. 6. GamRec: directed neighborhood search.

2 1 .

 21 points in the neighborhood are investigated. The complexity of the second step overall is O(( m i=1 size i ) • e 2 ), where size i is the size of the cluster belonging to the ith middle point. Hence, the complexity of the sweep is O(n + ( m i=1 size i ) • e 2 ). To specify O(( m i=1 size i ) • e 2 ) we investigate three cases: If all points would be middle points (m = n), then O(( m i=1 size i ) • e 2 ) = O(n•e 2

2 .

 2 If no point would be a middle point (m = 0), then O((m i=1 size i ) • e 2 ) = O(n • e 2 ) because then holds n = 0. 3. If 0 < m < n, then O(( m i=1 size i )•e 2 ) = O(n•e 2), as every point may belong to at most two clusters (just intersection points belong to two clusters) so that the sum of all cluster sizes size i has to be less than 2n.In every case the complexity of the sweep isO(n + n • e 2 ) = O(n • e 2 ).The noise removing step needs O(m + m • e 2 + m) time, as for every cluster h i two constant time operations are applied and e 2 points at four places in the rectangle around h i are investigated. Because of m ≤ n we assign a complexity of O(n + n • e 2 + n) = O(n • e 2 ) time to the noise removing. Altogether GamRec has a complexity of O(n • e 2 + n • e 2 ) = O(n • e 2 ) time. As one can see, the complexities O(n • K • t) and O(n • e 2 ) of K-means clustering and GamRec are similar efficient. In our experiments e.g. we have chosen K = 5, t = 20, e = 8, which results in O(n • 100) and O(n • 64).

  0008 m) (0.0008 m) (0.0015 m) (0.0031 m) (0.0033 m) (0.0037 m) K-means with K given 0.0678 m 0.1292 m 0.2413 m 0.5078 m 1.0707 m 0.2365 m (standard deviation) (0.0512 m) (0.0497 m) (0.1185 m) (0.3052 m) (0.5791 m) (0.4718 m)

Acknowledgments. This work is co-funded by the European Regional Development Fund project AcoGPR under the grant agreement no. WA3 80122470.


Najmeh Faridnia
email: faridnia@frec.co.ir

Haidar Samet
email: 2samet@shirazu.ac.ir

Babak Doostani Dezfuli

A New Approach to High Impedance Fault Detection Based on Correlation Functions

Keywords: high impedance fault detection (HIF), auto correlation function (ACF), partial auto correlation function (PACF), SACF, SPACF

Detection of High Impedance Fault (HIF) that is mainly occurred in MV power systems has been a challenge due to its low current magnitude. A new HIF detection approach based on correlation functions is applied on the voltage and current signals in this paper. Twelve indices based on correlation functions are implemented and tried on a lot of cases. Some of them which have clearly different values for no-fault and HIF conditions are chosen. Simulation results prove the efficiency of proposed method for HIF detection.

Introduction

The high impedance fault (HIF) is an abnormal electrical condition which may affect human lives wherever the feeders and the living areas have any overlaps. The HIF is usually occurred in the distribution feeders because of their low height in comparison with high voltage transmission lines. This fault is happened when an unexpected connection appears between a feeder and the ground by a high impedance object. Most of the time, it is generated by falling a broken conductor of an energized line on the ground or connection between tree branches to the line. In electrical point of view if a high impedance object produces a parallel loop with the load in a circuit while the line is loaded as much as nominal current, then the over-current which is imposed by fault would be low in contrast to load current. So, this phenomenon is called HIF. Since the basis of operation for conventional relays is current magnitude, HIF flow cannot be detected by mentioned relays.

A review of many approaches on HIF detection from classical to heuristically algorithms is done in [START_REF] Sedighizadeh | Approaches in High Impedance Fault Detection A Chronological Review[END_REF] that is a useful guideline.

Operating principle of HIF detection is reviewed and possible future technologies is discussed in [START_REF] Adamiak | High Impedance Fault Detection On Distribution Feeders[END_REF]. Discrete wavelet transformer technique in [2], [START_REF] Elkalashy | Modeling and Experimental Verification of a High Impedance Arcing Fault in MV Networks[END_REF] and [START_REF] Doaa Khalil | High-Impedance Fault Detection in EHV Transmission Lines[END_REF], artificial neural network in [START_REF] Khorashadi | An ANN-Based High Impedance Fault Detection Scheme: Design and Implementation[END_REF] is used to detection of HIF. Besides, some others such as [START_REF] Sharaf | High Impedance Fault Detection using Feature-Pattern Based Relaying[END_REF] and [START_REF] Sharaf | Residual Third Harmonic Detection of High Impedance Fault in Distribution System Using perception Neural Networks[END_REF] used low frequency harmonic analysis to detection.

In this paper, a new approach for HIF detection based on correlation functions is proposed. It is found that the values of correlation functions applied on voltage, current and also their derivations during the system healthy conditions are different from their values when high impedance fault occurs. To show the benefits of correlation functions in HIF detection, the paper will be continued with five sections. A brief summary of the correlation function theory and proposed indices for this purpose will be discussed in section 2. Then a case study system with using of two different models of high impedance fault is used for data collection in section 3. Simulation results of the case study system using PSCAD/EMTDC software are presented in section 4. The conclusion is given in section 5.

A brief theory of correlation functions

Auto correlation function (ACF) and partial auto correlation function (PACF) that are defined by equations 1 and 2 respectively are applied to detect any abnormal changes in a signal.

... , 3 , 2 , 1 ) ( ) )( ( 1 2 1            k Z Z Z r n t t k n t k t t k    (1)                        ... , 3 , 2 1 1 1 1 , 1 1 1 , 1 1 k if r r r r r k if r r k j j j k k j j k j k k kk (2) where 1 .., , 2 , 1 , 1 , 1 ,        k j r r r r j k k kk j k j k
r k and r kk symbolize ACF k and PACF k respectively. n is time series window length and µ is the mean of time series. To get a good performance it's better to use summation of ACF and PACF components. So, in the present approach two more functions are defined to distinguish the normal case from HIF. SACF: Summation of ACF components of a window including k number of samples that is called SACF. Related equation is as follows:

   max 1 k k k ACF r SUM ( 3
)
SPACF: This function is defined similar to SACF but the difference is that the PACF components have been used to make it. Equation (4) belongs to SPACF.

   max 1 k k kk PACF r SUM ( 4
)
Here six signals are used instead of Z in equations ( 1) and (2) which are as follows:

di n = i n -i n-1 ( 5
)
dv n = v n -v n-1 ( 6
)
d 2 i n = di n -di n-1 ( 7
)
d 2 v n = dv n -dv n-1 ( 8
)
where i and v are the current and voltage at relay placement. d and d 2 denote the first and second differences. Both SACF and SPACF functions are applied to all six above signals to obtain 12 indices as follow:

   max 1 ) ( k k k i r sacfi (9)    max 1 ) ( k k k v r sacfv ( 10
)    max 1 ) ( k k kk i r spacfi (11)    max 1 ) ( k k kk v r spacfv ( 12
)    max 1 ) ( k k k di r sacfdi (13)    max 1 ) ( k k k dv r sacfdv (14)    max 1 ) ( k k kk di r spacfdi (15)    max 1 ) ( k k kk dv r spacfdv (16)    max 1 2 2 ) ( k k k i d r i sacfd (18)    max 1 2 2 ) ( k k k v d r v sacfd (19)    max 1 2 2 ) ( k k kk i d r i spacfd (20)    max 1 2 2 ) ( k k kk v d r v spacfd ( 21
)
where r k (v) denotes kth component of voltage's ACF. Other equations are similarly defined.

Case Study system

One case study system with two different HIF models is implemented in figures 1 and 2. Case study system belongs to a radial cable feeder of a medium voltage (25 kv) system. Single line diagram of the case study system with the HIF model 1 is given in figure.1. V s and I st are the relay voltage and current which are used for sensing the HIF. The connected load to the feeder is assumed to be constant at the end of the line. HIF model 1 is explained through equation ( 22) and figure.1.

) ) ( 1 ( , 3 0 0   f f f f f i i R R mH L    ( 22
)
Where R f0 = 20, i f0 = 70, α= 0.6, β =2, V L-L = 25 kV , f=60 Hz As can be seen in figure .5, no sign of any events for model 2 is visible in both current and voltage of relay.

Implementation of new approach

After utilizing correlation functions on voltage, current and also their derivations samples, it is specified that both SACF and PSACF can be acceptable factors for HIF detection. For better perception, stem diagrams of some indices for both HIF models are shown in figure.6. k on the horizontal axis is plotted against ACF and PACF values. In this case, k max is equal to 20. As can be seen, the majority of pacfi (PACF i ) and pacfddv (PACF d 2 v ) values and all of acfdv (ACF dv ) values are positive in no fault condition while in HIF phenomenon there are both positive and negative. Briefly, feature pattern of three indices show differences for both conditions. If in each diagram the summation of all 20 data is calculated, the summation result will be a big positive number in no fault condition while it will be approximately zero in HIF event. So, it can be used for HIF detection. SACF and SPACF functions for all six of relay-side signals are introduced in section.2. Using a low-pass filter would be useful to make the curve flat to get a good and clear feature pattern for HIF detection. So, it is suggested that use a filter of degree 4 with equation number (23). Filtering of SPACF curve can be got similarly. Although voltage and current of relay don't show any problem in appearance, HIF can be detected easily by using of SACF and SPACF functions.

Conclusion

A new approach for HIF detection in distribution systems based on correlation functions is proposed in this paper. After simulating a sample system, SACF and SPACF functions were experimented on i, di, d 2 i, v, dv and d 2 v signals that resulted in twelve indices to HIF detection. After getting the results of both normal and HIF conditions for different values of variables, it was founded that SACF and SPACF operators can detect HIF phenomenon from voltage and current samples of relay.

Fig. 1 .

 1 Fig. 1. The case study system with HIF model 1 The model circuit with the description of second HIF model is given in figure.2. Figure.3 shows the case study system used in figure.1 which applied HIF model 2 in figure.2. Where i(t), u p , R and l p represent arc current, constant voltage parameter per arc length, resistive component per arc length and arc length respectively. It is considered that R (R=9Ω/cm) and l p (l p =10 cm) are constant. The fault is occurred in the middle of the feeder.

Fig. 2 .

 2 Fig. 2. Circuit and description of HIF model 2 [2]

Fig. 3 .

 3 Fig. 3. The case study system with HIF model 2The voltage and current curves of the case study system with both of HIF models which are simulated in PSCAD/EMTDC software are shown in figures 4 and 5. In these simulations, the time duration is 500 ms, the fault occurs in 202 ms and the sampling time is 50 µs.As can be seen in figure.5, no sign of any events for model 2 is visible in both current and voltage of relay.

Fig. 4 .Fig. 5

 45 Fig. 4. (a) current and (b) voltage results of case study system with HIF model 1

Fig. 6 .

 6 Fig. 6. comparing feature pattern of (a) pacfi (b) acfd2v (c) acfdv for both no-fault and HIF conditions (n=100, k=1,…,20) On the first step, the time window (including N samples with sampling time equal to 50 µs) is assumed. Two different windows are considered. The two windows are assumed to have the same length. The indices values are experimented for n=40, 80, 100, 200 and k max =10, 20, 40 in both HIF and no fault states. The values of all 12 mentioned indices for n=40 and k max =10 are given in tables I to IV.

7 .

 7 These features belong to spacfdv and spacfd 2 v for HIF model 1.

  like sacfd 2 i, sacfd 2 v and spacfi for HIF model 2 before and after filtering are brought in figure.8.

Fig. 7 .Fig. 8 .

 78 Fig. 7. Feature pattern of (a) spacfdv, (b) spacfd2v for HIF condition (time duration=3 s)

Table 1 .

 1 Current's ACF results for static window in no fault condition ( n=40 and kmax=10)

	Fault time (s)	sacfv	sacfdv	Current sacfd 2 v spacfv	spacfdv	spacfd 2 v
	0.200	6.5900	6.5085	6.3878	4.0011	4.5190	4.6588
	0.201	6.6192	6.4370	6.4153	4.0940	4.4195	4.6484
	0.202	6.6264	5.3030	6.4222	4.5049	4.9448	5.1292
	0.203	6.6248	6.2117	6.4206	4.2429	4.7193	5.0486
	0.204	6.6116	6.4814	6.4083	4.4522	4.6284	4.7855
	0.205	6.5568	6.5176	6.3564	4.4555	4.8640	5.0938
	0.206	5.9560	6.5267	5.7754	4.3321	4.6874	4.8726
	0.207	6.1057	6.5263	5.9225	4.1080	4.4797	4.7678
	0.208	6.5651	6.5162	6.3643	4.5088	4.6646	4.8468
	0.209	6.6134	6.4754	6.4099	3.8506	4.3387	4.6272
	0.210	6.6252	6.1363	6.4210	3.5469	4.3965	4.6285

Table 2 .

 2 voltage's acf results for static window in no fault condition ( n=40 and kmax=10)

	Fault time (s)	sacfv	sacfdv	Voltage sacfd2v spacfv	spacfdv	spacfd2v
	0.200	6.6228	6.3595	6.4187	5.0753	5.9513	7.0459

Table 3 .

 3 Current's acf results for static window in hif condition ( n=40 and kmax=10)

				Voltage
	Fault time (s)	sacfv	sacfdv	sacfd2v	spacfv	spacfdv	spacfd2v
	0.200	6.3466	3.0206	1.6699	0.3590	0.2429	0.6375
	0.201	2.6460	0.8672	-1.0995	-0.0993	1.0627	0.1453
	0.202	0.1420	0.1680	-0.4644	-0.8723	2.4235	0.1064
	0.203	-0.9982	0.0624	-0.1323	0.2219	1.8844	-0.3442
	0.204	-1.4295	-0.0277	-0.1417	0.2232	1.8580	-0.2611
	0.205	-1.5050	0.0226	-0.3410	-0.6238	1.7768	-0.1060
	0.206	-0.1160	0.2922	-1.0224	-0.3133	1.2205	0.5994
	0.207	4.8707	1.6456	-0.9479	-0.0176	0.6933	0.2696
	0.208	7.2593	2.4569	2.1355	-0.1606	0.4349	0.0973
	0.209	3.8737	1.6491	-1.0497	0.0396	0.8800	0.2729
	0.210	0.8931	0.2576	-1.0511	0.1691	0.8014	0.5908

Table 4 .

 4 voltage's acf results for static window in hif condition ( n=40 and kmax=10) These tables clearly show that all indices values are positive and almost large for both current and voltage in no-fault case, while most of them are around zero in HIF phenomenon. In this stage, sacfv, sacfi and sacfdi indices are put aside because their values in HIF are far from zero and they are not suitable for fault detection. Now, it's time to go one step ahead and analyze the results of moving time window for different values of n and k. Finally, feature patterns of two indices for HIF with clear changes are shown in figure.

				Voltage
	Fault time (s)	sacfv	sacfdv	sacfd2v	spacfv	spacfdv	spacfd2v
	0.200	6.9339	1.6671	-0.4695	0.2669	0.1739	0.0934
	0.201	5.4783	-0.9503	-0.2817	0.5827	0.0167	0.7039
	0.202	3.4873	-0.4204	-0.2167	0.6335	0.0104	0.2035
	0.203	0.9900	-0.0886	-0.0216	0.7565	-0.4900	-0.1591
	0.204	0.5949	-0.0708	-0.0063	0.6957	-0.5060	-0.8871
	0.205	3.2958	-0.1868	-0.0500	0.5217	-0.3864	-1.9170
	0.206	5.4064	-0.8062	-0.5091	-0.1666	0.6426	-0.8450
	0.207	6.5456	-1.2097	-0.8134	0.5620	0.1342	0.6332
	0.208	7.0033	0.8362	-0.0565	0.4719	0.1510	0.6231
	0.209	6.0259	-0.7936	0.2544	0.5499	0.0080	0.5214
	0.210	4.2548	-0.8951	-0.3937	-0.0565	0.5887	0.1744