Machine learning is constantly gaining popularity in real life applications. And one of them is prediction of various real-life events that depend on a huge number of factors that are hard to evaluate. In this article we describe the process of applying XGBoost - one of supervised machine learning methods - to help in prediction and localization of accidents in the district heating network of Kaunas region. We also investigate the importance of the diferent factors for these events.
cided to create a tool for the heating network supervising personal. And supervised machine learning alIn everyday operation centralised district heating com- gorithms seemed as an inexpensive and valid option pany “Kauno energija” is supervising more than 900 to help in the process of pipe breakage localization. kilometers of district heating networks that provide heating and hot water to 118 891 customers (as of the end of 2017) in Kaunas region. 2. Related works
Every year pipe breakages in the district heating network occur. In most of these times the district heat- We were able to find some similar studies where pipe ing services must be stopped for the customers. And break accidents were predicted using machine learndue to the aged infrastructure, it is dificult to deter- ing methods. In article [1] pipe breaks were predicted mine where did the accident happen. The only sign of for water distribution network using pipe attributes accident is often a critical pressure drop or a frequent and climatic data, since a monitoring network is not refill of the heating water in the system. There are a available [2]. The goal of the authors was to find pipes lot of cases when repair teams are excavating the area that can break soon to prioritize pipe replacements and but do not find the accident and sometimes small acci- repairs. Also diferent models are tested for best perdents cannot be found and are compensated by system formance: RankBoost.B, Cox proportional hazard morefill. del, Naive Bayes, Logistic Regression and Artificial Ne
When an accident happens any information that wo- ural Network. The provided results show that Rankuld help to determine its location is helpful. Due to the Boost.B is the most successful with AUC score of more amount and complexity of factors that cause accidents that 0.85. In article [3] an ensemble of models are used it is dificult to predict them. There are some complex to predict water utility pipeline condition. As input solutions with thermodynamic and hydro-mechanics the authors use physical pipe attributes, environmenmodels in the market, which allow to calculate pipe tal data, and operational factors and data obtained from breakage, but they are hard to use in everyday work physical models which are developed to understand and require a lot of investments and learning eforts. the physical process of pipe deterioration [4, 5, 6].
Also, a Web application for network accident man- Another article [7] describes pipe failure modelling agement TAVSIS was developed and it seemed like a for water distribution networks using boosted decigood idea to integrate accident localization algorithms sion trees. To predict pipe failures authors use Adwithin that system. All things considered, it was de- aBoost, RUSBoost, Random Forest, and Decision Tree models.
Although the mentioned articles describe similar meIKVTUUSS2a0n2t0a:kInafVoramlleayt,ioKnaSuoncaise,tyLiathnudaUnniaiversity Studies, 23 April 2020, thodology there are some major diferences to our ap" mantas.bukauskas@ktu.edu (M. Bukauskas); proach. First of all, we are targeting a district heating mantas.lukosevicius@ktu.edu (M. Lukoševičius) network. Secondly, we are using historic weather data as one of the inputs. And finally, we have a diferent goal - to locate pipe segment that have failed rather © 2020 Copyright for this paper by its authors. Use permitted under Creative CPWrEooUrckReshdoinpgs IhStpN:/c1e6u1r3-w-0s.o7r3g CCoEmUmoRns WLiceonrsekAsthtriobuptioPnr4o.0cIneteerdnaitniognasl ((CCC EBYU4R.0)-.WS.org) than analyse which pipes are most likely to fail in future.
As for a district heating network we were unable to ifnd any related work. In article [8] authors use machine learning approach to detect faults by analysing temperature readings and some additional data from district heating substations. In article [9] authors use a completely diferent methodology - a deterministic – probabilistic structural integrity analysis to predict pipeline lifetime and probability of failure.
Our approach is more similar to real-time trafic accident localization approach described in articles: trafifc accident prediction in the state of Utah (USA) [10]; predicting trafic accidents through heterogeneous urban data [11].
In this study we analyze pipe breakage accidents (see Fig. 1) that happened from the January 2013 to September 2019. There were 1 466 accidents in total that occurred during this time. We extracted these attributes: date and time of the accident occurrence (weekday, month, hour), geographic location. The accident data are collected using a GIS system by the company dispatchers that are supervising the district heating operations 24 hours a day in shifts. When an accident occurs a dispatcher marks its location on the pipe segment and this allowed us to use a spatial intersection with the pipe segments to determine how many accidents have occurred in each diferent segment.
There were more than 49 200 of pipe network segments with a total length of 904.17 kilometers that are used to provide district heating network services to customers • pipe line type (distribution in the data is displayed in Fig. 4); • length of the pipe; • pipe diameter.
Information about the weather conditions was extracted from archives on “Reliable prognosis” website[12]. We used one of Kaunas city weather stations. From the data set we used these attributes: • temperature; • atmospheric pressure; • humidity; • raining fact.
There were 19 724 records collected at 3-hour intervals during the analysed period January 2013 - October 2019. Data interpolation had to be used to get hourly records. Also, the raining fact was extracted from a human readable message rather than a numeric value.
The process of our study was: 1. Find and collect necessary data from the data sources; 2. Analyse, identify, and extract useful data; 3. Select machine learning algorithms that would best fit the case; 4. Prepare training data for the model; 5. Split the prepared data set into training, validation and testing data sets; 6. Set machine learning model parameters; 7. Execute the learning process; 8. Validate the results; 9. Repeat steps 6 – 8 until expected results are reached; 10. Export the prepared model and test it on the testing data; 11. Deploy the prepared model to production.
Probably the most complicated part of this study was data preparation. We had to connect diferent sources of information to one — training data set which can be used by the XGBoost machine learning framework. All the processing and data analysis was performed using these tools: • ArcGIS Pro software was used to manage geo
graphic data and perform initial analysis; • Jupyter Notebook software was used to develop
and run all the process and share the results; • Scikit-learn machine learning framework was u
sed to prepare training data; • XGBoost gradient boosting framework was used
for model preparation.
There were some attributes that were dropped as they were considered unimportant. These were attributes with none or very few values, or irrelevant ifelds: who edited the data, when the last edit was performed, etc. The quality of the remaining data was not perfect either, as some relevant attributes were missing. To mitigate the problem, we prepared the training data by using diferent techniques: • Data interpolation to increase the frequency of
the weather data to one-hour intervals; • Spatial intersection to connect accidents to pipe segments, to find the missing pipe attributes from the connected segments; • Calculation of mean values to fill the missing
continuous values; • Calculation of the most frequent values to fill the
Finally, after the training data was prepared, it was missing categorical values; • One-hot encoding to transform the categorical values to numeric as the model cannot handle non-numeric values. This method creates a bisplit into two parts: 70 % for training and 30 % for testing. To make sure that positive (accident happened) and negative (accident did not happened) samples are distributed in equal rates for both data sets. For this we had to set the stratify parameter in the scikit-learn nary column for each category and returns a spar- library for the output column indicating if the accident se matrix or a dense array (depending on the sparse parameter). This encoding is needed for feeding categorical data to many scikit-learn estimators, notably linear models and SVMs with the standard kernels[13].
Also we added some additional properties: month, day of week, hour of day, was it raining during the accident or not, the total count of accidents in the segment. We rounded time of accidents to hours to the lower side dropping any remaining minutes or seconds as there is always some delay between when the accident really happens and when it is noticed.
When all the data cleaning was finished, we joined the data set of pipe segments to the data set of accidents. We performed this operation by using spatial join method provided by Arcpy library which allowed us to join accident record to the nearest pipe segment within less than 10 meter distance. Next, due to data imbalance as described in section 4.5 we generated 5 times as many negative samples as we have accident records.
We used a negative sample selection technique [10]: • Randomly select an accident record from the pos- on gradient boosting algorithm called XGBoost. It was • Randomly alter: the pipe segment, the hour of popularity, won numerous Kaggle challenges, and is happened or not. 4.2. Model As we mentioned before, our goal is to find most vulnerable pipe segments. We had to create a model which would be able to predict the probability of an accident in all piping segments at given situation. And by sorting these probabilities from highest to lowest we would be able to provide district heating network dispatcher with information which pipe segments are most vulnerable to having accident with current conditions: time properties (month, day of week, hour), weather conditions (temperature, humidity, atmospheric pressure, raining conditions) and pipe segment properties. By identifying the vulnerable pipe segments after the occurrence of the accident, the dispatcher can send repair teams to investigate them. It is important to mention that we do not really care how high or low the probability is in absolute value, because we already know that an accident has happened but we do not know where exactly.
developed as a research project at the University of Washington. Since its release in 2016 it quickly gained used in real-life applications. It is available as Open Source project and is actively developed by a community of data scientists. The XGBoost algorithm is based on gradient tree boosting model with additional reguweights to avoid over-fitting [14]. The regularized objective, or loss function, can be described as (Φ) = ∑ ( ̂ , ) + ∑Ω( ),
where is a diferentiable convex loss function that measures the diference between the prediction ̂ and the target . The second term Ω penalizes the complexity of the model. Intuitively, the regularized objective will tend to select a model employing simple and predictive functions. This loss function can be integrated into the split criterion of decision trees leading to a pre-pruning strategy.
(1) itive examples; the day, or the day of the year; • If the new sample is not within the accident records, add it to the list of negative samples; samples as positive.
This allowed us to work with a relatively low amount of data (in total 6 619 pipe segment in time records).
After that we connected the weather data to all these records. The final data set contained these attributes: count of accidents in segment, segment length, year of segment installation, segment diameter, weather temperature, weather humidity, raining state, atmospheric pressure, segment material, type of segment insulation, type of segment installation, subtype of the segment, hour of sample, weekday of sample, month of sample.
• Repeat until we have 5 times as many negative larization term which helps to smooth the final learnt
Furthermore, randomization techniques are also im- 4.3. Training plemented in XGBoost both to reduce over-fitting and to increase the training speed.
There are multiple parameters that were tuned to get the best results (as described in XGBoost documentation[15]): When the data was prepared, machine training was an easy step. After few try-outs we were able to get a model with 88.88 % of AUC (area under the curve) rating. It took only 13 epochs to reach this value, but the training kept going 50 epochs to determine that • max_depth – the maximum depth of the tree. In- the value does not further increase. creasing this value will make the model more complex and more likely to over-fit.
4.4. Evaluation • min_child_weight – the minimum sum of instan- One of the advantages of decision tree-based models ce weight (Hessian) needed in a child. If the is that it allows us to trace the key factors of how the tree partition step results in a leaf node with the decision was made. There are multiple characterissum of instance weight less than the min_child- tics that determine how well the model performs. The _weight, then the building process will give up most important of them are: further partitioning. In a linear regression task, this simply corresponds to minimum number of instances needed to be in each node. The larger min_child_weight is, the more conservative the algorithm will be.
Accuracy Accuracy explicitly takes into account the classification of negatives, and is expressible both as a weighted average of Precision and Inverse Precision and as a weighted average of Recall and Inverse Recall [16, 17, 18]. For us it shows the rate of correctly predicted result (accident or no accident):
TP + TN Accuracy = TP + TN + FP + FN , (2) • eval_metric – type of evaluation metrics for validation data, a default metric will be assigned according to objective (RMSE for regression, and error for classification, mean average precision for ranking). • objective – learning task and the corresponding
learning objective. • eta - step size shrinkage used in update to pre- TN = True Negatives; vents over-fitting. After each boosting step, we can directly get the weights of the new features, FP = False Positives; and eta shrinks the feature weights to make the FN = False Negatives. boosting process more conservative.
where:
TP = True Positives; • max_depth = 6, • min_child_weight = 5.0, • eval_metric = ’auc’, • objective = ’binary:logistic’, • eta = 0.5.
• early_stopping_rounds – activates early stopping. Precision It denotes the proportion of predicted posValidation metric needs to improve at least once itive cases that are correctly real positives. This is what in every early_stopping_rounds round(s) to con- Machine Learning, Data Mining and Information Retinue training. trieval focus on, but it is totally ignored in ROC analysis. It can however analogously be called True Positive • num_boost_round – number of boosting itera- Accuracy, being a measure of accuracy of predicted tions. positives in contrast with the rate of discovery of real positives [17]. It shows how many accidents were predicted correctly compared to all predicted accidents
TP Precision = TP + FP .
(3) Recall It is the proportion of Real Positive cases that are correctly Predicted Positive. Recall has been shown to have a major weight in predicting the success of word alignment. In a medical context Recall is moreover regarded as primary, as the aim is to identify all
Real Positive cases, and it is also one of the legs on which ROC analysis stands. In this context it is referred to as True Positive Rate [17]. For us it shows how many accidents were predicted correctly from all the occured accidents.
TP Recall = TP + FN
(4) ROC curve (Receiver Operating Characteristic curve) It is a graph showing the performance of a classification model at all classification thresholds. The curve plots two parameters [17]: • True Positive Rate; • False Positive Rate. A ROC curve plots true positive rate vs. false positive rate at diferent classification thresholds. Lowering the classification threshold classifies more items as positive, thus increasing both False Positives and True Positives. A perfect classifier will score in the top left-hand corner (False Positive Rate=0, True Positive Rate=100 %). A worst-case classifier will score in the bottom right hand corner (False Positive Rate=100 %, True Positive Rate=0). A random classifier would be expected to score somewhere along the positive diagonal (True Positive Rate = False Positive Rate) since the model will throw up positive and negative examples at the same rate [17].
Missing or incomplete data After the initial run only 30 % of positive samples was valid for processing as data was incomplete. Important features were missing such as pipe diameter, year of installation, material, etc. And that was a problem because we already had data imbalance issue with only 1 466 positive samples (pipe segments with registered accidents).
One of the methods we used to calculate some missing values was spatial intersection. We have intersected pipes with each other and copied values from AUC (area under the curve) The area under such connected pipes assuming that connected pipes have a multipoint curve is thus of some value, but the op- the same parameters. Also we filled missing values timum in practice is the area under the simple trape- by replacing them with most frequent values for catezoid [17]. As shown in Fig. 5, the main diagonal rep- gorical values and mean values for continuous values. resents chance with parallel isocost lines represent- Finally, 1 104 out of 1 466 of registered accidents were ing equal cost-performance. Points above the diagonal successfully used in the model. represent performance better than chance and those below - worse than chance. For a single good (dotted green) system, AUC is area under the curve (trapezoid between the green line and x = [0,1]). The perverse (dashed red) system shown is the same (good) system with class labels reversed [17].
Data imbalance Every year has at least 8 760 hours and every hour we have more than 49 000 of pipe segments that would make more than 429 million records every year. And we have only 1 466 accidents registered since 2013 January to 2019 October. If we would use data with this rate of positive and negative samples, the model would not be able to predict any of accidents. As by predicting that accidents will not happen at all it would be almost always right. To mitigate (5) the problem we used negative sample selection technique as described in Section 4.1.
F1-score F-measure is defined as a harmonic mean of precision and recall [19] 1 = 2 +
.
During this research we met some issues that had to be overcome. Some of them are described below.
Non-linear factors of the accidents One of the problems of real-life event prediction is that there are many unknown and immeasurable factors that cause them to occur. There are many factors that might be important, but we cannot determine them, or it would require a lot of efort to collect and provide them to the model. Such factors can be faults during the pipe production, transportation or installation, human errors, pipe environment conditions like soil, humidity, temperature changes, materials used in pipe production, count electrical conductivity, chemical environment, etc.
• Test Precision: 44.92 % • Test Recall: 80.06 % • Test Accuracy: 80.31 % XGBoost library has a method to plot the most impor- • Test F1: 57.55 % tant features. Importance is a score that indicates how useful or valuable each feature was in the construc- But for our approach we care only about the hightion of the boosted decision trees within the model. est probability with given parameters. Because we alThe more an attribute is used to make key decisions ready know that an accident happened but we just do with decision trees, the higher its relative importance. not know where exactly. That makes even relatively As shown in Fig. 6, most important features for dis- low probabilities valuable to us as it is additional infortrict heating network accident prediction are weather mation that can help us to find the accident locations. pressure P, temperature T, humidity U, segment length And our user - piping network dispatcher can use this segment_length, diameter SalSkersmuo and year probability to decide where it most likely have hapof installation PaklojimoMetai. pened even if probability is relatively low. To compare
Also as shown in Fig. 7, according to ROC curve our our model performance we use AUC score. For XGmodel performs well compared to random guess. Of Boost we get score of 0.868 as linear regression model course, we have to keep in mind that these param- reaches AUC score of 0.857. It can also be seen in ROC eters are provided for our data set where we limited curve (see Fig. 7). our positive and negative sample rate to mitigate data Also, we plotted some histograms with the most imimbalance. XGBoost classifier predicts probability of portant features to see if they have any connection to accidents between 0 and 1. By default, it has a thresh- pipe breakage. old set to 0.5, meaning that probability higher than 0.5 As shown in Fig. 8, more accidents tend to happen will yield positive result and less than 0.5 - negative. when atmospheric pressure is between 755 and 765
With probability threshold set to 0.16 we get these mm Hg. results: As shown in Fig. 9, accidents are more likely to happen during warm temperatures when heating service [1] R. Wang, W. Dong, Y. Wang, K. Tang, X. Yao, Pipe failure prediction: A data mining method, in: 2013 IEEE 29th International Conference on Data
Engineering (ICDE), IEEE, 2013, pp. 1208–1218. [2] R. Giuliano, F. Mazzenga, A. Vizzarri, Satellitebased capillary 5g-mmtc networks for environmental applications, IEEE Aerospace and Electronic Systems Magazine 34 (2019) 40–48. [14] T. Chen, C. Guestrin, Xgboost: A scalable in: 2015 IEEE Symposium Series on Computatree boosting system, in: Proceedings of the tional Intelligence, IEEE, 2015, pp. 1476–1483. 22nd ACM SIGKDD International Conference on [17] D. Powers, Ailab, Evaluation: From precision, reKnowledge Discovery and Data Mining, KDD call and f-measure to roc, informedness, marked’16, Association for Computing Machinery, New ness and correlation, J. Mach. Learn. Technol 2 York, NY, USA, 2016, p. 785–794. URL: https:// (2011) 2229–3981. doi:10.9735/2229-3981. doi.org/10.1145/2939672.2939785. doi:10.1145/ [18] F. Beritelli, G. Capizzi, G. Lo Sciuto, C. Napoli, 2939672.2939785. M. Woźniak, A novel training method to preserve [15] Xgboost parameters, last accessed May 2020. generalization of rbpnn classifiers applied to ecg URL: https://xgboost.readthedocs.io/en/latest/ signals diagnosis, Neural Networks 108 (2018) parameter.html. 331–338. [16] C. Napoli, E. Tramontana, M. Wozniak, En- [19] Y. Sasaki, The truth of the f-measure, Teach Tutor hancing environmental surveillance against or- Mater (2007). ganised crime with radial basis neural networks,