Deepthi N Ratnayake
email: d.ratnayake@londonmet.ac.uk

Hassan B Kazemian
email: h.kazemian@londonmet.ac.uk

Syed A Yusuf
email: s.yusuf@londonmet.ac.uk

Azween B Abdullah
email: azweenabdullah@petronas.com.my

An intelligent approach to detect Probe Request attacks in IEEE 802.11 networks

Keywords: IEEE 802.11, DoS Attacks, Probe Request Flooding Attacks, Wireless. Supervised Feedforward Neural Network

In Wireless Local Area Networks (WLAN), beacon, probe request and response messages are unprotected, so the information is visible to sniffers. Probe requests can be sent by anyone with a legitimate Media Access Control (MAC) address, as association to the network is not required at this stage. Legitimate MAC addresses can be easily spoofed to bypass Access Point (AP) access lists. Attackers take advantage of these vulnerabilities and send a flood of probe request frames which can lead to a Denial-of-Service (DoS) to legitimate stations. This paper discusses an intelligent approach to recognise probe request attacks in WLANs. The research investigates and analyses WLAN traffic captured on a home wireless network, and uses supervised feedforward neural network with 4 input neurons, 2 hidden layers and an output neuron to determine the results. The computer simulation results demonstrate that this approach improves detection of MAC spoofing and probe request attacks considerably.

Introduction

The wireless technology today comes in several forms and in a multitude of solutions to provide availability and security. However, many risks remain unmanaged [START_REF] Goel | An Improved Method of Detecting Spoofed Attack in Wireless LAN[END_REF]. IEEE 802.11 is a set of standards specified by Institute of Electrical and Electronic Engineers (IEEE) for WLAN computer communication. IEEE 802.11 was first created in 1997 and improved over the years. IEEE 802.11w-2009 is currently the most powerful security standard available for WLAN users [START_REF] Sood | Discover how using the IEEE standards approach plugs vulnerabilities and thwarts attacks[END_REF][START_REF][END_REF]. The MAC layer of the 802.11 protocol is based on the exchange of request/response messages i.e. each request message sent by a station (STA) must be responded with a response message sent by the AP. Probe Request Flood (PRF) attacks are designed to take advantage of this request and respond design flaw [START_REF] Bernaschi | Access points vulnerabilities to DoS attacks in 802.11 networks[END_REF]. Flooding attacks cause serious performance degradation or prevent legitimate users from accessing network resources such as the bandwidth, access points, gateways, servers and target user systems. This vulnerability is increased due to the unprotected beacon or probe request and probe response frames which can be read by anyone to learn about the network.

We learned that before an attack, the attacker actively or passively monitors the network to learn vital network information. MAC address spoofing is the next step. Therefore, we recognised that any Wireless Intrusion Detection System (WIDS) should address these initial stages of an attack before moving on to more advance steps. After analysing the previous research work and the progress of IEEE 802.11 sub committees, it is understood that there is a gap of knowledge to develop a realistic WIDS that could detect MAC spoofing and probe request attacks on IEEE 802.11 networks. This research analyses real-time traffic captured on a wireless home network. This research works with real WLAN traffic as opposed to data from a sample database or synthetic traffic generated by a test bed used in many studies. Our work aims to detect an attack during an early stage of the communication. During our initial experiments, we observed that WLAN traffic pattern is usually unpredictable and also depends on the usage, operating system and applications of the user. Further, the monitoring STA can miss many frames due to its traffic load, or receive them out of order due to packet delay, packet jitter and lost packet or prioritisation services of network traffic such as Quality of Service (QoS) [START_REF] Rumín | VoIP over WLAN 802.11b simulations for infrastructure and adhoc networks[END_REF]. These inherent complexities and unpredictable nature of data made this research a good candidate for Artificial Neural Networks (ANN). Additionally, WLAN traffic and parallel processing nature of ANNs cause a considerable amount of overhead on the monitoring STA and therefore, can affect the performance of the monitoring STA. This research analyses only 4 parameters to detect an attack. This considerably reduces the overhead of the monitoring machine whilst producing the results expected.

The rest of the paper is organized as follows. Section 2 discusses the related work. Section 3 discusses the IEEE 802.11 security policy agreement and the basic concepts behind probe request attacks. Section 4 explains the philosophy behind our research. Section 5 presents WLAN environment of the experiment and computer simulation prototype and Section 6 discusses simulation results.

Related work

Intrusion detection is to identify an unauthorised user trying to gain access, or has already gained access or compromised the computer network [START_REF] Karygiannis | Wireless network security, 802.11, bluetooth and handheld devices: recommendations of the national institute of standards and technology[END_REF]. Many researchers have worked in this area looking for possible solutions. For example, [START_REF] Bicakci | Denial-of-Service attacks and countermeasures in IEEE 802.11 wireless networks[END_REF] presents a detailed review of most popular non-intelligent methods of detecting and preventing DoS attacks in MAC layer and [START_REF] Bansal | Non-cryptographic methods of MAC spoof detection in wireless LAN[END_REF] evaluates some commonly used non-cryptographic methods of MAC spoof detection in WLANs. They identify use of cryptography, sequence number analysis, Time Difference Of Arrivals (TDOA), decreasing re-try limits, and Network Interface Card (NIC) profiling: Signal Strength Indicator (SSI) and Radio Frequency (RF) finger printing for detecting and preventing PRF attacks. Cryptography may be the most reliable solution. But, it is expensive, may require a protocol repair and can easily be a DoS target itself. [START_REF] Malekzadeh | Security improvement for management frames in IEEE 802.11 wireless networks[END_REF] proposes security improvement in management frames by a shared key. However, this solution requires a change in the wireless card. A hardware upgrade is an unrealistic solution considering the number of wireless cards that will have to change.

Detection of spoofed frames plays a major role in detection of other attacks including probe request attacks. [START_REF] Guo | Sequence Number-Based MAC Address Spoof Detection[END_REF] introduces an algorithm to detect MAC spoofing based on sequence number gaps by leveraging the structure and behaviour of the sequence number field. [START_REF] Madory | New Methods of Spoof Detection in 802.11b Wireless Networking[END_REF] introduces time difference between consecutive frames and a sliding window of received signal strengths for spoof detection. [START_REF] Li | Detecting Spoofing and Anomalous Traffic in Wireless Networks via Forge-Resistant Relationships[END_REF] utilise a combination of window of sequence numbers and traffic inter-arrival statistics (FRR -Forge Resistance Relationship Method) to detect spoofing and anomalous traffic in wireless networks. [START_REF] Goel | An Improved Method of Detecting Spoofed Attack in Wireless LAN[END_REF] argue that these solutions work only when both the attacker and victim are transmitting and, also may be difficult to differentiate an attacker from a victim, when the victim is offline. They improved [START_REF] Li | Detecting Spoofing and Anomalous Traffic in Wireless Networks via Forge-Resistant Relationships[END_REF] solution by utilising transmission rate and by sending a probe request after every 9th frame. However, this solution generates an additional overhead on the network. [START_REF] Faria | Detecting identity-based attacks in wireless networks using signal prints[END_REF] proposes detecting identity-based attacks in wireless networks using only signal prints. However, this solution is ineffective when there is a single AP serving all STAs. Further, RSSI measurements by itself may not distinguish a genuine STA from an adversary if they are too close to each other [START_REF] Bicakci | Denial-of-Service attacks and countermeasures in IEEE 802.11 wireless networks[END_REF]. Above discussed non-intelligent WIDS methods use statistical, rule based, expert knowledge or pattern recognition approaches on known vulnerabilities or attack signatures and therefore consumes time, lacks flexibility to adaptation to environmental changes, and eventually becomes out-dated.

[14] presents a comparative analysis of IDS approaches. ANN is currently the most established approach for IDS considering the unpredictable behaviour of WLAN networks and attacks. Most intelligent IDSs for TCP/IP networks use Self Organizing Maps, Artificial Immune systems, Fuzzy Logic and Neural models, Adaptive Neural-Fuzzy Inference Systems and hybrid models. [START_REF] Lim | Wireless intrusion detection and response[END_REF] introduces a prototype of a standalone WID and response system based on NetStumbler attacks. This solution detects attacks only by calculating probe requests per second. It also responds to the attacker with a DoS attack in return, which can lead to attacking a own user. [START_REF] Pleskonjic | Wireless Intrusion Detection Systems (WIDS), 19th Annual Computer Security Applications Conference[END_REF] presents a corporative detection architecture based on intelligent agents with power of autolearning, incorporating NN and Fuzzy logics. This large and complex system was never implemented according to our knowledge. [START_REF] Yang | Intrusion detection solution to WLANs[END_REF] also propose a distributed and collaborative architecture using IDS agents but falls short of an implementation. [START_REF] Dasgupta | MMDS: Multilevel Monitoring and Detection System[END_REF] presents a multi-agent architecture using fuzzy decision support system which performs anomaly detection for ad-hoc and infra-structure networks. This solution does real time monitoring, but the solution is based only on sequence number anomalies. [START_REF] Ataide | An Architecture for Wireless Intrusion Detection Systems Using Artificial Neural Networks[END_REF] discusses a range of research architectures and open source and commercially available WIDSs. They propose a comparatively complex architecture for WIDS using ANNs based on real time data and tests virtual carrier sense, association flood, and de-authentication attacks. The solution focuses on the behaviour of the complete network, which can be challenging in a real network which has larger number of users.

The research observed that many of these solutions are designed based on non-real data sets and/or also identifies intrusive behaviours based on the exploration of known vulnerabilities [START_REF] Chavan | Adaptive neurofuzzy intrusion detection systems[END_REF][START_REF] Toosi | A Neuro-Fuzzy Classifier for Intrusion Detection Systems[END_REF][START_REF] Ahmad | Application of artificial neural network in detection of DOS attacks[END_REF]. Further, it is observed that some of these solutions are extremely complex and are simulated and tested without considering the practical implementation and computing power they may require. Therefore, these solutions are limited for academic research world as implementing is too complex or expensive.

WLAN security and probe request attacks

Fig. 1 shows the security policy agreement phase of IEEE 2007. A STA seeking to connect to a WLAN has a choice of passive scan or active scan. In passive scan, STA listens to successive channels, waiting to hear a beacon frame from an AP, while in an active scan, the STA broadcasts a probe request message on every channel its physical layer supports, until the STA finds an AP. These frames are unprotected and information passed between the frames can be read using freely available software like Wireshark [START_REF] He | Security analysis and improvements for IEEE 802.11i[END_REF]24,[START_REF] Ahmad | Application of artificial neural network in detection of probing attacks[END_REF]. APs keep a list of legitimate MAC addresses that can access its services to prevent unauthorised access. However, MAC addresses can easily be spoofed using ifconfig, macchanger (Linux) or using SMAC2 (Windows) to pretend it is a legitimate STA. Association to a network is not required to probe and receive a response. Hence, an adversary only requires a legitimate MAC address to send Probe Requests. Usually, probing is the initial phase of any other attack in computer networks [START_REF] Ataide | An Architecture for Wireless Intrusion Detection Systems Using Artificial Neural Networks[END_REF].

Security Policy Agreement

Access Point

[24] defines three frame types namely Management, Control and Data. The management frames establish and maintain communications. The control frames help in the delivery of data. The data frames encapsulate the Open System Interconnection (OSI) network layer packets. Each frame consists of a MAC header, frame body and a Frame Check Sequence (FCS). MAC header comprises of frame control, duration, address, sequence control information, and for QoS data frames, QoS control information. Frame body contains information specific to the frame type and subtype. FCS contains an IEEE 32-bit Cyclic Redundancy Check (CRC).

The below is a simple analysis of data captured during an attack on a test-bed.

Sequence number of the MAC frame is a 12-bit counter that starts from 0 when a NIC starts or resets, and wraps on 4095 at the overflow. Theoretically, a NIC can generate only one set of sequence numbers at a time [24]. However, Fig. 2 shows several parallel sequence number patterns generated from the same MAC address. The straight sharp lines (starting at 17.01) were formed due to a high frequency of sequence numbers generated during a spoofed attack whilst other fluctuating and scattered lines are from the genuine user (16.39-17.08). These parallel sequence number patterns may have generated due to QoS or packet delays [START_REF] Rumín | VoIP over WLAN 802.11b simulations for infrastructure and adhoc networks[END_REF]. Frame sub-type of the MAC frame identifies the type of the data frame [24]. Fig. 3 illustrates a high occurrence of probe request frames which complements the attack identified by the sequence number analysis in Fig. 2. Spoofed attacker cannot associate with the network without knowing the network key. Hence, other framesub-types have a lower frequency. SSI of a frame captured by Wireshark provides an indication of received transmission power of a NIC, which also gives an indication of the location [START_REF] Bansal | Non-cryptographic methods of MAC spoof detection in wireless LAN[END_REF]. Therefore, SSI patterns are useful in detecting spoofed attacks. In Fig. 4, from 16.36-17.00, shows the SSI pattern generated by the genuine user. Unusual SSI patterns were generated during the spoofed attack (from 17.01 to 17.08). The delta time value of a Wireshark captured frame indicates the time since the previous packet was captured. This is also useful in detecting attacks, as it gives an indication of the server response time, network round-trip time, and other delays.

Analysing frames of a WLAN test bed manually or statistically and detecting probing attacks are possible due to its controlled nature. However, a WIDS should be able to capture and analyse frames, and detect attacks automatically in a live environment that is unpredictable by nature. Therefore, after considering different models and their possible realistic and efficient application on detection of probe request attacks, the research utilised a supervised feed-forward NN architecture with 4 input neurons, 2 hidden layers with 20 sigmoid hidden neurons and an output layer with one linear output neuron which classify genuine frames from rogue frames.

WLAN environment and computer simulation prototype

This research designed to capture delta time value, sequence number, received signal strength and frame sub-type of the packets transmitted between an AP, users, and attackers of a wireless home network with 8 user stations. Note: Other wireless stations are not shown here Data capturing was performed for 75 minutes. User Test1-PC accessed internet to browse information, download software, watch a live TV channel, listen to a live radio channel and check/send emails. Attacker Test2-PC with its spoofed MAC address sent a flood of probe request frames to the AP. Both user and attacker performed start-up and shut down procedures, network scans, a network connect and disconnect procedures, and a NIC repair. The captured data set consisted of 190K frames. Frames generated from MAC address 00:0c:f1:5b: dd:b3 were filtered to prepare a training sample and labelled manually as rouge or genuine.

The prototype to detect probe request attacks using neural networks is designed and trained using MATLAB. The training sample consisted of approximately 175K frames that were generated from MAC address 00:0c:f1:5b:dd:b3. This was filtered to obtain delta time, sequence number, SSI and frame sub-type of each frame and fed into 4 input neurons. This sample was randomly divided and 70% of the data was used to train the network. 15% each was used for validation and testing respectively.

The network is trained using Levenberg-Marquardt back propagation algorithm. The mean squared errors of training, validation and testing were 3.86409e-3, 3.75033e-3 and 3.67129e-3 respectively. Training, validating and testing were converged outstandingly resulting overall regression value 0.98043.

Simulation results and discussion

The trained NN was tested using pre-defined scenarios (Table 1(. The results were generated based on 1000 frame samples. Target output of a user's frame considered as 1 (genuine), whilst a frame from an attacker considered as 0 (not-genuine). The results in Table 1 shows that the trained NN is capable of detecting known attacks 99%-100% (Cap24-28) and unknown attacks 89% (Cap30).The NN is also tested against the mobility effect of an attacker and a user. The user's mobility within the signal range does not affect the detection rate. However, when the attacker was at the same location as the user, the detection rate dropped to 75%. The detection rate of genuine frames of a user is about 98-100%. The detection rate drops when a genuine user scans a network excessively which generates unusually a large number of probe requests. This can occur due to an ill-configured WLAN card, a weak signal strength or user deliberately scanning the network which may require network administrator's attention. However, the issue can be solved within the system by setting a threshold value of warnings to be tolerated per second to suit to specific users or network. The research considered choosing only few variables to develop a WIDS considering the large data volume and computation power required for real-world implementation. Sequence number and frame sub-type of a MAC frame and signal attributes SSI and delta time values are the 4 independent variables that were carefully chosen by the research based on the following situations; Manipulation of these variables is nearly impossible. Some argue that attackers use sequence number synchronising software to generate sequence number patterns to match with the user STA. However, the precision and effectiveness of this technique is doubtful as they cannot predict the behaviour of the user STA, whether the STA is starting or resetting its NIC card or transmitting data or idling. Some argue that SSI can be manipulated by attackers controlling the NIC signal strength, moving close to the user STA or signal strength can simply fluctuate due to environmental factors. This is also more of a theoretical issue than a practical one as the attacker cannot perform all these activities in a WLAN without exposing itself. Some argue that excessive probe request frames can be generated by an ill-configured STA or a genuine user repeatedly attempting to log-in to the AP. In this instance, network administrator can correct if there are any problems with a genuine STA or the user. This solution also works when the genuine user is offline. Finally, each of the individual variables has the potential of indicating the possibility of a spoofed STA and a probe request attack. However, it is expected to be supported by other 3 variables if there is an uncertainty.

Conclusion:

This research has been carried out to identify an external attacker by analysing the traffic generated from a user MAC address in a single frequency band of a Wireless Local Area Network. A supervised feed-forward neural network with four distinct inputs, delta-time, sequence number, signal strength and frame sub-type is applied to identify and differentiate a genuine frame from a rogue frame. The experimental results show that the use of neural network can detect probe request attacks to a very high precision. This solution also allows WLAN users to be mobilised freely within the signal area. Further research will be conducted to enhance this experiment by monitoring the entire bandwidth, using coordinated multiple monitoring stations.

Fig. 1 :

 1 Fig. 1: Security policy agreement phase of IEEE 2007 [24].

Fig. 2 :

 2 Fig. 2: Analysis of sequence numbers generated by a single MAC address

Fig. 3 :Fig. 4 :

 34 Fig. 3: Analysis of frame types generated by a single MAC address

Fig. 5 :

 5 Fig. 5: WLAN including two attackers and a network monitor.

Fig. 5

 5 Fig.5shows the components of the wireless home network relevant to the research. Wireshark network monitoring software was used for data capturing. Ifconfig and SMAC2 were used to change the MAC address of the attackers Test2-PC and Test4-PC. Test3-PC was kept on promiscuous mode. Monitoring was restricted to IEEE 802.11 WLAN channel number 11 -2462 MHz due to heavy frame loss experienced when capturing on all channels. Therefore, monitoring statistics of the entire bandwidth for STA's behaviour on the entire bandwidth is unavailable.

Table 1 .

 1 Summary of tests conducted.

	Capture Code	Test Scenario	Detection Rate (genuine/rouge)	False +ve rate	False -ve rate
	Cap21	Unseen data set from user	96%	n/a	4%
	Cap29	Unseen data set from user far away from AP	98%	n/a	2%
	Cap24	Unseen data set from attacker using NetStumbler	99%	1%	n/a
	Cap25	Unseen data set with attacker closer to user	75%	25%	n/a
	Cap26	Unseen data set with attacker far away from user	100%	0%	n/a
	Cap27	Unseen data set with new attacker (attacker 2) using	99%	1%	n/a
		NetStumbler
	Cap28	Unseen data set with 2 attackers using NetStumbler	100%	0%	n/a
	Cap30	Unseen data set with an attacker using a Linux network	89%	11%	n/a
		scanning tool


Dimitris Voukantsis

Kostas Karatzas

Siegfried Jaeger

Uwe Berger

Personalized Information Services for Quality of Life: The Case of Airborne Pollen Induced Symptoms

Keywords: Allergy, Computational Intelligence, Personalized Health Services

Allergies due to airborne pollen affect approximately 15-20% of European citizens; therefore, the provision of health related services concerning pollen-induced symptoms can improve the overall quality of life. In this paper, we demonstrate the development of personalized quality of life services by adopting a data-driven approach. The data we use consist of allergic symptoms reported by citizens as well as detailed pollen concentrations of the most allergenic taxa. We apply computational intelligence methods in order to develop models that associate pollen concentration levels with allergic symptoms on a personal level. The results for the case of Austria, show that this approach can result to accurate and reliable models; we report a correlation coefficient up to r=0.70 (average of 102 citizens). We conclude that some of these models could serve as the basis for personalized health services.

Introduction

Allergy due to airborne pollen is a reaction of the human immune system to certain allergens carried by the pollen grains. It is estimated that approximately 15-20% of the European citizens suffer from pollen-related allergies [START_REF] Huynen | Phenology and Human Health: Allergic Disorders[END_REF]; additionally, there are well-established, yet not fully understood, associations with certain respiratory diseases, e.g., asthma [START_REF] Taylor | Links between Pollen, Atopy and the Asthma Epidemic[END_REF]. During the last years, there has been an increasing trend in pollen-induced allergic symptoms, while the World Allergy Organization [START_REF]World Allergy Organization[END_REF], reports an increase in severity of the symptoms.

In this context, there is an emerging need for quality of life information services that could assist sensitized citizens in making decisions concerning their daily life. This issue has already been partially addressed by existing information systems that provide information of pollen concentration levels of allergenic taxa [4]. However, the severity of allergic symptoms is strongly dependent on the citizen under consideration. Therefore, a more personalized approach that addresses symptoms rather than pollen concentration levels seems more appropriate.

The design and development of personalized information services addressing citizen-specific symptoms is a challenging task due to the lack of a generic (universal) mechanism responsible for the triggering of allergic symptoms. Moreover, airborne pollen concentration levels have not been regulated as in the case of chemical air pollutants [START_REF] Cafe -Derective | EC of the European Parliament and of the Council. 6. Pollen diary[END_REF], thus there is no regulatory or management framework, based on common limit values and threshold levels to serve as a reference for any health related decision. In order to avoid these difficulties, we adopt a data-driven modeling approach by utilizing the database of pollen-induced symptoms reported by citizens themselves through the Pollen Diary system [6]. Additionally, we use detailed pollen concentration levels, monitored in several sampling sites, for the area of interest. We demonstrate the development of such models for the case of Austria, using data collected during the year 2010. The resulting personalized models could serve as the main modeling tool of existing or future information systems, providing citizenspecific warnings concerning the occurrence and severity of airborne pollen-induced symptoms.

Materials and Methods

Data Presentation

The modeling approach adopted in this study has been based in two distinct databases. i) Pollen concentration data of 65 distinct taxa sampled at several monitoring sites across Austria during 2010, and ii) Pollen-induced symptoms reported by users of the Pollen Diary system. The total number of users was 716, with an average of 57 records per user. The symptoms were reported in detail (eye itching, eye watering, nose, etc.); however, in this case we have identified the "overall symptoms" as the main parameter of interest. The latter one ranges from 0 (no symptoms) up to 25 (strong symptoms in eyes, nose and lungs). The data were preprocessed excluding taxa with more than 10% missing values, resulting to a final of 20 taxa (Table 1). Furthermore, users with less than 100 records were not included in this study, resulting in a total of 102 users. The final preprocessing step was to normalize the data using variance scaling. The latter one results to variables with average value μ=0 and standard deviation σ=1. Fig. 1 presents the overall symptoms (average of all citizens) and average pollen concentrations (average of all taxa) as a function of time. It is evident that during spring (days 80-130) peak values of overall symptoms and pollen concentrations coincide.

Computational Intelligence Methods

In order to model the inter-relations between pollen concentrations and citizenspecific symptoms, we have applied a series of Computational Intelligence (CI) methods. The task under consideration has been addressed as a regression problem (the target was the numerical value of "overall symptoms"). Therefore the proper versions of the following CI methods were considered:

Multi-Layer Perceptron (MLP). MLP [START_REF] Haykin | NeuralNetworks: A Comprehensive Foundation[END_REF] models are probably the most popular artificial neural network models, capable of modeling complex and highly non-linear processes. In this case, we have evaluated network architectures consisting of one hidden layer, setting the number of neurons equal to N=(inputs+outputs)/2 [START_REF] Haykin | NeuralNetworks: A Comprehensive Foundation[END_REF]. Where inputs is the number of input variables of the model and outputs the number of output variables; in this case, outputs=1. Furthermore, we have use the Levenberg-Marquardt backpropagation algorithm during the training of the models, with a maximum number of epochs at 500. Support Vector Regression (SVR). SVR is a modification of the original algorithm introduced by Vapnik [8] capable of tackling regression problems. In this case we have chosen the SVR kernel to be the Radial Basis Function (RBF). Furthermore, we optimized the regularization parameter (C) and the ε-intensive zone of the model, as well as the γ parameter of the kernel. The SVR models have been developed using the Spider Toolbox for Matlab [START_REF]Spider Toolbox for Matlab[END_REF].

Least Squares Support Vector Regression (LS-SVR). LS-SVR is a modification of support vector machines [START_REF] Suykens | Least Squares Support Vector Machines[END_REF]. In this case, we have used the RBF kernel and optimized the regularization (C) and kernel parameter (γ) using leave-one-out cross-validation. The LS-SVR models were developed using the LS-SVMlab toolbox for Matlab [START_REF] Ls-Svmlab | [END_REF].

k-Nearest Neighbors (kNN). kNN is non-parametric CI algorithm that can be used to tackle regression problems. In this case, we have set the number of nearest neighbors to be k=5 and we applied the 1/r (r is the distance to the neighbor) weighting scheme.

Multiple Linear Regression (MLR). MLR is a traditional statistical method that provided the reference frame for comparison of the other CI methods.

Feature Selection

The selection of the appropriate features as input parameters of the models is an important modeling step as irrelevant or noisy features may result to poor performing models. Mutual Information (MI) as well as other information related criteria have been used in the past in order to identify relevant input variables. In this case, we have used MI following the approach presented in [START_REF] Peng | Feature selection based on mutual information: criteria of max-dependency, max-relevance, and min-redundancy[END_REF].

Training and Evaluation

Due to the limited number of available instances per citizen, the training of the models was based on the 10-fold cross validation, in order to utilize all instances. The Root Mean Square Error (RMSE), the Correlation Coefficient (r), as well as the Index of Agreement (d) have been used to validate the models. The latter index is defined by the following formula:

( ) 2 2 1 ∑ ∑ - + - - - = i i i i i i a a a p a p d (1)
where i p refers to predicted values and i a to observed ones, whereas with p and a are denoted the average of the predicted and observed values, respectively. The index of agreement ranges from 0.0 (theoretical minimum) to 1.0 (perfect agreement between observed and predicted values).

Results and Discussion

Table 2 summarizes the performance of the resulting 102 citizen-specific models. The results are presented as averages (and standard deviations) of the statistical indices presented in section 2.4. It is clear that the performance differs significantly depending on the algorithm used to develop the model. The non-parametric kNN models as well as the LS-SVR models indicate superior performance compared to the SVR and MLP models. This may be attributed to the limited number of training examples in the data set (at least for the MLP models). Furthermore, all CI methods perform on average better than the reference (MLR) method. This indicates that CI algorithms are more suitable to model the pollen-induced symptoms.

The results presented in Table 2 indicate that modeling pollen-induced symptoms on citizen level can result to models of acceptable performance. In contrast, a global approach, i.e., using the same data and algorithms to build a model for all citizens, does not result to acceptable performances (d<0.4, r<0.25, rmse>4). The prevalence of allergic symptoms in humans is a highly complex process that dependants on several factors such as pollen concentrations, meteorological and chemical (i.e. air quality) weather conditions, citizen habits (outdoor activity, travelling, medication). Therefore, in some cases the available data cannot result to accurate models. This can be demonstrated in more detail for the case of two specific citizens (Citizens with ID number 80 and 85, hereafter denoted as Cit.#80 and Cit.#85), with similar data records and maximum overall symptoms. Table 3 presents the performance of the two citizen-specific models (Cit.#80, Cit.#85). The results show that the allergic symptoms indicated by Cit.#85 can be successfully modeled; however, this is not the case for Cit.#80. Fig. 3 presents the overall symptoms for both citizens, as a function of time. It is evident that the pattern of symptoms for both citizens is different. In this paper, we have demonstrated the use of CI methods in order to develop models capable of estimating pollen induced symptoms on a personal level. The results show that the data collected by the Pollen Dairy system can provide a solid foundation to build accurate citizen-specific models in most of the cases. The latter ones can serve as the basic modeling tool to develop personalized health services. Future work in this field includes i) the use more databases and more advanced CI methods in order to improve the performance of the personalized models and ii) the use of the Pollen Diary database in order to categorize citizens into several sensitized groups and identify common characteristics. These goals will advance the accuracy and interpretability of the models, therefore providing a more solid background for the development of personalized health services.

Fig. 1 .

 1 Fig. 1. Overall symptoms (average of all citizens, solid line) and pollen concentration (average of all taxa, dash-dot line). Both variables are normalized, using variance scaling (μ=0, σ=1).

Fig. 2 .

 2 Fig. 2. Overall symptoms (solid line: predicted, dash-dot line: observed) for Cit.#80 and Cit.#85

Table 1 .

 1 Variables included in the preprocessed data set.

	Temporal Variables Pollen Concentration	Symptoms (target)
	Day of Year	Acer, Alnus, Ambrosia,	Overall Symptoms
		Betula, Carpinus, Corylus,
		Cupress, Fraxinus, Juglans,
		Pinus, Plantago, Platanus,
		Poaceae, Populus, Quercus,
		Rumex, Salix, Sambucus,
		Tilia, Urtica

Table 2 .

 2 Performance (average and standard deviation) of 102 citizen-specific models. The unit of RMSE is the number of symptoms.

	Model	d	r	RMSE
	MLR	0.74 ± 0.14 0.61 ± 0.20 2.26 ± 1.05
	MLP	0.76 ± 0.15 0.63 ± 0.22 2.17 ± 0.93
	kNN	0.80 ± 0.11 0.67 ± 0.17 2.07 ± 0.86
	SVR	0.74 ± 0.14 0.64 ± 0.19 2.11 ± 0.88
	LS-SVR 0.79 ± 0.14 0.70 ± 0.18 1.92 ± 0.80

Table 3 .

 3 Typical detailed performance (10-fold cross validation) of models for two citizens.

	ID	Model	MLR	MLP	kNN	SVR	LS-SVR
	Cit.#80 d	0.67 ± 0.16	0.73 ± 0.16	0.74 ± 0.15	0.71 ± 0.13	0.79 ± 0.12
		r	0.54 ± 0.23	0.65 ± 0.26	0.59 ± 0.28	0.61 ± 0.19	0.68 ± 0.21
		rmse	4.51 ± 0.95	4.16 ± 1.33	4.44 ± 1.18	4.48 ± 1.08	3.83 ± 1.20
	Cit.#85 d	0.92 ± 0.09	0.95 ± 0.06	0.97 ± 0.02	0.94 ± 0.02	0.95 ± 0.03
		r	0.88 ± 0.15	0.93 ± 0.06	0.94 ± 0.04	0.92 ± 0.04	0.93 ± 0.05
		rmse	2.87 ± 0.90	2.30 ± 1.02	2.08 ± 0.78	2.70 ± 0.35	2.34 ± 0.67

Acknowledgments. The authors greatly acknowledge the contribution of all colleagues involved in the pollen concentration data collection. Manfred Bobek, Sigmar Bortenschlager, Inez Bortenschlager, Uschi Brosch, Edith Bucher, Bernard Clot, Pramod Harvey, Veronika Kofler, Andreja Kofol-Seliger, Herta Koll, Sabine Kottik, Margit Langanger, Rudolf Litschauer, Karol Micieta, Anna Paldy, Ondrej Rybnicek, Hanna Schantl, Jutta Schmidt, Roland Schmidt, Ingrid Weichenmeier, Helmut Zwander.


Emmanouil Hourdakis

Panos Trahanias

Improving the classification performance of Liquid State Machines based on the separation property

Keywords: Liquid State Machines, Separation, Fisher Discriminant Ratio

Liquid State Machines constitute a powerful computational tool for carrying out complex real time computations on continuous input streams. Their performance is based on two properties, approximation and separation. While the former depends on the selection of class functions for the readout maps, the latter needs to be evaluated for a particular liquid architecture. In the current paper we show how the Fisher's Discriminant Ratio can be used to effectively measure the separation of a Liquid State Machine. This measure is then used as a fitness function in an evolutionary framework that searches for suitable liquid properties and architectures in order to optimize the performance of the trained readouts. Evaluation results demonstrate the effectiveness of the proposed approach.

Introduction

Liquid State Machines (LSMs) are based on the concept of using biologically realistic neural networks for processing continuous input channels of information [START_REF] Maass | Real-time computing without stable states: A new framework for neural computation based on perturbations[END_REF][START_REF] Natschlager | Computer models and analysis tools for neural microcircuits[END_REF]. Their powerful computational abilities can be attributed to the capacity of these networks to transform the temporal dynamics of an input signal into a high dimensional spatiotemporal pattern, which preserves recent and past information about the input. This information can be retrieved with a high degree of accuracy using certain types of classifiers. Following their introduction [START_REF] Maass | Real-time computing without stable states: A new framework for neural computation based on perturbations[END_REF], LSMs have been used in various pattern classification tasks, including speech recognition [START_REF] Verstraeten | Isolated word recognition with the liquid state machine: a case study, 13th European[END_REF] and movement prediction [START_REF] Burgsteiner | Movement prediction from real-world images using a liquid state machine[END_REF]. The notion behind LSMs has also been extended to problem domains outside computational modeling, where researchers use physical mediums for the implementation of the liquid, such as a bucket of water [START_REF] Fernando | Pattern recognition in a bucket[END_REF] or real cell assemblies [START_REF] Dockendorf | Liquid state machines and cultured cortical networks: The separation property[END_REF]. An LSM consists of three components: (i) the liquid, i.e. a pool 𝑀 of spiking neurons that accepts input from different sources and outputs a series of spike trains, (ii) a filter 𝐿 that is applied on the output of the liquid in order to create a state matrix 𝑆, and (iii) one or more memoryless readout maps that are trained to extract information from 𝑆. The main conception behind this setup is that the complex dynamics of the input are transformed by the liquid to a high dimensional space, in a way that preserves their recent and past properties. This can be compared to a pool of water with a stone thrown in it. The disturbances that are caused in the liquid could be used by a trained observer to deduce the properties of the motion of the stone before entering the water. To improve the classification performance of an LSM, one must ensure that for two different input histories, the liquid states produced are significantly different [START_REF] Legenstein | Edge of chaos and prediction of computational performance for neural circuit models[END_REF]. This property, known as separation, has recently received increased attention in the literature, due to its close correlation with the performance of the LSM. In [START_REF] Maass | Real-time computing without stable states: A new framework for neural computation based on perturbations[END_REF] the separation between two different liquid states is calculated by measuring the Euclidean distance of their state vector, i.e. the filtered neuron output sampled at one time instance. A similar geometric interpretation has been given by [START_REF] Goodman | Spatiotemporal pattern recognition via liquid state machines[END_REF], who measure the separation of the liquid as the Euclidean distance between the centroids of the states that belong to different classes. In [START_REF] Dockendorf | Liquid state machines and cultured cortical networks: The separation property[END_REF], the authors use spike train distance metrics instead of Euclidean distance. From the perspective of a classification system, it has been suggested that the rank of the state matrix 𝑆 can be used to measure the quality of the liquid [START_REF] Legenstein | Edge of chaos and prediction of computational performance for neural circuit models[END_REF]. According to this measure, the larger the number of linear independent variables produced by a liquid state, the better the classification that can be performed by the LSM. Attempts to improve the performance of an LSM in the literature have shown that it is very difficult to devise a proper measure or structural criterion to optimize the quality of the liquid. For example in [START_REF] Matser | Structured liquids in liquid state machines[END_REF] the authors have concluded that randomly generated liquids outperform any attempt to structurally modify the LSM. In [START_REF] Kok | Liquid State Machine Optimization[END_REF] the authors use both reinforcement learning and genetic algorithms in order to optimize the classification performance of the LSM. Finally in [START_REF] Norton | Improving the performance of liquid state machines through iterative refinement of the reservoir[END_REF] the authors use the centroid separation measure in order to drive the synaptic modification of the LSM. In the current paper we propose a criterion that measures the separation of the liquid states that correspond to different classes based on the class means and variances. The classification performance of an LSM is subsequently improved by employing an evolutionary framework to minimize the introduced criterion. Experimental results attest on the performance and accuracy of the proposed approach.

Separation Measure

Similarly to Support Vector Machines, the liquid of the LSM acts as a kernel that transforms the low-dimensional space of an input signal to the spatio-temporal space of the liquid. When used for classification, it is important that this transformation yields liquid states that are well separated across different classes. Thus an appropriate measure of the quality of the LSM is a criterion that incorporates the liquid means and variances in order to compute the separation of the data. For a classification task of 𝜔 𝑖 different classes we define three quantities.

(a) The between-class scatter matrix:

𝑆 𝑏 = � 𝑃 𝑖 𝑀 𝑖=1 (𝜇 𝑖 -𝜇 0 )(𝜇 𝑖 -𝜇 0 ) 𝛵 (1)
where

𝜇 0 = ∑ 𝑃 𝑖 𝜇 𝑖 𝛭 𝜄=1
is the global mean vector for all 𝑀 classes, 𝑃 𝑖 is the a priori probability of class 𝜔 𝑖 , and 𝜇 𝑖 is the mean of the liquid states that correspond to class 𝜔 𝑖 .

(b) The within-class scatter matrix:

𝑆 𝑤 = � 𝑃 𝑖 𝛴 𝜄 𝑀 𝑖=1 (2)
where 𝛴 𝑖 is the covariance matrix of the liquid states that correspond to class 𝜔 𝑖 , and 𝑃 𝑖 is as eq. 1.

(c) The covariance matrix 𝑆 𝑚 with respect to the 𝜇 0 global mean:

𝑆 𝑚 = 𝑆 𝑤 + 𝑆 𝑏 (3)
The trace of the 𝑆 𝑚 matrix adequately describes the sum of the variances of the features around the global mean, while the trace of the 𝑆 𝑤 matrix is the sum of the feature variances computed for all classes. Consequently, the separation between different classes can be determined using the following quantity:

𝐹𝐷𝑅 = 𝑡𝑟𝑎𝑐𝑒{𝑆 𝑤 -1 𝑆 𝑚 } (4)
which is the generalization of the Fisher Discriminant Ratio [FDR,[START_REF] Fisher | The Use of Multiple Measurements in Taxonomic Problems[END_REF] to more than two classes.

To apply the FDR measure on an LSM we sample the spatio-temporal patterns of the liquid's action potentials with respect to each input. These are used to construct the state matrix 𝑆:

𝑆 = � 𝑠 11 ⋯ 𝑠 1𝑗 ⋮ ⋱ ⋮ 𝑠 𝑖1 ⋯ 𝑠 𝑖𝑗 � (5)
where 𝒔 𝒊𝒋 is the filtered output of the j th spiking neuron in the liquid and i is the index of the time window in which the outputs of the neurons are sampled. For each task, we create 𝑛 matrices 𝑆, each corresponding to a different class. The FDR measure is then calculated based on eqs. 1-4.

Genetic Algorithm based Liquid Evolution

The separation of the liquid is positively correlated with the performance of the LSM [START_REF] Legenstein | Edge of chaos and prediction of computational performance for neural circuit models[END_REF]. Consequently, one can improve the classification performance of the trained readouts by designing a liquid that has a high separation measure. To accomplish this we minimize the FDR (eq. 4) of a liquid using genetic algorithms (GAs). GAs are a stochastic optimization method that can optimize complex functions by exploiting their parameter space [START_REF] Fogel | An introduction to simulated evolutionary optimization[END_REF]. Due to their ability to find good solutions in multi-modal and non-differentiable functions they have been extensively used to optimize the performance of neural networks. In these cases researchers encode properties such as the architecture, weights or neuronal models of the network and exploit them in order to minimize some objective optimization function (see [START_REF] Yao | Evolving artificial neural networks[END_REF] for a review). To evolve the LSM we utilize three types of properties: (i) the parameters of the neurons in the liquid, (ii) the architecture of the liquid, and (iii) the local properties of each architecture. The effect these parameters have on the liquid performance is discussed in the following. The first part of the chromosome encodes the firing threshold of all neurons in the liquid and the mean of the Gaussian noise added to the neurons' output on every step. These two parameters control the responsiveness and generalization properties of each neuron. In the first case, if the firing threshold of a neuron is low, then it will require to integrate more spikes before firing a post-synaptic potential, making the liquid less responsive to the perturbations of the low input signals. The second parameter controls the mean of the white noise added to each neuron, which affects the generalization properties of the training. The second part of the chromosome encodes three different architectures for the LSM. Each architecture specifies a different way for connecting the inputs to the liquid, and the inter-liquid connectivity (Fig. 1). In architecture 1, all input neurons are connected to all neurons within the liquid. Thus all task information is integrated in overlapping liquid locations. This is the original setup suggested by [START_REF] Maass | Real-time computing without stable states: A new framework for neural computation based on perturbations[END_REF]. In the second architecture the neurons that encode the input from different sources project to different locations in the liquid. In this case, the resulting LSM will produce a state vector whose entries correspond to particular input properties. The third architecture also incorporates the properties of the input, but discriminates it depending on whether they are temporally varying or static throughout the classification task. Temporally varying inputs project to different locations within the liquid, while the constant input signals are propagated to all neurons. The last part of the chromosome encodes the properties common to all architectures. These include the size of the liquid map and the locality of the connections in the liquid (i.e. the size of the neighborhood that each neuron is allowed to connect to). Because of the exponential descending output of the dynamic synapses of the LSM [START_REF] Natschlager | Computer models and analysis tools for neural microcircuits[END_REF], the last parameter affects the chaotic dynamics of the liquid, i.e. the period in which a certain input has an effect on the liquid state.

To evolve the chromosomes we use 3 different operators, mutation, crossover and selection. Mutation was implemented by adding a random number drawn from a Gaussian distribution, with zero mean and standard deviation that starts from 1 and decreases linearly until it reaches 0 in the final generation. To perform crossover, the GA selects (with probability 0.5) a bit from each parent chromosome in order to form a child. Selection was implemented using a roulette wheel function.

In the current section we evaluate the performance of the FDR measure and the GA optimization framework on a number of different classification tasks. For this reason, we focus on two issues: (i) the ability of the proposed measure to predict the quality of the liquid in an LSM, and (ii) whether the optimization framework can reduce the error of the readouts by minimizing the FDR.

Measure Evaluation

To evaluate the FDR measure we use three classification tasks, which incorporate different methods for encoding the input. This is important since diverse input encodings can have a different effect on the liquid dynamics. Population codes [START_REF] Rieke | Spikes: Exploring the neural code[END_REF] provide a consistent representation of the input by using distributed and partially overlapping neuron groups in order to encode the values of a variable. In contrast, rate codes [START_REF] Rieke | Spikes: Exploring the neural code[END_REF] produce a higher homogeneity when used as input because they employ the same neuron to represent different input values. Consequently in the three tasks discussed below, rate codes are used to encode the input in two cases, whereas population codes are used in the third case. The liquid in all the aforementioned classification tasks follows the initialization and topology settings discussed in [START_REF] Maass | Real-time computing without stable states: A new framework for neural computation based on perturbations[END_REF]. To evaluate the measure on different classification methods we employ four different readouts: the first readout is implemented with a multi-layer perceptron using the backpropagation rule to train the weights [START_REF] Rumelhart | Learning Internal Representations by Error Propagation[END_REF]. The second readout implements linear regression [START_REF] Duda | Pattern Classification[END_REF]. The third readout implements a classifier that uses least squares to find the regression coefficients [START_REF] Duda | Pattern Classification[END_REF], while the fourth the p-Delta rule on a parallel perceptron layer [START_REF] Auer | A learning rule for very simple universal approximators consisting of a single layer[END_REF]. For the first classification task, we compare the performance of the two most popular measures in the literature, namely the centroids and rank measures, against the FDR. For this reason we use an LSM with one linear regression readout to classify whether the rate of the input is above a particular value, in this case five Hertz. Input is encoded as a random Poisson rate and applied for 1000ms to a pool of Leaky Integrate and Fire [START_REF] Stein | Some models of neuronal variability[END_REF] spiking neurons. The liquid used for this purpose consists of a pool of 125 spiking neurons, arranged in a 3 dimensional grid. States are sampled every 10ms for 1 second and input to the 4 readout units. The FDR, Rank and Centroids measures were (a) calculated for the two state vectors 𝑆 1 and 𝑆 2 from eq. 5 that correspond to classes 𝜔 1 and 𝜔 2 , and (b) compared against the performance of a linear regression readout for 16 different simulations (Fig. 2). Fig. 2. Evaluation of the FDR, Rank and Centroids measures against their ability to predict the performance of the linear readout map of an LSM. In each subplot the x,y axes correspond to different configurations for an LSM. The color in each x,y entry corresponds to the value of the error (for subplot a) or the negative value of the measure (for subplots b,c,d).

A good measure should be positively correlated with the error of any trained readout that is used to extract information from the liquid. As Fig. 2 shows there is a clear correlation between the value of the error of a readout map and the value of the FDR measure (for both cases high values close to 1 are colored with red shades, while low values close to 0 are colored with blue). Furthermore, the results presented in Fig. 2 show that the proposed measure outperforms the Centroids measure and, at the same time, performs better than the Rank measure. By comparing Fig. 2a and Fig. 2b, it is evident that the FDR measure can predict with satisfying accuracy the performance of the linear regression readout and, therefore, the separation of the liquid in the LSM (readout error/FDR correlation was 0.86). Due to space limitations we do not provide the corresponding contour plots for the other three readouts, although we note that the results were similar to the ones presented in Fig. 2. The second task requires the LSM to classify two different motions of an object, based on the projection of its image on a 9x9 grid of receptive field neurons. The output of the retina field is encoded as a group of 81 neurons, each one corresponding to a different cell. These neurons fire random Poisson spikes of 30Hz when their corresponding cell in the retina field is occupied, and at a rate of 5Hz otherwise. This output is then projected to a liquid with 63 neurons, where we record the postsynaptic potentials of the neurons for 3 sec (3000ms). Information from the liquid response is used to classify whether the movement on the retina belongs to either one of the two behaviors for 9 different simulations. The error is calculated by subtracting the readout value from the actual behavior being performaned for each step of the simulaton, and normalized to 1. Due to space limitations, results are not presented for this task in the form of contour plots but rather as graphs of the errors of the readouts against the FDR (Fig. 3). As Fig. 3 illustrates, FDR follows quite closely the corresponding error in all cases. The correlation between the FDR measure and the readout error was in all cases above 0.8, indicating a close relationship between the two. For the third task we use an LSM that must classify the type of three different objects, a circle, a square and a hexagon. To encode the input we first sharpen each image using a Laplacian filter and consequently convolve it with 4 different Gabor filters with orientations π, respectively. The four convolved images from the input are projected into four neuronal grids of 25 neurons, where each neuron corresponds to a different location in the Gabor output. Information from the liquid response is classified by the above mentioned four readouts for 9 different simulations (Fig. 4). The results presented above for all three tasks, indicate that the FDR measure can describe the quality of the liquid over a broad range of tasks and input encodings.

Liquid optimization

In the current section we evaluate the extent to which the aforementioned GA framework can optimize the performance of an LSM. For this reason, we consider an additional classification task that requires the integration of temporal information in the liquid states. More specifically, we consider a binary classification task in which an LSM must classify whether the end point of a moving planar robotic arm is closer (or not) than a predefined distance to a given target location. The difficulty of the task lies in the fact that the time-varying control model of the arm must be combined with the static signal of the end point location and produce discrete liquid states, even in cases where the input dynamics are not so different. Input in this case consists of two different channels. The first encodes the spatial location of an end point position in (𝑥, 𝑦) space coordinates and the second the inverse kinematics of different arm trajectories. Input joint positions are generated by creating different trajectories using a two-link planar arm based on one start position (Fig. 5a), three speed profiles (Fig. 5b) and five random ending positions for the train and test sets (Fig. 5c,d). The training set consists of the trajectories between the initial position (Fig. 5a) and a random end position (Fig. 5c). The test set is generated using a different set of ending positions (Fig. 5d). To determine the trajectory between a starting and ending position, a random speed profile is chosen from the templates in Fig. 5b. The joint configurations of the robot across the pathway of a trajectory are obtained using an iterative solution to the inverse kinematics problem based on the pseudo-inverse of the robot's Jacobian. To encode the target position we use a population code with 10 neurons for each dimension (i.e. the 𝑥, 𝑦 coordinates). Thus for the two dimensional space 20 input neurons are used. The simulated robotic arm that is employed in the experiments consists of 2 joints, namely elbow and shoulder, which are also encoded using population codes. The classification task we consider requires the LSM to predict whether in the next location, the end point of the robot's arm will be closer than a predefined distance to the object in 𝑥, 𝑦 coordinates. To classify a given location correctly, the LSM must make a prediction on the speed of the arm at any given time. Hence, the liquid state must integrate information about the location of the robot's end-point effector position in previous time steps. To generate the different liquid states we have run 100 simulations for the train set kinematics and 30 simulations for the test set kinematics. To learn the classification task, these liquid states were input to 4 readout units, namely a feedforward neural network, a parallel perceptron layer, a linear regression and a linear classification readout. Even though the results reported here regard the optimal individual produced by the genetic algorithm, it should be noted that similar results were obtained for all chromosomes in the last generation. The GA was run for 18 generations in order to minimize the FDR criterion (Fig. 6). As the four rightmost plots in Fig. 6 show, while the genetic algorithm was used to minimize the FDR measure, it also reduced the error on all four readout maps. The subplots in Fig. 6 demonstrate how the GA was able to optimize the performance of the LSM by reducing the classification error of the readouts. The error was reduced from 0.3 to 0.1 in all four readouts maps (Fig. 6, right four subplots), simply by optimizing the FDR measure on the liquid (Fig. 6, left subplot) during the evolutionary process. The improvement in the liquid performance is also evident when we examine the output of the four classifiers in the optimal individual produced by the GA. As Fig. 7 shows, all the readouts are able to classify different movements with a high degree of accuracy. In the first example (plot a), the robot's arm never reaches the target location in a distance closer than required. In the second (plot b), it approximates the end location in the final 10 simulation steps. In both plots, we show the stimulus input to the liquid (top graph of plots a,b), the output of the four readouts (blue lines in bottom four graphs of plots a,b) and the target for each readout (red lines in bottom four graphs of plots a,b). Each graph is labeled with the corresponding readout map. The x-axis represents the simulation time in 100ms intervals for all graphs.

Discussion

In the current paper we presented a method for improving the computational capabilities of LSMs by focusing on the separation property of the liquid. This was measured using the Fisher's Discriminant Ratio, a measure that is maximized when the class means are far from each other, and the class variances are as small as possible. To evaluate the FDR measure against a broad range of classification tasks, we incorporated different types of neuron encodings for the input. As the results show, the FDR criterion accurately predicts the performance of the readouts without having any knowledge of the algorithm used to train them. Due to this fact, the GA, by minimizing the value of the FDR criterion, also improved the performance of all four classifiers. In contrast to other criterions, the FDR is a supervised measure, i.e. requires the class labels in order to compute the quantities in eqs. 1-3. Consequently, the evolutionary framework that was presented is also supervised. We consider this a benefit of the method, since it allows the design of liquid architectures that will be suited for the specific dynamics of the classification task. The proposed methodology will be used to develop a large-scale model of movement imitation and observational learning. In this context, the FDR measure can be employed to optimize the performance of the LSMs that underpin the agent's perception and motor control system. We also plan to investigate whether projecting the liquid states along the eigenvectors of the argument of the FDR measure, can further improve class separation.

Fig. 1 .

 1 Fig. 1. The three architectures used in the genetic algorithm, shown from a different perspective in order to highlight how the components are connected together. In each architecture the liquid component of the LSM is shaded with red, while the inputs with green.

Fig. 3 .

 3 Fig. 3. Graphs of the FDR measure (blue line) against the readout error (red line) of linear regression, backpropagation, linear classification and pDelta methods. The x-axis corresponds to the 9 different trials of the simulation, while the y-axis shows the output of the corresponding error and measure.

Fig. 4 .

 4 Fig. 4. Graphs of the FDR measure (blue line) against the readout error (red line) of linear regression, backpropagation, linear classification and pDelta methods. The x-axis corresponds to the 9 different trials of the simulation, while the y-axis shows the output of the corresponding error and measure.

Fig. 5 .

 5 Fig. 5. a. The initial position of the robot's arm. b. The three speed profiles used to generate random movements. c. The five different configurations of the arm of the robot for the five ending positions of the train set. d. The five different configurations of the arm for the test set.

Fig. 6 .

 6 Fig. 6. The results from the evolution of the LSM for 18 generations of the 100 individuals.The left plot shows the best and mean fitness (y-axis) of the population on each generation (xaxis). The right set of four plots shows the average error of the four readouts across the generations (y-axis) of the 100 individuals (x-axis) in the final population.

Fig. 7 .

 7 Fig. 7. The output of the four trained readouts for two sample movements (four bottom subplots) and the input projected to the liquid (top subplot).


Lazaros Iliadis
email: liliadis@fmenr.duth.gr2staxos@gmail.com3stavros.avramidis@ubc.ca

Stavros Tachos

Stavros Avramidis

Shawn Mansfield

Support Vector Machines versus Artificial Neural Networks for wood Dielectric Loss Factor estimation

Keywords: Wood Loss Factor, Gaussian Support Vector Machines, Polynomial Support Vector Machines, Artificial Neural Networks

This research effort aims in the estimation of Wood Loss Factor by employing Support Vector Machines. For this purpose experimental data for two different wood species were used. The estimation of the dielectric properties of wood was done by using various Kernel algorithms as a function of both ambient electro-thermal conditions applied during drying of wood and basic wood chemistry. Actually the best fit neural models that were developed in a previous effort of our research team were compared to the Kernels' approaches in order to determine the optimal ones.

Introduction

It has been shown in the literature that the knowledge of dielectric properties of wood can be used to determine its density and moisture content nondestructively. On the other hand this knowledge can lead to the detection of hand knots, spiral grain, and other defects [START_REF] Collet | Evaluation of wood characteristics: Internal scanning of the material by microwaves[END_REF]. From a theoretical perspective a better understanding of the molecular structure of wood and wood-water interactions can be obtained [START_REF] Kabir | Temperature dependence of the dielectric properties of rubber wood[END_REF]. Various approaches towards this direction have been reported in the literature, namely the radio frequency vacuum drying [START_REF] Koumoutsakos | a. Radio frequency vacuum drying Part I. Theoretical model[END_REF] [START_REF] Koumoutsakos | Radio frequency vacuum drying Part II[END_REF] [START_REF] Koumoutsakos | Radio frequency vacuum drying Part III. Two dimensional Model, Optimization and Validation[END_REF] and the high frequency electric field heating, such as the veneer and finger-joint gluing and parallam manufacturing [19]. Thus, estimation of the dielectric constant (ε'), the loss tangent (tanδ) and the loss factor (ε") plays a catalytic role in the process of design, control, optimization and simulation.

Recently Koumoutsakos [START_REF] Koumoutsakos | a. Radio frequency vacuum drying Part I. Theoretical model[END_REF] has shown that ε" is proportional to the thermal energy transferred to the wood, during radio frequency vacuum drying (RFV) process. RFV heating is a volumetric method where thermal energy is produced simultaneously through a pile of lumber which is placed inside an electromagnetic field. In this process the moisture starts its transportation from the center to the surface as soon as the wood is exposed to the field [START_REF] Avramidis | Radio Frequency Vacuum Drying of wood[END_REF]. In fact, the electric power converted to thermal is given by the following formula 1.

PD = 5.56x10 -11 E 2 f ε" (1)
It should be clarified that PD is the power density in W/m 3 , E is the field strength in V/m and f is the frequency measured in Hz Research efforts in the past have shown that there exists a strong relationship between the dielectric properties and wood attributes studied plus field frequency applied [18][19].

It is a fact that despite all conducted research so far, there is very little respective knowledge regarding the effect of wood chemical composition on the determination of ε". Actually, Norimoto [START_REF] Norimoto | Dielectric properties of wood[END_REF] [START_REF] Norimoto | Anisotropy of dielectric constants in coniferous wood[END_REF] have investigated the dielectric properties of some wood chemical constituents as a function of frequency and temperature, but no attempt was made towards the correlation of their percent content in the cell-wall composition to the gross wood ε" values.

Another point that motivated this study is the fact that the data that have been used so far to construct such models have originated from assorted species, under different thermo-physical conditions and variable frequencies and thus, they are not suitable to be employed in drying modeling [START_REF] Avramidis | Wood Dielectric Loss Factor Prediction with Artificial Neural Networks Wood Science and Technology[END_REF].

This research effort aims in the development of Support Vector Machines models (Gaussian, Fuzzy weighted and Polynomial) capable of determining reliably the value of ε" based on ambient electro-thermal conditions and also on basic wood chemistry. On the other hand a comparative study is performed with a previous study of our research team that used Artificial Neural Networks (ANNs) for the same purpose [START_REF] Avramidis | Wood Dielectric Loss Factor Prediction with Artificial Neural Networks Wood Science and Technology[END_REF]. This comparison aims in the determination of the optimal approaches that can be used to offer reliable and most of all cost and time effective approximation.

Materials and methods

Obtaining the experimental data

This study employees ε" values and macro-physical data that were already published by Zhou and Avramidis [22]. According to Zhou and Avramidis all-sapwood and allheartwood western hemlock [Tsuga heterophylla (Raf.) Sarg.], and all-heartwood western red cedar [Thuja plicata Donn] specimens were evaluated in the radial direction (thickness) and at various moisture contents and temperature levels were exposed to two levels of electric field voltage. After the exposure was terminated the ε" values were calculated indirectly by performing heating studies at a 13.56 MHz fixed frequency with a laboratory size RFV dryer [22] [START_REF] Avramidis | Wood Dielectric Loss Factor Prediction with Artificial Neural Networks Wood Science and Technology[END_REF]. Also the above described wood species were analyzed by Zhou and Avramidis [22] regarding their chemical composition as follows: air dried wood samples were ground in a Wiley Mill to pass a 40-mesh screen and extracted in a Soxhlet apparatus with acetone for 12 hours, and extractives determined gravimetrically [START_REF] Avramidis | Wood Dielectric Loss Factor Prediction with Artificial Neural Networks Wood Science and Technology[END_REF]. According to Zhou and Avramidis [22], the content of Lignin was determined using a modified Klason approach derived from the TAPPI standard method T222 om-98. For more details refer to [22].

5-Fold-Cross Validation

Cross-validation is a methods for estimating generalization error based on "resampling" [START_REF] Weiss | On the loss factor of two B.C. softwoods[END_REF] [7] [8] [START_REF] Shao | The Jackknife and Bootstrap[END_REF]. The resulting estimates of generalization error are often used for choosing among various models. In k-fold cross-validation, the data are divided into k subsets of (approximately) equal size. The model is trained k times, each time leaving out one of the subsets from training, but using only the omitted subset to compute whatever error criterion interests you. If k equals the sample size, this is called "leave-one-out" cross-validation. The division of the data set was done by the use of MATLAB's crossvalind function, which is included in the Bioinformatics Toolbox.

ε-SV Regression (ε-SVR)

Support Vector Machines (SVM) are keeping the training error foxed while at the same time they are minimizing the confidence interval [START_REF] Kecman | Learning and Soft Computing[END_REF]. Their primary target was pattern recognition [START_REF] Boser | A training algorithm for optimal margin classifiers In Fifth Annual Workshop on Computational Learning Theory[END_REF]. Now they are used for both regression and classification. Vapnik [START_REF] Vapnik | Support Vector method for function approximation regression estimation and signal processing[END_REF] [10] introduced the following loss function that ignores errors less than a predefined value

0         max 0 , y f x y f x         (2)
It is a fact that the ε-SVR algorithm offers in many cases the optimal function of the

form:     , , , N f x k w x b w x R b R       (3)
. The whole idea is based on the determination of the function with the minimum testing error. The problem in this case is that the minimization of the above function is not possible due to the fact that the probability distribution P is unknown. Consequently, the actual solution can be reached by minimizing the following normalized risk function 4. [START_REF] Davy | An Introduction to Support Vector Machines and other kernel algorithms[END_REF] and it ignores errors less than  .

  2 1 2 SVR emp w C R f    (4)   emp Rf  is the function of empirical risk     1 1 ,, p emp i i i R f L f x y p    (5) whereas the loss function is   ( , , ) i i i i L f x y y f x   (
Also it is clear that 2 w is related to the complexity of the model, whereas SVR C is a constant value determining the point that relates

  emp Rf  to 2 w [6].
Obviously, the objective which is the Minimization of the above risk function 4 can be phased as a specific optimization problem which phases the constraints described below: The aim is the minimization of the following function:

    2 1 1 1 ,, 2
p SVR i i i w w C p              (7)
The constraints in this case are: ii    (10). Obviously (and this can be seen in the following figure 1)

    , i i i k w x b y      (8) and     i i i y k w x b        (
i   , i
 are the distances of the training data set points from the e-zone. All of the errors that are smaller than  are ignored. As it is clearly shown, i  stands for the distance of a point that can be found above the e-tube zone and i   stands for the distance of a point that is located below the e-tube zone. It is well known that in the case of the application of Lagrange multipliers the problem would correspond to a double optimization one as follows: Maximize [START_REF] Koumoutsakos | Radio frequency vacuum drying Part III. Two dimensional Model, Optimization and Validation[END_REF].

           1 1 , 1 1 , , 2
p p p i i i i i i i j j i j i i i j W a a a a a a y a a a a k x x                     (11) subject to   1 0 p i i i aa      (12) and   0, SVR i C a p         (

Figure 1. ε-SVR for Regression

Based on Vapnik and Kecman [START_REF] Vapnik | Support Vector method for function approximation regression estimation and signal processing[END_REF][10], the optimization problem can be solved by the use of the following linear extension of the kernel functions ( 14) and ( 15)

  1 p i i i i w a a x      (14)       1 , , , p i i i i f x a a a a k x x b        ( 15
)
where b is estimated by the following function 16

      , i i i k j j j i j b average sign a a y a a k x x                 (16)
The training set vectors

    2 2 1 , , exp 2
p i ii i RBF xx f x a a a a b               (17)
The involved parameters   ,, RBF    are very significant and they determine to a great level of extend the good fit of the ε-SVR. It should be specified that RBF  is the RBF kernel's standard deviation,  has a constant value that determines the point where the empiric error is related to complexity and finally the parameter  is the width of the ε-zone. In statistics, when we consider the Gaussian probability density function σ RBF is called the standard deviation, and the square of it, the variance.

Fuzzy weighted Support Vector Regression

The fuzzy weighted SVR with a fuzzy partition first employs the fuzzy c-mean clustering algorithm to split training data into several training subsets. Then, the localregression models (LRMs) are independently obtained by the SVR approach for each training subset. Finally, those LRMs are combined by a fuzzy weighted mechanism to form the output. Experimental results show that the proposed approach needs less computational time than the local SVR approach and can have more accurate results than the local/global SVR approaches does [START_REF] Chuang | Fuzzy Weighted Support Vector Regression with a Fuzzy Partition Systems,Man, and cybernetics[END_REF].

Application

In this section of the paper, the application results of the Global SVR approach for all existing data vectors will be presented thoroughly. The regression was done for different kernels, namely the Gaussian (RBF) the Fuzzy weight SVM and the Polynomial one. In order to perform Global SVR a global regression model was developed.

The regression was performed by the use of the LIBSVM v2.9 (http://www.csie.ntu.edu.tw/~cjlin/libsvm/) [START_REF] Chang | LIBSVM A Library for Support Vector Machines[END_REF] which is encoded in C++ and offers a Matlab Interface. In this specific application the RBF-Kernel was applied.

Performing trial and error for the Gaussian Kernel

The initial data set comprised of one hundred forty-four (144) sets of experimental data, where the ε" was measured under different temperature and moisture conditions and for various types of chemical composition.

During Training and Testing, several hundreds of experiments were performed in order to determine the optimal parameters' values. More specifically in the case of the Gaussian (RBF) Kernel, 400 experimentation cases (applications of the 5-fold cross validation) were constructed with the integer value of σ RBF ranging from 1 to 20 and with the value of the parameter γ ranging from 1 to 14. In the same case the e parameter took real values from the following set {0.003, 0.004, 0.005, 0.006, 0.007, 0.008, 0.009, 0.1, 0.12, 0.14, 0.16, 0.18,0.20, 0.22} This means that totally 400*14=5600 trials were performed during the 5-fold cross validation process.

The Root Mean Square Errors (RMSE) produced in each Training or Testing cycle were stored and then they were averaged each time 5-fold cross validation was employed. In each trial the initial data set was divided into five subsets (4 training and one testing) where five training and testing cycles were performed (each time changing the testing set in a round robin manner). Totally 5600 Average Root Mean Square Errors were obtained (ARMSE) one for each 5-fold cross validation cycle and for each unique combination of the parameters, together with the corresponding % Mean Average Percent Error (MAPE).

In the case of the Gaussian Kernel the optimal SVM model was produced for the σ RBF = 2 and for the value of the parameter γ=12 in the case with e=0.04 as it is shown in the following table 1. The following figures 2 and 3 present the evolution of the Mean RMSE and of the MAPE respectively, according to the values of the e parameter for the Gaussian Kernel. It is clearly shown that the best model should have an e value equal to 0.04. For the case of the Polynomial Kernel, the σ RBF again took integer values from 1 to 20 and the parameter γ from 1 to 14. However the e parameter took real values from the following set {0.005, 0.1, 0.15, 0.20, 0.25}. The performance of the best fit Polynomial Kernel is shown in the following table 2 and it was achieved for σ RBF= 3 γ=18 and e=0.1 For the Fuzzy weight Support Vector Machine the σ RBF again took integer values from 1 to 20 and the parameter γ from 1 to 14. Also the e parameter took real values from the following set {0.005, 0.1, 0.15, 0.20, 0.25}. The performance of the best fit Fuzzy weight SVM is shown in the following table 3 and it was achieved for σ RBF= 3 γ=16 and e=0.05

Performing Gaussian SVM with 5-fold cross validation

Discussion -Comparative analysis

The results of the e-SVM regression were compared to the results of the Artificial Neural Network that was developed by our research team in a previous study [START_REF] Avramidis | Wood Dielectric Loss Factor Prediction with Artificial Neural Networks Wood Science and Technology[END_REF] and also the results of the three SVM kernels were compared to each other. In [START_REF] Avramidis | Wood Dielectric Loss Factor Prediction with Artificial Neural Networks Wood Science and Technology[END_REF] the optimal ANN that estimated wood loss factor for the same input variables was found to be a multi layer Back propagation one, with Sigmoid transfer function. It was found to have eleven [START_REF] Koumoutsakos | a. Radio frequency vacuum drying Part I. Theoretical model[END_REF] processing elements in the input layer, nine (9) neurons in the hidden layer and of course one neuron in the output layer [START_REF] Avramidis | Wood Dielectric Loss Factor Prediction with Artificial Neural Networks Wood Science and Technology[END_REF]. The above table 4 presents the performance of the best fit ANN. It is clearly shown that the Gaussian Kernel approach has by far much lower RMS error than the Polynomial and the Fuzzy weight SVM and also than the best fit ANN. Thus it can be considered the most reliable modeling approach.

9 )

 9 and of course , 0

  presents the Radial Basis function (RBF) kernel

Figure 2 .Figure 3 .

 23 Figure 2. Evolution of the Mean RMSE for different values of the parameter

Table 1 : Performance of the Optimal Gaussian Kernel Optimal Gaussian Kernel SVM (5-fold cross validation)

 1

				Mean	Mean
	Value of	Value	Value	RMSE	RMSE	MAPE	MAPE
	σ RBF	of γ	of e	Training	Testing	Training	Testing
	2	12	0.04	0.00009664	0.0001501	2.4824596	3.7328347

Table 2 : Performance of the Optimal Polynomial Kernel Optimal Polynomial Kernel SVM (5-fold cross validation)

 2

				Mean	Mean
	Value of	Value	Value	RMSE	RMSE	MAPE	MAPE
	σ RBF	of γ	of e	Training	Testing	Training	Testing
	3	18	0.1	0.000485	0.000549	11.74819	14.01062

Table 3 : Performance of the Optimal Fuzzy weight SVM Optimal Fuzzy weight SVM (5-fold cross validation)

 3

				Mean	Mean
	Value of	Value	Value	RMSE	RMSE	MAPE	MAPE
	σ RBF	of γ	of e	Training	Testing	Training	Testing
	3	16	0.05	0.000488	0.000542	12.13017	14.12399

Actual versus estimated values of Loss Factor for Fuzzy Weighted SVMs

Table 4 : Performance of the optimal ANN for the estimation of loss factor in wood

 4

			Optimal
	Optimal	Optimal	ANN	Optimal
	ANN	ANN R 2	RMSE	ANN R 2
	RMSE	Training	Testing	Testing
	Training
	0.017	R 2 =0.9989 0.0382	R 2 =0.9945


Primož Potočnik
email: primoz.potocnik@fs.uni-lj.si

Tomaž Berlec
email: tomaz.berlec@fs.uni-lj.si

Marko Starbek
email: marko.starbek@fs.uni-lj.si

Edvard Govekar
email: edvard.govekar@fs.uni-lj.si

SOM-Based Clustering and Optimization of Production

Keywords: production optimization, clustering, SOM neural network, k-means clustering, hierarchical clustering

des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Modern production strives toward the optimization of production costs by methods such as lean production, [START_REF] Hines | Going Lean, Lean Enterprise Research Centre[END_REF] which involves defining the internal value system and consequently analyzing and optimizing the internal production flow. Related methods known as "group technology" (GT) [START_REF] Kusiak | The generalized group technology concept[END_REF] have been proposed to implement cellular manufacturing systems (CMS) that group machines into production cells. A number of methods have been proposed for cell formation. [START_REF] Crama | Models for machine-part grouping in cellular manufacturing[END_REF][START_REF] Shanker | Models and solution methodologies for the generalized grouping problem in cellular manufacturing[END_REF][START_REF] Adenso-Díaza | Machine cell formation in generalized group technology[END_REF][START_REF] Fan | Decision support for proposal grouping: a hybrid approach using knowledge rules and genetic algorithms[END_REF] A unified approach that combines assigning parts to individual machines and forming machines into cells is proposed in [START_REF] Foulds | A network flow model of group technology[END_REF]. A clustering model is introduced in [START_REF] Andrés | Group technology in a hybrid flowshop environment: A case study[END_REF], exploiting similarities between products. The ART1 neural network-based cell formation in GT has been proposed by [START_REF] Yang | Machine-part cell formation in group technology using a modified ART1 method[END_REF] and group technology based clustering is applied by [START_REF] Morača | Application of Group Technology in Complex Cluster Type Organizational Systems[END_REF].

In this paper, we develop a framework for production optimization for small and medium enterprises (SMEs) that are characterized by individual and small batch production with many different products in their production range. This paper's findings are based on data obtained from a Slovenian manufacturing company. We propose a clustering approach for the segmentation of various groups of products that can be organized into small production cells. Based on the representatives of each production cell, an efficient floor plan can be designed that will support lean production objectives.

Data preprocessing

The data about required operations and features of the products need some preprocessing before being applied to various clustering algorithms. The following rules were applied to prepare the data: 1. Operations are encoded with a single value ('operation is not required' = -1, 'operation is required' = 1). 2. Materials are encoded with a single descriptor ('material not used' = -1, 'material used' = 1) for each possible material (Al/CuZn, SL, Fe).

3. Form is encoded with a single descriptor (-1,1) for each form applied (cast/forged, rod, platinum, band, profile/tube). 4. Shape is expressed as 'simple' = -1, 'combined' = 1. 5. Dimensional accuracy is originally expressed as ['0.1', '0.01', '0.001'] and we encode this property as '0. Finally, constant operations (packaging and progressive deep drawing) were eliminated. Without prior knowledge about the relative importance of various operations, we assumed equal importance for all operations. This assumption was encoded into data by using the same scaling in interval [-1,1] for all the attributes. The preprocessed data comprises 58 attributes: 37 operations, 3 materials, 5 forms, 1 shape, 1 dimensional accuracy, 1 appearance, 1 request for examination, 3 protections, and 6 for real valued categories.

Clustering algorithms

Several clustering algorithms were applied to the task of products clustering, as follows: hierarchical clustering, k-means clustering, and SOM neural network. An important ability of these algorithms is to also accept real-valued attributes, which is not possible with structures such as the ART1 neural network.

Hierarchical clustering

Hierarchical clustering algorithms produce a nested series of partitions based on a criterion for merging or splitting clusters based on similarity [START_REF] Jain | Data clustering: a review[END_REF]. In this paper, we apply agglomerative hierarchical clustering, with a Euclidean distance metric and several linkage methods: single linkage, complete linkage, average linkage, and Ward linkage.

K-means

The k-means is the simplest and most commonly used algorithm employing a squared error criterion [START_REF] Jain | Data clustering: a review[END_REF]. It starts with a random initial partition and keeps reassigning the patterns to clusters based on the similarity between the pattern and the cluster centres until a convergence criterion is met.

SOM

A self-organizing map (SOM) was proposed by Kohonen [START_REF] Kohonen | Self-Organizing Maps, Second Edition[END_REF]. In this paper, we examine a SOM algorithm with distances defined by the Euclidean distance metric, and two variants of two-dimensional grid (rectangular and hexagonal). We apply a two-dimensional grid with six clusters, arranged by two rows each with three elements. Such a topology seems to be well suited to production floor planning.

Evaluation of clustering results

Clustering validity measures fall broadly into three classes [START_REF] Brun | Model-based evaluation of clustering validation measures[END_REF]: a) internal validation (based on properties of the resulting clusters), b) relative validation (running the algorithm with different parameters), c) external validation (comparison with a given partition of the data). In our case study, there is no possibility of evaluating the correct clustering based on external validation measures; therefore, we have to rely in internal and relative validation that makes clustering essentially a subjective visualization tool. Following the recommendation proposed in [START_REF] Brun | Model-based evaluation of clustering validation measures[END_REF], we evaluate the clustering results based on silhouette values that seem to be a good internal validation measure and also provide good graphical representation of clustering quality. The silhouettes validation technique [START_REF] Rousseeuw | Silhouettes: a Graphical Aid to the Interpretation and Validation of Cluster Analysis[END_REF] calculates the silhouette width for each sample, the average silhouette width for each cluster and the overall average silhouette width for a total data set. In our study, we evaluate the clustering quality based on the average silhouette width for a total data set.

Clustering results

In this section, we present and compare various clustering results, obtained by the proposed clustering algorithms. Results are presented for hierarchical clustering, k-means, and SOM clustering. Within each clustering method, various parameters are optimized in order to obtain the best results.

From the perspective of production planning, the company would prefer a small number of clusters, i.e. condensed production cells sharing the necessary tools and operation within confined space. The number of clusters should be around K≈5 but this is only a recommendation and not a strictly defined condition. Consequently, in our study we fix the number of clusters to be K = 6, which supports the application of a two-dimensional clustering architecture (3 × 2).

Hierarchical clustering

Hierarchical clustering results based on Euclidean distance and several linkage methods are shown in Figure 1. Single linkage obviously yields the worst result, with an average silhouette S = 0.1937. Average and complete linkage give better clustering results, while the best result, S = 0.4698, is obtained by applying the Ward distance.

K-means

Due to the sensitivity of the k-means algorithm to become trapped in a local minimum, the algorithm was restarted 100 times from various random initial positions. This effectively converged into a unique solution, as presented in Figure 2.

The obtained average silhouette value amounts to S = 0.4742, which slightly exceeds the result obtained by hierarchical clustering.

SOM

Two variants of SOM grid organization are explored: rectangular and hexagonal grid.

For the distance metrics and for the neighbourhood distance functions, the Euclidean distance metric is applied. Figure 3 displays the clustering results for both topologies (rectangular and hexagonal). The best result is obtained by the hexagonal topology and yields an average silhouette of S = 0.4743. Figure 4 shows a hexagonal SOM topology with class labels and number of samples in each class.

Comparison of clustering results

The random initialization of the k-means and SOM clustering algorithm causes various final clustering arrangements in which class labels depend on initial conditions. In order to compare the silhouette plots of various clustering algorithms, the class labels should be aligned properly. We applied the method of class matching so that outcomes of various clustering methods were rearranged in a way that guarantees the highest cluster matching. The best results of hierarchical clustering, kmeans, and SOM clustering are shown in Figure 5 and presented in Table 2. Table 2 also presents the similarity measure between various clustering approaches. The methods are compared to the SOM result with respect to the number of equally classified samples. After rearrangement, all three methods exhibit very similar results, with k-means and SOM being almost identical and hierarchical clustering diverging below 5%. The overall best result is considered to be obtained by SOM. Its average silhouette, S = 0.4743, is the highest score obtained in this study. The result is almost the same as k-means and only slightly exceeds hierarchical clustering, but there are two more advantages to support the selection of SOM clustering: a) SOM result has no negative silhouettes, which means there are no products that are classified in a wrong cluster, b) SOM topology shown in Figure 4 can be directly interpreted as a production floor plan. SOM clusters maintain neighbourhood properties, which can be very helpful when designing a production floor plan. Operation clusters that are close to each other will probably share more equal operations and thus more material exchange than the clusters that are far apart. Therefore, we conclude that the SOM clustering method seems to be the most suitable approach for the task considered in this paper. The homogeneity of results obtained by various clustering methods only further supports the assumption that the obtained clustering result is meaningful and therefore applicable to production planning.

Application of clustering results

After obtaining a meaningful clustering result, the next step is to apply this result to the production environment. In this section, we propose an interpretation of clustering results that may yield an applicable industrial solution.

According to the SOM clustering result, six operation cells are arranged in a twodimensional hexagonal grid. This architecture (shown in Figure 4) can already be interpreted as an initial production floor layout. The next question is about which operations should be contained in the arranged production cells. This leads to the interpretation problem of clustering results as described below.

Interpretation of clustering results

The interpretation problem can be formulated as how to define a mapping from obtained SOM clusters into the real world production cells:

SOM clusters Production cells

As SOM clusters are represented by prototypes, an initial estimation could be to directly translate SOM prototypes into production cells, but this turns out to be inappropriate. SOM prototypes have continuous values in the interval [-1,1], which is acceptable for the product features (weight, volume, shape, etc.) but not for the operations that should be either included or not included in the production cell. Therefore, some kind of discretization of SOM prototypes should be performed for logical descriptors (such as operations, materials, form, shape, etc.).

We propose a percentile measure for the interpretation of SOM clusters into the production cells. For each operation in a particular cluster, we can provide a percentage of samples in this cluster, which should contain this operation (or property) in order for this operation (or property) to be included into the production cluster. Various percentile margins can be defined, such as: p = {50th, 75th, 100th}. For each percentile margin, the particular operation should be included into the production floor planning cell if at least (100-p)% of samples in a particular SOM cluster require this operation. The 100th percentile should be interpreted as a limit value: if at least one sample requires an operation, it should be included in the production cell.

Figure 6 presents SOM interpretation results based on various percentile margins p = {50th, 75th, 100th}. The result for p = 50th seems to be under-populated as there are a significant number of operations missing in a complete production scheme because they are simply not frequently required by the production process. On the other side, the result for p = 100th is probably over-populated. An optimal arrangement can be expected somewhere between the 50th and 100th percentile margins; therefore, a 75th percentile margin can be taken as a guideline to successfully interpret SOM clusters into the production cells.

Finalization of production clusters

If SOM interpretation based on the 75th percentile is applied to cell arrangement, some operations are still not assigned to any production cell. This is a consequence of their rare application, but in practice the production process can not be completed unless all the required operations are available. Therefore, as a final stage, we propose assigning each missing operation into the cell that has the highest requirement for this operation among all cells. The final result is shown in Figure 7, where filled missing discrete operations are displayed in blue. The following two properties support the assumption that the resulting production cell arrangement is suitable for production planning: a) Organization of production cells supports production of similar products in closed production units, which optimizes material and work flow, and reduces production costs. Clustering evaluation based on silhouette widths confirms good clustering quality, which means that the proposed clustering is meaningful. b) Production cells are arranged according to SOM topology. This guarantees that neighbourhood properties of clusters are maintained and consequently, it can be expected that this will lead to the minimization of material and work piece exchange between production cells that are not close to each other. The results reported can be considered as a recommendation to the production planning managers. We hope the proposed results will be a useful guidance to production planning in the company.

Fig. 1 .

 1 Fig. 1. Hierarchical clustering results compared by various linkage measures. Ward linkage results in best average silhouette value S = 0.4698.

Fig. 2 .

 2 Fig. 2. K-means clustering, average silhouette value S = 0.4742.

Fig. 3 .Fig. 4 .

 34 Fig. 3. SOM clustering results for rectangular and hexagonal 2-dimensional topology

Fig. 5 .

 5 Fig. 5. Comparison of clustering results (hierarchical clustering, k-means, SOM).

Fig. 6 .Fig. 7 . 6 Conclusions

 676 Fig. 6. Interpretation of SOM clusters into production cells based on percentile margins

Table 1 .

 1 List of operations applied during the production.

	No	Operation	No	Operation	No	Operation
	1 Band cutting	14 Thermal cutting	27 Cutting on shears
	2 Service-casting	15 Service-forging	28 Deploy profile
	3 Broaching	16 Drilling	29 Single-spindle thread.
	4 Multi-spindle thread.	17 CNC turning	30 3-axis CNC machin.
	5 Brushing	18 Service-blasting	31 Service-deburring
	6 Honing	19 Powder Coating	32 Dip galvanizing
	7 Galvanic coating	20 Artificial aging	33 Carbonitrirating
	8 Service-hardening	21 Washing	34 Precision polishing
	9 Before assembly	22 Assembly 1	35 Hand deburring
	10 Viewing area	23 Testing	36 Packaging
	11 Chamfering machine	24 Tumble deburring	37 Compression-cutting
	12 Remodeling	25 Compression	38 Progressive compress.
	13 Progr. deep drawing	26 Deep Draw-transfer	39 Welding

Table 2 .

 2 Comparison of clustering results by average silhouette values.

	Method	Average silhouette S	Similarity to SOM
	Hierarchical clustering	0.4698	95.2%
	K-means	0.4742	99.6%
	SOM	0.4743	100%