Petr Olivka
email: petr.olivka@vsb.cz

Michal Krumnikl
email: michal.krumnikl@vsb.cz

Pavel Moravec
email: pavel.moravec@vsb.cz

David Seidl
email: david.seidl@vsb.cz

Process of Point Clouds Merging for Mapping of a Robot's Working Environment

Keywords: Lidar, Laser Range Finder, Point Clouds, Mapping, Robot

The lidar is nowadays increasingly used in many robotic applications. Nevertheless the 3D lidars are still very expensive and their use on small robots is not economical. This article briefly introduces a construction of cheap 3D lidar for indoor usage based on the 2D laser range finder. Subsequently, this article introduces process of merging acquired point clouds. The every pair of neighboring point clouds is oriented in space to fit together in the best possible way. The result of this process is a 3D map of the robot working environment. This map can be segmented and further used for a navigation. In this way, it is also possible to map inaccessible and dangerous areas.

Introduction

Nowadays there are many types of laser devices for a distance measurement. Probably the best known application of this type of devices is the use of lidar in the autonomous car driving systems. For example, Google cars use lidar Velodyne HDL-64E which can be seen in Fig. 1.

It is not possible to use this type of lidar on small robots in industrial applications. There are many reasons for not using Velodyne lidar. At first, it is very expensive. The Velodyne HDL-64E price is around $8000. For small robotic applications such price is absolutely unacceptable. The next reason is the narrow measuring range. Despite this lidar measures 360 degree around, the vertical field of view is only 28 degrees [START_REF] Lidar | HDL-64E[END_REF]. The visualization of the raw data from HDL-64E is in Fig. 2. The small viewing angle is suitable in traffic environment, but it is insufficient for indoor application like recognition of nearby objects and obstacles. But what is for internal use completely inappropriate is the blind space above the lidar. For the indoor robot navigation it is very important to detect space above the robot.

In contrast to this high-end devices there are several cheap devices available on the market suitable for home or office use. The first example is Kinect as Fig. 1. Velodyne lidar HDL-64E [START_REF] Lidar | HDL-64E[END_REF] Fig. 2. Raw data from HDL-64E [START_REF] Lidar | HDL-64E[END_REF] depicted in Fig. 3 and a similar device Xtion PRO can be seen in Fig. 4. The disadvantage of these devices is that they are not designed for industrial applications and they have very limited measurement field of view in vertical and horizontal axes. These devices are also not prepared for measurement in 360 degrees range.

Because 3D lidars are nowadays still very expensive and cheap customer devices are not suitable for industrial usage, it is necessary to focus to a different device design for a reliable measurement of distances in 3D space. There is a large variety of 2D Laser Range Finders (LRF) available on the market. These devices are produced for measurement of short or long distances with different measurement precision, they are equipped with different interfaces and they are designed for indoor or outdoor use. These 2D devices can be easily used for the measurement in 3D. It is only necessary to add the third dimension. The probably first design of 3D LRF device was introduced in 2009 and it can be seen in Fig. 5. In this figure there is 2D LRF Hokuyo URG-04LX mounted on a servo in order to achieve the ability to measure in 3D.

This design has two disadvantages. The servo has a problem with precision of positioning. It is impossible to achieve the precise angular position and the repeatability of the positioning is also problematic. The second problem of this design lies in the mechanical features of the servo. The design of the servo does Fig. 6. Zebedee design of 3D LRF with 2D LRF mounted on spring [START_REF] Bosse | Zebedee: Design of a spring-mounted 3-D range sensor with application to mobile mapping[END_REF] not handle axial and radial forces correctly. The use of this design on the robot will let to a quickly damage of the servo. The inertial forces of mounted 2D LRF will very quickly damage gears and their mounting.

Alternative design of 3D LRF was introduced by SCIRO [START_REF] Bosse | Zebedee: Design of a spring-mounted 3-D range sensor with application to mobile mapping[END_REF]. In this design the 2D LRF is mounted on a spring and the third dimension is added by a swinging. The overall design is shown in Fig. 6. Such design has also some disadvantages. It is necessary to manually swing with the mounted LRF and maintain very precise synchronization between the gyroscope mounted under the LRF and the LRF scanning.

3D LRF Desing

Based on the introduced constructions and their disadvantages we have developed a new design of the 3D LRF. The diagram of this design is in Fig. 7. The 2D LRF Hokuyo URG-04LX is mounted on a step motor by a carrier. This construction guarantees the sufficient strength of the design. All inertial forces are captured in ball bearings of the step motor and the step motor also guarantees the precise positioning of the mounted LRF. In the diagram we can see that the rotation axis of the step motor is identical with the internal axis x i of the LRF. It simplifies the future computation. The current position of the step motor is given by angle α r and the position and measurement range of laser beam is defined by angle α s . The mounting position of LRF on the step motor is marked in the diagram as the angle α m . Now the principle of 3D LRF measurement will be described. The 2D LRF measures distances of obstacles by a laser beam. The 2D measurement is performed in the plane x i y i . The current position of the laser beam is marked in the diagram as an angle α s . The resolution of used 2D LRF is 360 • /1024. The angular range of the measurement is from -120 • to 120 • and in this range 628 points are measured.

The carrier with 2D LRF is mounted on the step motor and is rotated around axis x i . The rotation of the whole LRF is performed in the range of a half turn (180 • ). This rotating movement guarantees the whole coverage of the environment. The current position of the 2D LRF is marked as α r . The resolution of the step motor is 360 • /800.

The overall design of the 3D LRF is depicted in Fig. 8. The figure shows a complete device with the control computer and the microstepping unit [START_REF] Krumnikl | PWM Nonlinearity Reduction in Microstepping Unit Firmware[END_REF]. For the following experiment the 3D LRF was mounted on the small robotic undercarriage Roomba.

The 2D LRF measures distances in all positions and measured distances will be marked as l r,s . This measured distances can be transformed to a position vector l r,s = [l r,s , 0, 0, 1] T . All measured distances can be then transformed to point cloud, where every point will be represented by its position vector p r,s . The computation is performed by following formula:

p r,s = R z (α r ) • R y (α m ) • R z (α s ) • l r,s , (1)
where R y and R z are transformation matrices for rotations.

During the single measurement the designed 3D LRF captures point clouds composed of 628 × 401 points. It is in total 252828 points.

Practically every manufactured LRF Hokuyo URG-04LX unit has some manufacturing inaccuracies. For the precise measurement it is necessary to calibrate every LRF. The process of calibration is described in details in [START_REF] Olivka | Calibration of Short Range 2D Laser Range Finder for 3D SLAM Usage[END_REF]. Acquired results are very inaccurate without this calibration.

The set of point clouds were acquired for the following experiment. Inside the building 26 meters long corridor was selected with variety of surfaces and irregular shapes. This corridor were measured by 3D LRF four times in both directions and measurements were performed approximately at a distance of one meter. Thus 25 point clouds were acquired in every measurement. Overall 100 point clouds were acquired. These points clouds will be marked as a set of point clouds M 1 to M 100 . Point clouds in range from M 1 to M 25 are from the first measurement, point clouds in range from M 26 to M 50 are from the second measurement, etc. An example of single point cloud, concretely M 1 , is visible in Fig. 9. The following task is to merge all point clouds from every measurement into one whole.

Point Clouds Preprocessing for Merging

Nowadays there are many preprocessing methods, algorithms, filters and analysis known. But for point clouds merging only a few of them are suitable. The main focus should be concentrated mainly on these three areas:

features, keypoints, registration.

Features

Features is general term for many characteristics of point clouds. It is possible to compute density of points, normal vector of points, collect neighbours, compute curvature of surface etc. From this list the most suitable features for preprocessing before merging are the computation of normal vectors and curvature.

In this experiment the Point Cloud Library [START_REF]Point Cloud Library[END_REF] (PCL) is used. The normal vectors are in this library computed by eigenvalues and eigenvectors. The algorithm is described in more details in [START_REF] Rusu | Semantic 3D Object Maps for Everyday Manipulation in Human Living Environments[END_REF]. The knowledge of this algorithm is very important for the computation of curvature and segmentation. Without this knowledge it is not possible to understand used terms and parameters in PCL implementation.

Keypoints

The aim of keypoints is to select a subset of points from the point cloud which will represent the whole point cloud in the processing steps. The subset selection should represent the point cloud, but the reduced amount of points will increase speed of the following processing. This is similar procedure as searching for correspondence in stereometry. The keypoints detection can be performed by many algorithms. The description and the comparison and recommendation of the best algorithms after testing were introduced by Alexandre L.A. [START_REF] Alexandre | 3D descriptors for object and category recognition: a comparative evaluation[END_REF] and Filipe S. [START_REF] Filipe | A comparative evaluation of 3D keypoint detectors in a RGB-D object dataset[END_REF].

Registration

The previous algorithms are designed for usage with single point cloud and they detect important attributes in point cloud. However, the registration is directly designed for two point clouds merging. The principles of registration process based on planar surfaces are described in [START_REF] Xiao | 3D point cloud registration based on planar surfaces[END_REF] and advanced registration for outdoor environment was introduced in [START_REF] Xiao | Planar Segment Based Threedimensional Point Cloud Registration in Outdoor Environments[END_REF]. The implementation of registration is based on features and keypoints.

Preprocessing Recommendation

Even there are many algorithms for point clouds preprocessing, all of them mentioned above are designed and tested for point clouds representing the small area in front of a robotic device. The point cloud acquired by the 3D LRF introduced above captures the whole environment around robotic device, not only in front. Thus problem of merging the point clouds is more complex. The merging process should be based on main principles introduced in this section, but must be designed to accept different type of point clouds.

Point Clouds Segmentation

The point cloud acquired by the presented 3D LRF covers big part of the environment around the device. In this situation the point cloud contains much more information than from devices which acquire only a small range of the environment. In this case it is possible to expand small details -keypoints -to larger object. Every point cloud can be divided to separated segments. The number of segments is significantly smaller than number of points in cloud. Moreover all segments can be represented by some typical features, e.g. center of gravity, normal vector, size, etc. So the usage of segmentation for point cloud preprocessing is promising and will be consequently verified.

Segmentation Algorithms

The PCL library offers a few segmentation algorithms. Most of them are described by R. B. Rusu [START_REF] Rusu | Semantic 3D Object Maps for Everyday Manipulation in Human Living Environments[END_REF].

The algorithm Fitting simplified geometric model is based on searching of certain structures or patterns of points. The Basic clustering techniques is based on space decomposition and clustering by proximity. The Finding edges algorithm is complementary algorithm of the next one. It is segmentation by Region growing.

This algorithm is well known from image processing. But the implementation for point clouds must be based on different attributes. In the PCL it is based on normal vectors and curvature.

A few additional segmentation algorithms are implemented in the PCL library. The testing of all algorithms showed that the best method for segmentation is the Region growing. The comparison is not presented here.

Point Clouds Preprocessing by Segmentation

The example of the segmentation of the point cloud from Fig. 9 is visible in Fig. 10. The figure contains only large segments. The threshold for the large segment was set to 500 points and every segment is distinguished by a color.

Verification of Region Growing Algorithm

Fig. 10 shows that the segmentation method Region growing is working well for the point cloud M 1 . To be sure that this segmentation can be used for the whole data set it is necessary to apply segmentation to every point cloud from data set and create overall statistic of all segmentation results.

All segments of point clouds will be separated into six groups: left, right, front, rear, top and bottom. This separation can be based on normal vector of segments. This vector can be computed in the centroid of all segments, where the centroid can be computed by the following formula:

c s = 1 N s • Ns i=1 p s,i , (2)
where c s is the centroid of the segment s and N s is a number of points of that segment. Vectors p s,i are position vectors of all points of segment s computed by [START_REF] Bosse | Zebedee: Design of a spring-mounted 3-D range sensor with application to mobile mapping[END_REF].

The computation of normal vector of segments in PCL library automatically orients the vector direction to half-space with the origin of coordinate system. So the coordinate of normal vector can divide segments to their proper group. The largest positive coordinate x of normal vector specifies rear segment, the largest positive y specifies right segment, the positive z specifies bottom segment, etc.

The segmentation of all point clouds from the data set was performed and for all large segments their proper group was determined. The overall statistic is visible in Table 1. In this table the column Minimum indicates that in all point clouds at least 1 large top, left and right segment were recognized. The median for these segments is 3. One bottom segment was also recognized in all point clouds. Only front and rear segments are not recognized in all point clouds. This is an expected result, because the measurement was performed in the long narrow corridor. This statistic proved that the segmentation method Region growing is suitable for preprocessing the whole data set of our experiment.

Point Clouds Merging Process

The process of point clouds merging will be performed in four stages. All of them are illustrated for clarity in Fig. 11 The main advantage of proposed 3D LRF design will be used in following process. The advantage is the possibility of scanner to measure the whole environment including a space even above the scanner. Just the ceiling is the most visible part of the environment. It is flat, without obstacle and in long term in the unchanging state. Thus the ceiling will be in center of attention in the most stages of point clouds merging process.

Horizon Leveling

The measuring device does not stay during the measurement always in ideal horizontal position. The floor has minor irregularities which cause small inclination of the device. The tilt in some cases is greater than few degrees.

In this stage the horizon will be balanced according to the ceiling plane. The segment of ceiling is in all point clouds the largest segment of all. It covers the area of approximately 10 to 16 square meters. So the informative value of these segments is very high and it can be used in a simple way. The normal vector of the ceiling segment can be easily computed in centroid of the segment ceiling. This normal vector will consequently define a new direction of the vertical axis z.

The balancing of the whole point cloud will be realized by rotations around axes x and y. The necessary angles of rotations can be easy computed from the angle between the normal vector n and the axis z. All points in every point cloud will be transformed by following formula:

p r,s = R y (α y ) • R x (α x ) • p r,s . (3)

Direction Alignment

In this stage the ceiling segment will be reused again. The direction of the point cloud will be based on left and right edges of the ceiling segment. The top view of this situation is depicted in Fig. 12. The measuring device is in general position in the origin of the coordinate system O. The virtual lines, parallel with the axis y, are regularly marked at distances d on both sides of the axis y. On the left and right side these lines leaves ceiling segment in edge points L i and R j . The connecting lines of these points creates a set of vectors l i and r j . Now these vectors will be used to determine the direction in two steps: estimation of direction and improving its accuracy.

d d d d d d x L 1 y l 1 L 2 L 3 L 4 L 5 L 6 L 7 l 2 l 3 l 4 l 5 l 6 R 1 R 2 R 3 R 4 R 5 R 6 R 7 r 1 r 2 r 3 r 4 r 5 r 6
O Fig. 12. The segment representing the ceiling used for the direction alignment

The estimation of the direction will be computed separately for the left and right side of the ceiling segment using the following formulas:

l = M i=1 l i M i=1 |l i | , r = N j=1 r j N j=1 |r j | , (4)
where M and N is a number of left and right vectors. The estimation will be used for a simple filtering. Vectors from the left and right set, which divert from the vector l and r more than the given angular threshold, will be removed. In this experiment this threshold is set to 10 degrees. This simple filtering removes bounces and niches in wall from both ceiling edges. Removed vectors are marked in Fig. 12 by dotted line.

The number of vectors will be reduced by this filtering to number M and N and final direction d is computed from all remaining vectors by the formula:

d = M i=1 l i + N j=1 r j M i=1 |l i | + N j=1 |r j | . ( 5
)
Now it remains only to calculate the angle α z between the vector d and the axis x. Every point cloud will be rotated around the axis z by this angle:

p r,s = R z (α z ) • p r,s . (6)

Centering

This is the first stage where a pair of neighbouring point clouds will be used. Furthermore, the odometric value of the shift between two measurements is needed. As was mentioned above, this distance is approximately 1 meter. This distance will be marked as d o .

In this stage, two neighbouring point clouds will be placed side by side, on its joined axes x at the distance d o . This situation is partially visible in Fig. 13 Now it is necessary to compute distances between all points L i of cloud a and edges l j of cloud b and vice versa. The same procedure is applied for right edges of both point clouds. The distances between points and removed edges are not computed.

The result is a set of distances y k and the number of distances can be marked as K. The average value of all y k can be computed by the following formula:

y a = K k=1 y k K . (7)
The average distance will be used for translation of the whole point cloud b by distance y a :

p r,s = T y (-y a ) • p r,s . (8)

Insertion

The last stage of point clouds merging will use odometric distance d o and corresponding front and back segments of neighbouring point clouds.

p r,s = T x (d o ) • p r,s . (9)
Now it is possible to merge all points from point cloud a and b into single point cloud.

Merged Point Cloud

The merging process described in the previous section must be used for all point clouds in the data set as well as for all neighbouring point clouds. The result of this process is a single point cloud that merges 25 point clouds. The result is visible in Fig. 14. This point cloud is composed of 5.4 × 10 6 points. Thus in this figure it seems solid, not as a point cloud.

Fig. 14 shows that the merging process did not create an exact straight object. The result is a little bit curved. But this result was achieved by using point clouds only without any additional position or incline sensors. For further usage it is necessary to reduce the amount of points of the merged point cloud. For the reduction of points a voxelization algorithm is used. It this process the whole space is divided into small cubes, for this experiment the cube size is 100 mm, and all points inside every cube is reduced to a single point. Obtained results directly correspond to the Octomaps used for the robot navigation [START_REF] Wurm | OctoMap: A probabilistic, flexible, and compact 3D map representation for robotic systems[END_REF]. The result of this process is visible in Fig. 15. This point cloud is composed of only 46 × 10 3 points. It is less than 1% of the original merged point cloud size. The segmentation algorithm can be repeatedly applied on the merged point cloud with reduced number of points.

The result of this segmentation is visible in Fig. 16. There is clearly visible that the segmentation recognizes the ceiling very well as the single (red) segment. The same result is for the floor which is depicted as a single (blue) segment. In Fig. 17 it is also visible that the segmentation can very well recognize niches and doors. The resulting Based on this result, it is possible to say that the process of point clouds merging was successful.

Conclusion

The paper introduces a novel concept of using cheap 2D LRF as a 3D scanning device. In contrast to more advanced devices, additional procedures of calibration and point cloud merging is necessary, however the final results are comparable with the more expensive devices. This article focuses on the process of merging point clouds and describes it in more details. Special properties of indoor measurements are used as a basis for establishing merging parameters. Merged point clouds represent whole environment and can be used after the segmentation as a source for the robot navigation and environment mapping. The experiments showed that for the successful merging only a small part of raw measurements is needed (in the described example less than 1% of acquired points were used).

Fig. 3 .Fig. 4 .Fig. 5 .

 345 Fig. 3. Kinect [http://www.microsoft.com]

Fig. 7 .

 7 Fig. 7. 3D LFR schema Fig. 8. 3D LRF overall design

Fig. 9 .

 9 Fig. 9. Acquired point cloud M1

Fig. 10 .

 10 Fig. 10. Point cloud M1 after segmentation

  using a matchbox. In sub-picture 11 a) depicts the general startup position of two neighbouring point clouds. Subsequent sub-pictures depicts following four stages of merging process: 1. Fig. 11 b) -the horizon leveling. 2. Fig. 11 c) -the direction alignment. 3. Fig. 11 d) -the centering. 4. Fig. 11 e) -the insertion.

Fig. 11 .

 11 Fig. 11. Illustration of point clouds merging stages on the process of closing a matchbox

Fig. 13 .

 13 Fig. 13. The left edges of two point clouds

Fig. 14 .

 14 Fig. 14. The top view and the side view of the whole corridor -25 merged point clouds

Fig. 15 .

 15 Fig. 15. The merged point cloud after voxelization (100 mm)

Fig. 16 .

 16 Fig. 16. Side view of the resulting segmentation of the reduced point cloud

Fig. 17 .

 17 Fig. 17. Top view of the resulting segmentation of the reduced point cloud

Table 1 .

 1 Statistic of the whole data set segmentation

	Segments	Number of segments Average Median Minimum Maximum
	Left	3.1	3	1	7
	Right	3.3	3	1	8
	Top	3.2	3	1	7
	Bottom	1.1	1	1	3
	Front	1.2	1	0	4
	Rear	1.2	1	0	4

  . The illustration depicts only left edges of two neighbouring point clouds. Elements of the first point cloud are marked with additional index a and elements of the second point cloud with index b.

								y a				y b
									L 5a
	L 1a	l 1a	L 2a	l 2a	L 3a	l 3a	L 4a	l 4a			l 5a	L 6a	l 6a	L 7a
				y 1 y 2 y 3					y 4 y 5 y 6
				L 1b	l 1b	L 2b	l 2b	L 3b d o	l 3b	L 4b	l 4b	L 5b	l 5b	L 6b	l 6b	x L 7b

This work was partially supported by the Grant of SGS No. SP2016/58, V ŠB -Technical University of Ostrava, Czech Republic.


Tomáš Martinovič
email: tomas.martinovic@vsb.cz

Kateřina Janurová
email: katerina.janurova@vsb.cz

Kateřina Slaninová
email: katerina.slaninova@vsb.cz

Jan Martinovič
email: jan.martinovic@vsb.cz

Automated Application of Inventory Optimization

Keywords: inventory optimization, sales forecasting, model selection, parallelization

We present automated application of inventory optimization based on sales forecast. Inventory stock optimization is very required issue by companies recent years, however inventory models are based on the sales expectation. Therefore, the problem of optimizing inventory stock is divided into two parts, sales forecast and setting optimal inventory. We describe an automated solution to model selection for sales forecast and the inventory setting based on those predictions. In the end, we present our validation of the system through historical simulation. We compare simulations results against real inventory levels. Due to the large number of different length time series, this simulation was run in parallel on cluster and was parallelized in R. The algorithms were developed and tested on inventory time series from real data sets of the K2 atmitec company 1 .

Introduction

The goal of our work was the design and implementation of automated application of inventory control for optimization and simplification of the decision making process of sales managers. Models for inventory control are based on the products sales expectations. Therefore, as the first step of inventory control, the sales forecast needs to be estimated. In our application, the sales forecast can be determined by four different models applied on several modifications of sales history. The result of the model with the best prediction according to the selected measure of accuracy is used as the input for the inventory control model.

Because of the lack of information needed for the standard inventory control models like Economic Order Quantity model or newsvendor model, we had to base our approach for the second step of inventory control on a fixed time review period model with variable demand and modify it to better suit our needs. This modification was verified by the simulation, which imitates behavior of sales manager on the assumption of complete trust in our application. The results of such simulation were compared with the real historical data of inventory time series.

The simulations were parallelized in R language and run in parallel on cluster because of the large number of products used for verification and repeated model estimations.

Inventory Stock Optimization

As described above, the inventory stock optimization consists of the two steps: sales forecasting and optimal inventory setting.

Sales Forecasting

Sales forecasting represents a typical application of time series analysis and prediction. Classical time series analysis is based on the assumption that the time series y t for t = 1, 2, . . . , n can be decomposed into four components: trend, cyclic, seasonal and irregular. The trend component characterizes a long-term tendency of the examined phenomenon that is influenced by factors operating in the same direction eg. production technology, demographic conditions, market conditions, etc. The first approach we used in time series modeling is based on the idea that if the long-term tendency of time series corresponds to a certain function of time eg. a linear, quadratic, exponential, S-curve etc., the trend component can be modeled using this function [START_REF] Hyndman | Forecasting: principles and practice[END_REF].

Assume that the time series y t can be written as the sum of trend and residual process

y t = T t + ε t t = 1, 2, . . . , n, (1)
where T t can be described by the appropriate mathematical function of the time variable t, and ε t is the residual component with the character of white noise process.

We used three different mathematical functions for modeling the trend components in our simulation:

1. linear trend function in the form

T t = β 0 + β 1 t, t = 1, 2, . . . , n, (2)
2. quadratic trend function in the form

T t = β 0 + β 1 t + β 2 t 2 , t = 1, 2, . . . , n, (3)
3. and cubic trend function in the form

T t = β 0 + β 1 t + β 2 t 2 + β 3 t 3 , t = 1, 2, . . . , n. (4)
The regression parameters of the trend functions β 0 , β 1 , β 2 and β 3 were estimated by the least squares estimation method. Although this approach may appear to be too simple, in many practical cases the best solution is the simplest one.

The second approach we used in time series modeling is provided by the ARIMA model [START_REF] Asteriou | Applied Econometrics[END_REF]. ARIMA -Autoregressive integrated moving average -is part of the Box-Jenkins methodology, which was developed for the analysis of time series, where we assume that the residual component ε t does not have the character of white noise. Unlike the first approach, where the regression model is built only as the function of time, the Box-Jenkings models are built as the function of linear combination of p past values of the modeled variable itself, hence the term autoregressive, and regression-like model of q + 1 past forecast errors. In such way the ARMA(p, q) -autoregressive moving averages -model is obtained. ARIMA(p, d, q) model is then essentially an ARMA model applied to the dth difference of the original series. Such approach is used in the cases when the studied series is not stationary (it's properties depend on the time at which the series is observed). The stationarity can be achieved by applying the difference of the first or the second order, difference of the higher order are usually not used in practice.

The final sales forecast is made in four steps. At first, all the models described above are fitted to the training part of the data, y t for t = 1, 2, . . . , m, where m < n, then the sales forecasts ŷt for t = m + 1, m + 2, . . . , n are made on the basis of fitted models. The different forecasts are compared with the test part of the data y t for t = m + 1, m + 2, . . . , n using appropriate measure of accuracy, and finally the best forecast is chosen as the input for inventory control model. Additionally, the models are fitted simultaneously on the original daily data and data with weekly and monthly aggregation.

For the best forecast selection, several measures of accuracy eg. RMSE, MAE, MAPE, etc. were considered [START_REF] Palit | Computational Intelligence in Time Series Forecasting: Theory and Engineering Applications[END_REF]. However, because of the different time aggregation of the compared time series, the accuracy measure, that is independent of the scale of the data, seems to be the most convenient one. This kind of accuracy measure, Mean Absolute Scaled Error, was proposed by Hyndman and Koehler [START_REF] Hyndman | Another look at measures of forecast accuracy[END_REF] as

MASE = mean (|q t |) , (5)
with the scaled error q t defined as

q t = e t 1 n-m-1 n i=m+2 |y i -y i-1 | , (6)
where e t is the forecast error defined as e t = y t -ŷt , t = m + 1, m + 2, . . . , n.

Inventory Control

Numerous models for inventory control were developed such as Economic Order Quantity model (EOQ) or newsvendor model. Those are the basic deterministic and stochastic models for the inventory control. Based on the stochastic nature of our problem and available data we used fixed time review period model with variable demand as a basis for our solution.

It was difficult to estimate ordering costs and holding costs due to the character of the available data. Therefore, we use just quantities and we focus on providing inventory level measures based on the user specified service level in our approach. To achieve that, three measures need to be computed. The first is the safety stock, which should cover most of the uncertain demand. This uncertainty means that it would be too costly, and often impossible, to hold enough inventory to cover all the variations in demand. Therefore, the user of the software may set requested level of service. Safety stock then depends on standard deviation of sales, lead time, service level and probability distribution of sales. Assuming standard normal distribution of sales, constant lead time and constant variation of sales, we can define the safety stock [START_REF] Russell | Operations Management[END_REF] as

S = ϕ(sl) t l σ 2 d , (7)
where ϕ(sl) is the quantile function for normal distribution and its quantile is given by service level. Variable t l stands for lead time and σ d is standard deviation of demand. By leaving the assumption of constant lead time and adding uncertainty to lead time we get

S = ϕ(sl) t l σ 2 s + E(s)σ 2 l , (8)
where E(s) is expected mean sales and σ l is lead time standard deviation. Now we can define reorder point, which is the level of inventory needed to cover the demand until the arrival of the next batch of inventory plus safety stock. If we assume the constant demand and constant lead time we can define reorder point as follows

R = S + E(s)t l . (9)
The last step is to compute the order quantity. This is defined as the quantity necessary to cover the demand for the user specified time. Hence, the order quantity is affected by mean expected sales, order time, incoming goods and goods reserved in advance. Additionally, the ordered quantity should be adjusted by the difference between actual stock and reorder point. We can write this as

Q = E(s)t o -G i + G r + R -A, (10)
where t o is the order time, G i is the incoming goods quantity, G r is reserved goods quantity and A is actual level of inventory. This definition of order quantity covers the demand for the defined time with set service level for covering the stochastic behaviour of customers.

Automated Application of Inventory Optimization

The theory described in sections 2.1 and 2.2 was used for development of an approach for inventory optimization that was implemented as part of ERP software by K2 atmitec company. This software can simplify the decision making process of sales managers. It includes three main parts, which are "best model selection", "sales prediction" and "inventory control".

Verification of Automated Application of Inventory Optimization

To evaluate our approach, we created a simulation which imitates behavior of sales manager on the assumption of complete trust in our approach. The simulation includes all three parts mentioned in previous section, which are complemented by an automatically controlled decision making process for running specific parts. We used historical time series of sales and inventory stock levels for verification of our simulation, the results are given in section 3.3.

Data Description

The real historical data for our software verification was provided by the K2 atmitec company. The dataset contained the information about approximately 33 000 products in several warehouses in Czech Republic. We were not able to simulate the behavior for all of the provided products due to the small number of observations in the followed history. The software was finally verified on about 950 and 750 products in two selected warehouses, what was sufficiently large enough number.

Simulations

On the beginning of the simulation, some parameters must be set. The lead time for the given goods, order interval and maximum level of inventory. The lead time and maximum level of inventory are variable for every type of stock. The review period was set to 30 days. In addition, we had to select the service level, which was set to 95%. The starting point for the simulation is a year and a half. That means we take the first year and a half of sales time series and the inventory level at that time. The "best model selection" part begins with running the model selection algorithm, when the original time series is preprocessed in several steps. In order to find the best model for given data, the time series with originally daily sales is aggregated in weekly and monthly sales. There is a possibility of removing outliers for each aggregation and also filtering by Kalman [START_REF] Durbin | Time Series Analysis by State Space Methods[END_REF] or Baxter-King [START_REF] Baxter | Measuring Business Cycles Approximate Band-Pass Filters for Economic Time Series[END_REF] filter for weekly and monthly aggregation. At the end of preprocessing step, each original time series leads to 10 final modifications on which the four models given in section 2.1 are fitted and compared by the MASE accuracy measure. This chooses a predictive model used for the next three months. Then, the simulation runs the "sales prediction" part, that calculates the predicted sales quantity on the basis of previously selected predictive model. The result of this step is the predicted time series, that is used as the input for the "inventory control" After the first day, the simulation iterates by day, reducing the inventory level by the same amount as it is reduced in the real inventory time series. When it reaches the reorder point, the "sales prediction" and "inventory control" parts are run again. The reorder point and safety stock are updated according to the new results and order for the new inventory is made. The ordered inventory is added to the simulated inventory after the lead time period. Every three months new predictive model is chosen based on the "best model selection" part. This process is repeated until the end of the time series.

The main output of this simulation is simulated inventory time series. This may be compared with the actual inventory time series, as is shown in Fig. 2. In general, if the mean level of simulated inventory is less than the actual inventory and the stockout time is within boundaries of chosen service level, we can say, that the simulation is successful. In addition, to the simulated inventory time series, we have other outputs. Those are chosen predictive models, history of inventory control settings, orders history and graphical output. These outputs allow us to make a quick analysis of reasons for the success or failure of our simulation. Because of a large number of products at each warehouse and repeated model estimations, we had to parallelize this algorithm. Due to the independency of simulations for specific goods, the computation was easily parallelized through Rmpi2 package in R software. This package uses MPI for communications between nodes. Our algorithm uses distributed apply function, where a master node sends time series data to free slave nodes, which execute the simulation on given data. This method allows us to run this algorithm both single node or on the cluster and scales well, because there is minimal need for the communication. The computation was made on Salomon cluster equipped with Intel Xeon E5-2680v3 processors. We used 8 compute nodes, each with 24 cores and 128GB RAM.

DateX

In the future, we would like to make improvements in the computations of reorder point and safety stock, since they are very susceptible to sudden changes in demand as is shown in Fig. 3. The most sensible seems to be the solution by moving averages since this will allow for changes to happen, but it would smooth them at the same time. However, this approach needs to be tested.

Comparison of Simulation Results and Real Inventory Data

Simulations were completed on two warehouses, one with 950 products and the second with 761 products. For each simulation, the mean level of real inventory and the mean level of simulated inventory were determined, and the relative percentage difference with respect to mean level of real inventory was calculated. The higher relative percentage difference, the more successful simulation was. Stockout was calculated for the sake of real information as the number of days for which the product was unavailable. The mean and median values for the relative percentage difference as well as the mean values of stockout days for both warehouses are given in Tab.1 As can be seen from the summary statistics, the mean stouckout is within the required level of 5 % for both warehouses. The mean level of relative percentage difference was 17.4 % and 20.0 %, with the meaning that if the sales manager trusted completely in our software, the mean inventory level would be lower by 17.4 % and 20.0 % The values of median relative percentage difference suggest, that the decrease could be even more significant, and that the mean relative percentage difference value is influenced by several unsuccessful simulations, as can be seen from the distribution of relative percentage difference with removed outliers in Fig. 4 and Fig. 5. Those were mainly simulations on products with the lead time greater than 40 days. Improvement of the simulation on such products is one of the tasks for our future research.

Conclusion

Our paper describes design, implementation and verification of the automated application of inventory control that is based on predictions of real historical sales data. The sales prediction is estimated by several approaches of time series modeling of which the best is chosen by the selected measure of accuracy. Modification of fixed time review period model with variable demand was proposed for inventory control model and verified by the described simulation process which compares inventory time series simulated by our application with the real inventory time series. The results of the simulation showed, that if the sales manager acted completely in accord with our application recommendation, the mean inventory level would decrease by 17.4 % in one tested warehouse and by 20.0 % in the second tested warehouse.

The future work in this area involves, besides the tasks mentioned in previous section, the research and the implementation of other suitable models for sales forecasting and their ensembles, an improved model selection, and alternative inventory control models.

Fig. 1 .

 1 Fig. 1. Prediction of sales quantity of real inventory time series

Fig. 2 .

 2 Fig. 2. Comparison of simulated inventory time series and real inventory time series

Fig. 3 .

 3 Fig. 3. Inventory quantity, safety stock, reorder point and order quantity history for inventory control simulation

Fig. 4 . 1 2 RelativeFig. 5 .

 4125 Fig. 4. Distribution of relative percentage difference with respect to real inventory in Warehouse 1

Table 1 .

 1 Summary statistics for the comparisons of the mean level of real inventory and the mean level of simulated inventory.

	Warehouse 1 Warehouse 2

Sales forecasting and inventory optimization are parts of ERP solution K2, which is provided by K2 atmitec s.r.o. company, http://www.k2.cz/en/.

Yu, H. (2002). Rmpi: Parallel Statistical Computing in R. R News 2/2:10-14.

Acknowledgment

This work was supported by The Ministry of Education, Youth and Sports from the National Programme of Sustainability (NPU II) project 'IT4Innovations excellence in science -LQ1602' and from the Large Infrastructures for Research, Experimental Development and Innovations project 'IT4Innovations National Supercomputing Center -LM2015070'.


Wladyslaw Homenda

Tomasz Penza

Mapping points back from the concept space with minimum mean squared error

In this article we present a method to map points from the concept space, associated with the fuzzy c-means algorithm, back to the feature space. We assume that we have a probability density function f defined on the feature space (e.g. a normalized density of a data set). For a given point w of concept space, we give explicitly a set of points in feature space that are mapped onto w and we give a formula for a reverse mapping to the feature space which results in minimum mean squared error, with respect to density f , of the operation of mapping a point of feature space into the concept space and back. We characterize the circumstances under which points can be mapped back into the feature space unambiguously and provide a formula for the inverse mapping.

Introduction

There is an effort among researchers today, to capture the mechanisms of human cognition that allow the mind to process abstract information. Their goal is to enable the machine to discern patterns, spot relations and make connections between objects and ideas, to reason and predict on the basis of knowledge and observation. To make this possible, computers must be able to form and process concepts. A cluster in a data set can be thought of as representing an abstract concept. We can thus use data clustering methods to extract concepts. With each clustering method, comes a way to measure to what degree a given data point belongs to a given cluster. If we regard these degrees of membership in different clusters as coordinates, we obtain a concept space of dimension equal to the number of clusters. We can then map numeric data from the feature space into the concept space and process it at the level of concepts. Many points of the feature space may end up being mapped onto the same point of concept space. Consequently, we run into a problem when we try to put the result of conceptlevel computations back into the feature space. In certain circumstances, if we want to have a reverse mapping from the concept space into the feature space, we have to accept that it will make errors, i.e. when we map an observation into the concept space, this reverse mapping will not map it back onto itself.

In this text we work with a transformation from the feature space to the concept space associated with the fuzzy c-means algorithm which we denote with symbol M. We will describe sets of points that M maps onto the same point of concept space. We will also determine under what conditions M is invertible and give a formula for its inverse. We will assume that there is a continuous probability density function defined on the feature space and we will present a formula for a reverse mapping G that allows minimum mean squared error (MSE), with respect to that probability density, of an operation of mapping an observation to the concept space via M and back via G. In practice such probability density might arise as a density of data points of some data set in the feature space. We propose an application of this reverse mapping as a defuzzification-like transformation of the results of concept-level processing, which is the last stage of the process outlined at figure 1.1. There is a long record of defuzzification methods and we refer to selected results in References without explicit discussion on them. Our approach is different in that we study a general reverse mapping and defuzzification discussed in this paper is a special case of such general transformation.

Preliminaries

Fuzzy c-means algorithm

An analogue of k-means algorithm within fuzzy clustering is the fuzzy c-means algorithm (FCM) which lets us find c fuzzy clusters in a data set D = {x 1 , x 2 , . . . , x n } -a subset of feature space R h . We denote the i-th cluster by C i and its centroid as µ i . For a given observation x i , its degree of membership in the cluster C j is denoted by w ij and given by the formula

w ij = ||x i -µ j || -2 m-1 c k=1 ||x i -µ k || -2 m-1 = 1 c k=1 ||x i -µ j || ||x i -µ k || 2 m-1
where m > 1 is a parameter called the fuzzifier. For every data point, its membership in different clusters sums up to 1. Each fuzzy cluster is represented by its weighted mean µ i = n j=1 w ji x j / n j=1 w ji . These two formulas are mutually referential. At each step of the algorithm only one of them will be satisfied, but the clustering that the algorithm will output will approximately satisfy them both.

The fuzzifier m controls how fuzzy the resulting clustering will be. The higher the fuzzifier, the less distinct the values of membership and thus the fuzzier the clustering. When m converges to 1, for each observation one of its cluster memberships converges to 1 and the rest converge to 0, so in the limit the clustering becomes hard.

Before running the algorithm, we must choose the number c of clusters that we want to form, the fuzzifier m and the threshold ε. The threshold will let us terminate the algorithm when upon iteration all the memberships w ij changed by less than ε. The algorithm is as follows.

1. Randomly initialize all the degrees of membership w ij with numbers from the unit interval. 2. Compute all the weighted means µ i using current values w ij . 3. Compute all the membership degrees w ij using current values µ i . 4. If any of values w ij changed by more than ε, go to step 2.

FCM endeavors to minimize the loss function

n i=1 c j=1 w ij ||x i -µ j || 2
which never increases after an iteration of the loop. Thanks to the threshold ε the algorithm always converges. The final clustering depends on the initial values w ij , which are assigned randomly, thus the clusters will vary each time we run FCM. It is advised to run the algorithm several times and choose the best result. We can evaluate clusterings using cluster validity indices. They can also help us pick an optimal number of clusters c.

Mapping M

Suppose that in the feature space R h , we used FCM with fuzzifier set to m to find c fuzzy clusters and we obtained the set of centroids C = {µ 1 , µ 2 , ... , µ c }. Mapping M C : R h -→ [0, 1] c , assigns to x ∈ R h a vector of its degrees of membership in the respective clusters, which is a point of the concept space [0, 1] c . For x ∈ R h \C, we have

M C (x) = 1 c k=1 ||x -µ i || ||x -µ k || 2 m-1 i=1, 2, ..., c
We extend M C continuously onto R h by defining M C (µ i ) to be a vector whose i-th coordinate is 1 and all its other coordinates are 0. For now we will write simply M, but we will invoke the subscript in section 3.3, when we will need to refer to mapping M C for a subset C of the set of centroids C. M is continuously differentiable on R h \C. We will invert M in circumstances when it is possible and otherwise we will provide a formula for the most accurate reverse mapping to the feature space.

If we exclude the trivial case of just one centroid, the values of M on set R h \C are all contained in the subset W = (w 1 , w 2 , ..., w c ) ∈ (0, 1) c c i=1 w i = 1 of concept space. For c = 3 this subset would be an interior of a triangle inside the cube [0, 1] 3 with vertices (1, 0, 0), (0, 1, 0) and (0, 0, 1), i.e. points of the concept space corresponding to the centroids. Even though elements of W have c coordinates, they have only c -1 degrees of freedom, because the last coordinate is given by the first c -1 coordinates:

w c = 1 - c-1 i=1 w i .
Let P denote the projection onto the first c -1 coordinates, i.e. P (w) = P (w 1 , w 2 , ..., w c ) = (w 1 , w 2 , ..., w c-1 ). For any w ∈ W and X ⊂ W , we will denote their images P (w) and P (X) under P by w and X respectively. Set

W = (w 1 , w 2 , ..., w c-1 ) ∈ (0, 1) c-1 c-1 i=1 w i < 1 is an open subset of R c-1 .
For any point w in W , w c can be reconstructed using the equation above, so P is a bijection. Identity w = (w, w c ) relates the elements of these two sets and whenever symbols w and w appear in the same context, they are always bound by this identity. P is actually a homeomorphism between W and W which makes W a (c -1)-dimensional manifold.

3 Sets mapped onto the same point of concept space Preimages of points of W under M, when they are nonempty, are spheres of some dimension embedded in R h , except for a special case when the preimage is a linear manifold. We define R 0 to be a space containing only a zero 1 x 1 vector and B 0 (0, R) and S -1 (0, 0) to be equal to R 0 .

Preimage of a point of concept space

For a set containing at least two elements, C = {µ 1 , µ 2 , ... , µ c } ⊂ R h , a point in W , w = (w 1 , w 2 , ..., w c ) and m > 1, we define the following values

1. M is a (c -1) × h matrix whose i-th row is (µ i -µ c ) T 2. u is a column vector of length c -1 such that u i = 1 2 ||µ i -µ c || 2 3. v is a column vector of length c -1 such that v i = -1 2 (w 1-m i -w 1-m c ) 4. M + is the Moore-Penrose pseudoinverse of M 5. Q = I -M M + 6. a = M + u, b = M + v 7. r = rankM , d = h -r 8. If d > 0, U is a h × d matrix whose columns form an orthonormal basis of the linear subspace {(I -M + M )x | x ∈ R h } of R h of dimension d. If d = 0, U is a zero vector of R h 9. If Qu = 0, matrix H is U b ||b|| for d > 0 and b ||b|| for d = 0. If Qu = 0, it is U .
Matrix U can be obtained by selecting a maximal set of linearly independent columns of matrix I -M + M , orthogonalizing it, normalizing each vector and taking the resulting set of vectors as columns of U . Values of M , u, a, Q, d and U depend only on C, while values of v, b and H depend on both C and w. We will treat v, b and H as functions of variable w = (w 1 , w 2 , ..., w c-1 ) ∈ W . These functions are continuous.

We need above definitions to provide a description of the preimage M -1 (w) of w ∈ W , which is split into two cases depending on if the set of centroids lies on a sphere or not. In each case if the conditions are not met, then the preimage is empty (no point of R h is mapped onto w). We define w eq = ( 1c , 1 c , . . . , 1 c ). Symbol S n stands for the unit n-sphere

S n = y ∈ R n+1 ||y|| = 1 which is a subset of R n+1 .
Set of centroids lies on a sphere (equivalent to Qu = 0)

If Qv = 0 and w 1-m c ≥ 2(a • b + ||a||||b||), then the preimage is nonempty. For w = w eq , it is a d-sphere of radius R embedded in R h centered at point p M -1 (w) = p + RHy y ∈ S d with R = (2a • b -w 1-m c ) 2 -4||a|| 2 ||b|| 2 2||b|| p = µ c + a + w 1-m c -2a • b 2||b|| 2 b
For w = w eq the preimage is a linear manifold of dimension d

M -1 (w eq ) = µ c + a + U y y ∈ R d
Set of centroids does not lie on a sphere (equivalent to Qu = 0)

If Qv = 0, Qu • Qv = -||Qu||||Qv|| and
||b|| 2 ||Qu|| 2 + (2a • b -w 1-m c )||Qu||||Qv|| + ||a|| 2 ||Qv|| 2 ≤ 0 ||b|| 2 ||Qu|| 2 + (2a • b -w 1-m c )||Qu||||Qv|| + ||a|| 2 ||Qv|| 2 = 0, when d = 0 then the preimage is nonempty. It is a (d -1)-sphere of radius R embedded in R h centered at point p. If d > 0 M -1 (w) = p + RU y y ∈ S d-1
and if d = 0

M -1 (w) = p with R = -||b|| 2 ||Qu|| ||Qv|| 2 -(2a • b -w 1-m c ) ||Qu|| ||Qv|| -||a|| 2 p = µ c + a + ||Qu|| ||Qv|| b
Let n represent a dimension of the spherical preimages for a given C. If preimage of w = w eq is nonempty, then it is given by M -1 (w) = p + RHy y ∈ S n . Columns of h × (n + 1) matrix H form an orthonormal set of vectors which entails that a linear transformation defined by this matrix is distance-preserving. It maps objects from R n+1 into R h without distortion, so the formula for a preimage describes embedding of S n into R h , scaling it to have radius R and the translating it so that it is centered at poin p.

Functions p and R are defined on set W \{w eq }. They are continuous which together with continuity of H implies that a slight change in w results in only a slight change of position, radius and spatial orientation of sphere M -1 (w). If d > 0, then the set of centroids lies on some linear manifold of dimension smaller than h. Mapping M is not invertible if and only if Qu = 0 or d > 0, in other words if and only if the set of centroids lies on a sphere or on a linear manifold of dimension smaller than h. Thus, if there are at least h + 1 centroids, then M is typically invertible (unless they are unfortunately positioned) and we can map points from the concept space back into the feature space precisely. On the other hand, if there are at most h centroids, then they always lie on some sphere, so M is not invertible. If M is invertible, then M -1 (w) = p(w).

Derivatives of b, p and R.

Functions b and p are continuously differentiable everywhere on W \w eq and R on its subset where it is nonzero. We introduce symbols for three (c -1)

× 1 vectors 1 = 1, 1, . . . , 1 T , v I = (1 -m) w -m 1 , w -m 2 , . . . , w -m c-1
T and

v c = (1 -m) w -m c , w -m c , . . . , w -m c
T . We let diag(x) denote a diagonal matrix with coordinates of vector x on its diagonal and we put s = R 2 + ||a|| 2 ||b|| 2 . We can now write the differentials

Dv = -1 2 v c 1 T + diag(v I ) Db = M + Dv
In case of Qu = 0 we have

∇R = 1 R Db T √ sa -sb + √ s 2 v c Dp = 1 ||b|| 2 b Db T 2 √ sb -a - 1 2 v c T - √ sDb
In case of Qu = 0 we define t 0 = ||Qu|| ||Qv|| and we have

∇R = t 0 2R 1 ||Qv|| 2 2t 0 ||b|| 2 + 2a • b -w 1-m c Dv T Q T Qv -2Db T a + t 0 b -v c Dp = t 0 M + I - 1 ||Qv|| 2 vv T Q T Q Dv 3.

Reduced set of centroids

In this section, we will specify a c element subset C of the set of centroids C, such that for every w ∈ W , preimage M -1 (w) under M is the same as preimage of some w ∈ W under mapping M C : R h -→ [0, 1] c , which corresponds to the set of centroids C. We will use the hat symbol to indicate that a given object corresponds to this reduced set of centroids, e.g. b is b computed for the set of centroids C.

The i-th row of M is of the form (µ i -µ c ) T , so let's say that i-th row of M and centroid µ i correspond to each other. Let S be a set containing µ c and all the centroids that correspond to some set of r linearly independent rows of M . If Qu = 0, then we set C = S. This set lies on a sphere. If Qu = 0, we take another centroid µ k such that S ∪ {µ k } does not lie on a sphere and we set C = S ∪ {µ k }. We can find µ k by checking for each centroid if amending S with it results in Q u = 0.

Reverse mapping with minimal MSE

Let's assume that there is a continuous probability distribution, with a continuous probability density function f , defined on the feature space R h . Function f can be interpreted as density of points of some data set D ⊂ R h . For a real data set, such density can be approximated using statistical density estimation methods. Figure 2 shows an example of a simple probability density function derived from a data set. In this chapter, we will look for a reverse mapping G from the concept space to the feature space which allows minimum mean squared error (MSE), with respect to f , of an operation of mapping points of feature space to concept space via M and then back to feature space via G, cf. figure 3. Let X be a random vector distributed according to density f . MSE of this operation is given by expression E X -G(M(X)) 2 , so the function that minimizes it is equal to conditional expectation G(w) = E X M(X) = w almost everywhere.

Parameterization of feature space

For n > 0, we parametrize almost all of the unit sphere S n ⊂ R n+1 by n angular coordinates θ = (θ 1 , θ 2 , ..., θ n ) that range over set Θ = (0, π) n-1 × (0, 2π) ⊂ R n .

y(θ) =          cos(θ 1 ) sin(θ 1 ) cos(θ 2 ) sin(θ 1 ) sin(θ 2 ) cos(θ 3 ) . . . sin(θ 1 ) • • • sin(θ n-1 ) cos(θ n ) sin(θ 1 ) • • • sin(θ n-1 ) sin(θ n )          y -1 (x) =           arccos x1 √ x 2 n+1 +x 2 n +•••+x 2 1 arccos x2 √ x 2 n+1 +x 2 n +•••+x 2 2 . . . arccos xn-1 √ x 2 n+1 +x 2 n +x 2 n-1 2arccot xn+ √ x 2 n+1 +x 2 n xn+1          
Vector y has length n + 1 and y -1 length n. The (i, j)-th element of (n + 1) × n matrix Dy is

∂y i ∂θ j =          0 if i < j -sin(θ 1 ) . . . sin(θ i-1 ) sin(θ i )
if i = j sin(θ 1 ) . . . sin(θ j-1 ) cos(θ j ) sin(θ j+1 ) . . . sin(θ i-1 ) cos(θ i ) if j < i ≤ n sin(θ 1 ) . . . sin(θ j-1 ) cos(θ j ) sin(θ j+1 ) . . . sin(θ n ) if i = n + 1

For n = 0, we set Θ = {-1, 1} and we parametrize S 0 with an identity y : Θ -→ {-1, 1}. Let I s be a subset of W such that preimages of its points are nondegenerate spheres. We define transformation ψ :

I s × Θ -→ R h as follows ψ( w, θ) = ψ w (θ) = p(w) + R(w)H(w)y(θ)
For a fixed w, it is a bijection between Θ and M -1 (w). Mapping ψ is a bijection onto

M -1 (I s ) whose complement in R h is of measure zero if M is not invertible. ψ is continuously differentiable. (if n = 0, it is continuously differentiable for a fixed θ). If Qu = 0, then D w ψ = D p + Hy∇ R T + Ry n+1 ||b|| I - 1 ||b|| 2 bb T D b and if Qu = 0, then D w ψ = D p + U y∇ R T If n > 0, then also D θ ψ = RHDy
For n > 0, we define J (w, θ) as the Jacobian of ψ, J (w, θ) = det[D w ψ | D θ ψ] and for n = 0, as J (w, θ) = det D w ψ .

Reverse mapping

We define a subset of W ,

I f = w ∈ W ∃x ∈ M -1 (w) f (x) > 0}.
The following mapping G : I f -→ R h , allows the minimum MSE of the operation we discussed above. For w ∈ I f ∩ I s ,

G(w) = Θ ψ w (θ)f (ψ w (θ))J ( w, θ)dλ n Θ f (ψ w (θ))J ( w, θ)dλ n
and for all other w = w eq in its domain G(w) = p(w). Symbol λ n stands for Lebesgue measure on R n . Value of an integral of a vector-valued function is a vector of integrals of its component functions. For w ∈ I f ∩ I s , we can rewrite the above formula as

G(w) = p(w)+ Θ f (ψ w (θ))J ( w, θ)dλ n -1 R(w)H(w) Θ y(θ)f (ψ w (θ))J ( w, θ)dλ n
If w eq ∈ I f , then G(w eq ) can be arbitrary as this single value has no bearing on the value of MSE. G is continuous on I f \{w eq }.

For w ∈ I f ∩ I s , θ ∈ Θ and x ∈ M -1 (I s ), we define φ(w, θ) = det D w ψ( w, θ) T D w ψ( w, θ)

Φ w (x) = φ w, ψ -1 w (x) = φ w, y -1 1 R U T x -p
Integrals below are over the embedded n-sphere that is the preimage of w and is an n-dimensional manifold. If Qu = 0 and n > 0, then for w ∈ I f ∩ I s ,

G(w) = M -1 (w) xf (x)Φ w (x)dS n M -1 (w) f (x)Φ w (x)dS n
This alternative form enables approximation of this value by summing over sets of points spaced regularly and finely on the n-sphere M -1 (w) embedded in R h .

Algorithm

Based on this work we may consider the following algorithm for mapping a point from the concept space to the feature space. We are given set C ⊂ R h of centroids of c clusters and a probability density function f defined on R h or a data set D ⊂ R h . If we are given a data set, we will find values of f with density estimation methods. We compute M , M + , Q, u, a, r, d and U . Based on Qu we ascertain which case we are dealing with and we compute H and n. We determine a reduced set of clusters C. The input of the algorithm is w ∈ [0, 1] c and the output is x ∈ R h . When value of x is set, the algorithm stops. Let µ be a centroid of the cluster in which w has the highest membership. We tested above algorithm against assigning to w ∈ W a mean of M -1 (w) (which in case of FCM is p). Each data set consisted of c clouds of points generated using multivariate normal probability distributions of different means and covariance matrices. Results are in table 1.

Application

The reverse mapping from the concept space to the feature space that we derived, is best fitted -in terms of MSE -to a given probability density function defined on the feature space which can be a density of some data set in the feature space. We will discuss an application of this mapping to the problem of defuzzification of results of computation on fuzzy data. Suppose we have a training set T whose each element consists of k input points x i 1 , x i 2 , . . . , x i k of feature space R h and one output point x i 0 of feature space. All the feature space points from T put together, constitute a data set D on which we use FCM to find c fuzzy clusters and a set of their centroids C. We transform every x i j into M(x i j ), to obtain a concept-level training set T c . Using some machine learning method, we train on T c an operation F c , working on the level of concept space, that returns a point of concept space based on k input points of concept space. Let G be a minimum MSE reverse mapping associated with C, found for f -a density of data set D out ⊂ D consisting only of output points. G has minimum MSE when defuzzifying output points of the training set T c , so we may expect it to defuzzify other outputs well.

Furthermore, for all sets of input points from T c , F c is fitted to return a point close to a correct output point M(x i 0 ) and G is fitted to map to a point of feature space close to x i 0 . Thus, the composition F = G • F c • M may be a good fit for training set T as an operation that returns a point of feature space based on k input points of feature space. This approach may be useful in case of data which is easier modeled at the concept level than at the numeric level.

Conclusion

We presented a reverse mapping from concept space to feature space, best fitted -in terms of MSE -to a given probability density function defined on feature space or a given data set in feature space. Similar reverse mappings could be obtained for different transformations into the concept space. Maybe transformation that allows more precise reverse mapping could be constructed. Many aspects of this work leave room for further research. In the described application the concept-level operation may return points that do not lie in set W . Effectiveness of various approaches to amend this situation can be investigated, such as orthogonal projection onto W or projection onto W in the direction of the origin. Also there may be better ways to treat the points of W that drop out at steps 3 and 6 of the algorithm. For a point of W that has a nonempty preimage, but density f is zero on all its points, it might be a good treatment to map it onto the point of its preimage that is closest to the mean of D weighted by f . The second important direction of research is approximation of f in a situation when we do not have access to any data set in the feature space.

Fig. 1 .

 1 Fig. 1. Processing of data

Fig. 2 .

 2 Fig. 2. Preimages of points (w1, 1 -w1) for h = 2 and c = 2.

Fig. 3 . 1 CT

 31 Fig. 3. Contour plot of a simple probability density derived from a data set in a twodimensional feature space.

Fig. 4 .

 4 Fig. 4. Point x of R h is mapped by M onto a point w of W which in turn is mapped back into R h by G.

1 .

 1 Project w onto W orthogonally/in the direction of the origin. 2. Compute v, b and Qv. 3. If conditions for M -1 (w) to be nonempty are not met, set x = µ. 4. If w = w eq , compute p and R. 5. If d = 0 and Qu = 0 or if R = 0, set x = p. 6. Check if f is positive on any point of M -1 (w). If not, set x = p. 7. Compute an approximation of G(w) and set x to its value.

Table 1 .

 1 Comparison of MSE

	h	c	p	G(w)
	3	3	22.58	15.35
	4	4	25.37	19.19
	7	7	8.21	6.75
	4	3	48.65	25.95
	5	4	23.72	20.22
	6	5	25.71	23.23


Bożejko Wojciech

Grymin Rados

Jagie L Lo Szymon

Pempera Jaros

Robust tabu search algorithm for planning rail-truck intermodal freight transport

Keywords: intermodal transport, optimization algorithm, tabu search

In this paper a new efficient tabu search algorithm for assigning freight to the intermodal transport connections was developed. There were also formulated properties of the problem that can be used to design robust heuristic algorithms based on the local search methods. The quality of solutions produced by the tabu search algorithm and by often recommended greedy approach were also compared.

Introduction

Road transportation is very expensive. Because of that, carrying goods over a long distance is achieved as a combination of different types of transportation (e.g.with the use of rail, ships and trucks). Such means of transportation is called intermodal freight transport. The problems related to the intermodal transport are intensely studied by the Operations Research. Thanks to such activiites, transport infrastructure can be better designed, vehicle routes and delivery schedules can be planned more consciously which results in savings and consequently may lead to lower prices of transported goods.

A widely used definition of intermodal freight transport was introduced at the European Conference of Ministers of Transport. It was defined as "the movement of goods in one and the same loading unit or vehicle by successive modes of transport without handling of the goods themselves when changing modes". Second commonly used definition was introduced by C. Macharis and Y. Bontekoning [START_REF] Macharis | Opportunities for or in intermodal freight transport research: A review[END_REF]. They define it as "the combination of at least two modes of transport in a single transport chain, without the change of container for the goods, with most of the route travelled by rail, inland waterway or ocean-going vessel and with the shortest possible initial and final journeys by road". More definitions of intermodal freight transport can be found in [START_REF] Bontekoning | Is a new applied transportation research field emerging?--a review of intermodal rail-truck freight transport literature[END_REF].

Problems related to the intermodal transport are much more complex to solve than problems which take into consideration a fixed type of transportation (aka unimodal problems). Moreover, models used to solve unimodal problems are also applied to intermodal transport problems. For this reason, there is a group of people convinced that intermodal transportation research is emerging as a new transportation research field and it is still in the pre-paradigmatic phase and proper models are not known at the moment (see [START_REF] Bontekoning | Is a new applied transportation research field emerging?--a review of intermodal rail-truck freight transport literature[END_REF][START_REF] Macharis | Opportunities for or in intermodal freight transport research: A review[END_REF]). Y. Bontekoning, C. Macharis and J. Trip [START_REF] Bontekoning | Is a new applied transportation research field emerging?--a review of intermodal rail-truck freight transport literature[END_REF] described the actual state of intermodal transport research field and marked actions which should be done to make it similar to "normal science". They proved that the intermodal transport research is in the pre-paradigmatic phase. They found plenty of small research communities that address intermodal transport problems. They proposed transformation of these small communities to one or two large research communities.

C. Macharis and Y. Bontekoning [START_REF] Macharis | Opportunities for or in intermodal freight transport research: A review[END_REF] underlined that the intermodal transport is a very complex process that involves many decision-makers. They distinguished four types of them and emphasized that some decisions can have longterm effects (e.g. planning and building a railway infrastructure) when other can have short-term effects (temporary changes in timetable). Decision-makers should work in a close collaboration to achieve supreme results. Many decisions lead to a variety of areas where optimization can be used, whereas many authors introduced overview of articles and methods used in the intermodal transport research and classified them by the type of decision maker and by the time horizon of operations problem.

A. Caris, C. Macharis and G. Janssens [START_REF] Caris | Decision support in intermodal transport: a new research agenda[END_REF] proposed new research fields regarded to decision support in intermodal freight transport. They made an overview of applied applications that support decisions of policy makers, terminal network design, intermodal service network design, intermodal routing, ICT (Information and Communication Technologies) innovations and drayage operations. They pointed out that there is no link between models for terminal network design and intermodal service network design. They recognized the need for solution methods solving intermodal freight transport optimization problems that can accept multiple objective functions, transportation mode schedules, economies of scale and demanded times of delivery.

Intermodal freight transport, due to its big complexity and plenty of constraints imposed on solutions, constitutes a challenge for many types of heuristics and metaheuristics. A. Caris and G. Janssens [START_REF] Caris | A local search heuristic for the pre-and end-haulage of intermodal container terminals[END_REF] optimized pre-and endhaulage of intermodal container terminals using heuristic approach. They modeled problem as a Full Truckload Pickup and Delivery Problem with Time Windows (FTPDPTW) where vehicles carry full truckloads to and from an intermodal terminal. Time windows were used to represent the time interval in which the service at a customer must start. They proposed two-phase insertion heuristic and the improvement heuristic with three types of neighborhood. The solution is obtained by two-phase insertion heuristic and afterwards improved by the improvement heuristic.

There are some articles devoted to comercial decision support systems (DSS). G. Kelleher, A. El-Rhalibi and F. Arshad [START_REF] Kelleher | Scheduling for intermodal transport[END_REF] described features of PISCES, the integrated system for planning intermodal transport. They presented methods used in PISCES for dealing with the triangulation in the pick-up and drop scenario. Another research on decision support system was published by A. Rizzoli, N. Fornara and L. Gamberdella [START_REF] Rizzoli | A simulation tool for combined rail/road transport in intermodal terminals[END_REF]. They presented the terminal simulator component of the Platform project, funded by the Directorate General VII of the European Community. Presented software can model processes taking place in an intermodal road or rail terminal. It was designed on the basis of discreteevent simulation paradigm. Software user can define the structure of terminal and different input data. It allows us to check how changing terminal structure may have an influence on its performance.

Noteworthy are also researches on systems that support decision-making by more than one policymaker. A. Febraro, N. Sacco and M. Saeednia [START_REF] Di Febbraro | An agent-based framework for cooperative planning of intermodal freight transport chains[END_REF] proposed an agent-based framework for cooperative planning of intermodal freight transport chains. In this system, many actors can work together and negotiate their decisions to achieve a common goal.

Transport companies are not willing to share their tariffs on their websites. This information may be hidden due to many factors, eg. cost changes. Costs can change overnight because of fluctuations of exchange rates, political situation etc. Moreover, prices vary from one company to another. Researchers need such data to develop better algorithms. Special price models come forward. T. Hanssen, T. Mathisen, F. Jorgensen [START_REF] Hanssen | Generalized transport costs in intermodal freight transport[END_REF] proposed a generalized transport costs model that can be used to assess mean prices of different types of transport on the given distance. Intermodal transport research field demands knowledge base that will ease scientists to conduct their researches on the real data. Experimental results will be more reliable. It should be in the interest of governments and all shipping companies to support building of such database.

Problem formulation

In a planning phase, transportation management company has to realize a certain number of transport tasks. Transport task is to carry some amount of goods from suppliers to customers. Transport can be organized in two ways: (i) by a single truck, (ii) by the intermodal transport (truck-train-truck). Goods are transported in containers or semitrailers customized to rail transport.

There are known locations of customers, suppliers and intermodal terminals, distances between: (i) intermodal terminals, (ii) intermodal terminals and customers, (iii) suppliers and customers. Of course in the first case this is the length of the rail route and in the other cases it is the length of the road route. Cargo trains implementing intermodal freight transport follow the schedule of courses. Each course specifies initial and final intermodal terminal, time of delivery, unit cost of the course and the amount of free wagons. The number of free wagons is updated online based on reservations of transportation management companies. Attachment of wagons on intermediate stations is forbidden. The objective of optimization is the assignment of task to the train course simultaneously minimizing the overall costs of transport.

Let J = {1, . . . , n} be a set consisting of n transport tasks and let T = {1, . . . , t} be the set of railway courses. For each task j ∈ J and each course i ∈ T the distance achieved by the traffic transport d j,i is given. Note, that the distance d j,i is the result of summing distances between supplier and the initial intermodal terminal and between final intermodal terminal and the customer. The distance achieved by the traffic transport between supplier and customer for the fixed j ∈ J is marked as d j,0 . The railway distance for the course i ∈ T is r i . The course i ∈ T has l i ≥ 0 free cargo wagons to load. The overall price for the carriage of freight from the supplier to the customer using cargo train i ∈ T is c j,i . The cost of direct road transport from the supplier to the customer specified in the task j is c j,0 . Let the assignment of the course to the task j be marked as a j , a j ∈ {0} ∪ T (a j = 0 if the transport is carried only by the truck). Vector α = [a j ] denotes the assignment of all tasks to the courses. The total cost of transport for assignment α is

Cost(α) = n j=1 c j,aj . (1)
We would like to find such assignment α * that the total cost of transport is as small as possible

Cost(α * ) = min α∈Λ Cost(α), (2)
where Λ is the set of all possible assignments, |Λ| = n m+1 .

Properties of the problem

In the current subsection we formulate certain properties of the problem, which can be used in the design of efficient heuristic algorithms based on local search methods. The first one relates to the method that allows us to determine the lower bound of the objective function value for the optimal solution, the second one allows us to reduce the number of solutions in the neighborhood by eliminating the subset with worse solutions.

Proposition 1. Let u i = min k=0,...,m c j,k , then the lower bound

LB = n j=1 u j . (3)
Proposition 2. Let α and β be the assignments of tasks to the trains such that Cost(β) < Cost(α). Then, at least one task has lower cost of transporting in the assignment β than in the assignment α.

Proposition 3. Let β be the assignment of tasks to the trains resulting from α by the assignment of the task j to train i, then

Cost(β) = Cost(α) -c j,aj + c j,i . (4)
Proposition 4. Let β be the assignment of tasks to the trains resulting from α by interchanging assignment of tasks j and k, then

Cost(β) = Cost(α) -c j,aj -c k,a k + c j,a k + c k,aj . (5)
Note, if Cost(α) is known, expressions (4) and ( 5) can be determined in time O(1).

Example

A certain logistic company has to realize n = 5 transportation tasks. Transportation may be achieved with t = 3 courses of cargo train. All trains have l i = 2 free cargo wagons to load. The transportation costs (traffic and intermodal) are given in Table 1.

It is easy to see that the price of road transport is 9623 whereas, the lowest price of transport is LB = 5472 (marked with bold). Transport with the lowest price requires use of 5 wagons of the course 3. It is not possible because rail connection 3 has only two free wagons. In Table 1 a feasible connection (that takes into account number of free wagons) with total price 5893 was marked in bold.

Let us consider a greedy strategy that assigns tasks to the cheapest intermodal freight transport connection. In the first row, tasks with the biggest difference in price between road and rail transport will be assigned. Tasks that cannot be achieved by the rail transport with the cheapest cost will be achieved by truck transport. The described strategy is used in many logistic companies, where assignment of transport tasks is realized with the use of forwarding agents. A forwarding agent is concentrated on finding the cheapest solution of transport problem, because his salary depends (directly or indirectly) on the income that is the difference between the price of transport negotiated with a customer and the real price.

Let us assume that road transport prices were negotiated with customers, profits from the intermodal freight transport are formed respectively: 905, 707, 739, 875, 925. Tasks 1 and 5 generate the biggest profits, therefore they will be realized by the intermodal freight transport. The rest of tasks will be realized by trucks.

An approximation algorithm

In order to solve the stated problem we propose a local search algorithm based on the tabu search (TS) approach [START_REF] Glover | Tabu search part i[END_REF][START_REF] Glover | Tabu search part i[END_REF]. The tabu search is one of the best methods of constructing heuristic algorithms. This is confirmed for many optimizing problems, as a main method (for scheduling of tasks [START_REF] Bożejko | Parallel tabu search algorithm for the hybrid flow shop problem[END_REF][START_REF] Nowicki | A fast taboo search algorithm for the job shop problem[END_REF][START_REF] Bożejko | A fast hybrid tabu search algorithm for the no-wait job shop problem[END_REF][START_REF] Bożejko | Parallel genetic algorithm for the flow shop scheduling problem[END_REF], vehicle routing [START_REF] Renaud | A tabu search heuristic for the multi-depot vehicle routing problem[END_REF][START_REF] Bożejko | Solving permutational routing problems by populationbased metaheuristics[END_REF], packing [START_REF] Leung | Extended guided tabu search and a new packing algorithm for the two-dimensional loading vehicle routing problem[END_REF], container loading [START_REF] Bortfeldt | A parallel tabu search algorithm for solving the container loading problem[END_REF]), as well as a key element of higher metaheuristics, e.g. golf method [START_REF] Bożejko | The new golf neighborhood for the flexible job shop problem[END_REF]. Neighborhood determination in the tabu search metaheuristics is also frequently parallelized (see [START_REF] Bożejko | Parallel path relinking method for the single machine total weighted tardiness problem with sequence-dependent setups[END_REF][START_REF] Bożejko | On single-walk parallelization of the job shop problem solving algorithms[END_REF]).

An algorithm based on this method, in each iteration searches neighborhood of the basic solution for the solution with the best objective function value. In every iteration the best solution replaces the basic solution. To prevent searching loops, the tabu mechanism is used. It is usually implemented as a list with the limited length. In each iteration, selected attributes of subsequently visited solutions are stored. Contents of the list divide the neighborhood into two subsets: a set of forbidden and a set of feasible solutions. Forbidden solutions are not searched except the case when forbidden solution is better than the best solution found so far.

The search stops when the given number of iterations without improving the criterion value has been reached or the algorithm has performed a given number of iterations.

Moves and neighborhood

The neighborhood of a solution is generated by moves. In our problem, the solution is represented by a vector of task assignments α. The neighborhood of the solution can be created due to exchanges of assignment of two tasks or change of assignment of a single task.

Let v = (a, b) be a pair of exchanged tasks. We define a new assignment α (v) obtained from α = (α 1 , . . . , α n ) by exchanging assignment of tasks a and b as follows:

α (v) = (α 1 , . . . , α a-1 , α b , α a+1 , . . . , α b-1 , α a , α b+1 , . . . , α n ) f or a < b (6)
and

α (v) = (α 1 , . . . , α b-1 , α a , α b+1 , . . . , α a-1 , α b , α a+1 , . . . , α n ) f or a > b. (7)
Let EX be a set of some such moves and N (EX, α) = {α (v) : v ∈ EX} be a neighborhood of solution α generated by a move set EX. For a feasible solution α every move v ∈ EX generates a feasible solution.

Let v = (a, k) be the move that changes the assignment of the task a to the course k. We define the new assignment α v obtained from α = (α 1 , . . . , α n ) by execution move v in α as follows: We propose a reduction of the neighborhood size to the set of promising moves. The move v is promising if its execution gives the chance to receive a better solution. From Proposition 2 we have simple conditions for obtaining a better solution:

α (v) = (α 1 , . . . , α a-1 , k, α a+1 , . . . , α n ). ( 8
-c a,αa > c a,k for move v = (a, k) ∈ IN S, -c a,αa > c a,α b for move v = (a, b) ∈ EX, -c b,α b > c b,αa for move v = (a, b) ∈ EX.
We will mark the reduced set of moves as V .

Computation results

The main objective of experimental studies was to evaluate the usefulness of advanced heuristics in assigning transportation tasks to intermodal transport. Experimental test was carried out on the randomly generated data. The set of 120 instances is divided into 12 groups. Each group consists of 10 instances with the same number of tasks n and freight trains t. The study was conducted for groups, where the number of tasks n ∈ {50, 100, 200, 500} and the number of cargo trains t ∈ {10, 20, 30}.

Railway distances r i for course i ∈ T were generated from the uniform distribution on [2000,2500], in intermodal transport the distances achieved by traffic transport d j,i were generated from uniform distribution on [30,200]. A traffic transport distance is usually shorter than distance of intermodal transport thus we determine this distance from the expression d j,0 = min i=1,...,t (r i -d j,i ).

We assumed unit costs of traffic transport c j,0 = d j,0 . The cost of task j ∈ J carried by freight train i ∈ T includes the cost of road transport (d j,i ), the cost of handling in intermodal terminal (h j ), the cost of transporting freight train depending on the distance (r j ) and is expressed by the formula c j,i = d j,i + h j + γr j , where γ is the factor of the cost of rail transport to the cost of traffic transport. The research was carried out for three values of factor γ (0.8, 0.65, 0.5) and h j = 60. Note that for γ = 0.8 and transport distance 300 the costs of traffic and intermodal transport are comparable.

The number of free wagons was the same for each cargo train. We considered two levels of wagon availability for loading: (i) a few free wagons: l i = n/t + 1, (ii) many free wagons:

l i = 2 • n/t .
The algorithm TS was implemented with the reduced neighborhood V , written in C++ and ran on Lenovo T540p personal computer with processor i7-4710 2.5 GHz. Further, we wrote a greedy algorithm G (see subsection 2.2 for details) and an algorithm R which compute the total cost of traffic transportation. The algorithm TS performed 1000 iterations and started from the solution in which all tasks were assigned to road transport.

Since there are no algorithms for solving the considered problem in the literature, we made a comparison of TS, G and R with lower bound LB (3). For each instance, we defined the following values:

Cost(α A ) -the total cost of transportation of tasks forms set J found by the algorithm A, A ∈ {T S, G, R}, P RD(A) -the mean value of the relative cost of solution found by algorithm A with respect to the lower bound LB i.e.

P RD(A) = Cost(α A ) -LB LB • 100[%], (9)
CPU -the mean computation time (in seconds).

The results of computer computations are summarized in Table 2. The first column contains the number of tasks and freight trains in each instance of the group, the second contains the average relative cost of traffic transport, the next three columns refer to instances of the small number of available wagons and include: the average relative cost of the solution generated by greedy algorithm R, the average relative cost of the solution generated by tabu search algorithm TS and the average number of transports carried out exclusively by road in the solution generated by the algorithm TS. The other three columns refer to instances with many free cargo wagons. The table shows the results for different γ values.

At the beginning of the analysis of the results collected in Table 2, it should be noted that the proposed tabu search algorithm successfully finds the task assignments for intermodal transport. The solutions generated by TS algorithm for the intermodal transport with a limited number of free cargo wagons are only a few percent worse than transport with minimum cost i.e. with the unlimited number of free cargo wagons and trucks. It is easy to notice that the algorithm TS finds significantly better solutions for instance with a large number of wagons to be loaded.

According to Table 2, tabu search heuristic performs significantly better than the greedy heuristic. The average relative cost does not exceed 3.2% for instances with large number of free wagons and 8.3% for instances with small number of free wagons. In the case of greedy algorithm this cost varies accordingly from 3.2% to 38.3% and from 6.0% to 52.6%. The superiority of the algorithm TS over the greedy algorithm R increases with decreasing γ (with increasing attractiveness of intermodal transport).

While comparing the cost of road transport and intermodal transport, it can be noted that with decreasing γ values increases the difference between the cost of road transport and the intermodal one. For γ = 0.5 it is close to 70%. The experiment shows that for the highest γ value, the profit of using intermodal transport is admittedly less (approximately 10%) however, in our opinion, it is important from the standpoint of business activity. In addition, with decreasing Table 3 shows the average computational time for 9 groups of instances (for groups with n = 50 the computation time was less than 0.1). The calculations were performed for two versions of the algorithm TS: T S(V ) with reduced neighborhood V and T S(U ) with the full neighborhood U . It is easy to see that the computation time increases with the increasing number of tasks n and the number of cargo trains t. While comparing computation time T S(V ) and T S(U ), it should be highlighted that T S(V ) runs faster than T S(U ) from 30% to 50%. The computation time of algorithm T S(V ) does not exceed 6 seconds for instance with the biggest number of tasks and cargo trains.

Conclusion

In this paper, we developed the tabu search algorithm to minimize the total cost of transport. We have considered the intermodal transport as an alternative to the most used traffic transport. We have formulated several properties of the problem, which were used not only to increase the efficiency but also to reduce the computation time of TS algorithm. We experimentally proved that the global cost optimization of intermodal transport allows us to achieve significantly higher profits than using in practice of the greedy approach.

The results obtained and our research experience [START_REF] Smutnicki | A new approach for multi-criteria scheduling[END_REF] encourage us to extend the ideas proposed to multi-criteria problems generated by intermodal transport.

  )Let IN S be a set of some such moves and N (IN S, α) = {α(v) : v ∈ IN S} be a neighborhood of solution α generated by a move set INS. For the feasible solution α, there are infeasible moves v = (a, k) ∈ IN S such that α a = k and cargo train k ∈ T has not free wagons.LetU = {v : α (v) is f easible ∧ v ∈ EX ∪ INS} be a set of feasible moves and N (U, α) be a neighborhood of solution α. The neighborhood N (EX, α) has n 2 neighbors, while N (IN S, α) has nm neighbors. The determination of the total cost explicitly from the formula (1) requires O(n) time for each solution, therefore the search of neighborhood can require a great computational effort. The use of expression (4) for moves from IN S and expression (5) for moves from EX reduces computation time to O(1) (n times) for each generated solution.

Table 1 .

 1 Costs of intermodal and traffic transport

		Task 1	Task 2	Task 3	Task 4	Task 5
	0-traffic	1962	1863	1879	1947	1972
	1-course	1323	1276	1375	1289	1410
	2-course	1328	1283	1304	1327	1184
	3-course	1057	1156	1140	1072	1047

Table 2 .

 2 Relative cost of traffic and intermodal costs , the percentage share of road transport in the solutions generated by the algorithm TS is reduced. For γ = 0.8 it varies from 25% to 45% and from 2.4% to 15.4%, for the remaining values of coefficient γ it equals 0.

		G	a few free wagons	many free wagons
	n × t	traffic	Greedy TS use of trucks	Greedy TS use of trucks
				γ = 0.80
	50×10	10.5	6.0	3.7	25.0	3.2	1.5	8.0
	50×20	10.5	7.4	5.3	37.0	5.4	2.5	8.8
	50×30	10.5	7.7	5.4	32.4	5.8	2.3	2.4
	100×10	10.7	6.4	4.3	32.0	3.6	1.3	3.5
	100×20	10.5	7.8	5.6	38.8	5.9	3.0	10.4
	100×30	10.4	7.7	5.5	37.8	5.9	2.7	7.3
	200×10	10.6	6.3	4.4	33.4	3.7	1.4	6.7
	200×20	10.7	7.5	5.1	33.7	5.4	2.4	4.7
	200×30	10.4	7.9	6.0	45.1	6.1	3.2	12.8
	500×10	10.4	6.9	5.2	42.7	4.6	2.5	15.4
	500×20	10.6	7.5	5.7	42.4	5.3	2.7	12.5
	500×30	10.7	8.1	5.9	41.6	6.2	3.0	7.8
				γ = 0.65
	50×10	34.7	19.5	4.5	0.0	10.4	1.5	0.0
	50×20	34.8	24.0	6.7	0.0	17.0	2.4	0.0
	50×30	34.9	25.6	6.5	0.0	19.1	2.1	0.0
	100×10	35.0	20.5	5.2	0.0	11.6	1.1	0.0
	100×20	34.8	25.1	7.2	0.0	18.9	2.8	0.0
	100×30	34.8	24.7	6.9	0.0	18.8	2.5	0.0
	200×10	34.8	20.5	5.7	0.0	11.3	1.3	0.0
	200×20	35.0	24.2	6.1	0.0	17.0	2.1	0.0
	200×30	34.8	26.1	8.3	0.0	19.9	3.2	0.0
	500×10	34.4	22.9	6.8	0.0	15.5	2.5	0.0
	500×20	34.9	24.4	7.6	0.0	16.8	2.7	0.0
	500×30	35.1	26.3	8.2	0.0	19.9	2.8	0.0
				γ = 0.50
	50×10	72.4	38.1	3.9	0.0	18.7	1.1	0.0
	50×20	72.8	47.5	6.2	0.0	32.4	2.1	0.0
	50×30	73.3	51.1	6.0	0.0	36.6	1.8	0.0
	100×10	73.0	39.8	4.7	0.0	20.0	0.9	0.0
	100×20	72.8	49.2	6.6	0.0	35.4	2.4	0.0
	100×30	73.1	49.1	6.4	0.0	35.0	2.2	0.0
	200×10	72.8	40.7	5.2	0.0	20.7	1.0	0.0
	200×20	73.3	47.7	5.6	0.0	31.8	1.8	0.0
	200×30	73.2	53.2	7.8	0.0	38.8	2.9	0.0
	500×10	71.9	44.9	6.0	0.0	28.5	1.9	0.0
	500×20	73.2	49.0	7.1	0.0	31.7	2.3	0.0
	500×30	73.6	52.6	7.7	0.0	38.3	2.5	0.0

γ

Table 3 .

 3 Computation time in seconds

		a few free wagons	many free wagons
	n × t	T S(V )	T S(U )	T S(V )	T S(U )
	100×10	0.1	0.1	0.0	0.1
	100×20	0.1	0.2	0.1	0.2
	100×30	0.2	0.3	0.1	0.4
	200×10	0.3	0.5	0.2	0.5
	200×20	0.4	0.8	0.3	1.0
	200×30	0.7	1.2	0.4	1.2
	500×10	1.9	3.2	1.4	3.2
	500×20	3.9	6.2	2.9	6.0
	500×30	5.8	9.3	4.6	9.7