Zbornik 21. mednarodne multikonference INFORMACIJSKA DRUŽBA - IS 2018 Zvezek D Proceedings of the 21st International Multiconference INFORMATION SOCIETY - IS 2018 Volume D Mednarodna konferenca o visokozmogljivi optimizaciji v industriji, HPOI 2018 International Conference on High-Performance Optimization in Industry, HPOI 2018 Uredila / Edited by Bogdan Filipič, Thomas Bartz-Beielstein http://is.ijs.si 8.–12. oktober 2018 / 8–12 October 2018 Ljubljana, Slovenia Zbornik 21. mednarodne multikonference INFORMACIJSKA DRUŽBA – IS 2018 Zvezek D Proceedings of the 21st International Multiconference INFORMATION SOCIETY – IS 2018 Volume D Mednarodna konferenca o visokozmogljivi optimizaciji v industriji, HPOI 2018 International Conference on High-Performance Optimization in Industry, HPOI 2018 Uredila / Edited by Bogdan Filipič, Thomas Bartz-Beielstein http://is.ijs.si 8.–12. oktober 2018 / 8–12 October 2018 Ljubljana, Slovenia Urednika: Bogdan Filipič Odsek za inteligentne sisteme Institut »Jožef Stefan«, Ljubljana Thomas Bartz-Beielstein Institute for Data Science, Engineering and Analytics, TH Köln Založnik: Institut »Jožef Stefan«, Ljubljana Priprava zbornika: Mitja Lasič, Vesna Lasič, Lana Zemljak Oblikovanje naslovnice: Vesna Lasič Dostop do e-publikacije: http://library.ijs.si/Stacks/Proceedings/InformationSociety Ljubljana, oktober 2018 Informacijska družba ISSN 2630-371X Kataložni zapis o publikaciji (CIP) pripravili v Narodni in univerzitetni knjižnici v Ljubljani COBISS.SI-ID=31870503 ISBN 978-961-264-138-2 (pdf) PREDGOVOR MULTIKONFERENCI INFORMACIJSKA DRUŽBA 2018 Multikonferenca Informacijska družba (http://is.ijs.si) je z enaindvajseto zaporedno prireditvijo osrednji srednjeevropski dogodek na področju informacijske družbe, računalništva in informatike. Letošnja prireditev se ponovno odvija na več lokacijah, osrednji dogodki pa so na Institutu »Jožef Stefan«. Informacijska družba, znanje in umetna inteligenca so še naprej nosilni koncepti človeške civilizacije. Se bo neverjetna rast nadaljevala in nas ponesla v novo civilizacijsko obdobje ali pa se bo rast upočasnila in začela stagnirati? Bosta IKT in zlasti umetna inteligenca omogočila nadaljnji razcvet civilizacije ali pa bodo demografske, družbene, medčloveške in okoljske težave povzročile zadušitev rasti? Čedalje več pokazateljev kaže v oba ekstrema – da prehajamo v naslednje civilizacijsko obdobje, hkrati pa so notranji in zunanji konflikti sodobne družbe čedalje težje obvladljivi. Letos smo v multikonferenco povezali 11 odličnih neodvisnih konferenc. Predstavljenih bo 215 predstavitev, povzetkov in referatov v okviru samostojnih konferenc in delavnic. Prireditev bodo spremljale okrogle mize in razprave ter posebni dogodki, kot je svečana podelitev nagrad. Izbrani prispevki bodo izšli tudi v posebni številki revije Informatica, ki se ponaša z 42-letno tradicijo odlične znanstvene revije. Multikonferenco Informacijska družba 2018 sestavljajo naslednje samostojne konference:  Slovenska konferenca o umetni inteligenci  Kognitivna znanost  Odkrivanje znanja in podatkovna skladišča – SiKDD  Mednarodna konferenca o visokozmogljivi optimizaciji v industriji, HPOI  Delavnica AS-IT-IC  Soočanje z demografskimi izzivi  Sodelovanje, programska oprema in storitve v informacijski družbi  Delavnica za elektronsko in mobilno zdravje ter pametna mesta  Vzgoja in izobraževanje v informacijski družbi  5. študentska računalniška konferenca  Mednarodna konferenca o prenosu tehnologij (ITTC) Soorganizatorji in podporniki konference so različne raziskovalne institucije in združenja, med njimi tudi ACM Slovenija, Slovensko društvo za umetno inteligenco (SLAIS), Slovensko društvo za kognitivne znanosti (DKZ) in druga slovenska nacionalna akademija, Inženirska akademija Slovenije (IAS). V imenu organizatorjev konference se zahvaljujemo združenjem in institucijam, še posebej pa udeležencem za njihove dragocene prispevke in priložnost, da z nami delijo svoje izkušnje o informacijski družbi. Zahvaljujemo se tudi recenzentom za njihovo pomoč pri recenziranju. V letu 2018 bomo šestič podelili nagrado za življenjske dosežke v čast Donalda Michieja in Alana Turinga. Nagrado Michie-Turing za izjemen življenjski prispevek k razvoju in promociji informacijske družbe bo prejel prof. dr. Saša Divjak. Priznanje za dosežek leta bo pripadlo doc. dr. Marinki Žitnik. Že sedmič podeljujemo nagradi »informacijska limona« in »informacijska jagoda« za najbolj (ne)uspešne poteze v zvezi z informacijsko družbo. Limono letos prejme padanje državnih sredstev za raziskovalno dejavnost, jagodo pa Yaskawina tovarna robotov v Kočevju. Čestitke nagrajencem! Mojca Ciglarič, predsednik programskega odbora Matjaž Gams, predsednik organizacijskega odbora i FOREWORD - INFORMATION SOCIETY 2018 In its 21st year, the Information Society Multiconference (http://is.ijs.si) remains one of the leading conferences in Central Europe devoted to information society, computer science and informatics. In 2018, it is organized at various locations, with the main events taking place at the Jožef Stefan Institute. Information society, knowledge and artificial intelligence continue to represent the central pillars of human civilization. Will the pace of progress of information society, knowledge and artificial intelligence continue, thus enabling unseen progress of human civilization, or will the progress stall and even stagnate? Will ICT and AI continue to foster human progress, or will the growth of human, demographic, social and environmental problems stall global progress? Both extremes seem to be playing out to a certain degree – we seem to be transitioning into the next civilization period, while the internal and external conflicts of the contemporary society seem to be on the rise. The Multiconference runs in parallel sessions with 215 presentations of scientific papers at eleven conferences, many round tables, workshops and award ceremonies. Selected papers will be published in the Informatica journal, which boasts of its 42-year tradition of excellent research publishing. The Information Society 2018 Multiconference consists of the following conferences:  Slovenian Conference on Artificial Intelligence  Cognitive Science  Data Mining and Data Warehouses - SiKDD  International Conference on High-Performance Optimization in Industry, HPOI  AS-IT-IC Workshop  Facing demographic challenges  Collaboration, Software and Services in Information Society  Workshop Electronic and Mobile Health and Smart Cities  Education in Information Society  5th Student Computer Science Research Conference  International Technology Transfer Conference (ITTC) The Multiconference is co-organized and supported by several major research institutions and societies, among them ACM Slovenia, i.e. the Slovenian chapter of the ACM, Slovenian Artificial Intelligence Society (SLAIS), Slovenian Society for Cognitive Sciences (DKZ) and the second national engineering academy, the Slovenian Engineering Academy (IAS). On behalf of the conference organizers, we thank all the societies and institutions, and particularly all the participants for their valuable contribution and their interest in this event, and the reviewers for their thorough reviews. For the sixth year, the award for life-long outstanding contributions will be presented in memory of Donald Michie and Alan Turing. The Michie-Turing award will be given to Prof. Saša Divjak for his life-long outstanding contribution to the development and promotion of information society in our country. In addition, an award for current achievements will be given to Assist. Prof. Marinka Žitnik. The information lemon goes to decreased national funding of research. The information strawberry is awarded to the Yaskawa robot factory in Kočevje. Congratulations! Mojca Ciglarič, Programme Committee Chair Matjaž Gams, Organizing Committee Chair ii KONFERENČNI ODBORI CONFERENCE COMMITTEES International Programme Committee Organizing Committee Vladimir Bajic, South Africa Matjaž Gams, chair Heiner Benking, Germany Mitja Luštrek Se Woo Cheon, South Korea Lana Zemljak Howie Firth, UK Vesna Koricki Olga Fomichova, Russia Mitja Lasič Vladimir Fomichov, Russia Blaž Mahnič Vesna Hljuz Dobric, Croatia Jani Bizjak Alfred Inselberg, Israel Tine Kolenik Jay Liebowitz, USA Huan Liu, Singapore Henz Martin, Germany Marcin Paprzycki, USA Karl Pribram, USA Claude Sammut, Australia Jiri Wiedermann, Czech Republic Xindong Wu, USA Yiming Ye, USA Ning Zhong, USA Wray Buntine, Australia Bezalel Gavish, USA Gal A. Kaminka, Israel Mike Bain, Australia Michela Milano, Italy Derong Liu, USA Toby Walsh, Australia Programme Committee Franc Solina, co-chair Matjaž Gams Vladislav Rajkovič Viljan Mahnič, co-chair Marko Grobelnik Grega Repovš Cene Bavec, co-chair Nikola Guid Ivan Rozman Tomaž Kalin, co-chair Marjan Heričko Niko Schlamberger Jozsef Györkös, co-chair Borka Jerman Blažič Džonova Stanko Strmčnik Tadej Bajd Gorazd Kandus Jurij Šilc Jaroslav Berce Urban Kordeš Jurij Tasič Mojca Bernik Marjan Krisper Denis Trček Marko Bohanec Andrej Kuščer Andrej Ule Ivan Bratko Jadran Lenarčič Tanja Urbančič Andrej Brodnik Borut Likar Boštjan Vilfan Dušan Caf Mitja Luštrek Baldomir Zajc Saša Divjak Janez Malačič Blaž Zupan Tomaž Erjavec Olga Markič Boris Žemva Bogdan Filipič Dunja Mladenič Leon Žlajpah Andrej Gams Franc Novak iii iv KAZALO / TABLE OF CONTENTS Mednarodna konferenca o visokozmogljivi optimizaciji v industriji, HPOI 2018 / International Conference on High-Performance Optimization in Industry, HPOI 2018 ............................................... 1 PREDGOVOR / FOREWORD ....................................................................................................................... 3 PROGRAMSKI ODBORI / PROGRAMME COMMITTEES ........................................................................... 5 On Using Real-World Problems for Benchmarking Multiobjective Optimization Algorithms / Tušar Tea .... 7 Bridging Theory and Practice Through Modular Graphical User Interfaces / Rehbach Frederik, Stork Jörg, Bartz-Beielstein Thomas ...............................................................................................................11 Expensive Optimisation Exemplified by ECG Simulator Parameter Tuning / Breiderhoff Beate, Naujoks Boris, Bartz-Beielstein Thomas, Filipič Bogdan .....................................................................................15 A Hybrid Optimization Strategy with Low Resource Usage for Large Scale Multi-objective Problems / Monteiro Wellington Rodrigo, Reynoso-Meza Gilberto ..........................................................................19 Electric Vehicle Routing Problem: State of the Art / Serrar Jihane, Ellaia Rachid, Talbi El-Ghazali .........23 Optimization of End-to-End Deep Learning for Obtaining Human-Like Driving Models / Dovgan Erik, Sodnik Jaka, Filipič Bogdan ...................................................................................................................27 A Bi-Objective Maintenance-Routing Problem: Service Level Consideration / Rahimi Mohammad, Talbi El-Ghazali ...............................................................................................................................................31 Study on Reducing Turn-Around Time of Multi-Objective Evolutionary Algorithm on an Industrial Problem / Fukumoto Hiroaki, Oyama Akira .........................................................................................................35 Evolution of Electric Motor Design Approaches: The Domel Case / Papa Gregor, Petelin Gašper, Korošec Peter .........................................................................................................................................39 Model-Based Multiobjective Optimization of Elevator Group Control / Vodopija Aljoša, Stork Jörg, Bartz- Beielstein Thomas, Filipič Bogdan .........................................................................................................43 From a Production Scheduling Simulation to a Digital Twin / Papa Gregor, Korošec Peter ......................47 Indeks avtorjev / Author index ......................................................................................................................51 v vi Zbornik 21. mednarodne multikonference INFORMACIJSKA DRUŽBA – IS 2018 Zvezek D Proceedings of the 21st International Multiconference INFORMATION SOCIETY – IS 2018 Volume D Mednarodna konferenca o visokozmogljivi optimizaciji v industriji, HPOI 2018 International Conference on High-Performance Optimization in Industry, HPOI 2018 Uredila / Edited by Bogdan Filipič, Thomas Bartz-Beielstein http://is.ijs.si 8. oktober 2018 / 8 October 2018 Ljubljana, Slovenia 1 2 PREDGOVOR Z optimizacijskimi problemi se v realnem svetu, zlasti pa v industriji, srečujemo vsakodnevno. Visokozmogljiva optimizacija temelji na združevanju računske moči in naprednih optimizacijskih algoritmov in se je pojavila kot odgovor na izzive, ki jih predstavljajo zahtevni optimizacijski problemi, ki so lahko visokodimenzionalni, multimodalni, šumni, dinamični, večkriterijski ali pa njihovo reševanje vključuje časovno zahtevne simulacije. Mednarodna konferenca o visokozmogljivi optimizaciji v industriji ( High-Performance Optimization in Industry, HPOI 2018) je mišljena kot forum za predstavitev primerov uporabe in izmenjavo izkušenj med akademskimi in industrijskimi partnerji o uvajanju visokozmogljive optimizacije. Poleg tega spodbuja nadaljnje širjenje metodologije in neposredno sodelovanje med akademskimi ustanovami in industrijo. Konferenca je aktivnost projekta Synergy for Smart Multiobjective Optimization (SYNERGY, http://synergy-twinning.eu) iz programa Twinning v Obzorju 2020. Eden od ciljev tega projekta je prenesti znanje, ki so ga pridobili partnerji v konzorciju, na druge raziskovalne ustanove in v industrijo, zlasti podjetja, ki sodelujejo v Slovenski strategiji pametne specializacije (S4). Pri doseganju tega cilja so člani projekta že predstavili svoje dosežke v visokozmogljivi optimizaciji na specializirani delavnici na Gospodarski zbornici Slovenije, nekatere pa predstavljajo tudi na tej konferenci. Program konference obsega 11 predstavitev, vsi prispevki pa so objavljeni v konferenčnem zborniku. Prispevalo jih je 21 (so)avtorjev, od katerih je večina sodelavcev projekta SYNERGY. Obravnavane teme vključujejo optimizacijsko metodologijo, pristope k premoščanju vrzeli med akademskimi ustanovami in industrijo ter študije primerov s področij transporta, avtomobilske industrije, inženirstva in proizvodnje. Zahvaljujemo se avtorjem za oddajo in predstavitve njihovih del, članom programskega odbora za ocenjevanje prispevkov, Institutu »Jožef Stefan« kot gostitelju srečanja in organizatorjem 21. Mednarodne multikonference Informacijska družba (IS 2018), katere del je tudi HPOI 2018, za organizacijsko podporo. Bogdan Filipič, Thomas Bartz-Beielstein 3 FOREWORD Optimization problems are met in the real world, and particularly in industry, on a daily basis. High-performance optimization (HPO) is founded on the coupling of high computing power and advanced optimization algorithms, and has emerged in response to the challenges posed by hard optimization problems that can be high-dimensional, multimodal, noisy, dynamic, multiobjective or involve time-consuming simulations in order to be solved. The International Conference on High-Performance Optimization in Industry (HPOI 2018) is meant as a forum for presenting use cases and exchanging experience among academic and industrial partners on deploying HPO. Apart from that, it stimulates further proliferation of the methodology and direct collaboration between academia and industry. The conference is an activity of the Horizon 2020 Twinning project “Synergy for Smart Multiobjective Optimization” (SYNERGY, http://synergy-twinning.eu). One of the objectives of this project is to spread the knowledge gained by the consortium partners to other research institutions and the industry, in particular to the companies participating in the Slovenian Smart Specialization Strategy (S4). Pursuing this goal, the project members have already presented their achievements in HPO at a specialized workshop at the Chamber of Commerce and Industry of Slovenia, and some of them are also being presented at this conference. The conference program consists of 11 presentations and the related papers are published in the proceedings. They were contributed by 21 (co)authors, most of them being the SYNERGY project members. The topics discussed include the optimization methodology, approaches to bridging the gap between academia and industry, and case studies from the domains of transportation, automotive industry, engineering, and manufacturing. We are grateful to the authors for submitting and presenting their work, the program committee members for reviewing the papers, the Jožef Stefan Institute for hosting the event, and the staff of the 21st International Multiconference on Information Society (IS 2018) that HPOI 2018 is part of for organizational support. Bogdan Filipič, Thomas Bartz-Beielstein 4 PROGRAMSKI ODBOR / PROGRAMME COMMITTEE Bogdan Filipič (chair) Thomas Bartz-Beielstein (chair) Erik Dovgan Peter Korošec Nouredine Melab Marjan Mernik Boris Naujoks Gregor Papa Jörg Stork El-Ghazali Talbi Tea Tušar Martin Zaefferer Jernej Zupančič This event has received funding from the European Union’s Horizon 2020 research and innovation programme under grant agreement No 692286. 5 6 On Using Real-World Problems for Benchmarking Multiobjective Optimization Algorithms Tea Tušar Department of Intelligent Systems Jožef Stefan Institute Ljubljana, Slovenia tea.tusar@ijs.si ABSTRACT ment. The artificial test problems that are being consis- Although the motivation to study multiobjective optimiza- tently used for benchmarking EMO algorithms have charac- tion algorithms comes from practice, there are only a few teristics that are not representative of real-world problems. challenging real-world problems freely available to the re- They also fail to incorporate the peculiarities of real-world search community. Because of this, algorithm benchmarking problems, which means that the algorithms need additional is performed primarily on artificial test problems. The most adjustments before they can be applied to real-world prob- popular artificial test problems have characteristics that are lems [8]. Furthermore, most studies do not investigate the not well-represented in real-world problems. This and the influence of the problem dimension on the performance of predominant inadequate performance assessment methodol- the algorithms and the performance assessment is often done ogy widen the gap between theory and practice in the field only at a predefined number of evaluations. This makes it of multiobjective optimization. The paper suggests to in- hard to predict which algorithm will perform best on a par- stead compare the algorithms with the anytime performance ticular real-world problem when less evaluations are allowed benchmarking approach of COCO (the Comparing Continu- than the (high) numbers usually used in the studies. ous Optimizers platform) on more realistic artificial problem suites as well as suites with diverse real-world problems. By The COCO platform [2, 10] resolves many of these issues listing the benefits of sharing the real-world problems with by providing an alternative to the overused test suites and a the community, the paper hopes to encourage domain ex- more rigorous approach to algorithm benchmarking. How- perts to embrace this practice. ever, in order to bridge the gap between theory and prac- tice, multiobjective optimization algorithms should be stud- Keywords ied and compared not only on well-understood and easy-to- multiobjective optimization, real-world problems, algorithm compute artificial functions, but also on real-world problems benchmarking with various characteristics. Currently, only a small num- ber of challenging real-world problems are freely available 1. INTRODUCTION to the EMO community, which hinders the development of algorithms that could be used ‘off the shelf’. Most real-world optimization problems found in science and engineering are inherently multiobjective. For example, the The purpose of this paper is to show the advantages of task of many engineering design problems is to find solutions benchmarking algorithms on real-world problems and to en- of high quality and low cost. Such problems seldom have a courage domain experts to share their hardest problems with single solution (called the ideal solution) that would opti- the researchers to their mutual benefit. mize all objective simultaneously. Rather, they have (possi- bly infinitely) many Pareto-optimal solutions that represent In the remainder of the paper, we first recall the purpose of different trade-offs among the objectives. These solutions algorithm benchmarking (Section 2). Then, we review the form the so-called Pareto set in the decision space and Pareto existing practice of benchmarking multiobjective optimiza- front in the objective space. tion algorithms on artificial test problems and remind of an available alternative in the form of the COCO platform (Sec- Evolutionary Multiobjective Optimization (EMO) [4] is one tion 3). Next, we mention some real-world problems that of the most active research areas that deal with multiobjec- have been made publicly available, discuss the benefits of tive problems. It studies algorithms that make no assump- sharing real-world problems and give recommendations for tions on the properties of the optimization problems, such proposing new real-world problems and performing bench- as linearity, continuity and unimodality, and are therefore marking with them (Section 4). We conclude with some applicable to a variety of problems, including black-box op- closing remarks (Section 5). timization ones. EMO algorithms have successfully solved numerous challenging real-world optimization problems [3]. 2. THE PURPOSE OF ALGORITHM Nevertheless, there is a large gap between theory and prac- BENCHMARKING tice in the EMO field (stemming from the one in Evolution- The no free lunch theorem implies that no optimization algo- ary Computation [18]), which is widened by the dominat- rithm performs best for all possible problems [22]. The ob- ing (inadequate) paradigm of algorithm performance assess- served differences in performance are due to the (more/less) 7 successful adaptation of the algorithms to the problem land- ample the C-DTLZ test suite [15], there is no established scapes [12]. It is therefore crucial that the test problems used test suite containing mixed-integer problems with multiple in comparison studies have characteristics that are represen- objectives. tative of real-world problems. Furthermore, although the problems from the mentioned Algorithm benchmarking, either when comparing variants suites are scalable in the number of variables (the problem of the same algorithm or a novel algorithm to an established dimension) and the number of objectives, performance stud- one, can be used to gain an understanding of the algorithms ies rarely investigate the scaling of the algorithms with the at hand. However, the ultimate purpose of algorithm bench- problem dimension. This is usually simply fixed to a value marking is to find the algorithm that is expected to perform (often 30), while the number of objectives is being changed. best for a specific target problem—a real-world problem of Such an approach to performance assessment is problematic interest. This entails that we have as it disregards one of the most defining characteristics of a problem—its dimension. (a) some knowledge about the characteristics of the target problem, Finally, most studies compare the performance of the al- (b) information on the performance of a number of algo- gorithms only at a specific point in time, determined by rithms on test problems with similar characteristics as the number of function evaluations. Because they provide those of the target problem, and no data on the performance of the algorithms prior to that (c) an understanding of what best is, i.e., we can define moment, the findings of such studies cannot be used to in- and measure the desired algorithm performance. fer algorithm performance when less evaluations are avail- able, making them effectively useless for the main purpose Then, machine learning methods can be used to select the of benchmarking mentioned earlier. most appropriate algorithm for the given target problem [16]. 3.2 Benchmarking with the COCO Platform 3. USING ARTIFICIAL PROBLEMS FOR COCO (Comparing Continuous Optimizers) [2, 10] is an ALGORITHM BENCHMARKING open-source platform for benchmarking black-box optimiza- Benchmarking multiobjective algorithms on artificial opti- tion algorithms. It implements different test problem suites mization problems has several advantages. The evaluations and provides an anytime performance assessment method- are cheap (computed instantaneously), the characteristics ology that is in line with the purpose of benchmarking as of the problems can be controlled, and the problems can be described in Section 2. Furthermore, COCO incorporates implemented in any programming language. If constructed the results of various optimization algorithms on its tests with care, the artificial problems can be scaled in the num- suites that are regularly being collected at BBOB (Black- ber of decision variables, constraints and objectives, and the Box Optimization Benchmarking) workshops [1] and can be Pareto sets and fronts can be known, which considerably readily used for comparisons with new algorithms. facilitates performance assessment. In addition to singleobjective test suites, such as the estab- The main question when using artificial test problems for lished bbob suite [11], COCO currently provides two test benchmarking algorithms is whether they are good repre- suites with biobjective problems, bbob-biobj with 55 func- sentatives of real-world problems. tions and its extended version bbob-biobj-ext with 92 func- tions [21], each instantiated in six dimensions (n ∈ {2, 3, 5, 10, 20, 40}) and ten instances (small alterations of the func- 3.1 Issues with the Prevailing Benchmarking tion, such as shifts, etc.). Every biobjective function is con- Methodology structed using two separate bbob functions—one for each Since the introduction of the DTLZ [6] and WFG [13] test objective. This approach is motivated by the nature of real- suites in 2001 and 2006, respectively, the vast majority of world multiobjective problems, where each objective corre- studies in EMO have been comparing algorithms on one or sponds to a separate singleobjective function. It is there- both of these two suites. In fact, they have been overused fore closer to real-world conditions than the constructions to such a degree that we can speculate on overfitting of op- with distance and position variables used by the DTLZ and timization algorithms to these problems. This is especially WFG test suites. However, this approach results in unknown concerning because they have some properties that are ben- Pareto sets and fronts, which is not convenient for perfor- eficial when designing test suites, but are not likely to be mance assessment purposes. In order to alleviate this issue, found in real-world problems. For example, in order have COCO provides approximations of the Pareto fronts for all a known Pareto set and a controllable shape of the Pareto problems, collected during several runs of various EMO al- front, the problems are parameterized by two sorts of vari- gorithms. These can be used in plots to showcase the charac- ables: distance variables, which indicate the distance of a so- teristics of the Pareto fronts and to compute the best known lution from the Pareto front, and position variables, which hypervolume [23] values for these problems. indicate the position of a solution along the Pareto front. The resulting Pareto sets and fronts are much easier to work The anytime performance assessment approach from COCO with than the irregularly shaped real-world ones. is based on the notion of runtime, i.e., the number of function evaluations needed to achieve a target hypervolume (see [9] Many real-world problems have additional difficulties, such and [21] for more details). This makes it possible to study as constraints or a mixed-integer decision space. While there the results for each problem separately as well as aggregate are some multiobjective test suites with constraints, for ex- them over all problems in a suite. For example, the plot 8 1.0 bbob-biobj f1-f55, 5-D HMO-CMA-ES Loshc tion three that are of different nature, but are very demand- 58 targets: 1..-1.0e-4 best 2016 10 instances ing and therefore suitable for algorithm benchmarking: UP-MO-CMA-ES Kra 0.8 RM-MEDA Auger bb • The Radar Waveform problem has an integer decision DEMO Tusar bbob- SMS-EMOA-SA Wess space that can be scaled from four to 12 decision vari- 0.6 NSGA-II-MATLAB A ables, and nine objectives [14]. SMS-EMOA-DE Auge SMS-EMOA-PM Auge • The HBV Benchmark Problem consists of calibrating RANDOMSEARCH-5 A the HBV rainfall–runoff model [19]. It has 14 real- 0.4 MO-DIRECT-HV-Ran valued decision variables and four objectives. RANDOMSEARCH-4 A MO-DIRECT-ND Won • The recently proposed Mazda Benchmark Problem [17] 0.2 MO-DIRECT-Rank W is a car structure design optimization problem with RANDOMSEARCH-100 222 integer decision variables, two objectives and 54 Fraction of function,target pairs MAT-DIRECT Al-Du constraint functions that make it hard to find a feasible 0.0 v2.2.1.417, hv-hash=ff0e71e8cd978373 MAT-SMS Al-Dujai 0 2 4 6 solution. log10(# f-evals / dimension) There are multiple reasons why only a few black-box real- world problems are being publicly shared. Sometimes, the Figure 1: Bootstrapped empirical cumulative dis- companies that have such problems hide them to protect tribution of the number of objective function their trade secrets. Other times, the reasons are of an im- evaluations divided by dimension for 58 targets plementation nature, for example because some proprietary with target precision in {100, 10−0.1, . . . , 10−4.9, 10−5, 0, software is needed to perform the evaluations. It is also pos- −10−5, −10−4.8, . . . , −10−4.2, −10−4} for 16 algorithms sible that people do not make their problems public simply on all 5-D functions of the bbob-biobj test suite. because they see no benefit in doing so. Most of these issues can be amended. If the domain ex- in Figure 1 shows the proportion of targets (on the y axis) perts wish to keep the details of the problem hidden, this that an algorithm is expected to achieve given the number can be achieved by sharing an executable program without of function evaluations (divided by the problem dimension, the source code. If the companies fear that their competitors on the x axis). The plot presents the results aggregated over could retrieve useful information already from how the prob- all instances of the 5-D functions of the bbob-biobj suite. lem is defined, a simple linear transformation can be used Note that such plots allow to compare the performance of to transform a box-constrained continuous decision space to algorithms that were run using a different budget of function [0, 1]n without affecting the nature of the problem landscape evaluations (up to the minimal common budget). (an integer or mixed-integer decision space can be handled in a similar way). Although the least noteworthy, some im- The COCO platform could similarly be used to benchmark plementation issues can be hardest to bypass. The best way real-world problems. might be to use freely available software instead of the pro- prietary one (this, of course, might not always be possible). 4. USING REAL-WORLD PROBLEMS FOR If conceivable, time-consuming evaluations using specialized ALGORITHM BENCHMARKING software can be replaced by surrogate models as was done, for example, in [17]. 4.1 Availability of Real-World Problems Real-world problems can be separated into those whose ob- jectives and constraints can be given in an analytic form and 4.2 Benefits of Sharing Real-World Problems others that are truly black-box problems, for example those Suppose a real-world problem is interfaced with the COCO that require complex computations or simulations to evalu- platform and used in the BBOB workshops to benchmark ate the functions and constraints of the problem. Note that multiobjective algorithms. This means that the researchers as soon as one function or constraint behaves like a black not only run their algorithms on the problem, but also sub- box, the entire problem is considered to be a black box. mit their results to COCO for use in future comparisons. The first and most obvious benefit of such a setting is that There are quite a few multiobjective real-world problems the interested EMO community would most likely find bet- of the first type, i.e., with a known analytic form. See for ter solutions to the problem in question than a single team example the problems from [5], [7] and [20]. Similarly to of researchers. Next, if the problem has some characteris- the artificial problems, they can be evaluated quickly and tics that are not well-represented in artificial test problems, implemented in any programming language. However, as such as a mixed-integer decision space, sharing such a prob- recently shown in [20], many such problems are not chal- lem will motivate the researchers to adapt their algorithms lenging enough to distinguish between algorithms and can to its characteristics. This means that in time, there will therefore be useful for benchmarking purposes only in test be more versatile algorithms for these kinds of problems to suites containing other, harder problems. choose from. Finally, it is likely that in the future, the same experts who shared this problem, will face another problem On the other hand, there are also many black-box real-world of similar nature. Then, the algorithms that performed best problems from various domains, but only a few of them are on the original problem might be readily used on the future freely available to EMO researchers. Here, we briefly men- alternative versions of this problem. 9 4.3 Recommendations [9] N. Hansen, A. Auger, D. Brockhoff, D. Tušar, and When proposing real-world benchmark problems, domain T. Tušar. COCO: Performance assessment. ArXiv experts should try to make them as flexible as possible. Ide- e-prints, arXiv:1605.03560, 2016. ally, it should be possible to instantiate them in a few differ- [10] N. Hansen, A. Auger, O. Mersmann, T. Tušar, and ent dimensions and also to create some instances of the same D. Brockhoff. COCO: A platform for comparing problem (minor modifications that do not change the nature continuous optimizers in a black-box setting. ArXiv of the problems). In addition to providing better grounds e-prints, arXiv:1603.08785, 2016. for performance assessment, this might also help to better [11] N. Hansen, S. Finck, R. Ros, and A. Auger. understand the problems in question. Real-parameter black-box optimization benchmarking 2009: Noiseless functions definitions. Research Report When benchmarking EMO algorithms, artificial test suites RR-6829, INRIA, 2009. with properties reflective of the real-world problems should [12] Y. Ho and D. Pepyne. Simple explanation of the be used in order to gain understanding about the algorithms. no-free-lunch theorem and its implications. Journal of In addition, the algorithms should also be tested on real- Optimization Theory and Applications, world problems to show their applicability in practice. Since 115(3):549–570, 2002. real-world problems come from various domains and might [13] S. Huband, P. Hingston, L. Barone, and R. L. While. have particular characteristics, the algorithms should be run A review of multiobjective test problems and a on suites of real-world problems from different domains. scalable test problem toolkit. IEEE Transactions on Evolutionary Computation, 10(5):477–506, 2006. 5. CONCLUSIONS [14] E. J. Hughes. Radar waveform optimisation as a This paper reviewed the many drawbacks of the existing many-objective application benchmark. In Proceedings practice of benchmarking multiobjective algorithms with the of EMO 2007, pages 700–714, 2007. over-used DTLZ and WFG test suites. Using the COCO [15] H. Jain and K. Deb. An evolutionary many-objective platform most can be amended, but the performance assess- optimization algorithm using reference-point based ment is still being done solely on artificial problem functions. nondominated sorting approach, part II: Handling The paper proposes to benchmark algorithms using COCO’s constraints and extending to an adaptive approach. anytime performance assessment on suites of real-world al- IEEE Transactions on Evolutionary Computation, gorithms in addition to the artificial ones. Some benefits of 18(4):602–622, 2014. sharing real-world problems with the EMO community are [16] P. Kerschke and H. Trautmann. Automated algorithm presented in hope to encourage greater exchange of knowl- selection on continuous black-box problems by edge between academia and industry. combining exploratory landscape analysis and machine learning. ArXiv e-prints, arXiv:1711.08921, 2017. 6. ACKNOWLEDGMENTS [17] T. Kohira, H. Kemmotsu, A. Oyama, and The author acknowledges the financial support from the T. Tatsukawa. Proposal of benchmark problem based Slovenian Research Agency (project No. Z2-8177). on real-world car structure design optimization. In Companion Proceedings of GECCO 2018, pages 7. REFERENCES 183–184, 2018. [1] The BBOB workshop series. [18] Z. Michalewicz. Quo vadis, evolutionary computation? https://numbbo.github.io/workshops/index.html, On a growing gap between theory and practice. In 2018. Proceedings of WCCI 2012, pages 98–121, 2012. [2] COCO: Comparing Continuous Optimizers. [19] P. M. Reed, D. Hadka, J. D. Herman, J. R. Kasprzyk, https://github.com/numbbo/coco, 2018. and J. B. Kollat. Evolutionary multiobjective [3] C. A. C. Coello and G. B. Lamont, editors. optimization in water resources: The past, present, Applications of Multi-Objective Evolutionary and future. Advances in Water Resources, 51:438–456, Algorithms. World Scientific, Singapore, 2004. 2013. [4] K. Deb. Multi-Objective Optimization Using [20] R. Tanabe and A. Oyama. A note on constrained Evolutionary Algorithms. John Wiley & Sons, multi-objective optimization benchmark problems. In Chichester, UK, 2001. Proceedings of CEC 2017, pages 1127–1134, 2017. [5] K. Deb and A. Srinivasan. Innovization: Innovative [21] T. Tušar, D. Brockhoff, N. Hansen, and A. Auger. design principles through optimization. KanGAL COCO: The bi-objective black-box optimization report number 2005007, Indian Institute of benchmarking (bbob-biobj) test suite. ArXiv e-prints, Technology, Kanpur, 2005. arXiv:1604.00359, 2016. [6] K. Deb, L. Thiele, M. Laumanns, and E. Zitzler. [22] D. H. Wolpert and W. G. Macready. No free lunch Scalable test problems for evolutionary multi-objective theorems for optimization. IEEE Transactions on optimization. KanGAL report number 2001001, Indian Evolutionary Computation, 1(1):67–82, 1997. Institute of Technology, Kanpur, 2001. [23] E. Zitzler, L. Thiele, M. Laumanns, C. M. Fonseca, [7] EMO’2017 Real-World Problems, Black Box and V. Grunert da Fonseca. Performance assessment Optimization Competition. of multiobjective optimizers: An analysis and review. https://bbcomp.ini.rub.de/, 2017. IEEE Transactions on Evolutionary Computation, [8] B. Filipič and T. Tušar. Challenges of applying 7(2):117–132, 2003. optimization methodology in industry. In Companion Proceedings of GECCO 2013, pages 1103–1104, 2013. 10 Bridging Theory and Practice Through Modular Graphical User Interfaces Frederik Rehbach Jörg Stork Thomas Bartz-Beielstein TH Köln TH Köln TH Köln Institute for Data Science, Institute for Data Science, Institute for Data Science, Engineering and Analytics Engineering and Analytics Engineering and Analytics Steinmüllerallee 1 Steinmüllerallee 1 Steinmüllerallee 1 51643 Gummersbach, 51643 Gummersbach, 51643 Gummersbach, Germany Germany Germany frederik.rehbach@th-koeln.de joerg.stork@th-koeln.de thomas.bartz-beielstein@th- koeln.de ABSTRACT rion (e.g. best model function value) is evaluated on the State-of-the-art evolutionary algorithms and related search expensive objective function and further used to update the heuristics are well suited to solve problems from industry. model. The process is repeated in an iterative fashion. A Unfortunately, easy to use graphical user interfaces (GUI) more in-depth explanation of SMBO and its applications can are not available for many algorithms. We claim that the be found in [5] and [2]. availability of well-designed GUIs might increase the accep- tance of these algorithms in the real-world domain. The SPOT has been further improved and developed for many spotGUI R-package, which is introduced in this paper, pro- years. Today the package provides a vast set of different vides a GUI for the already well-established SPOT package. models, optimizers, and sampling schemes, each of which It includes state-of-the-art algorithms and modeling tech- can be configured to user specific requirements. The system niques that can be used without the requirement of opti- was initially targeted to parameter optimization tasks, but mization or programming knowledge. Using the spotGUI is well suited to any costly to evaluate optimization prob- in industry, as well as education, delivered first promising lem. The availability of these methods together with their results. respective documentation in the R-package is a first step towards an easy to use modular optimization tool. How- Keywords ever, SPOT remains a high-level toolbox, which requires SPOT, Graphical User Interface, Real-World Applications user experience and some R programming skills. Further- more, since R is rarely used by engineers in industry, this 1. INTRODUCTION again leads to problems (a) and (b) as previously discussed. The presented spotGUI tries to address these problems by Industrial problems are highly complex and challenging for making the tools included in SPOT accessible to everyone even the most advanced state-of-the-art algorithms. How- through an easy to use graphical interface. ever, the difficulty in solving such problems is often not their high complexity, but rather the challenge for a non-expert The rest of this paper is structured as follows: Section 2 user to apply a suitable algorithm. For a significant subset gives an overview of the basic functionality and some con- of the existing optimization problems in industry, suitable ceptual ideas of the spotGUI. In Section3 two practical ex- state-of-the-art algorithms already exist. Yet, they are often ample applications for the spotGUI applied in industry are still not applied because they are presented. One of which is the Electrostatic Precipitator (ESP) Problem, a current, costly-to-evaluate, discrete opti- a) not known to the field specialist or mization problem from industry. Lastly, the software, future b) no simple implementation is available. opportunities, and room for improvements are discussed in Section 4. This paper presents a simple to use GUI that bridges the gap between existing algorithms and real-world problems. 2. WORKFLOW The core of the new package relies on the Sequential Param- eter Optimization Toolbox (SPOT) [1]. SPOT provides a 2.1 Availability modular structure for combining sampling methods, mod- The spotGUI package shall give more users easy access to eling techniques and optimizers for an all-in-one Surrogate SPOT. All stable versions are available on CRAN. Develop- Model-Based Optimization (SMBO) toolbox. In SMBO, a ment versions are published on GitHub. One of the primary data-driven surrogate model is fitted to the data of an ex- goals of the spotGUI is to allow non-R-users and even non- pensive to evaluate objective function, e.g., a complex sim- programmers to use SPOTs model-based optimization tech- ulation or a real-world experiment. Under the assumption niques. Additionally, it can benefit experienced SPOT users that the surrogate is cheap to evaluate, an extensive search by enabling a faster setup and even code generation which on the model becomes feasible. The predicted candidate will be covered in more detail in Section 2.5. The spotGUI solution, which best fulfills some user-specified infill crite- is developed in the R extension Shiny [4]. It is divided into 11 I: Setup II: Parameters III: Experiment IV: Save Objective Function Spot Config Run Spot Exports Variables and Model and Run, Evaluate R-Code Dimensions Optimizer Setup and Visualize Log Figure 1: Typical optimization workflow for SMBO in the spotGUI four separate tabs, arranged in a typical workflow order as presented in Figure 1 and Algorithm 2.1. Each of the tabs is explained in more detail in the following. Algorithm 2.1: Surrogate Model-based Optimization 1 step I: setup 2 select and parametrize objective function 3 begin Figure 2: Screenshot illustrating the objective func- 4 step II: parameters tion setup in the spotGUI. The user has to define the 5 select and parametrize surrogate model function as well as it’s dimensionality and variable 6 select and parametrize experimental design types. 7 step III: experiment 8 generate design points 9 evaluate design points with objective function some real-world experiments by entering / importing the 10 build initial surrogate model experiment results back into the spotGUI. The only con- 11 while not termination-condition do figuration required in this scenario is to insert information 12 search for optimum on surrogate model on the problem dimensions. Each dimension is configured 13 evaluate new point on the objective function with a type (numeric/integer/factorial), as well as upper 14 update surrogate model and lower bounds. If there are multiple dimensions with the 15 end same upper and lower bounds, the convenience option ”am- 16 step IV: save ntDimensions” can be used to specify that the same bounds 17 end are required multiple times. 2.2 Setup 2.3 Parameters The objective function is specified and parametrized on the One of the main benefits of the spotGUI becomes evident first setup tab. A screenshot of the configuration window is during the setup of SPOT itself. As previously mentioned shown in Figure 2. Additionally to having an option to insert SPOT features a wide variety of different models and op- any function through the R-Environment and supporting timizers, each of which again provides a variety of config- manual result input, the spotGUI provides a broad set of uration options. In the spotGUI, these are conveniently preconfigured test functions. selectable from drop-down menus. Showing each available option together with simple explanations through tooltips, The set of provided test functions is loaded from the ’smoof’ tackles the requirement of any documentation reading for R-package [3], which provides an interface to many single- the user. The settings are arranged in four categories cov- and also multi-objective test functions. Of these, the spot- ering a general setup, modeling setup, optimizer setup and GUI only includes the current set of single-objective func- lastly design setup. Skipping the ’Spot Config’ tab alto- tions, totaling in 76 test functions. Each of these functions gether results in a robust default setup for SPOT. is loaded with its respective bounds as well as dimensional- ity. Scalable functions are loaded as 2-dimensional functions 2.4 Experiment and can then be adapted by the user to any desired dimen- The previously configured processes are executed in the ”Run sionality. The ’smoof’ package also allows the user to filter Spot” tab. The available options include creating a DOE, the functions by specific tags such as ”separable”, ”differen- fitting a model, running a model-based optimization, and tiable” or ”weak-global-structure”. This makes it possible to more. In the following, these methodologies will be briefly test a given optimizer on a particular type of test function explained. In many expensive real-world applications, an that should behave somewhat similar to a real-world prob- initial screening for variable importance and interaction is lem that shall be solved. Different settings for SPOT and desired. The spotGUI provides the option to do so with its tools can quickly be tested by using the spotGUI with a configured sampling method to build a design of experi- the given set of test functions. ments. Depending on the objective function configuration, the generated experiments can be evaluated automatically or The possibility to manually input evaluation results enables manually, e.g. a real-world experiment. Such manual results non-programmers to use the spotGUI without any require- can either be imported into the spotGUI or directly entered ments for an objective function definition in code. Thus into the result table. A surrogate-model is fitted to the for example making it possible to use SPOT to optimize given data making interactive 3D-visualizations available. 12 3. EXAMPLE APPLICATIONS 3.1 Applying the Manual Mode The spotGUI offers a couple of functionalities to be easily usable and applicable to problems where real-world experi- ments are required. We can imagine the following example where the user is not too affine with software programming: A machine engineer who needs to set up a new metal hard- ening machine to deliver good performance. Through the machine’s interface, he is allowed to control two temperature parameters which define a temperature curve that the machine runs through in the hardening process. Additionally, he can change two time parameters which de- fine the duration of the heating as well as the cooling phase in the hardening process. He is looking for the set of opti- mized parameters which result in the hardest end product. However, each test requires to run the hardening machine for a few hours and involves material costs. In this scenario, the manual mode of the spotGUI could help the engineer in this parameter optimization problem. First of all, by using Figure 3: Auto generated plot showing the fitted the spotGUI in the manual mode, no coded fitness function surrogate model. Red dots indicate evaluated can- is used. Instead, parameter settings are proposed by SPOT, didate solutions. Hovering the mouse over the plot manually evaluated on the hardening machine and inserted results in the black info box showing more detailed into the results table by the engineer. information for the given plot location. The tool- bar above the plots provides features for easy plot The detailed workflow is as follows: After an initial setup exports. in the spotGUI, defining the bounds and types of the in- put parameters, a DOE (Design of Experiments) is built. This is quickly done via the ’createDOE’ button in the ’run- The graphics are generated through plotly, an R-library for Mode’ tab. A model can be fitted, and a visualization of it creating web-based graphs [7]. The availability of interactive is available. With the now to him available information, the 3D plots enables the user to learn more about the landscape engineer could continue in a few different ways. He could of their objective function intuitively and gives a deeper in- straightforward accept the best solution found in the DOE. sight into variable behavior. After a model fit, it is easily However, this should not be done if resources for more ma- possible to run an optimizer on the model to propose a sin- chine tests exist. Continuing with a more in-depth DOE, gle next candidate solution, thus enabling SMBO even to a he could increase the DOE budget and optionally shrink the manual user / non-programmer. parameter bounds to an area that is considered as promis- ing by the fitted model. The second option to spend the Further options are again aimed at enhancing the automatic remaining test budget is to run an optimizer on the fitted evaluation and optimization of a configured objective func- model via the ’propose new point’ functionality. This ad- tion. As sometimes even just a few objective function evalu- ditional point is the model optimum for some configured ations might take a long time, the spotGUI execution can be infill-criterion. This criterion might be the best-predicted interrupted and restarted from the last completed function point, but depending on the model, it could for example evaluation. For users who only want to use the spotGUI as also be the point with the highest expected improvement as a quick setup tool for their code, another option exists. By utilized in EGO [6]. After evaluating the proposed point on entering the ’Log Only’ mode, all computations that would the machine, the model can be refitted to include the new usually be applied to the objective function are skipped. In- data point. After that, the ’propose new point’ functionality stead, the actions are only written to the code log. From is usable again. Therefore, by using this feature, surrogate there they can be exported and used in any R-Script, en- model-based optimization is available in a manual use case, abling an extra fast setup for new SPOT projects. making a well-known and powerful optimization technique available to a broader audience. Lastly, the configuration 2.5 Save of the spotGUI can easily be changed during the optimiza- Each action that is executed in the spotGUI is written into tion process, allowing for a more interactive optimization an exportable R-Code log. The log is accessible on the ’Ex- approach. port’ tab of the GUI, it can easily be exported or copied to the clipboard through the provided button. The resulting R-Code can be run standalone (given the spotGUI library 3.2 The Electrostatic Precipitator Problem is installed) and generates the same results as previously Electrostatic precipitators (ESP)s are large scale electrical shown in the experimentation tab. This also ensures re- filtering/separation devices. They are used to remove solid producibility of any work that was done with the help of particles from gas streams, such as from the exhaust gases spotGUI. of coal-burning power plants. An overview of the struc- ture of an ESP can be seen in Figure 4. The illustrated separator has three central separation zones in which the 13 4. SUMMARY The SPOT package has been available for many years. It has been continuously updated and grew to a very large and useful platform. However, through the growing amount of possible configurations and use cases it simultaneously be- came more complex to dig through all settings and find the best ones for each problem. The here introduced spotGUI package reduces the configuration complexity back down to a level where any beginner can use the package. It was suc- cessfully applied to industry use cases as well as in student courses. Thus, demonstrating its ease of use and capability to provide easy to access visual information. The playful style with which different optimization methods can be ap- plied makes the software a useful tool in education. One of the most significant drawbacks of the current ver- sion of the spotGUI is its dependency on R. Till now, the Figure 4: Electrostatic precipitator with 3 separa- spotGUI can only be published as a web application avail- tion zones. This figure was kindly provided by Stein- able through a browser or started directly in R. Future work müller Babcock Environment GmbH. on the spotGUI will, therefore, concentrate on making the software available as a standalone executable without the requirement of starting it through R. Additionally, more fea- particles are separated from the gas flow by the precipita- tures are planned or even already are under construction, in- tor. Gas streams in from piping through the inlet hood and cluding: Parallelization support for SPOT, more DOE and exits through an outlet hood. The entrance and exit pip- analysis functionality, additional exports, and report gener- ing of the separator has a much smaller cross-section and ation. therefore a higher gas velocity than desired in the separa- tor. Without additional measures the fast gas stream would 5. ACKNOWLEDGEMENT rush through the center of the precipitator, resulting in very Parts of this work were supported by the ”Ministerium für low separation efficiency. The primary optimization target Kultur und Wissenschaft des Landes Nordrhein-Westfalen” is the so-called gas distribution system (GDS). The GDS is (FKZ: 005-1703-0011). mounted directly behind the flue gas inlet of the precipita- tor. It is used to distribute the gas flow from the small inlet 6. REFERENCES cross-section to the much larger cross-section of the precip- [1] T. Bartz-Beielstein and M. Zaefferer. A Gentle itation zones. The GDS in the given application consists of Introduction to Sequential Parameter Optimization. 49 configurable slots. Each of these slots can be filled with CIplus 1/2012, Fakultät 10 / Institut für Informatik, different types of metal plates, porous plates, angled plates, 2012. or be left completely empty. Increasing the separators effi- [2] T. Bartz-Beielstein and M. Zaefferer. Model-based ciency by achieving a more evenly distributed gas flow al- methods for continuous and discrete global lows a smaller overall separator. A reduced separator size, optimization. Applied Soft Computing, 55:154 – 167, together with lowered operating costs would accumulate to 2017. multiple millions of euro in cost reduction. [3] J. Bossek. smoof: Single- and Multi-Objective Optimization Test Functions. The R Journal, Two central factors reveal a complex to solve optimization 9(1):103–113, 2017. problem: [4] W. Chang, J. Cheng, J. Allaire, Y. Xie, and J. McPherson. shiny: Web Application Framework for a) The amount of configurable slots together with the R, 2018. R package version 1.1.0. amount of available configurations per slot leads to ≈ 1041 possible configurations for the overall system [5] A. Forrester, A. Keane, et al. Engineering design via b) Each objective function evaluation requires a costly surrogate modelling: a practical guide. John Wiley & CFD-simulation in order to judge the gas flow through Sons, 2008. the system [6] D. R. Jones, M. Schonlau, and W. J. Welch. Efficient global optimization of expensive black-box functions. The ESP optimization was approached with a combination Journal of Global optimization, 13(4):455–492, 1998. of a parallelized model-based evolutionary algorithm that [7] C. Sievert, C. Parmer, T. Hocking, S. Chamberlain, was equipped with newly created task-specific mutation and K. Ram, M. Corvellec, and P. Despouy. plotly: Create recombination operators. Tuning these operators was re- Interactive Web Graphics via ’plotly.js’, 2017. R quired in order to be able to reduce the overall runtime of package version 4.7.1. each optimization to fit into standard project run times. In this industry project, the spotGUI was successfully applied to set up parameter tuning for the evolutionary algorithm and its operators. 14 Expensive Optimisation Exemplified by ECG Simulator Parameter Tuning Beate Breiderhoff, Boris Naujoks, Bogdan Filipič Thomas Bartz-Beielstein Jožef Stefan Institute Institute of Data Science, Jamova cesta 39 Engineering and Analytics Ljubljana, Slovenia Steinmüllerallee 6 bogdan.filipic@ijs.si Gummersbach, Germany {firstname.lastname}@th-koeln.de ABSTRACT model is constructed using a three-dimensional grid [5, 6]. This article describes the tuning of an Electrocardiogram For the simulation, the APs are described mathematically (ECG) simulator as a benchmark problem to show the ap- and represent voltage as a function of time for an individual plication of surrogate modelling in complex global optimi- cell. The function AP(t) is parameterised with nine param- sation. After presenting the background on ECG, its sim- eters and is approximated by a combination of exponential ulation and the optimisation task, the main concepts and functions [20]. Out of nine AP function parameters, two methods of surrogate modelling and Efficient Global Opti- have predefined values, while the remaining seven are sub- misation (EGO) are presented. Here, next to the standard ject to optimisation. As three layers of heart muscle cells techniques regularly involved in the algorithm, alternative are considered in the model, the total number of optimisa- approaches are discussed briefly. Finally, first results apply- tion variables is 21. The optimisation goal is to find the ing the depicted algorithm on the ECG simulator optimisa- best set of parameters to produce properly shaped APs and tion problem are presented. approximate simulated ECG waveforms to a measured ECG waveform by perfecting the shape of the APs. Keywords Pearson’s correlation coefficient (PC) is the covariance of ECG, evolutionary algorithms, surrogate assisted optimisa- the two variables divided by the product of their standard tion deviations. The coefficient PC1 between the measured ECG waveforms and the simulator output builds the objective 1. INTRODUCTION function, which is required to be maximized to obtain a good The heart muscle pumps blood in a specific rhythm through- match between the two waveforms. out the entire body. In order to do this, the heart muscle requires an electrical impulse to contract. This electrical 3. THE EGO APPROACH impulse acts as a natural pacemaker. The electric current As each run of the ECG simulator takes around 15 minutes, is then transmitted via specific pathways throughout the finding the best solution is a time consuming process that heart, enabling regular contraction and relaxation. ECG is can take days or weeks. One way of relieving the burden of the result of recording this electrical activity of the heart expensive simulation runs is by constructing approximation over a period of time using electrodes placed on the skin. models that mimic the behavior of the simulator as closely It provides information about the heart’s rhythm and rate. as possible while being computationally cheaper to evaluate. The normal ECG shape and some typical defects are well The basic idea of using surrogate models in optimisation can known, but the transfer function that maps the ECG mea- be quite simple. First, the surrogate models for the objective sured on the skin to individual cells of the middle layer of the function with sufficient accuracy are built; second, the opti- heart wall is unknown. Gathering additional knowledge on mum is found by an optimizer, with the objective function the transfer function would help to improve ECG-based di- evaluated by surrogate models, rather than by the expensive agnostics and enable better prediction of health condition, simulation runs. Since prediction with the surrogate models based on the ECG reading. Simulation models may help, is much more efficient than that by the expensive simulation but the simulation of a human ECG signal is a complex runs, the optimisation efficiency can be largely improved. optimisation problem. Although the framework of the surrogate-based optimisation 2. THE ECG SIMULATOR is very intuitive and simple, questions may arise, e.g.: Is the The ECG simulator1 considered within the SYNERGY2 pro- surrogate model accurate enough and has the true optimum ject is a tool which tries to mimic the activity of the left been reached? The solution gained by the surrogate model ventricle of the heart, by producing ECG waveforms for a is only an approximation to the true optimum. One has to given set of Action Potential (AP) parameters. The heart refine the surrogate models by adding new sample points, which are observed by running the ECG simulator. The 1https://github.com/synergytwinning/ekgsim flowchart of the surrogate-based optimisation is sketched in 2http://synergy-twinning.eu/ Figure 1. 15 the Expected Improvement (EI, cf. [11]) infill criterion. EI Start not only considers the objective function provided by the no model but also the model quality to suggest new points for time-consuming evaluations. DoE EGO starts by building an initial Kriging model using some initial design points which are often produced by an exper- iment design method. Then, in each iteration, the point with the highest EI value is selected by using a traditional Build Surrogate Model optimisation algorithm. The selected point is evaluated us- (e.g. Kriging) ing the real expensive objective function and used to update the Kriging model. In such a way, the EI criterion guides the search toward the optimum of the real problem. Main - Optimise on Surrogate Sub - 3.2 Design of Experiments optimisation (e.g. GA) optimisation The first mandatory step in surrogate modelling is the col- lection of data to set up an initial model. This is normally done by a DoE approach [2], which results in an initial sam- Converge pling plan. By choosing an initial sampling plan the chal- No Add New lenge is to limit the number of samples but nevertheless get a Yes Data Point good and suitable design. There are various sampling tech- niques available such as Uniform Random Sampling, Latin Hypercube Sampling, and Orthogonal Array Sampling. A Run ECG Simulator common choice is Latin Hypercube Sampling (LHS, [16]) a statistical method for generating a near-random sample of Yes parameter values from a multidimensional distribution. No Converge 3.3 Modelling Approaches An important issue is the huge number of surrogate models available in the literature. Here we limit our discussion to Stop three popular techniques that are shortly described below. 3.3.1 Kriging Figure 1: Flowchart of the surrogate-based optimi- Kriging is a popular choice of surrogate models. It under- sation stands observations as realisations of a Gaussian process. The popularity of this technique is due to the fact that it not only produces accurate predictions, but also provides an The steps from the figure like Design of Experiments (DoE), estimate of the prediction uncertainty [14, 18]. building the surrogate model, optimising on the surrogate etc. selection are explained in a bit more detail in the follow- 3.3.2 Random Forests ing. Here, we mention common techniques next to promi- Random Forests [3] are ensembles of prediction trees such nent alternatives. A focus is put on the techniques that are that each tree depends on the values of a random vector used for addressing the ECG simulator optimisation prob- sampled independently and with the same distribution for lem. all trees in the forest. The generalisation error of a forest of tree classifiers depends on the strength of individual trees 3.1 Single-objective Surrogate Modelling in the forest and the correlation between them. Internal In single-objective surrogate-assisted optimisation, there ex- estimates monitor the error, strength, and correlation, and ists only one objective function, which is the fitted surrogate these are used to show the response to increasing the num- model using the acquired data points. The most straightfor- ber of features. Internal estimates are also used to measure ward approach is to find the global optimum of this model. variable importance [3]. The major problem is that the search may stall at a local op- timum. Solving this problem implies that the search needs 3.3.3 Support Vector Regression to combine exploration and exploitation; i.e., the search ex- SVR is a modelling technique based on the theory of sup- plores the total experimental area and zooms in on the local port vector machines [7, 19]. SVR models produce a pretty area with the apparent global optimum. accurate estimate of the objective function, provided that a suitable kernel is selected and parameters are appropri- Efficient Global Optimisation [11, 10] is a popular search ately tuned. This tuning process is expensive, especially for heuristic that tries to realize this exploration and exploita- models with higher dimensions and a high amount of sample tion. Many alternatives exist, one is the pre-selection ap- points [12]. proach described in [9]. EGO is a widely used surrogate- based optimisation algorithm for expensive single-objective optimisation specialised on utilising Kriging modelling and 16 3.4 Optimisers Choosing a suitable search strategy which can perform effec- tive global optimisation is the most difficult part in surrogate- assisted optimisation. In our work we use two well-known optimisers, namely an Evolutionary Algorithm and Simu- lated Annealing. 3.4.1 Evolutionary Algorithms (EAs) EAs are metaheuristics inspired by the process of natural selection that belong to the larger class of evolutionary al- gorithms. They are commonly used to generate high-quality solutions to optimisation and search problems by relying on bioinspired operators such as mutation, crossover and selec- tion [4, 8]. A subclass of EAs are Genetic algorithms (GAs). 3.4.2 Simulated Annealing (SA) SA is a method to solve complex optimisation problems [13]. This method models the physical process of heating a ma- terial and then slowly lowering the temperature to decrease Figure 2: Box-plot of RMSE obtained for cross- defects, thus minimizing the system energy. At each iter- validation of PC1 ation of the simulated annealing algorithm, a new point is randomly generated. The distance of the new point from the current point, or the extent of the search, is based on worst (as its weight was optimized to a value of zero). a probability distribution with a scale proportional to the temperature. The algorithm accepts all new points that improve the objective value, but also, with a certain proba- 5. CONCLUSIONS AND OUTLOOK bility, points that worsen the objective value. By accepting The efficient global optimisation approach is presented. Next points that raise the objective, the algorithm avoids being to standard settings from the methods involved in the algo- trapped in local optima in early iterations and is able to rithm, some alternatives were discussed briefly. The result- explore better solutions globally. ing algorithm is applied to the ECG simulator optimisation problem and first results are presented. 3.5 Model Selection and Validation However, the parameters used for optimisation (i.e., models K-fold cross-validation is an improved scheme which allows invoked, evolutionary algorithm, simulated annealing etc.) us to use most of the data for constructing the surrogates. were not tuned for the best performance. Parameter tuning In general, the final quality of the surrogate model is judged of optimisers might further enhance the surrogate-assisted using the mean and the standard deviation of the root-mean- optimisation process. Here, SPOT [1] might be invoked, square error (RMSE) for each cross-validation set [17]. which provides a set of tools for model-based optimisation and tuning of algorithms. It also includes surrogate models, 4. RESULTS optimisers and design of experiment approaches. The purpose of this article is to summarize approaches to surrogate modelling which are applicable to the ECG simu- The simulator provides two simulated ECG signals at differ- lator. The various surrogate models selected were Kriging, ent positions on the body surface. A second coefficient PC2 SVR, RF and a convex combination ensemble of the former could be used for multi-objective optimisation, also known three models. as multi-criteria optimisation or Pareto optimisation. It is a special case of solving optimisation problems with more The ensemble model performed best in K-fold cross-validation than one objective function to be optimised simultaneously. tests, while SVR performed worst. This provided an insight The final result is a set of solutions known as Pareto opti- of how the models would actually perform during the op- mal solutions. The Pareto front is a set of non-dominated timisation process as shown in Figure 2. Single-objective solutions, being chosen as optimal, if no objective can be optimisation was carried out to investigate the performance improved without sacrificing at least one other objective. of the surrogates in a practical scenario. The single-objective surrogate-assisted optimisation yielded some pretty interest- 6. ACKNOWLEDGMENTS ing facts about the behavior of the ECG simulator. Firstly, This work is funded by the European Commission’s H2020 pro- the maximum value that is achieved for the objective func- gramme, through the Twinning project SYNERGY under Grant tion is 0.31. The EGO algorithm based on Kriging (Ex- Agreement No. 692286 as well as through the UTOPIAE Marie pected Improvement) using simulated annealing performed Curie Innovative Training Network, H2020-MSCA-ITN-2016, Grant superiorly relative to other strategies when comparing the Agreement No. 722734. mean and standard deviation of the best obtained values as shown in Figure 3 [15]. 7. REFERENCES The optimisation of the weight vector for building the en- [1] T. Bartz-Beielstein, C. W. G. Lasarczyk, and semble model revealed that the SVR model did perform M. Preuss. Sequential parameter optimization. In 2005 17 Figure 3: B-plot of best observed values for optimisation of PC1 IEEE Congress on Evolutionary Computation, [11] D. R. Jones, M. Schonlau, and W. J. Welch. Efficient volume 1, pages 773–780 Vol.1, 2005. global optimization of expensive black-box functions. [2] G. Box, J. Hunter, and W. Hunter. Statistics for Journal of Global Optimization, 13(4):455–492, 1998. Experimenters: Design, Innovation, and Discovery. [12] S. Karamizadeh, S. M. Abdullah, M. Halimi, Wiley, New York, 2nd edition, 2005. J. Shayan, and M. j. Rajabi. Advantage and drawback [3] L. Breiman. Random forests. Machine Learning, of support vector machine functionality. In Computer, 45(1):5–32, 2001. Communications, and Control Technology (I4CT), [4] K. A. DeJong. Evolutionary Computation: A Unified pages 63–65, 2014. Approach. MIT Press, Cambridge, MA, 2006. [13] S. Kirkpatrick, C. D. Gelatt, and M. P. Vecchi. [5] M. Depolli, V. Avbelj, and R. Trobec. Optimization by simulated annealing. Science, Computer-simulated alternative modes of U-wave 220(4598):671–680, 1983. genesis. Journal of Cardiovascular Electrophysiology, [14] G. Matheron. Principles of geostatistics. Economic 19(1):84–89, 2008. Geology, 58(8):1246–1266, 1963. [6] M. Depolli, R. Trobec, and B. Filipič. Asynchronous [15] A. Mazumdar. Performance analysis of different master-slave parallelization of differential evolution for models in surrogate assisted optimization. Master’s multiobjective optimization. Evolutionary thesis, TH Köln, 2018. Computation, 21(2):261–291, 2013. [16] M. D. McKay, R. J. Beckman, and W. J. Conover. A [7] H. Drucker, C. J. C. Burges, L. Kaufman, A. Smola, comparison of three methods for selecting values of and V. Vapnik. Support vector regression machines. In input variables in the analysis of output from a International Conference on Neural Information computer code. Technometrics, 21(2):239–245, 1979. Processing Systems (NIPS96), pages 155–161. MIT [17] N. Queipo, R. Haftka, W. Shyy, T. Goel, Press Cambridge, MA, 1996. R. Vaidyanathan, and P. Kevin Tucker. [8] A. E. Eiben and J. E. Smith. Introduction to Surrogate-based analysis and optimization. Progress in Evolutionary Computing. Natural Computing Series. Aerospace Sciences, 41:1–28, 2005. Springer, Berlin, 2 edition, 2015. [18] J. Sacks, W. J. Welch, T. J. Mitchell, and H. P. [9] M. T. Emmerich, K. C. Giannakoglou, and Wynn. Design and analysis of computer experiments. B. Naujoks. Single- and multiobjective evolutionary Statistical Science, 4(4):409–423, 1989. optimization assisted by Gaussian random field [19] A. J. Smola and B. Schölkopf. A tutorial on support metamodels. IEEE Transactions on Evolutionary vector regression. Statistics and Computing, Computation, 10(4):421–439, 2006. 14(3):199–222, 2004. [10] A. Forrester, A. Sobester, and A. Keane. Engineering [20] B. Wohlfart. A simple model for demonstration of design via surrogate modelling: A practical guide. John STT-changes in ECG. European Heart Journal, Wiley & Sons, 2008. 8(4):409–416, 1987. 18 A Hybrid Optimization Strategy with Low Resource Usage for Large Scale Multi-objective Problems Wellington Rodrigo Monteiro Gilberto Reynoso-Meza Pontifical Catholic University of Paraná Pontifical Catholic University of Paraná R. Imaculada Conceição, 1155 R. Imaculada Conceição, 1155 Curitiba, Paraná, Brazil Curitiba, Paraná, Brazil wellington.monteiro@pucpr.edu.br g.reynosomeza@pucpr.edu.br ABSTRACT might be more interested to purchase the thighs or other The use of multi-objective approaches to solve problems in cuts, all the other parts must also be processed and sold industry grew in the last years. Nevertheless, these strate- somehow. gies are still unused in many fields where their performance is suboptimal or when they are too complex to be implemented As a result, the production plan is by itself a problem com- or even are simply unknown. One example is in the poultry posed of a large number of variables (at least 2000) [7]. This industry with its particularly complex chain. In this paper, plan is usually executed as a single-objective problem (the we will discuss a hybrid multi-objective approach with low objective being the overall profits) since it is currently able computational resource usage intended for this scenario as to provide results within the same day. However, it is well well as other similar ones. known that these profits differ from the real values since the reliability of the plants vary. By reliability impacts we Keywords mean both internal (e.g. different production costs between multi-objective, optimization, many variables, low resource the plants, worker strikes, unscheduled maintenance, stock usage. issues) or external (suppliers, weather) causes that are re- sponsible to reduce the projected profits. Therefore, the 1. INTRODUCTION company could benefit if the reliability of the production plans was known beforehand—if it was converted to a multi- In corporate environments, the use of simpler tools in some objective problem (MOP) with the objectives being the ex- situations might be favored against other tools that would pected profits and the reliability of said plans, depending generate better results. This business decision might be on the case the analyst could choose a production plan that caused due to the lack of technical understanding of these has less expected profits, but higher reliability rates. On tools or due to cost and performance reasons. One practical the other hand, he or she could also be more aggressive and example shown in a previous work of the author [7] is the attempt higher profits, but also with higher chances of not Production Plan algorithm used by one of the largest meat achieving the expected value. companies in the world, specifically its poultry business line. The work presented in [7] proved it was possible to convert By definition, the supply chain in industry is a complex set the production plan of the company into a multi-objective of operations and resources that must be extremely opti- problem. However, two problems were found: 1) there was mized in order to achieve its maximum potential which does an issue with the input data (the reliability rates of each include the management of upstream and downstream rela- plant), which resulted in very similar grades for all the plants tionships in order to achieve an outcome which is more prof- and 2) the multi-objective optimization algorithm took too itable to all the parties in the chain [2]. One of its parts is long (more than 24 hours in an i7 desktop with 16 GB RAM) the Production Plan, which defines what should one or more to generate the production plans. Nevertheless, before its plants build considering a myriad of variables—i.e. market implementation can be greenlit by the company, additional demand, production line capacity, logistics and stock limits, work on both sides is required. As such, while the data issue suppliers constraints, raw material limits, etc. Therefore, an was corrected by the plants themselves so that it can be accurate Production Plan is a key component to maximize usable, the multi-objective optimization algorithm needed the potential profits. to be greatly improved in order to be executed faster and with lower resource usage so that it could be used by an In the meat industry the challenges are greater. Since it is off-the-shelf corporate laptop and able to generate a Pareto livestock, at the same time the supply chain is very large and front approximation in under one hour. very tight [3]. The former because its production involves genetics, feeding, breeding, growth control from the current Considering this background, the proposal is to generate a animal up to its grandparents and the latter because there multi-objective optimization (MOO) algorithm intended for are very strict sanitary (including, but not limited to the problems with a large number of variables (more than 2000 health, safety and environment) controls with the ration, since the original problem is expected to grow in complex- water, effluents, temperature and vaccination, for example. ity). As such, the main objective here is to have an algorithm Also, the whole animal must be pushed to the market—in that balances both the performance (i.e. low computational the case of the chicken, for instance, even though the market 19 resource usage) and the scalability (i.e. capable of processing and the Ackley’s function and, as shown in [1], implements problems with thousands or tens of thousands of variables). features such as a chaos-based pseudo random number gen- On that end, a test set with similar characteristics of the erator (to address the nonuniform grouping), a correlation real-world problem will be used to evaluate the algorithm. matrix to keep track of the correlations between variable groups and the objective groups, different basic fitness land- The remainder of this document is structured as follows. scape functions (to achieve the mixed separability) and also The second section explains the background of the test prob- making use of a linkage function to define the variable link- lems that will be used as well as their rationale. The third ages. section shows the proposal to modify and test the presented case into a MOP. Then, the following section specifically 3. PROPOSAL shows the technical details of the created MOP as well as As mentioned in the Section 2, this article proposes to imple- its results after optimization. The fifth section presents the ment an algorithm capable of resolving multi-objective prob- conclusions and the future work to be done from this docu- lems with many variables with low resource usage. For this ment. reason, the approach chosen was a hybrid multi-objective genetic algorithm. By hybrid, as shown in [4], it is consid- 2. BACKGROUND ered to be an algorithm (in this case a genetic algorithm) Currently, in the meat industry some optimization solutions enhanced with an additional local search step shortly after are both single-objective and using specialized algorithms the selection of the individuals. The algorithm, shown in built from scratch with the profitability in mind such as the Algorithm 1 and implemented in MATLAB c , heavily OtimixTM. However, said algorithms might not enable the focuses on the parallelization for performance purposes. On industries to use two or more objectives or provide greater the other hand, it attempts to overuse large matrices in or- parameter tuning possibilities. Since such algorithms are der to reduce the memory footprint. targeted towards only one objective, the problem designer usually has only one solution as the result of the minimiza- Algorithm 1: The proposed algorithm tion/maximization, eliminating the possibility to analyze the tradeoffs between different production plans consider- Data: design space, objective vector Result: the Pareto front approximation ing two or more objectives. As a result, the production generate n random solutions; planners are required to empirically consider the differences rank the solutions by a dominance filter; between the plants from a reliability standpoint, leaving no for each generation until it reaches the stopping criteria do possibility to compare the solutions based on this factor. store previous solutions and their objective values; generate offspring by tournament, recombination and The scenario that originated this algorithm had two objec- mutation; join the offspring to the other solutions; tives: profitability and reliability [7]. 2032 variables were partially rank all solutions with a dominance filter; employed from which all were integers—however, this sce- locally improve the best solutions; nario was known to be a test—therefore, more variables were replace the best solutions with locally improved expected. As such, considering these characteristics and the solutions; other business requirements, the new algorithm had to meet rank all the joined solutions with a dominance filter; the following objectives: prune some solutions by their crowding distances. end • Be a multi-objective optimization algorithm; • Be able to resolve problems with many variables (more The local improvement algorithm (shown in Algorithm 2) than 1000, ideally with more than 15000), all integers; has a new, optional parameter named number of random • Low computer resource usage (preferably less than 1 neighbors for the local search (NumberRandomNeighbors). GB RAM per 5000 variables); It is used to improve the performance in scenarios where • Be able to generate a Pareto front approximation (even there are too many variables and this algorithm would take if there is still room for improvements) in less than one too long to process all of the variables. hour. For example, if a solution has 1000 variables and this pa- Since the original data needed additional work from the rameter is not used, this algorithm will create 1000 new teams responsible for it, the alternative was to choose eas- solutions based on the original solution where the first so- ily configurable and reliable mathematical problems. The lution assigned a new random value for the first variable choice was the test problems in [1]—nine large-scale, multi- while keeping all the other variables intact; the second so- objective problems were considered to evaluate the proposed lution assigned a new random value for the second variable algorithm, each configured with 1000, 5000, 15000, 30000 while keeping all the other variables intact and so on. If 100 and 50000 variables. All the variables are integers similar to solutions are involved in the local search, at least 100 thou- the production plan problem—all of the problems were set sand new solutions would be created as a result. On the to have 2 objectives, analogous to the production plan prob- other hand, if a solution had 15000 variables, considering lem. These test problems are henceforth named LSMOPn, the same 100 solutions as a result 1.5 million new solutions where n is the problem number ranging from 1 to 9 and form- would be created. Considering the local search would hap- ing a different problem. All the tests use a combination of pen more than once during the algorithm execution, this six basic single-objective functions. These functions are the method—based on the exhaustive neighborhood exploration sphere function, the Schwefel’s function, the Rosenbrock’s [5] would take too long. This results in greater performance function, the Rastrigin’s function, the Griewank’s function gains between the different generations. 20 Algorithm 2: Local improvement Problem NVar gamultiobj sp-Mode II new algorithm Data: a solution LSMOP1 1000 0.998933 0.916951 0.966355 Result: N eighbors LSMOP2 1000 0.531942 0.472204 0.497616 LSMOP3 1000 0.999881 0.916340 0.960246 initialize the list of neighbors N eighbors; LSMOP4 1000 0.562390 0.494451 0.521077 if N umberRandomN eighbors was given then LSMOP5 1000 0.997482 0.427428 0.628847 select N umberRandomN eighbors random variables; LSMOP6 1000 0.999599 0.626870 0.916751 else LSMOP7 1000 0.999644 0.615992 0.947478 select all the variables; LSMOP8 1000 0.995152 0.508990 0.686416 end LSMOP9 1000 0.991079 0.483426 0.695966 foreach one of the variables selected do LSMOP1 5000 0.990506 0.931265 0.958904 copy the original solution; LSMOP2 5000 0.512829 0.470381 0.488354 randomly modify its value according to its bounds; LSMOP3 5000 0.996004 0.925004 0.954693 evaluate the new solution; LSMOP4 5000 0.545665 0.496839 0.516004 add it to N eighbors; LSMOP5 5000 0.995279 0.561639 0.609145 LSMOP6 5000 0.947114 0.575706 0.677455 end LSMOP7 5000 0.999411 0.752861 0.882027 replace N eighbors with only its anchors. LSMOP8 5000 0.993857 0.599694 0.609580 LSMOP9 5000 0.958439 0.507450 0.593416 LSMOP1 15000 0.979783 - 0.945874 4. TEST LSMOP2 15000 0.513371 - 0.484434 LSMOP3 15000 0.998651 - 0.941039 The new, proposed algorithm and the other algorithms were LSMOP4 15000 0.535847 - 0.510049 tested on an i7 desktop equipped with 16 GB RAM and LSMOP5 15000 0.974157 - 0.688771 LSMOP6 15000 0.771705 - 0.585666 a dedicated video card running Windows 10 Pro running LSMOP7 15000 0.999059 - 0.845907 MATLAB c R2015a. This desktop ran all the tests over the LSMOP8 15000 0.989229 - 0.680188 course of three weeks with one weekly system restart. In LSMOP9 15000 0.990760 - 0.540451 MATLAB c , the start and end time of each algorithm were LSMOP1 30000 - - 0.840306 tracked (therefore storing the time taken for each execution) LSMOP2 30000 - - 0.466373 LSMOP3 30000 - - 0.885842 as well as the Pareto front approximation found for each ex- LSMOP4 30000 - - 0.486888 ecution, each limited to 300 seconds. Outside MATLAB c , LSMOP5 30000 - - 0.192226 the Windows Performance Monitor (perfmon.exe, a known, LSMOP6 30000 - - 0.604263 LSMOP7 30000 - - 0.223234 built-in performance monitor tool in Windows) was used to LSMOP8 30000 - - 0.184088 track performance and memory usage in MATLAB c . As LSMOP9 30000 - - 0.668866 such, it was able to properly track resource usage by each al- LSMOP1 50000 - - 0.831989 gorithm. The new algorithm was tested against MATLAB c ’s LSMOP2 50000 - - 0.476713 own gamultiobj, a built-in, optimized algorithm based on LSMOP3 50000 - - 0.887729 LSMOP4 50000 - - 0.485269 NSGA-II and a Multi-objective Differential Evolution with LSMOP5 50000 - - 0.220266 Spherical Pruning algorithm (sp-MODE II) [8]. LSMOP6 50000 - - 0.573581 LSMOP7 50000 - - 0.081000 LSMOP8 50000 - - 0.208853 4.1 Evaluation Methods LSMOP9 50000 - - 0.669033 With the aforementioned tools both the hypervolume found for each algorithm considering the best Pareto front approxi- Table 1: Median hypervolume values found for each mations found for them and the memory and processor usage algorithm and problem. The values are relative to for each algorithm were measured. the utopia and nadir determined from all solutions in the Pareto front approximations found for all runs For each of the nine LSMOP problems and number of vari- and all algorithms. The best values are in bold. ables (1000, 5000, 15000, 30000 and 50000 variables), 51 runs were executed for each one of the three algorithm. Num- berRandomNeighbors was set to 50 in the new algorithm On the performance side, the memory allocation behavior and the initial number of random solutions set to 20 for all was similar for the problems with the same size for all the algorithms, limited to 200 * number of generations and a three algorithms. The Figure 1 represents the memory usage maximum of 5 generations without improvement (where an for all the runs with 5000 variables for a given problem. improvement is determined when the utopia had improved First, from 11:34:49 PM to around 03:48:00 AM all the 51 in at least 0.0001% for at least one objective) per run. runs for the new algorithm took place, followed by the 51 runs of the gamultiobj from that time to around 7:00:00 AM and later by the 51 runs of sp-MODE II. The new algorithm, 4.2 Findings as specially shown in Figure 1, used a little more than 1 GB The values shown in Table 1 represent the median values in memory in order to achieve the results. The processor for each one of the runs for each algorithm and for each usage registered the same pattern with around 15% in overall problem based on the implementation in [6]. From the hy- usage for the new algorithm. pervolume in itself the new algorithm proved incrementally better performance when the problem has more variables— Furthermore, the memory usage grew considerably when in fact, starting with 15000 variables the other algorithms changing the problem size to 15000 variables, as seen in are unable to run these problems due to out of memory er- an attempt recorded in the Figure 2. While the new al- rors (as shown by their lack of median values depending on gorithm, with its runs recorded between 11:28:23 PM to the case). around 3:30:00 AM kept the memory usage around 1 GB in 21 Figure 3: Hypervolume distribution for all the prob- lems for 15000 variables. The y -axis represents the Figure 1: Results for the memory usage for three hypervolume values, where ga represents the runs algorithms with 5000 variables. The y-axis represent with gamultiobj and new with the new algorithm. the memory usage in hundreds of megabytes. variables since each industrial boiler, production line and memory, gamultiobj (executed from that time up to around work shift could also be accounted. 8:00:00 AM) registered around 10 GB in usage. sp-MODE II, on the other hand, attempted to allocate more than 16 6. ACKNOWLEDGMENTS GB (as registered by the last spike to the right), culminating The authors would like to thank the Pontifical Catholic Uni- in an out-of-memory error. This behavior happened again versity of Paraná for supporting the research. This study with 30000 and 50000 variables, but with gamultiobj as well. was financed in part by the Coordenaç˜ ao de Aperfeiçoa- The processor usage also showed the same behavior from an mento de Pessoal de N´ıvel Superior - Brasil (CAPES) - Fi- usage standpoint. With 15000 variables the new algorithm nance Code 001 and the PQ2|304066/2016-8 grant. used around 25% from the processor, peaking in around 35% with 50000 variables, depending on the problem. 7. REFERENCES [1] R. Cheng, Y. Jin, M. Olhofer, and B. Sendhoff. Test problems for large-scale multiobjective and many-objective optimization. IEEE transactions on cybernetics, 2016. [2] M. Christopher. Logistics & supply chain management. Pearson UK, 2016. [3] F. Flanders and J. R. Gillespie. Modern livestock & poultry production. Cengage Learning, 2015. [4] V. Kelner, F. Capitanescu, O. Léonard, and L. Wehenkel. A hybrid optimization technique coupling an evolutionary and a local search algorithm. Journal of Computational and Applied Mathematics, 215(2):448–456, 2008. [5] A. Liefooghe, J. Humeau, S. Mesmoudi, L. Jourdan, and E.-G. Talbi. On dominance-based multiobjective local search: design, implementation and experimental Figure 2: Results for the memory usage for three al- analysis on scheduling and traveling salesman gorithms with 15000 variables. The y-axis represent problems. Journal of Heuristics, 18(2):317–352, 2012. the memory usage in hundreds of megabytes. [6] MATHWORKS. Hypervolume approximation - matlab central, 2015. 5. CONCLUSIONS [7] W. R. Monteiro and G. Reynoso-Meza. A multi-criteria The proposed algorithm showed that it is better suited for based approach for the production distribution in the problems with a large number of variables. As of the time of poultry industry. In 24th ABCM International writing this paper the company still could not make a new, Congress of Mechanical Engineering. Brazilian Society real-world case available for tests. On the other hand, the of Mechanical Sciences and Engineering, 2017. new algorithm can also be implemented in other complex [8] G. Reynoso-Meza, J. Sanchis, X. Blasco, and scenarios where, for example, the objectives are both the M. Mart´ınez. Design of continuous controllers using a maximum profits for a given distributed production plan and multiobjective differential evolution algorithm with the energy and steam consumption reduction in the plants. spherical pruning. In European Conference on the The energy and steam consumption variables would result Applications of Evolutionary Computation, pages in an estimated multi-objective problem with around 40000 532–541. Springer, 2010. 22 Electric Vehicle Routing Problem: State of the Art Jihane Serrar Rachid ELLAIA El-Ghazali Talbi LERMA Laboratory, E3S center LERMA Laboratory, E3S center BONUS, INRIA Lille-Nord Europe Mohammed V university of Rabat Mohammed V university of Rabat University of Lille, France BP765, Ibn Sina av, Rabat, Morocco BP765, Ibn Sina av, Rabat, Morocco 59655 - Villeneuve d’Ascq cedex serrar.jihan@gmail.com ellaia@emi.ac.ma el-ghazali.talbi@univ-lille1.fr ABSTRACT we assign to each arc the cost considering that the dis- Electric vehicles (EVs) represent a clean alternative but have tance is known and given for each arc. Time windows some limitations especially in terms of autonomy. Therefore, associated with nodes or arcs may also be defined in efficient routing of EVs is crucial to encourage their use. some problems. This article surveys the existing research related to electric vehicle routing problems (EVRP) and their variants. It ex- • Demand: The demands are either given for each node amines EVRP in terms of their definitions, their objectives, and are known in advance in the case of deterministic and algorithms proposed for solving them. problem or given by probabilistic formulas. Keywords • Fleet: The fleet refers to a set or a group of EVs. In optimization, electric vehicles, routing problem fact, it is associated to the electric vehicles available to the routing problem, hence we can either have a 1. INTRODUCTION homogeneous or a heterogeneous fleet. Although, electric vehicles face several challenges such as: the low energy density of batteries compared to the fuel of • Electric vehicles recharging technologies: Un- combustion engined vehicles, the long recharge times com- like the combustible vehicles, the electric vehicles are pared to the relatively fast process of refueling a tank and charged by plugging the car to the electric grid. There the scarcity of public charging stations, they contribute to are four main technologies: a sustainable and environmental friendly freight transporta- tion by reducing the air pollution. – Household charging: EVs can be charged by a conventional household plug using a cable and a Electric vehicle routing problem and variants are considered connector in the vehicle. This technology is slow. as optimization problems and, more specifically, they be- – Fast charging: This technology is a conductive long to the combinatorial optimization problem that can be charging method. It’s faster than the previous solved by two types of solution methods: exact methods and one. approximate ones. – Wireless charging systems: also known as induc- This paper presents an overview of different problems re- tive charging is an emerging technology that al- lated to the electric vehicle routing, different variants and lows EV recharging without the use of a cabled solution methods found in the scientific literature. The rest connection. of the paper is structured as follows. First, Section 2 gives the main characteristics of EVRP. Then Section 3 enumer- – Battery swap: It’s a a high-speed method. ates the electric vehicles drawbacks. Section 4 describes the variants of EVRP presented in the literature. Finally, Sec- • Cost: The cost is a term that depends on different tion 5 reports on different solution methods for EVRP. parameters. It depends on the distance traveled, the energy consumed and the time of the travel. In ad- 2. CHARACTERISTICS OF EVRP dition, in the case of the time window variant, there are some penalties that could be added if the window The electric vehicle routing problem aims at routing a fleet isn’t respected. Moreover, the cost changes from one of EVs on a given network, or a subset of a network, to recharging technology to another. serve a set of customers under specified constraints in order to optimize one or several fixed objective(s). So, an EV routing problem can be defined in terms of the following • Objectives: It could be a single-objective problem or components: a multi-objective problem according to the number of objectives considered. The objectives are very diver- sified because the EVRP has a lot of components in • Network: The network can be represented as a graph its definition, for instance, minimizing the total trav- composed of nodes referring to cities, customers and elling distance, the delay time and the waiting time, depots and arcs standing for connections. Sometimes, the total cost, etc. 23 3. ELECTRIC VEHICLES CHALLENGES 4.3 Electric Vehicle Routing Problem To combat environmental and energy challenges, electric ve- Lin et al. [7] presents a general Electric Vehicle Routing hicles may provide a clean and safe alternative to the inter- Problem (EVRP) that seeks to optimize the routing problem nal combustion engine vehicles. However, electric vehicles while minimizing the total cost related to the distance as well are still facing several weaknesses: as to the energy consumption by the battery. The proposed EVRP finds the optimal routing strategy in which the total 3.1 Autonomy Limitations cost is minimized such that each customer is visited exactly once by one vehicle on its route, the total demand of the The vehicles have a much smaller driving range due to the customers served on a route does not exceed the vehicle limited battery capacity. The range of an electric vehicle capacity. depends on the number and type of batteries used but gen- erally the driving range varies between 80 and 130 km for light duty EVs according to [13]. 4.4 Electric Vehicle Routing Problem with Non- Linear Charging Functions (EVRP-NL) 3.2 Long Charging Times Montoya [12] extended current EVRP models to consider EVs often have long recharge times compared to the rel- partial charging and non-linear charging functions which is atively fast process of refueling a tank which takes just a more realistic for the charging process. In EVRP-NL, the couple of minutes. Its charging time ranges between 0.5 and task consists of minimizing the total traveling distance as 12 hours as mentioned in [12]. Hence, the user must think well as the charging time since it does not depend on the about refueling at night for example. total tour distance. 3.3 Scarce Charging Infrastructure 4.5 Electric Vehicle Routing Problem with Time The number of electric recharging stations is still very small Windows (EVRP-TW) compared with that of conventional fuel stations as the elec- This variant seeks to satisfy the order of customers within tric fuelling points are still in the development stages. So, certain time window. Many researches have been interested the driver must do a research about the plug-in stations in studying this variant. Some of works found in the litera- localisation to know where and when he will have the op- ture are outlined below. portunity the recharge his EV. 4.5.1 EVRP-TW with recharging stations 4. EVRP AND VARIANTS In fact the time window variant of EVRP was first intro- Several versions and extensions of the basic electric vehicle duced by Schneider et al. [14]. They studied the electric routing problem have been presented in the literature. vehicle routing problem with time windows and recharging stations (E-VRPTW) which incorporates the possibility of recharging at any of the available stations considering that 4.1 Green Vehicle Routing Problem (GVRP) the required recharging time depends on the state of the Erdo˘ gan and Miller-Hooks [3] are the first to introduce the charge. Hence, electric vehicles, which have a restricted ca- Green VRP which consists of alternative fuel-powered ve- pacity, must reach cutomers whithin a time window while hicle fleets with limited driving range and limited refueling minimizing the number of vehicles used and the total travel infrastructure. The objective is to minimize the total dis- distance. tance traveled by the vehicles while allowing them to visit stations when necessary. 4.5.2 Electric vehicle routing problems with time win- In [10], Koç et al. proposed the same problem as Erdo˘ gan dows and Miller-Hooks with the motivation of saving the ecosys- Desaulniers et al. [2] tackled the routing problem in which tem and the health of humans while serving and executing route planning has to take into account the limited driving the transportation and good distribution process. range of EVs and the customer time window. The authors studied four variants of this problem. The first one allows More recently, J. Andelmin et al. [8] also studied the green a single recharge per route knowing that batteries must de- vehicle routing problems taking into account the several par- part fully recharged from the station, the second one permits ticularities of autonomy and charging process of this type multiple recharges but only full rechargement are allowed of vehicles. Hence, the refueling stops are allowed. Their unlike the next one where partial battery recharges are al- model aims to find optimal routes while minimizing the to- lowed but just one time and the last one with partial but tal distance and by using a homogeneous fleet of vehicles. multiple recharges permitted. Contrary to Erdo˘ gan and Miller-Hooks, they didn’t put the restriction on the number of vehicles that must be used. 4.5.3 The electric fleet size and mix vehicle routing problem with time windows and recharging sta- 4.2 The Green Vehicle Routing Problem With tions Multiple Technologies And Partial Recharges Hiermann et al. [6] aim to optimize the fleet and the vehicle Felipe et al. [4] presents a variant in which different charg-routes including the choice of recharging times and recharg- ing technologies are considered and partial EV charging is ing stations as the refuelling operation is assumed necessary allowed in recharging stations when needed in order to en- for EVs because of the limited capacity storage of electricity sure the continuity of the route. by batteries. They considered that the fleet is heterogeneous 24 which adds complexity to the problem. Furthermore, they modified Clarke and Wright savings (MCWS) heuristic as incorporate the time windows constraint where customers the original Clarke and Wright algorithm was developed to have to be reached within a specified time interval. tackle the classical vehicle routing problem and its variants, thus it was modified to take into consideration the need to 4.5.4 The recharging vehicle routing problem with visit stations that have to be inserted in the routes while avoiding redundant. Meanwhile, the second heuristic is the time window density-based clustering algorithm (DBCA) that consists of Conrad and Figliozzi [1] introduced the recharging vehicle forming clusters in a clustering step dedicated for that and routing problem wherein vehicles with a limited range must then the MCWS algorithm is applied for each single cluster. service a set of customers, but may recharge at certain cus- tomer locations instead of using only dedicated recharging 5.2 Exact Algorithms stations while operating whithin customer time window. In Desaulniers et al. [2] decided to solve the different variants other words, the battery of a vehicle can be recharged while of EVRP-TW presented in their paper using exact meth- servicing the customer if needed. Also, the authors showed ods. They used the exact branch-price-and-cut algorithms the impact of the customer time windows on the tour dis- adapted to each variant. Hence, for each variant a set of tance taking into account that the driving range is limited routes is generated and for that monodirectional and bidi- and the recharging time is long. rectional labeling algorithms are presented. 4.5.5 Electric vehicle routing problem with time win- Branch-and-price is a metaheuristic that was used by Hier- dows and mixed fleet mann et al. [6] to solve the E-FSMFTW which is formulated Goeke et al. [5] proposed to study a mixed fleet of electric as a mixed-integer linear program (MILP). In fact, the algo- vehicles and internal combustion vehicles. They consider rithm has to insert the charging constraints in its procedure. that the energy consumption function isn’t linear and follows a realistic model depending on multiple parameters like the Exact methods were also used by J. Andelmin et al. [8] speed of the vehicle and the load distribution. Hence, EVs to solve set partitioning (SP) formulation of the green ve- can recharge anytime en route to enhance the driving range. hicle routing problem where each variable corresponds to a simple circuit of a route, thus each SP contains a limited subset of routes. The authors proposed an exact method 4.5.6 Partial recharge strategies for the electric ve- composed from two phases: Phase I computes the lower and hicle routing problem with time windows upper bounds, while Phase II executes the set partitioning In their work, M. Keskin et al. [9] relax the full recharge heuristic and the dynamic programming algorithm. restriction and allow partial recharging in order to minimize time. Therefore, shorter recharging durations are allowed Koç & Karaoglan [10] implemented the B & C (branch and especially when the customer time window is set. The ob- cut) algorithm for the exact solution of the GVRP where jective of the model proposed is to minimize the total dis- the initial solution is generated using the classical simulated tance while respecting the time constraints. Concerning the annealing. In addition, the authors adapted the simulated partial recharge scheme, the charging process is identified annealing to the problems related to the electric vehicle rout- by a continuous decision variable. ing problem by adding the GVRP constraints to improve the results. At each step of the method the new solution is 4.5.7 Heterogenous electric vehicle routing problem compared with the current one so that the best solutions is accepted. with time dependent charging costs and a mixed fleet 5.3 Local Search Heuristics Sassi et al. [13] studied a new real-life routing problem in In [4], some constructive and local search heuristics have which they consider a number of realistic features such as: been proposed by Felipe et al. to find feasible routes while different charging technologies, coupling constraints between considering the recharge constraints as well as the real-world vehicles and charging technologies, charging station avail- size problems. In addition, the authors used the 48A algo- ability time windows, and charging costs depending on the rithm in which they consider 48 combinations of improving time of the day. Also, partial charging is allowed and the algorithms with different neighborhood structures. cost of vehicles as well as the total travel and charging costs. In their study, Sassi et al. [13] formulated the Heteroge- 5. SOLUTION APPROACHES TO EVRP FROM nous Electric Vehicle Routing Problem with Time Depen- LITERATURE dent Charging Costs and a Mixed Fleet (HEVRP-TDMF) using a Mixed Integer Programming Model. And to solve it, In the literature, many studies work on finding sophisticated they worked with a Charging Routing Heuristic (CRH) in and efficient solution methods that can be applied to EVRP. order to find feasible routes. This algorithm is based on two main steps: the first one manages the charging of EVs in 5.1 MCWS and DBCA depot and the second one solves the problem starting from Two heuristics were proposed by Erdo˘ gan and Miller-Hooks the depart of EVs from the depot. Moreover, a Local Search [3] with the goal of finding a set of routes that represents Heuristic based on the Inject-Eject routine with three dif- the feasible solution of the green vehicle routing problem ferent insertion strategies has been introduced. knowing that the authors have formulated it as a mixed- integer linear program (MILP). Actually the first one is the 25 5.4 Hybrid Heuristics [2] G. Desaulniers, F. Errico, S. Irnich, and M. Schneider. Exact algorithms for electric vehicle-routing problems 5.4.1 Hybrid VNS/TS heuristic with time windows. Operations Research, To solve the E-VRPTW, Schneider et al. [14] used a combi- 64(6):1388–1405, 2016. nation of Variable Neighbourhood Search (VNS) and Tabu [3] S. Erdo˘ gan and E. Miller-Hooks. A green vehicle Search (TS) heuristics in order to make use of the diversi- routing problem. Transportation Research Part E: fication of the search provided by the VNS algorithm and Logistics and Transportation Review, 48(1):100–114, the efficiency of TS as many combinatorial problems have 2012. proved that this last heuristic is very strong. This combina- [4] Á. Felipe, M. T. Ortu˜ no, G. Righini, and G. Tirado. A tion has the aim to find feasible solutions while respecting heuristic approach for the green vehicle routing all the constraints. problem with multiple technologies and partial recharges. Transportation Research Part E: Logistics 5.4.2 Multi-space sampling heuristic + Hybrid meta- and Transportation Review, 71:111–128, 2014. heuristic [5] D. Goeke and M. Schneider. Routing a mixed fleet of Montoya [12] adapted the multi-space sampling heuristic electric and conventional vehicles. European Journal of (MSH) used before to tackle the VRP with stochastic de- Operational Research, 245:81 – 99, 2015. mand [11] to the green vehicle routing problem by designing [6] G. Hiermann, J. Puchinger, S. Ropke, and R. F. Hartl. a tailored route extraction procedure. MSH is a heuristic The electric fleet size and mix vehicle routing problem that consists of two main phases: the sampling phase and with time windows and recharging stations. European the assembling phase. Furthermore, a hybrid metaheuristc Journal of Operational Research, 252(3):995–1018, is proposed to tackle the EVRP with non linear charging 2016. function. The metaheuristic combines two heuristics: the [7] O. W. Jane Lin, Wei Zhou. Electric vehicle routing iterated local search and heuristic concentration. problem. The organization, Transportation Research Procedia, 2016. The 9th International Conference on 5.5 Adaptive Large Neighborhood Search Al- City Logistics, Tenerife, Canary Islands (Spain). [8] E. B. Juho Andelmin. An exact algorithm for the gorithm green vehicle routing problem. Transportation Science, The ALNS algorithm was also used by Hiermann et al. [6]. 71:1288 – 1303, 2017. In order to optimize the location of the refueling stations [9] M. Keskin and B. Catay. Partial recharge strategies during the routing process, a hybrid heuristic has been pro- for the electric vehicle routing problem with time posed. This heuristic is a combination of the Adaptive windows. Transportation Research Part C: Emerging Large Neighbourhood Search (ALNS) and an embedded lo- Technologies, 65:111 – 127, 2016. cal search procedure that uses different neighbourhoods. In- [10] C. Koc and I. Karaoglan. The green vehicle routing deed, the local search was used to itensify and strengthen problem: A heuristic based exact solution approach. the search operation guided by the ALNS. Applied Soft Computing, 39:154 – 164, 2016. [11] J. E. Mendoza and J. G. Villegas. A multi-space Moreover, like Hiermann et al., Goeke et al. [5] developed sampling heuristic for the vehicle routing problem the Adaptive Large Neighborhood Search algorithm to ad- with stochastic demands. pages 1503–1516, 2013. dress the Electric Vehicle Routing Problem with Time Win- [12] J.-A. Montoya. Electric Vehicle Routing Problems : dows and Mixed Fleet. They also enhanced the algorithm models and solution approaches. Computation and by a local search for intensification. Language. PhD thesis, Université de Nantes Angers Le Mans, 2016. Also, the ALNS algorithm was proposed by M. Keskin et al. [9] to solve the EVRP with time window. The authors [13] O. Sassi and A. Oulamara. Electric vehicle scheduling formulated the problem as a mixed integer linear program. and optimal charging problem: complexity, exact and heuristic approaches. International Journal of Production Research, 55(2):519–535, 2017. 6. CONCLUSIONS [14] M. Schneider, A. Stenger, and D. Goeke. The electric Over the last several years, the green vehicle routing problem vehicle-routing problem with time windows and has been widely studied. This survey lists the main works recharging stations. Transportation Science, that exist in the scientific literature since its appearance in 48(4):500–520, 2014. 2011 by Erdo˘ gan and Miller-Hooks. Based on this paper, the models that have been proposed are single-objective. Yet, most of real problems in industry are multi-objective by nature, so a multi-objective variant of EVRP must be proposed. 7. REFERENCES [1] R. G. Conrad and M. A. Figliozzi. The recharging vehicle routing problem. In IIE Annual Conference. Proceedings, pages 1–8. Institute of Industrial Engineers-Publisher, 2011. 26 Optimization of End-to-End Deep Learning for Obtaining Human-Like Driving Models Erik Dovgan Jaka Sodnik Bogdan Filipič Department of Intelligent Faculty of Electrical Department of Intelligent Systems Engineering Systems Jožef Stefan Institute University of Ljubljana Jožef Stefan Institute Jamova cesta 39 Tržaška cesta 25 Jamova cesta 39 SI-1000 Ljubljana, Slovenia SI-1000 Ljubljana, Slovenia SI-1000 Ljubljana, Slovenia erik.dovgan@ijs.si NERVteh, raziskave in razvoj, bogdan.filipic@ijs.si d.o.o. Kidričeva ulica 118 SI-1236 Trzin, Slovenia jaka.sodnik@fe.uni-lj.si ABSTRACT data, i.e., driving data of only a small subset of driving sit- Modeling human driving with human-like driving models uations. Consequently, the time to collect the driving data can help companies in the evaluation of human drivers. Whi- is reduced, while the driver or more precisely his/her clone le a human-like driving model can be tested in various sce- is still evaluated in a large number of situations. narios, this is not feasible for driver evaluation due to time constraints. During the evaluation, only a small set of driv- Existing work has demonstrated that end-to-end approach ing data can be typically collected for each driver, which for learning to drive is appropriate when large sets of learn- represents an issue for advanced modeling approaches such ing data are available [5, 6, 14]. On the other hand, the as deep learning. To overcome this issue, an optimization problems with small sets of learning data have not been ad- approach is proposed, which tunes deep learning when a dressed appropriately. This paper aims at tackling this issue small learning dataset is available. by enhancing end-to-end deep learning approach with opti- mization in order to obtain human-like driving models from Keywords small sets of learning data. optimization, deep learning, human-like driving models The paper is further organized as follows. Section 2 presents 1. INTRODUCTION the optimization approach for end-to-end deep learning. Ex- periments and results are described in Section 3. Finally, Human-like driving models have been learned with several Section 4 concludes the paper with ideas for future work. methods, such as ARX models [8], Gaussian processes [11], Gaussian mixture models [1], artificial neural networks [15], 2. OPTIMIZATION OF END-TO-END DEEP support vector regression [13], etc. Recently, Deep Neural Networks (DNN) are being effectively used in learning tasks LEARNING from various application fields. For example, when driving End-to-end deep learning approach applies deep neural net- a vehicle, DNN can be used to recognize the road, other works to learn the transformation between the input and the vehicles, pedestrians, etc. from video data [7]. Moreover, output data. The main property of this approach is that a DNN has been also applied to directly learn the control single model is used to obtain this transformation. There ex- actions from video data without firstly reconstructing the ist also other approaches that decompose the problem and scene. This approach is called end-to-end learning and its apply specific models for each subproblem. For example, examples aim to learn steering, throttle and braking control one model can be used to recognize the objects, while an- actions, etc. [5, 6, 14]. other model can be used for higher-level reasoning [7]. The end-to-end approach aims at solving all the subproblems at Unfortunately, deep learning has a significant drawback: it once with a single model [5]. requires a lot of learning data. Existing driving datasets used for training DNN models vary from about 10 hours Existing work in the field of end-to-end deep learning for to up to 10,000 hours [14]. However, in some cases such a obtaining human-like driving models has shown that the se- large set of driving data is not available. For example, the lection of deep learning model and its parameter values is deep learning approach can be used to assess a driver, e.g., not straightforward [9]. In addition, the data need to be if he/she drives safely, is able to avoid critical situations, augmented to learn how to recover from poor positions or etc. [12]. This can be done by building a human-like driver orientations [3]. We propose to automate the selection of model, i.e., a clone of the driver, and test it in a large number appropriate parameter values and data augmentation func- of driving situations. A similar approach has been applied in tions with an evolutionary algorithm. Evolutionary algo- related domains where the goal was to learn human behav- rithms are search and optimization algorithms inspired by ior [10]. This procedure requires only a small set of driving the principles of biological evolution. They work with a set 27 • Batch size: Parameter of the deep learning algorithm. Defines the number of training examples utilized in one learning iteration. • Number of epochs: Parameter of the deep learning al- gorithm. Defines the number of passes through the training dataset during learning. • Image multiplier: Data augmentation parameter. De- fines how many times an image is multiplied. If it is multiplied, it is divided into overlapping subimages. For example, if the image is multiplied by 3, three im- ages are created containing: 1) left 80 % of the original image; 2) central 80 % of the original image; 3) right 80 % of the original image. The control actions are also appropriately adapted. For the left images steer- ing is added to simulate turning right, while for the right images steering is subtracted to simulate turning left. • Noise added to output: Data augmentation parame- ter. Defines the amount of noise an to be added to the control actions. The amount of noise is randomly selected at each time step with a uniform distribution Figure 1: Overview of the algorithm for obtaining U (−an, an). human-like driving models. • Flip image: Data augmentation parameter. Defines whether randomly selected images should be vertically of solutions that are improved through several generations flipped. If the image is flipped, the control action is by applying genetic operators, i.e., selection, crossover and also appropriately adapted. mutation [4]. • Activation function: Parameter of the neural network We propose to discover human-like driving models in two model. Defines the activation function of the neural steps. In the initial step, driving models that are able to network layers. drive the vehicle along a route are built, while in the final • Kernel regularizer: Parameter of the neural network step, these models are enhanced to imitate human driving. model. Defines the regularization of the neural net- The approach presented in this paper focuses on the initial work layers, which applies penalties on layer weights. step by applying an evolutionary algorithm to maximize the The penalties try to keep the weights small, which re- length of the route that has been traveled by the driving duces the possibility of overweighting a small subset of model during the simulation. Each solution (consisting of layer’s input data and prevents overfitting. parameters of model construction) is evaluated by applying the following steps: 3. EXPERIMENTS AND RESULTS The developed approach was tested on two scenarios. Both 1. The learning data are augmented to enable recovery scenarios did not contain traffic vehicles or pedestrians. For from poor situations or orientations. both scenarios, the same architecture of the neural network 2. The deep learning algorithm is used to learn a human was used. This architecture is shown in Figure 2 and is based driving model. on the architecture presented in [2]. It contains five convo- lutional layers and three fully connected layers. The con- 3. The driving model is evaluated on a route to measure volutional layers extract features, from simple features such the route length of feasible driving. as lines to complex features such as road contour. The fully connected layers implement the vehicle controller, which cal- The driving simulation stops if the driving becomes infeasi- culates the control action based on the extracted features. ble (e.g., the vehicle goes offroad) or when the entire route is traveled. The evolutionary algorithm applies tournament 3.1 First scenario selection (tournament size = 2), two-point crossover (prob- The first scenario consisted of a circular route of around 2 ability = 0.9) and uniform mutation (probability = 0.1) to km, which is shown in Figure 3a. An example of a route improve the solutions over generations. An overview of the image as input to the neural network is shown in Figure 4a. developed algorithm and its steps, i.e., evolutionary algo- The learning data were obtained from one driving along the rithm steps (selection, crossover and mutation) and solution route. evaluation steps (data augmentation, model building and model evaluation), is shown in Figure 1. The proposed approach was evaluated by tuning only a sub- set of the parameters listed in Section 2, which already en- The evolutionary algorithm optimizes the following deep abled us to obtain models that drove along the entire route learning and data augmentation parameters: for this scenario and consequently no additional parameters 28 (a) (b) Figure 4: Examples of the input images: (a) first scenario, (b) second scenario. Table 1: Parameter values for the first scenario Parameter Values image multiplier {1, 3, 5} noise added to output {0, 0.1} flip image {true, f alse} batch size 40 (not tuned) number of epochs 100 (not tuned) activation function elu (not tuned) kernel regulizer none (not tuned) driven. Figure 2: Architecture of the neural network. The proposed approach was evaluated with the parameter values shown in Table 3. The results show that the built models were able to drive only short routes (see Figure 5). However, it should be noted that due to high time complex- ity of deep learning, only a small number of generations were executed. More precisely, it took more than 17 days to ex- ecute 30 generations on a 3.6 GHz desktop computer with 16 GB RAM. The analysis of the results also shows that the activation function had the most significant effect on the re- sults. It turned out that the majority of models that were able to drive more than 450 m, contained the relu activa- tion function. In addition, the models were able to drive on (a) (b) straight segments, but had issues with crossroads. This is Figure 3: Maps of the testing routes: (a) first sce- probably due to the relatively simple architecture of neural nario, (b) second scenario. network. For example, images of the first route (see Figure 4a) are significantly less complex in comparison to the im- ages of the second route (see Figure 4b), since they do not needed to be tuned. The values of tuned and not tuned contain any buildings, sidewalks, crossroads, etc. Conse- parameters are shown in Table 1. quently, more complex architectures of the neural network are needed for the city roads. These can be obtained by The feasible solutions, i.e., those solutions that drove the optimizing also the topology of the neural network. entire route, are shown in Table 2. These results show that feasible solutions multiply the images by 3 or 5 and flip 4. CONCLUSIONS images, while the noise added to output does not influence This paper presented an optimization approach for tuning the results. In addition, the results also show that a lower end-to-end deep learning that builds human-like driving mod- number of epochs is needed if the images are multiplied more times. Table 2: Tuned parameter values of feasible solu- 3.2 Second scenario tions for the first scenario The second scenario was related to a city whose map is Noise added Image Flip Epochs to feasible shown in Figure 3b. Figure 4b shows an example of the to output multiplier image solution city image, which was given as input to the neural network. 0 3 true 30 The learning data were obtained from one driving through 0.1 3 true 30 several crossroads. In contrast to the first scenario, the sec- 0 5 true 18 ond scenario does not predefine the route. Nevertheless, the simulation stops if a distance of more than 2 km has been 0.1 5 true 16 29 Table 3: Parameter values for the second scenario of the IEEE Intelligent Vehicles Symposium IV, pages 814–819, 2011. Parameter Values [2] M. Bojarski, D. Del Testa, D. Dworakowski, B. Firner, batch size {20, 40, ..., 200} B. Flepp, P. Goyal, L. D. Jackel, M. Monfort, number of epochs {10, 20, ..., 50} U. Muller, J. Zhang, X. Zhang, J. Zhao, and K. Zieba. image multiplier {1, 3, 5, 7} End to end learning for self-driving cars. arXiv, noise added to output {0, 0.05, ..., 0.20} 1604.07316, 2016. flip image {true, f alse} [3] Z. Chen and X. Huang. End-to-end learning for lane activation function {linear, elu, relu} keeping of self-driving cars. In Proceedings of the kernel regulizer {none, l2(0.001)} IEEE Intelligent Vehicles Symposium IV, pages 1856–1860, 2017. [4] A. E. Eiben and J. E. Smith. Introduction to Evolutionary Computing (2nd ed.). Springer, Berlin, 2015. [5] H. M. Eraqi, M. N. Moustafa, and J. Honer. End-to-end deep learning for steering autonomous vehicles considering temporal dependencies. CoRR, abs/1710.03804, 2017. [6] T. Fernando, S. Denman, S. Sridharan, and C. Fookes. Going deeper: Autonomous steering with neural memory networks. In Proceedings of the IEEE International Conference on Computer Vision Workshops, pages 214–221, 2017. [7] A. Garcia-Garcia, S. Orts-Escolano, S. Oprea, Figure 5: Length of the feasible route through gen- V. Villena-Martinez, P. Martinez-Gonzalez, and erations for the second scenario. J. Garcia-Rodriguez. A survey on deep learning techniques for image and video semantic segmentation. Applied Soft Computing, 70:41–65, 2018. els. This approach aims at learning good driving models [8] K. Mikami, H. Okuda, S. Taguchi, Y. Tazaki, and when a low quantity of learning data is available. It was eval- T. Suzuki. Model predictive assisting control of vehicle uated with one neural network architecture on two routes: following task based on driver model. In Proceedings of a circular route and a city route. The results show that the IEEE Conference on Control Technology and this approach was able to find driving models for the cir- Applications, pages 890–895, 2010. cular route, but did not manage to find driving models for [9] J. Patterson and A. Gibson. Deep Learning: A handling crossroads inside the city. Practitioner’s Approach. O’Reilly, Sebastopol, 2017. [10] A. Tavčar, B. Kaluža, M. Kvassay, B. Schneider, and Future work will focus on determining the most appropri- M. Gams. Surrogate-agent modeling for improved ate neural network architecture for urban environments. In training. In Proceedings of the Twenty-first European addition, the efficiency of the evolutionary process needs to Conference on Artificial Intelligence (ECAI), pages be increased by, for example, introducing parallelism in the 1103–1104, 2014. model learning. Furthermore, the behavior of the obtained driving models will be compared to human driving behavior [11] Q. Tran and J. Firl. Modelling of traffic situations at to determine how well the models reproduce human driving. urban intersections with probabilistic non-parametric In case of unacceptable reproduction, these models with be regression. In Proceedings of the IEEE Intelligent enhanced to obtain driving models that are able to imitate Vehicles Symposium IV, pages 334–339, 2013. human driving. [12] K. Trontelj, T. Čegovnik, E. Dovgan, and J. Sodnik. Evaluating safe driving behavior in a driving simulator. In Proceedings of the 7th International 5. ACKNOWLEDGMENTS Conference on Information Society and Technology, The authors acknowledge the financial support from the pages 299–302, 2017. Slovenian Research Agency (research core funding No. P2- [13] D. Wei and H. Liu. Analysis of asymmetric driving 0209, and research project “Multiobjective discovery of driv- behavior using a self-learning approach. ing strategies for autonomous vehicles”, funding No. Z2- Transportation Research Part B: Methodological, 7581). This work was also partially funded by NERVteh, 47:1–14, 2013. raziskave in razvoj, d.o.o., and is part of a project that has [14] H. Xu, Y. Gao, F. Yu, and T. Darrell. End-to-end received funding from the European Union’s Horizon 2020 learning of driving models from large-scale video research and innovation programme under grant agreement datasets. In Proceedings of the IEEE Conference on no. 692286. Computer Vision and Pattern Recognition, pages 3530–3538, 2017. 6. REFERENCES [15] L. Xu, J. Hu, H. Jiang, and W. Meng. Establishing [1] P. Angkititrakul, C. Miyajima, and K. Takeda. style-oriented driver models by imitating human Modeling and adaptation of stochastic driver-behavior driving behaviors. IEEE Transactions on Intelligent model with application to car following. In Proceedings Transportation Systems, 16(5):2522–2530, 2015. 30 A Bi-Objective Maintenance-Routing Problem: Service Level Consideration Mohammad Rahimi El-Ghazali Talbi INRIA Laboratory, CRISTAL-CNRS, INRIA Laboratory, CRISTAL-CNRS, 59650 Lille, France 59650 Lille, France Mohammad.Rahimi@inria.fr El-ghazali.Talbi@inria.fr ABSTRACT To the best of our knowledge, few studies attempted to investigate We study a joint maintenance and routing problem and investigate the simultaneous maintenance scheduling and vehicle routing the impact of service level on the optimization of the total expected problem and consider two described assumptions. López-Santana cost. We propose a new bi-objective mathematical model to et al. [6] combine maintenance and routing problems to schedule determine an optimized maintenance-routing policy, maintenance operations for a set of geographically distributed simultaneously. In this model, the first objective function machines and plan to assign a set of technicians to perform minimizes the total costs due to traveling and a delay in start time preventive maintenance at the customer sites. The authors use a of a Preventive Maintenance (PM)/Corrective Maintenance (CM) distribution function for taking into account failures of machines as operation. The second objective function considers the service level an uncertain parameter. In this study, they use two-step iterative which is measured based on waiting times before beginning of the approach to solve the model which causes minimizing the total CM operations. In the proposed model, we consider time windows maintenance and routing cost, waiting time at each customer and in repairing the machines and skill-based technician assignment in failure probability. performing PM/CM operations. The proposed framework is In this study, we propose a new framework to model and to modelled as a mixed-integer linear program and is solved by using establish the trade-off between the service level (measured based the software GAMS. on waiting times before beginning of the CM operations) and Keywords different maintenance costs by taking into account the presented issues. preventive and unforeseen maintenance, vehicle routing problem, scheduling, service level, multi-objective mathematical model 2. PROBLEM DESCRIPTION 1. INTRODUCTION In this section, a bi-objective mathematical model is proposed to determine optimized routing-maintenance policy. In this model, Regularly planned and scheduled maintenance is a critical first objective function minimizes the total costs due to traveling, requirement to reduce the occurrence of an unforeseen failure and delay in start time of a Preventive Maintenance (PM)/Corrective keeping the equipment running at peak efficiency. Maintenance Maintenance (CM) operation at customer while second objective scheduling becomes complex when the machines are function attempts to minimize the waiting times before beginning geographically distributed. In this case, in addition to assigning the of the CM operations. maintenance operations to technicians, it is needed to find the best set of routes for technicians’ visits. In fact, it is necessary to study In this study, we consider a system with geographically distributed the maintenance and the routing scheduling decisions customers, where each customer has one machine that should be simultaneously. Such a joint decision problem is known as the visited and repaired by technician in different cycles. The PM maintenance-routing problems. operations are scheduled with a certain frequency to reduce the occurrence of unforeseen failure in the long term. Regarding the In the literature there are various studies which investigate previous experiences, the time of unforeseen failure occurrence is combination of maintenance and routing problem [1]–[5]. In the known for each machine at each customer, but its repairing can be most of these studies, authors have two initial assumptions: postponed until defined period. The time interval between  The replacement would be done immediately, if an occurrence of unforeseen failure and its repairing named waiting unforeseen failure occurs for the machines. In fact the time. The set of technicians, who need to visit the set of machines authors do not consider waiting time for performing a CM to perform the PM/CM operations to prevent the system failure. operations. While considering the waiting time is important The technician are different in duration time of doing a PM/CM especially where the machines are geographically operation which causes different in salary. A central depot is distributed and the number of technicians and machines are concerned as the point of departure and final destination. Since each limited. technician should travel to perform PM/CM operation at the customer location, the distance between each two customer is  defined. The main aim of this study consist of determining a joint The scheduling is predefined and authors try to assign the routing-maintenance policy for all machines taking into account technicians to machines considering skill of technician, making a balance between the waiting time and total cost of system. time windows and etc. An unforeseen failure causes The optimized maintenance policy determines in which periods the changes of the maintenance scheduling. In this case, PM and CM operations should be performed at each customer. The maintenance scheduling and routing should be done optimized routing policy determines that which technician is simultaneously. assigned to which customers and in which sequence should visit and perform PM and CM operations at each period. 31 The detailed conditions of system are summarized as follows: 2.1 Mathematical Formulation  The time required to perform a maintenance operation The following notations are used in the proposed model. depends on the skill of the assigned technician. Sets M set of customers, index for customers (1,2,…, m)  More skilled technicians receive more salary. M’ set of customers and central depot, (0,1,2,…, m+1) K index for technicians (1,2,…, k)  All technicians are able to perform any PM/CM operation. t, t’, t’’ index for period (1,2,…, T) Parameters  The technicians start in the central depot in the beginning of ck one unit time cost of a PM/CM operation by technician k each period and should return to the central depot by the end pm of the period. k time required to perform a PM operation by technician k time required to perform a CM operation by technician cmk k  Each machine should be repaired by only one technician at duration of the interval between two consecutive PM each period. It means if the machine should be repaired in λ operations the specific period, only one technician should be assigned L allowed duration to repair occurred unforeseen failure to the machine. zit a binary parameter which determines occurrence of unforeseen failure in customer i at period t  The PM operation should be performed on all the machines tjj traveling time between customer i and j at the first period. r transportation cost per unit time earliest and latest possible start time of a PM/CM [ai, bi]  If no unforeseen failure occurs on the machine at planning operation at customer i horizon, the PM operations will be performed regarding the pi penalty cost of one unit time delay due to start time of a defined frequency. The frequency is defined regarding PM/CM operation at customer i after latest possible planning horizon and the duration of the interval between arrival time two consecutive PM operations. G a large value number Variables  In the case of unforeseen failure occurrence on the machine, 1 if customer j is visited exactly after customer i by x no predictive maintenance can be scheduled and performed ijkt technician k at period t, otherwise 0 before performing CM operation. In this case, CM 1 if PM operation is planned in customer i at period t, y operation should be scheduled to assign a technician on the it otherwise 0 machine until maximum L period. Moreover, next PM 1 if a CM operation is planned in customer i at period t operation will be scheduled and performed after λ period. uit’t for the an occurred unforeseen failure at period t’, otherwise 0  After performing a CM operation, the machine returns to 1 if delay occurred in visiting customer i at period t, βitt’ the good condition and no unforeseen failure occurs until otherwise 0 the next repairing that will be a PM operation in λ period. It 1 if customer i is visited by technician k at period t to µikt means two unforeseen failure cannot occur consequently. perform a PM operation, otherwise 0 1 if customer i is visited by technician k at period t to πikt  The time required to perform a CM operation is longer than perform a CM operation, otherwise 0 the time required to perform a PM operation on each start time of an operation by technician k in customer I, Tikt machine. period t delay in start time of a PM/CM operation in customer i at dit period t  The CM cost is larger than the PM cost.  The mathematical model associated with the presented framework The machines impose time windows on the system which is provided in this section. Each equation in this model is detailed means the technician should start maintenance operation below. before the latest possible start time. In cases where this time windows is not respected, a delay penalty applies if the Min f   x .t .r d .p  .c .pm technician starts after the latest allowed time. 1 ijkt ij it i ikt k k i, j ,k ,t i,t i,k,t  The travel time between two customers depends on the    (1) .c .cm ikt k k speed of the vehicle in the rout at each period. i,k,t Min f    2 (2) , , itt i t t   S. t.  T 1  y   1 i   M (3) t it    32 y  1 z i  M,t performing a CM operation return the machine to as good as new (4) it it condition again and no PM operation is needed until next λ periods. Constraint (7) ensures that when a PM operation is performed at the y  1 i  M i1 (5) period t and no unforeseen failure occurs on the machine until the next λ periods, then a PM operation should be scheduled and t  1   y    i  M,t,t performed at the period of t+λ. Constraint (8) checks that when a   1 (y u  ) (6) t t 1  it it it t CM operation is performed, then a PM operation can be scheduled at the interval of λ periods or an unforeseen failure can be occurred t y  y    z i  M,t (7) until next λ periods. Equation (9) determines in which period a CM it i (t  ) tt it operation should be scheduled and performed to repair the occurred unforeseen maintenance. Moreover, this equation checks that CM t u    i  M,t  y  z (8) operation should be scheduled in a way to assign a technician on it t i(t  ) tt 1  it the machine until maximum L period after the failure. Equation t L  (10) calculates the waiting time until performing a CM operation in u  i  M ,t (9)   z  it t it t t  the case where an unforeseen failure occurs. Equation (11) ensures in case of unforeseen failure occurrence, the CM operation should     be performed once. Constraint (12) guarantees that CM operation  (t t).u i  M,t,t itt itt (10) and PM operation cannot be scheduled and performed for the same  u  period, simultaneously. Equations (13) and (14) determine that  1 i  M ,t (11) t itt visiting the customer is related to a PM operation or a CM operation. y  u   1 i  M ,t,t (12) it it t Equation (15) makes a connection between routing and  maintenance variables. This equation checks when a PM/CM   y i  M,t (13) should be performed on the machine, a technician should be k ikt it assigned to the machine.     u i  M,t (14) Constraint (16) guarantees that only one technician should be k ikt t it t assigned on the each machine at each period. Constraint (17) M  1   x     j  M, j  ,ik,t ensures that technician leave the current customer to the next one, (15) i M  , i0 ijkt jkt jkt after finish the PM/CM operation. Equation (18) determines the start time on the machine, which is calculated as the start time of M 1   x  1 j   M , k,t the immediate previous customer, increased by the PM/CM (16) j M  , j 1  0 jkt operation time and the traveling time between the two customers. Equation (19) checks the time windows constraint and calculates M  1 M    x   x  0 j  M,k,t (17) the delay. Finally, (20) and (21) impose bounds on the variables. i M  , i0 ijkt i M  , i2 jikt 3. RESOLUTION METHOD T   .pm  .cm  t In this section we firstly introduce the instance generation method ikt ikt    k ikt k ij i, j M (18) and solution procedure, briefly. Then, a numerical analysis is  T  . G (1 x ) k , t presented which derives managerial results. jkt ijkt Problem instances have been generated by a random generator. In M  1 M   1  i  M  this way, parameters of the problem are generated using random a .  x  T  b .  x  d (19) numbers by a discrete uniform distribution. Then, to solve the i jikt ikt i jikt it  problem, we use the weighted sum method [7]. Under this method, j M  , j 0 j M     , j 0 k , t the problem is solved by considering each objective function x , y , u     separately in both the maximization and the minimization for  ,   ,  , 0,  1 i, j M , k , t (20) ijkt it it t it t ikt ikt finding extreme points of each objective function. Then, a new single objective is considered that aims to minimize the weighted T , d  0 i  M , k , t (21) sum over the normalized and non-dimensional objective function. ikt it In order to show the feasibility and applicability of proposed model, The first objective function (1) minimizes the total cost which a small size problem is generated and it is solved based on consist of traveling cost between customers, penalty cost due to generated instance problem. It is assumed that there are 6 customers start time out of time windows and the wages of technicians for (m=6) where 3 periods are defined as duration of the interval PM/CM operations. The second objective function (2) optimizes between two consecutive PM operations (λ=3) and 2 periods are the customer satisfaction level by minimizing the waiting times considered as allowed duration to repair occurred unforeseen until performing a CM operation in the case where an unforeseen failure (L=2) by using 3 types of technicians (k=3) during 10 failure occurs. periods. To solve this problem, the “GAMS v22.2” optimization Constraint (3) checks number of PM operations on the machine of software using solver CPLEX v10.1 is used. each customer should not be exceeded. Constraint (4) guarantees At the beginning, the problem is solved without considering the that if the unforeseen failure occurred, then the PM cannot be second objective. In this case the total cost is optimised. The results scheduled and performed for the same period. Equation (5) show that minimum value of total cost is 402 while waiting time in determines that at the first period, PM operation should be this situation is 14. In the next step, the first objective function is performed on the all the machines. Equation (6) guarantees that relaxed and model is solved by minimizing the second objective 33 function. The obtained results shows when the minimum value of 4. CONCLUSION waiting time (second objective function) is 6, the value of total cost In this paper the integration of maintenance and routing problem is is 1,053. Table 1, shows the minimum and maximum value of investigated by taking into account waiting time for performing a objective functions. CM operation when unforeseen failure occurs. For this Purpose, a Table 1: Min. and Max. value of objective functions bi-objective mathematical model is proposed to find the optimized Minimum value Maximum Value policy of maintenance and routing problems and make a trade-off First objective 402 1,053 between maintenance costs and service level which is measured by Second objective 6 14 waiting time for performing a CM operation. In the proposed model The bi-objective model can be converted to a MILP model with one the time windows is considered for starting maintenance operation objective function using the equation (22). on the machine by technicians. Moreover, the technician’s skill regarding required time to perform a maintenance operation is min min f  f f  f 1 1 f    1  2 2 considered. Our results for a small size instance show that to  max min max min f  f f  f decrease by 57% of the waiting time, we have to increase the costs 1 1 2 2 by 62%. In this equation, α presents the importance degree of each objective function and varies between 0 and 1. Our future research in this area includes the consideration of stochastic parameters and proposing an efficient solution approach. Furthermore, the Objective Functions Value (OFV) by changing α value is introduced in Table 2. 5. ACKNOWLEDGMENTS Table 2: OFV against α This project is supported by The ELSAT2020 project which is co- α financed by the European Union with the European Regional 0 0.3 0.5 0.7 1 Development Fund, the French state and the Hauts de France First OFV 1,053 886 689 602 402 Region Council. Second OFV 6 8 9 12 14 6. REFERENCES Run time (second) 157 166 147 133 158 [1] R. Macedo, R. Benmansour, A. Artiba, N. Mladenović, and D. Urošević, “Scheduling preventive railway maintenance According Table 2, the total cost from 1,053 to 402 causes activities with resource constraints,” Electron. Notes increasing 57% in waiting time (from 6 to 14). It means the best Discret. Math., vol. 58, pp. 215–222, Apr. 2017. value of waiting time can be reached by increasing 62% in total [2] F. Blakeley, B. Argüello, B. Cao, W. Hall, and J. cost. Knolmajer, “Optimizing Periodic Maintenance Operations for Schindler Elevator Corporation,” Interfaces (Providence)., vol. 33, no. 1, pp. 67–79, Feb. 2003. [3] A. A. Kovacs, S. N. Parragh, K. F. Doerner, and R. F. Hartl, “Adaptive large neighborhood search for service technician routing and scheduling problems,” J. Sched., vol. 15, no. 5, pp. 579–600, Oct. 2012. [4] A. Goel and F. Meisel, “Workforce routing and scheduling for electricity network maintenance with downtime minimization,” Eur. J. Oper. Res., vol. 231, no. 1, pp. 210– 228, Nov. 2013. [5] M. D. Berrade, P. A. Scarf, C. A. V. Cavalcante, and R. A. Dwight, “Imperfect inspection and replacement of a system  with a defective state: A cost and reliability analysis,” Figure 1: Variation of total cost against waiting time by Reliab. Eng. Syst. Saf., vol. 120, pp. 80–87, Dec. 2013. changing value of α [6] E. López-Santana, R. Akhavan-Tabatabaei, L. Dieulle, N. Labadie, and A. L. Medaglia, “On the combined The variation of objective functions value by changing of α value maintenance and routing optimization problem,” Reliab. is presented in Figure 1. In this figure, X-axe shows value of total Eng. Syst. Saf., vol. 145, pp. 199–214, Jan. 2016. cost and waiting time while Y-axe presents different value of α. By [7] M. Rahimi, A. Baboli, and Y. Rekik, “A bi-objective this figure changes of total costs and waiting time is visualized inventory routing problem by considering customer against variation of α. satisfaction level in context of perishable product,” in 2014 IEEE Symposium on Computational Intelligence in Production and Logistics Systems (CIPLS), 2014, pp. 91– 97. 34 Study on Reducing Turn-Around Time of Multi-Objective Evolutionary Algorithm on an Industrial Problem Hiroaki Fukumoto Akira Oyama Institute of Space and Astronautical Science, Institute of Space and Astronautical Science, JAXA JAXA 3-1-1, Yoshinodai, Chuo, Sagamihara 3-1-1, Yoshinodai, Chuo, Sagamihara Kanagawa, Japan Kanagawa, Japan fukumoto@flab.isas.jaxa.jp oyama@flab.isas.jaxa.jp ABSTRACT scale many-core architectures [2] and the computational al- Multi-objective evolutionary algorithms (MOEAs) are pop- gorithms, say MOEAs, should utilize the large-scale compu- ulation based global optimization algorithms and it is said tational resources efficiently. One of the simple yet effective that the performance of the MOEAs depends on the pop- ways of MOEAs for utilizing the large-scale computational ulation size. Considering that the recent trends of com- resources would be to increase the number of concurrent so- puter development is in large-scale many-core architectures, lution evaluations, i.e., the population size. Note that the and massive parallel computation is getting feasible in more increased number of objectives of multi-objective optimiza- companies and laboratories, the available population size tion problems also gives a reason to increase the population is increasing and the efficiency of MOEA with large pop- size: the necessary number of solutions to cover the entire ulation size should be enhanced. This study examines the Pareto front exponentially increases as the number of objec- effect of the population size on MOEAs’ performance on tives increases [9, 18]. Therefore, the increase in the popu- a real-world-derived benchmarking optimization problem, lation size would be the right direction for recent MOEAs. with large population size. In this paper, three mate se- lection schemes with different degree of elitist strategy are This study aims to reduce the turn-around time of MOEAs adapted to NSGA-II-M2M. The experimental results show when large population size is used. This paper demonstrates that the elitist strategy can efficiently make use of the ef- the population size effect on the performance of an MOEA fect of the large population size, therefore can reduce the on a real-world-derived benchmarking problem and the re- turn-around time. duction of the turn-around time by making use of the pop- ulation size effect is attempted. This paper is organized Keywords as follows. In Section 2, the experimental settings are ex- multi-objective optimization, large population size, mate se- plained first, and the results demonstrating the impact of lection, real-world problem the population size on the performance of the MOEAs is presented. Then the method to reduce the turn-around time 1. INTRODUCTION is described and the experimental results are provided. Sec- tion 3 concludes this paper. Many of industrial design problems involve multiple objec- tives and constraints and they are so-called constrained multi- objective optimization problems. Considering that creating 2. REDUCTION OF TURN-AROUND TIME high value-added products in industries is getting more and more important along with the increase of the sophistica- 2.1 Experimental setting tion and diversity of social needs, it is very important to • Problem: The Mazda CdMOBP problem [11]. This catch up to the changes in customer demands and so short problem has two objectives, 54 constraints, and 222 development time of each product is highly appreciated. variables. The problem originates from an actual de- sign optimization of car models and the constraints For multi-objective optimization problems, multi-objective comprise the requirements for crashworthiness, body evolutionary algorithms (MOEAs) have been regarded as torsional stiffness, and low frequency vibration modes. a promising approach. With respect to the application of These constraints are evaluated by finite element sim- MOEAs to industrial design problems, the development time ulations on a supercomputer in actual design process, of the products corresponds to the turn-around time for however, in the benchmark problem these simulation MOEAs. The turn-around time of MOEAs corresponds to results are modeled with radial basis functions so as to the number of generations in MOEA, supposing that the shorten the evaluation time while retaining the nonlin- runtime for MOEA itself is negligible compared with the earity as much as possible. The details are presented runtime for solution evaluations. Here the turn-around time in [11] and the problem is available from the website is the time from the beginning of the optimization to the [12]. end of the optimization when a desired quality of solution set is obtained. • MOEA: NSGA-II-M2M [15] with the subproblem size of 10. The probability that the parents are chosen from One of the recent trends of computer development is in large- the corresponding subproblem δ is set to 0 . 9. 35 • Constraint handling technic: Multiple constraint rank- 0.40 ing (MCR) [5], which generally performs well on con- 0.35 strained optimization benchmarking problems [7]. The constraint handling technic is incorporated into NSGA- 0.30 II-M2M with the MOEA-CHT incorporation frame- 0.25 work [7]. N = 100 0.20 • Mate selection schemes: Random selection, binary tour- N = 316 nament (BT) [1], or Elitist BT (EBT, explained in next 0.15 N = 1000 subsection) [8]. The random selection scheme is the 0.10 N = 3162 default mate selection scheme for NSGA-II-M2M [15] Mean hypervolume values and its modified version of Jain et al. [10] is employed 0.05 so as to handle constraints. 0 0.00 60 120 180 240 300 • Reproduction operators: The crossover and mutation Number of generations operators with the same control parameter values as (a) Random selection in [14, 15]. • Direction vector generation method: Das and Dennis’s Figure 1: Convergence history of the HV values with systematic approach [4]. various population sizes for the case with random se- • Stopping criterion: The number of generations of 300. lection. The mean and the standard deviation values The number of fitness evaluations differs according to are plotted. the population size at a given generation, but the focus is in this study is on the reduction of the required generations, and so the differences in the number of 2.4 Reduction of turn-around time by enhanc- fitness evaluations is not considered in this study. ing the population size effect • Independent runs: Each case run for 31 times indepen- The population size affects the diversity of the solutions and dently. the convergence speed, and now it is commonly accepted • The population size N: N is set to be the numbers in that the population size should be large enough to guarantee a geometric progression with a scale factor of 100, and √ the diversity of the solutions while the large population size a common ratio of 10 is used to see the population makes the convergence slow [16, 3, 17, 19]. size effect. Specifically, the population sizes of 100, 316, 1000, and 3162 are used, especially for drawing Considering that the phenomenon of the population size ef- Figure 3. fect can be explained by a term “selection pressure” [19], we attempt to mitigate the slow convergence with large popu- 2.2 Performance Metric lation by somehow strengthening the selection pressure. In this study, a standard and popular mate selection of BT and The hypervolume (HV) indicator [20] is used as the perfor- a recently proposed mate selection scheme with a strong eli- mance indicator. In this study, the solution set used for the tist strategy named EBT [8], both of which have stronger calculation of the HV value is the solutions not only in the selection pressure than the random selection, are employed. final population but also in the external unbounded archive In EBT, i) the usual BT selection is conducted at first for [13], considering that the designers in actual industries use all the solutions in each subproblem then the indices of the MOEAs as design support tools for decision making and so selected solutions are sorted according the number of times use of the unbounded external archive is more practical than the each index is selected. Apart from that, ii) the indices of use of the solutions obtained only at the final generation. For the solutions are also sorted according the solutions’ fitness. calculating the HV value for a generation, non-dominated so- Finally, every sorted indices i) is replaced by the index in ii) lutions are extracted from all the feasible solutions obtained with the same rank order with i), so that the solution with by the generation and are used to calculate the HV value. higher rank is selected more. For further details of EBT, For the details of the formulation for the HV calculation, please refer to [8]. please refer to [11]. The larger the HV value, the better the approximation to the Pareto front. The most elitist is EBT, followed in order by BT and random selection. 2.3 Impact of the population size on the per- formance It must be noted that the strong elitist strategy tends to Figure 1 presents the convergence history of the mean HV deteriorate the diversity of the solutions, and the negative values with various population sizes for NSGA-II-M2M with effect of the strong elitist strategy should be compensated random selection. It is observed that the cases with higher by using some diversity-enhancing method. In this study, we population size show generally higher mean and smaller stan- enhance the MOEAs’ capability of keeping diversity by em- dard deviation values. This result supports the motivation ploying M2M, and this is the reason why the base algorithm for increasing the population size, however, the effect of the in this study is not NSGA-II [6] but NSGA-II-M2M. increased population size is not clearly observed until around the number of generations of 200, between the cases with the Figure 2 shows the convergence history of the mean HV val- population size of 1000 and 3162. ues with various population sizes for NSGA-II-M2M with 36 0.40 0.40 0.35 0.35 0.30 0.30 0.25 0.25 N = 100 0.20 0.20 N = 316 0.15 0.15 N = 1000 0.10 0.10 N = 3162 Mean hypervolume values 0.05 Mean hypervolume values 0.05 0 0.00 60 120 180 240 300 0.00 Number of generations 0 60 120 180 240 300 Number of generations (a) BT (b) EBT Figure 2: Convergence history of the HV values with various population sizes for the cases with BT and EBT selection. The mean and the standard deviation values are plotted. BT and EBT. Comparing Figure 1 and 2, it can be observed will include further improvement of the population size ef- that the strong elitist mate selection enhances the large pop- fect, even with much smaller population size. ulation size effect and the differences in the mean HV values can be observed more clearly and from earlier generations. Acknowledgment This research is supported by the HPCI System Research With regard to the reduction of the turn-around time, Fig- Project “Research and development of multi-objective de- ure 3 shows the generation that is required to attain a HV sign exploration and high-performance computing technolo- value against the population size. For example, in Figure gies for design innovation”. 3, the HV value of 0 . 2 can be attained with the number of (Project ID: hp160203 and hp170238) generations of approximately 300 with the population size of approximately 300, and with the number of generations 4. REFERENCES of approximately 160 with the population size of approx- [1] T. Blickle and L. Thiele. A mathematical analysis of imately 1000. The subfigures in Figure 3 show that the tournament selection. In Proceedings of the Sixth required generation to attain a certain HV value is reduced International Conference on Genetic Algorithms, with stronger mate selection scheme. pages 9–16. Morgan Kaufmann, 1995. Compared with the case with BT, the results for EBT shows [2] S. Borkar. Thousand core chips: A technology relatively poor performance with small population sizes, and perspective. In Proceedings of the 44th Annual Design so further development of more robust and more efficient al- Automation Conference, DAC ’07, pages 746–749, gorithm for reducing the turn-around time will be required. New York, NY, USA, 2007. ACM. [3] T. Chen, K. Tang, G. Chen, and X. Yao. A large population size can be unhelpful in evolutionary 3. CONCLUSIONS algorithms. Theoretical and Computational Science, 436:54–70, 2012. This paper demonstrates the population size effect on the performance of an MOEA on a real-world-derived bench- [4] I. Das and J. E. Dennis. Normal-boundary marking problem (Mazda CdMOBP) and the reduction of intersection: A new method for generating the pareto the turn-around time by making use of the population size surface in nonlinear multicriteria optimization effect is attempted. problems. SIAM Journal on Optimization, 8(3):631–657, 1998. By the demonstration of the population size effect, it is [5] R. de Paula Garcia, B. S. L. P. de Lima, A. C. shown that the large population size can improve the perfor- de Castro Lemonge, and B. P. Jacob. A rank-based mance of an MOEA, and it is also shown that the population constraint handling technique for engineering design size effect is not clearly shown until late stage of the evolu- optimization problems solved by genetic algorithms. tion with random mate selection scheme. Computers & Structures, 187(Supplement C):77 – 87, 2017. The late-appearing population size effect is then improved [6] K. Deb, S. Agrawal, A. Pratap, and T. Meyarivan. A by employing two techniques: a strong mate selection scheme fast and elitist multiobjective genetic algorithm: and its complementary scheme of M2M. The results show NSGA-II. IEEE Transactions on Evolutionary that the case with stronger elitist strategy exhibits relatively Computation, 6(2):182–197, 2002. faster large population size effect and the aim of reducing [7] H. Fukumoto and A. Oyama. A generic framework for turn-around time is achieved in some degree. Future work incorporating constraint handling techniques into 37 (a) Random selection (b) BT (c) EBT Figure 3: Plot for the generation that is required to attain a HV value against population size. multi-objective evolutionary algorithms. In K. Sim sub-regional search. In 2009 IEEE Congress on and P. Kaufmann, editors, Applications of Evolutionary Computation, pages 1928–1934, 2009. Evolutionary Computation, pages 634–649, Cham, [15] H. L. Liu, F. Gu, and Q. Zhang. Decomposition of a 2018. Springer International Publishing. multiobjective optimization problem into a number of [8] H. Fukumoto and A. Oyama. Study on improving simple multiobjective subproblems. IEEE efficiency of multi-objective evolutionary algorithm Transactions on Evolutionary Computation, with large population by m2m decomposition and 18(3):450–455, 2014. elitist mate selection scheme. In 2018 IEEE [16] R. Mallipeddi and P. N. Suganthan. Empirical study Symposium Series on Computational Intelligence on the effect of population size on differential (SSCI), 2018, to appear. evolution algorithm. In 2008 IEEE Congress on [9] H. Ishibuchi, Y. Sakane, N. Tsukamoto, and Evolutionary Computation (IEEE World Congress on Y. Nojima. Evolutionary Many-Objective Computational Intelligence), pages 3663–3670, 2008. Optimization by NSGA-II and MOEA/D with Large [17] D. Mora-Meliá, F. J. Martinez-Solano, P. L. Populations. In IEEE International Conference on Iglesias-Rey, and J. H. Gutiérrez-Bahamondes. Systems, Man and Cybernetics, pages 1758–1763, Population size influence on the efficiency of 2009. evolutionary algorithms to design water networks. [10] H. Jain and K. Deb. An Evolutionary Many-Objective Procedia Engineering, 186:341 – 348, 2017. XVIII Optimization Algorithm Using Reference-Point Based International Conference on Water Distribution Nondominated Sorting Approach, Part II: Handling Systems, WDSA2016. Constraints and Extending to an Adaptive Approach. [18] T. Tatsukawa, T. Watanabe, and A. Oyama. IEEE Transactions on Evolutionary Computation, Evolutionary computation for many-objective 18(4):602–622, 2014. optimization problems using massive population sizes [11] T. Kohira, H. Kemmotsu, O. Akira, and on the k supercomputer. In 2016 IEEE Congress on T. Tatsukawa. Proposal of benchmark problem based Evolutionary Computation (CEC), pages 1139–1148, on real-world car structure design optimization. In 2016. Proceedings of the Genetic and Evolutionary [19] T. Weise, Y. Wu, R. Chiong, K. Tang, and J. Lässig. Computation Conference Companion, GECCO ’18, Global versus local search: The impact of population pages 183–184. ACM, 2018. sizes on evolutionary algorithm performance. Journal [12] T. Kohira, H. Kemmotsu, A. Oyama, and of Global Optimization, 66(3):511–534, 2016. T. Tatsukawa. Mazda/JAXA/TUS simultaneous car [20] E. Zitzler, L. Thiele, M. Laumanns, C. M. Fonseca, structure design problem, and V. G. da Fonseca. Performance assessment of http://ladse.eng.isas.jaxa.jp/benchmark/index.html, multiobjective optimizers: an analysis and review. 2018. IEEE Transactions on Evolutionary Computation, [13] O. Krause, T. Glasmachers, N. Hansen, and C. Igel. 7(2):117–132, 2003. Unbounded population mo-cma-es for the bi-objective bbob test suite. In Proceedings of the 2016 on Genetic and Evolutionary Computation Conference Companion, GECCO ’16 Companion, pages 1177–1184. ACM, 2016. [14] H. l. Liu and X. Li. The multiobjective evolutionary algorithm based on determined weight and 38 Evolution of Electric Motor Design Approaches: The Domel Case Gregor Papa Gašper Petelin Peter Korošec Computer Systems Faculty of Computer and Computer Systems Department Information Science Department Jožef Stefan Institute University of Ljubljana Jožef Stefan Institute Ljubljana, Slovenia Ljubljana, Slovenia Ljubljana, Slovenia gregor.papa@ijs.si peter.korosec@ijs.si ABSTRACT Slovenian smart specialisation strategy [2], it is planned to The paper presents the evolution of geometry design ap- transfer this solution into Slovenian industry. proaches in the optimisation of an electric motor, more specif- ically its rotor and stator. It starts with the initial manual The rest of the paper is organized as follows: Section 2 approach, which was replaced with the automatic approach briefly describes the geometry elements of an electric mo- that introduced evolutionary algorithms to allow the intel- tor and the optimisation goal; Section 3 presents the con- ligent search in collaboration with evaluation tools. Next, ventional manual approach to the motor design; in Section 4 the new platform for remote optimisation was recently in- the use of evalutionary algorithms in electric motor design is troduced that allows remote optimisation with various algo- outlined; Section 5 introduces the new developed platform rithmic approaches, including multi-objective optimisation. for remote optimisation; and Section 6 draws conclusions At the end we propose further solutions that will improve and proposes possible future work. high performance of the design process. Keywords electric motor, design, evolution, high-performance 1. INTRODUCTION Many widely-used home appliances (e.g., mixers, vacuum cleaners, drills, etc.) use electric motors. These small motors are required to have high power and provide high starting and running torques, despite their small sizes. While having sufficient output power they should be energy efficient and inexpensive to manufacture [12]. There is a number of past works addressing the geometry optimisation design of rotor and stator parts [6], [10], [12], electric motor casing [7] and impeller [4]. These works, per- formed on various products of Domel company [1], intro- duced various artificial intelligence methods to implement automatic search of an optimal design. The reported optimi- sation approaches were mostly single objective. Still, there were some initial steps identified towards multi-objective handling of the design process. This paper focuses on the approaches for automatic optimi- sation of the electric motor geometry. The main parts of the electric motor, i.e., stator and rotor, are presented in Figure 1. Figure 1: Rotor and stator of an electric motor [10]. While improving the applicability of the multi-objective op- timisation, supported by parallelisation and surrogate mod- 2. PROBLEM DESCRIPTION elling through the support of the Horizon 2020 Twinning The rotor and the stator of an electric motor are constructed project SYNERGY - Synergy for smart multi-objective op- by stacking the iron laminations. The shape of these (rotor timisation [3], we implemented a platform for an efficient and stator) laminations is described by several geometry pa- optimisation with different methods and approaches. The rameters that define the rotor and stator in two dimensions platform is briefly presented in this paper. In line with (2D). 39 The whole set of geometry parameters consists of invariable and variable ones. Invariable parameters are fixed, as they cannot be altered, either for technical reasons (e.g., the air rer gap) or because of the physical constraints on the motor (e.g., the radius of the rotor’s shaft). Variable parameters, on the other side, do not have predefined optimal values. ryt Among these parameters, some are dependent (upon oth- ers variables), while some variable parameters are mutually independent and without any constraints. The mutually in- rpw dependent set of variable parameters of the rotor and stator a) geometry (see details in Figure 2) can be subject to optimi- sw sation: syh ssr • rotor yoke thickness (ryt), • rotor external radius (rer), • rotor pole width (rpw). • stator width (sw), sml • stator yoke horizontal thickness (syh), str • stator yoke vertical thickness (syv), • stator middle part length (sml), • stator internal edge radius (sie), • stator teeth radius (str), sie • stator slot radius (ssr). syv b) One of the optimisation tasks is to find the values of geom- etry parameters that would generate the rotor and stator Figure 2: Geometry parameters of a) rotor and b) geometry with minimum power losses. stator [12]. 2.1 Mathematical formulation of the problem The efficiency of an electric motor is defined as the ratio of two components: the eddy-current losses and the hysteresis the output power to the input power. It depends on various losses: power losses (see details in [9]), which include: • Copper losses: the joule losses in the windings of the PF e = keB2f 2 rotmrot+keB2f 2 statmstat+khB2f 2 statmstat (2) stator and the rotor. • Iron losses: including the hysteresis losses and the where ke is an eddy-current material constant of 50 Hz, kh is eddy-current losses, which are primarily in the arma- a hysteresis material constant of 50 Hz, B is the maximum ture core and in the saturated parts of the stator core. magnetic flux density, f is the frequency, and m is the mass. • Other losses: brush losses, ventilation losses and fric- Three additional types of losses also occur, i.e., brush losses tion losses. PBrush, ventilation losses PV ent, and friction losses PF rict. The output power P The overall copper losses (in all stator and rotor slots) are 2 of the motor is a product of the elec- tromagnetic torque T , and the angular velocity ω, as follows: X P P 2 = T ω (3) Cu = (J 2Aρlturn)i (1) i where ω is set by the motor’s speed, and T is a vector prod- where i stands for each slot, J is the current density, A is uct of the distance from the origin r, and the electromagnetic the slot area, ρ is the copper’s specific resistance and l force F . turn is the length of the winding turn. The overall efficiency of an electric motor is defined as: Due of the non-linear magnetic characteristic, the calcula- P2 tion of the iron losses is less exact; they are separated into η = (4) P2 + PCu + PF e + PBrush + PV ent + PF rict 40 2.2 Fitness evaluation 3. For evaluation of each solution (i.e., their fitness) each Each solution candidate of the population was decoded into geometrical configuration is analyzed using some FEM a set of the rotor and stator parameters. The fitness was program (e.g., ANSYS). This step requires a decoding estimated by performing a finite-element numerical simula- of the encoded parameters into a set of geometrical tion to calculate the iron and the copper power losses (using parameters that define the rotor and the stator. the above mentioned equations). The sum of power loses 4. After the fitness calculation, the reproduction of the corresponds to the solution’s fitness. individual solutions is performed and the application of various recombination operators to a new population For multi-objective version we can also introduce additional are done. objective like material costs, making it a typical price/ per- formance optimisation. The cost is calculated by taking into 5. The evolutionary algorithm repeats the above proce- account the amount of materials (i.e., iron and copper), that dure until some predefined number of iterations have are used to produce the electric motor, and their correspond- been accomplished or some other stopping criteria is ing prices. met. 3. MANUAL OPTIMISATION A manual design procedure of an electric motor consists of Some evident advantages of this approach are: the geometry estimation of the rotor and the stator of an electric motor by an experienced engineer. The suitability • There is no need for an experienced engineer to be of the proposed geometry is usually analyzed by means of present during the whole process. He is required only numerical simulation (e.g., FEM with an automatic finite- at the beginning to decide on the initial design. element-mesh generation) of the electromagnetic field of each proposed solution separately. • There is no need to know the mechanical and physical details of the problem. The problem can be solved, by The manual procedure can be repeated until the satisfied the use of optimisation algorithm, irrespective of any evaluation results is obtained. Similarly, the conventional knowledge about the problem. approach in most new designs starts with manual design, as there exist no prior design. Some possible drawbacks of this approach can appear: The advantage of the manual approach is that the engineers can significantly influence the progress of the design pro- • The improper use of recombination operators leads to cess with their experiences and react intelligently to any slow search progress. noticeable electromagnetic response with proper geometry redesign. • An initial solution set that is not divergent enough, can lead to a longer convergence time. The drawback of this approach is that an experienced en- gineer and a large amount of time (that is mostly spent on 5. REMOTE OPTIMISATION PLATFORM computation) are needed. The multi-objective optimisation is a natural approach to solve difficult real-world problems. As the presented elec- 4. AUTOMATIC OPTIMISATION tric motor geometry design can have several contradictory The above-described manual design approach can be up- constraints, it is useful to introduce the multi-objective al- graded with one of the stochastic optimisation techniques gorithms (e.g., NSGA-II, IBEA) into this process [11]. which, in connection with reliable numerical simulators, al- low for highly automated design process where the need for Within the project SYNERGY, we developed and imple- an experienced engineer to navigate the process is signifi- mented a platform for an efficient optimisation with different cantly reduced. methods and approaches. Its main role is to allow compari- son and testing of an effective optimisation methods for the So far, several evolutionary approaches have already been optimisation of electric motor geometry. The platform al- proven to be efficient in the process of the electric motor ge- lows comparison of single objective as well as multi-objective ometry optimisation; e.g., electromagnetism-like algorithm algorthms. [5], multi-level ant-stigmergy algorithm [6], adaptive evolu- tionary search algorithm [8], genetic algorithm [9], particle The platform is based on web-based services to allow remote swarm optimization, and differential evolution [12]. work of different experts, while keeping some important, se- cret features and characteristics hidden. The remote tool The automatic approach with the use of an evolutionary also allows for parallel processing, which allows for fast cal- algorithm can be summarized into the following steps: culations, without any intervention from the expert. Remote access enables experts to use the evaluation of the 1. The initial set of solutions is defined according to an proposed solution regardless of his location. The platform initial electric motor. allows remote access towards any simulation tools (e.g., FEM 2. It provides a set of problem solutions (i.e., different analysis). Furthermore, all evaluations are being stored in configurations of the mutually independent geometri- database and in case the same solution is being put to eval- cal parameters of the rotor and the stator). uation, the result is immediately returned without the need 41 to wait for it to be actually evaluated again, which furthers [2] Slovenia’s Smart Specialisation Strategy S4. speed up the evaluation process. http://www.svrk.gov.si/fileadmin/svrk.gov.si/ pageuploads/Dokumenti za objavo na vstopni strani/ Since actual parameter values are not relevant for optimisa- S4 strategija V Dec17.pdf. Accessed: 2018-08-26. tion process and to ensure that no secrets about the problem [3] SYNERGY – Synergy for Smart Multi-Objective are being shared, the platform hides the important prop- Optimisation. http://synergy-twinning.eu/. Accessed: erties of the solutions. Meaning all parameter values and 2018-08-26. evaluation results are being normalised within the interval [4] P. Korošec, J. Šilc, K. Oblak, and F. Kosel. [0.0, 1.0]. This way, the problem can be tackled by any op- Optimizing the shape of an impeller using the timisation expert without acquiring any relevant knowledge differential ant-stigmergy algorithm. In (e.g., actual dimensions, problem specifications) about the R. Wyrzykowski, J. Dongarra, K. Karczewski, and problem. J. Wasniewski, editors, Parallel Processing and Applied Mathematics, pages 520–529, Berlin, Parallelisation within the platform is considered on the level Heidelberg, 2008. Springer. of solution evaluation. Any other parallelisation on the level [5] P. Korošec, G. Papa, and J. Šilc. Optimization of optimisation algorithm is left to the optimisation expert. algorithms inspired by electromagnetism and stigmergy in electro-technical engineering. WSEAS 6. CONCLUSION Transactions on Information Science and This paper presented the evolution of approaches to the Applications, 5(2):587–591, May 2005. optimisation of the geometry design of the electric motor. [6] P. Korošec and J. Šilc. The distributed multilevel From the initial manual approach, through the automatic ant-stigmergy algorithm used at the electric-motor approach that uses some evolutionary algorithm combined design. Engineering Applications of Artificial with evaluation tools, towards the platform that allows re- Intelligence, 21(6):941–951, 2008. mote optimisation with various algorithms. The latter al- [7] K. Oblak, P. Korošec, J. Šilc, and F. Kosel. Stigmergic lows simple comparison and study of different methodologies optimization of plane constructions. Electrotehniški and algorithms. vestnik, 74(5):279–284, 2007. In Slovenian. [8] G. Papa. Parameter-less algorithm for In the future version of the optimisation platform we plan to evolutionary-based optimization. Computational introduce some surrogate models as well as some multi-level Optimization and Applications, 56(1):209–229, 2013. approaches, which would allow for additional speed up of [9] G. Papa and B. Koroušić-Seljak. An artificial the evaluation process, since most of the real-world problems intelligence approach to the efficiency improvement of have time-complex evaluations. a universal motor. Engineering Applications of Artificial Intelligence, 18(1):47–55, February 2005. 7. ACKNOWLEDGMENTS [10] G. Papa, B. Koroušić-Seljak, B. Benedičič, and The authors acknowledge the financial support from the T. Kmecl. Universal motor efficiency improvement Slovenian Research Agency (research core funding No. P2- using evolutionary optimization. IEEE Transactions 0098). This work is also part of a project SYNERGY that on Industrial Electronics, 50(3):602–611, June 2003. has received funding from the European Union’s Horizon [11] G. Petelin and G. Papa. Application for multi-criteria 2020 research and innovation programme under grant agree- optimization and visualization of solutions. Technical ment No 692286. Report 12380, Jožef Stefan Institute, 2017. [12] T. Tušar, P. Korošec, G. Papa, B. Filipič, and J. Šilc. 8. REFERENCES A comparative study of stochastic optimization [1] Domel company. http://www.domel.com/. Accessed: methods in electric motor design. Applied Intelligence, 2018-09-17. 27(2):101–111, 2007. 42 Model-Based Multiobjective Optimization of Elevator Group Control Aljoša Vodopija1,2, Jörg Stork3, Thomas Bartz-Beielstein3, Bogdan Filipič1,2 1Department of Intelligent Systems, Jožef Stefan Institute, Jamova cesta 39, SI-1000 Ljubljana, Slovenia 2Jožef Stefan International Postgraduate School, Jamova cesta 39, SI-1000 Ljubljana, Slovenia 3TH Köln, Institute of Data Science, Engineering, and Analytics, Steinmüllerallee 1, D-51643 Gummersbach, Germany {aljosa.vodopija, bogdan.filipic}@ijs.si, {joerg.stork, thomas.bartz-beielstein}@th-koeln.de ABSTRACT Sahin et al. [6], a real-time monitoring system is installed to Finding a suitable control strategy for the elevator group reduce the number of redundant stops, and improve passen- controller (EGC) is a complex optimization problem with ger comfort and energy consumption. In [1], an approxima- several objectives. We utilize the sequential-ring (S-Ring) tion model for EGC systems, the so-called sequential ring model of EGC systems and propose a biobjective formula- (S-Ring) [4], is used to benchmark single-objective heuris- tion of the EGC optimization problem. Unlike the previous tics. Using the S-Ring model, it is possible to retain a high work, we use true multiobjective optimizers in solving this level of complexity and optimize an EGC control strategy problem. Their results on three real-world elevator systems using modern heuristics with a high number of strategy eval- reveal the possible trade-offs between the objectives and of- uations, while keeping a feasible computational load. fer a valuable insight into the problem. In this paper, we utilize the S-Ring model of EGC systems Keywords and propose a biobjective formulation of the EGC optimiza- tion problem. In this formulation, the objectives are normal- elevator group control, S-Ring, perceptron, multiobjective ized to allow for comparison of results for elevator systems of optimization, NSGA-II, DEMO various configurations. As opposed to previous work, we ap- ply true multiobjective optimizers capable of finding approx- 1. INTRODUCTION imations for Pareto-optimal solutions that represent trade- With larger number of people living in urban areas and offs between the objectives. Specifically, we use two multiob- modern barrier-free building design, elevator systems are jective evolutionary algorithms (MOEAs) and demonstrate becoming more and more important. Modern multi-car el- their performance in optimizing EGC for three real-world evator systems are controlled by elevator group controllers elevator systems. (EGC) that assign elevator cars to their destinations based on the customer service calls. The control strategy strongly The paper is further organized as follows. Section 2 intro- affects the desired service quality, customer satisfaction, en- duces the S-Ring model, explains its elements and illustrates ergy consumption, and material attrition. Thus, finding an it with an example. Section 3 provides the optimization adequate control strategy depicts a complex optimization problem formulation. In Section 4, numerical experiments problem with several objectives, which is further dependent on the three test elevator systems and the results are pre- on the building structure and the passenger traffic situa- sented. Section 5 concludes the paper by summarizing the tion. Optimization of EGC imposes challenges, such as be- study and planning future work. ing nonlinear and multimodal, as well as highly dynamic and stochastic due to the stochasticity of customer arrivals. 2. S-RING MODEL OF ELEVATOR GROUP This renders classic gradient-based optimizers as not appli- CONTROL cable to these problems [1]. Moreover, EGC simulators are The S-Ring is a discrete, nontrivial event system to optimize computationally expensive and limit the number of control and benchmark control strategies without the need to use strategy evaluations. expensive EGC simulators [4]. It focuses on modeling the operation of an elevator system, i.e., handling the passenger While EGC optimization problems are widely discussed and traffic and serving passengers in the fastest and most com- known for involving conflicting objectives, they are seldom fortable way. We adapted the S-Ring model to feature two solved with true multiobjective optimization. Hakonen et service quality related objectives as described in Section 3. al. [3] utilize a set of objectives, such as the customer waiting time, the ride time, and the total number of elevator stops, In general, the S-Ring consists of three key elements: but combine them linearly into a single objective. Tyni and Ylinen [7] use a weighted aggregation method to optimize the landing call waiting time and energy consumption with • The deterministic state-space representation of the el- an evolutionary algorithm in a real-time environment. In evator control inputs for the customers ci and elevator 43 cars si, i = {1, ..., Ns}, where Ns = 2n − 2 is the the average number of states with waiting customers, and number of states, while n is the number of floors. Fig- ii) the total number of elevator stops [3, 6, 7]. In contrast to ure 1 shows an example of this state-space representa- previous publications, we do not combine the objectives into tion. The size of the S-Ring depends on the number a single function, but adopt the multiobjective perspective. of floors n, and the number of active elevator states is Moreover, to make it possible to compare the performance equal to the number of elevator cars m. The number of elevator systems of various configurations (determined by of currently active customer states is influenced by the the number of floors n and the number of elevator cars m), probability of a new arriving customer, p. we consider normalized objective function values. • The state transition table, which is explained in detail The first objective (h1) is the proportion of states with wait- by Markon [4], defines fixed and dynamic rules for a ing customers. It is expressed as the average number of transition in the current position of the S-Ring. If no states with waiting customers, Mw, divided by the number fixed rule is triggered, the dynamic rules decide how of all states, Ns: the state transition is performed. They are established M by a control policy. w h1 = . (1) Ns • The control policy π can be realized by a lookup table, The second objective (h2) is the proportion of elevator stops. but as its size grows exponentially with n, it is main- It is equal to the total number of elevator stops, Mt, divided tained by a perceptron with a weight vector of length by the maximum possible number of elevator stops. The |w| = 2Ns. The perceptron represents the most ele- latter can be calculated as the number of elevator cars m mentary implementation of neural network (NN). For multiplied by the number of EGC simulation cycles, which a given setup of n, m and p the objectives are only is in turn equal to the number of state transition steps, Nt, influenced by the weight vector w of the NN controller divided by the number of states, Ns, therefore and the number of state transition steps, Nt. At each M state, it is first checked whether a new customer ar- t h2 = . (2) rived. Next, if the current state is an active elevator mNt/Ns state, the controller determines whether the elevator Intuitively, the customers’ discomfort with long waiting times car stops or continues to the next state. Finally, the and long riding times due to many elevator stops does not indication of the customer active state is updated de- increase linearly with time, but more drastically. To model pending on whether or not the customers were served. this effect, we have additionally modified the original objec- tives as f1 = hα 1 and f2 = hβ, (3) 2 3rd floor where α, β ∈ [1, 2] are the objective function coefficients. The choice of their values is subjective, but the idea is to reflect the elevator system characteristics and the custumer 2nd floor preferences. customer ci The control policy π is represented by a perceptron as π(x) = θ(wTx), where x is a binary input vector, i.e., a concate- elevator si nation of the waiting customer and the elevator car state 1st floor vectors of total length equal to 2(2n − 2) = 4(n − 1), θ is the Heaviside function, and w a vector of perceptron weights from W = [−1, 1]4(n−1). In this framework, the policy π is ground floor defined by the weight vector w only. Therefore, the decision space of the EGC optimization problem as defined here is Figure 1: S-Ring: No waiting customer at the equal to W . ground and floor (“0”), two customers who want to go up on the first and second floor (“1,1”), and no 4. EXPERIMENTS AND RESULTS customers who want to go down on the third, sec- The multiobjective optimization of EGC was experimentally ond and first floor (“0,0,0”). By combining these evaluated on three test problems reflecting the characteris- information we obtain the following state vector for tics of real-world elevator systems operating in various build- waiting customers: (0,1,1,0,0,0). The state vector ings in Ljubljana, Slovenia. They are as follows. for the elevator is obtained in a similar manner. • S1: This system operates in a parking building (“Park- Due to its low computational costs, the S-Ring can quickly ing garage Šentpeter”) situated in the city center. In- evaluate a broad variety of EGC instances as benchmarks tensive passenger traffic can be observed in the build- for the proposed multiobjective optimization approach. ing on workdays. • S2: This is an elevator system installed in a typical 3. OPTIMIZATION PROBLEM FORMULA- residential building in the densely populated neighbor- TION hood (“Nove Fužine”) in the eastern part of Ljubljana. In this work, we deal with two EGC objectives that are often Here the traffic intensity alternates between high (e.g., studied in the literature and both need to be minimized: i) early in the morning) and low (e.g., at midday). 44 • S3: This is the elevator system in the “Crystal Palace”, robust and repeatable algorithm behavior on all three ele- a skyscraper situated in the north-western area of the vator systems. Similarly, small deviations are present for city. With its 89 meters it is currently the tallest build- execution times no matter which MOEA is used to produce ing in Slovenia. As an office building it has low pas- approximations for Pareto fronts. senger traffic. Figures 2 and 3 show Pareto front approximations for the test elevator systems resulting from typical runs of NSGA- The characteristics of these elevator systems are summarized II and DEMO, respectively (there were negligible differences in Table 1. between the results of different runs). As we can see, both MOEAs obtain well-distributed and very similar sets of so- Table 1: Characteristics of the test elevator systems: lutions. The best solutions with respect to both objectives number of floors n, number of elevator cars m, prob- were found for system S3. This was expected since S3 has ability of new arriving customer p, objective func- more elevator cars and a lower probability of new customer tion coefficients α and β, number of states in the arrivals than S1 and S2. S-Ring representation Ns. System n m p α β Ns 0.6 System S1 7 2 0.6 1.0 1.5 12 S1 S2 13 2 0.3 1.4 1.8 24 S2 S3 21 4 0.2 1.5 1.5 40 ator stops S3 v 0.4 Based on the multiobjective formulation of the EGC op- timization problem, the experimental evaluation aimed at tion of ele finding sets of trade-off solutions representing approxima- 0.2 tions for Pareto fronts. For this purpose we used two well- known MOEAs: Nondominated Sorting Genetic Algorithm f2: Propor II (NSGA-II) [2] and Differential Evolution for Multiobjec- tive Optimization (DEMO) [5]. The algorithms were as- 0.0 sessed from the point of view of both effectiveness (quality 0.7 0.8 0.9 1.0 of results) and efficiency (spent computational resources). f1: Proportion of states with waiting customers The experimental setup was defined in the following way. Figure 2: Pareto front approximations for the test Both algorithms were run with populations of 100 solutions elevator systems produced by NSGA-II. for 100 generations. Specifically, in NSGA-II, the crossover probability was set to 0.7 and the mutation probability to 0.2, while DEMO was run using the SPEA selection proce- 0.6 System dure, the crossover probability of 0.3 and the scaling factor S1 of 0.5. On each test problem every MOEA was run 30 times, S2 each time with a new randomly initialized population. ator stops S3 0.4 v Population members were the perceptron weight vectors of length 2Ns = 4(n − 1). Each solution was evaluated through a computer simulation of the perceptron EGC during which tion of ele the values of objectives f 0.2 1 and f2 were calculated. The sim- ulation was performed for a predefined number of simulation cycles which was 100.000 for all test problems. As a conse- f2: Propor quence, the number of state transition steps was equal to 0.0 Nt = 100.000Ns. 0.7 0.8 0.9 1.0 f1: Proportion of states with waiting customers The quality of results of an algorithm run was measured with the hypervolume of the Pareto front approximation found in Figure 3: Pareto front approximations for the test that run. Given f1, f2 ∈ [0, 1], the reference point for hyper- elevator systems produced by DEMO. volume calculations was set to (1.1, 1.1)T. As the compu- tational efficiency measure the execution time of algorithm runs was recorded. The experiments were run on a 3.50 GHz An additional experiment was devoted to the analysis of Intel(R) Xeon(R) E5-2637V4 CPU with 64 GB RAM. hypothetical variants of system S3 with various numbers of elevator cars. While S3 has its fixed configuration, such The hypervolume and execution time results are shown in a study is relevant for designing elevator systems for new Table 2, both averaged over 30 runs of every MOEA on each buildings and assessing potential configurations. test problem. From these results it is evident that regard- less of the elevator system, the hypervolumes obtained with Pareto front approximations obtained with NSGA-II for vari- NSGA-II and DEMO are very similar. Standard deviations ants of S3 with 2, 3, 4, 5 and 6 elevator cars are presented for both optimizers are small (less than 10−3), indicating in Figure 4. The figure clearly shows how the number of 45 Table 2: Average hypervolume and average execution time for both optimizers on the test elevator systems. NSGA-II DEMO Elevator system Hypervolume Time [min] Hypervolume Time [min] S1 0.28066 ± 0.00005 38 ± 1 0.28069 ± 0.00005 28 ± 1 S2 0.32455 ± 0.00016 147 ± 1 0.32450 ± 0.00014 128 ± 1 S3 0.46506 ± 0.00081 398 ± 2 0.46543 ± 0.00037 401 ± 2 In the future we plan to further assess the resulting elevator 0.4 No. of control policies through a comparisson with the results of elevators single-objective optimization and investigate the scalability 2 of the applied optimization methodology. We will also an- 0.3 3 ator stops alyze the produced trade-off solutions in the design space, v 4 and deal with alternative, potentially more transparent pol- 5 icy implementations. 0.2 6 tion of ele 6. ACKNOWLEDGMENTS 0.1 This work is part of a project that has received funding from the European Union’s Horizon 2020 research and innovation f2: Propor programme under grant agreement no. 692286. We acknowl- 0.0 edge financial support from the Slovenian Research Agency 0.5 0.6 0.7 0.8 0.9 1.0 (young researcher programme and research core funding no. f1: Proportion of states with waiting customers P2-0209). Figure 4: Pareto front approximations for variants 7. REFERENCES of S3 with different numbers of elevator cars (2, 3, [1] T. Bartz-Beielstein, M. Preuss, and S. Markon. 4, 5, 6) found by NSGA-II. Validation and optimization of an elevator simulation model with modern search heuristics. In Metaheuristics: Progress as Real Problem Solvers, cars affects the trade-off EGC policies. Higher number of pages 109–128. Springer, 2005. cars implies policies that can reduce the proportion of states [2] K. Deb, A. Pratap, S. Agarwal, and T. Meyarivan. A with waiting customers and the proportion of elevator stops fast and elitist multiobjective genetic algorithm: simultaneously. For example, in the case of only 2 elevator NSGA-II. IEEE Transactions on Evolutionary cars the lowest value of objective f1 is about 0.8, while with Computation, 6(2):182–197, 2002. 6 elevator cars it can be reduced to 0.5. However, one should [3] H. M. Hakonen, A. Rong, and R. Lahdelma. be careful in comparing the results with respect to f2, since Multiobjective optimization in elevator group control. the maximum possible number of elevator stops increases In P. Neittaanm with the number of elevator cars. Nevertheless, these results äki, T. Rossi, S. Korotov, E. O˜ nate, J. Périaux, and D. Kn allow for better problem understanding and are insightful to örzer, editors, European Congress on Computational Methods in Applied Sciences and various stakeholders involved in deciding on elevator system Engineering, ECCOMAS 2004, Jyv configurations, ranging from EGC designers to investors. äskylä, 2004. 17 pages. [4] S. Markon. A solvable simplified model for elevator 5. CONCLUSIONS group control studies. In 2015 IEEE 4th Global We studied the optimization of EGC needed in the design Conference on Consumer Electronics (GCCE), pages and operation of multi-car elevator systems. Utilizing the S- 56–60, 2015. Ring model of EGC systems, we proposed a biobjective for- [5] T. Robič and B. Filipič. DEMO: Differential evolution mulation of the EGC optimization problem that takes into for multiobjective optimization. In C. A. Coello Coello, account the proportion of states with waiting customers and A. Hernández Aguirre, and E. Zitzler, editors, the proportion of elevator stops, both subject to minimiza- Proceedings of the Third International Conference on tion. In this formulation, the objectives are normalized to Evolutionary Multi-Criterion Optimization, EMO 2005, support comparative empirical studies on elevator systems Guanajuato, Mexico, volume 3410 of Lecture Notes in with various numbers of floors and elevator cars. Computer Science, pages 520–533. Springer, 2005. [6] Y. G. Sahin, S. Uzunbayir, B. Akcay, and E. Yildiz. As opposed to previous work, we applied true multiobjec- Real-time monitoring of elevators to reduce redundant tive optimizers capable of finding approximations of Pareto- stops and energy wastage. In Proceedings of the 2013 optimal solutions. The results of two MOEAs for three real- International Conference on Systems, Control and world elevator systems were comparable regarding both the Informatics, pages 264–269, 2013. quality and execution time. They revealed the nondom- [7] T. Tyni and J. Ylinen. Evolutionary bi-objective inated sets of trade-off control policies for the considered optimisation in the elevator car routing problem. elevator systems. Moreover, we demonstrated how the ap- European Journal of Operational Research, proach can be used to support the elevator system configur- 169(3):960–977, 2006. ing at the design stage. 46 From a Production Scheduling Simulation to a Digital Twin Gregor Papa Peter Korošec Jožef Stefan Institute Jožef Stefan Institute Jamova cesta 39 Jamova cesta 39 Ljubljana Slovenia Ljubljana Slovenia gregor.papa@ijs.si peter.korosec@ijs.si ABSTRACT If the approach would be left as is, it would be considered Digital twins are becoming ever more important in smart only as multi-objective approach using a simulation tool to specialisation of factories of the future. Transition from find an approximation set of non-dominated solutions. But using current state in industry to using digital twins is a since it can be introduced into the actual production, mean- big step. We propose an initial step to upgrade simulations ing that the current information of the state of production, to digital twins to enhance the productivity even further. with regard to standing orders and orders which have al- The multi-objective optimisation approach is important in ready been processed so far, we can consider such an en- achieving high efficiency of production scheduling. The goal hanced simulation model to be a digital twin of the produc- of the optimisation is to find a production schedule that tion. With it, we could not only simulate theoretical future satisfies different, contradictory production constraints. We capacities, but also include actual production and its daily take a simulation tool that was used by a memetic version specifics to predict future events with higher accuracy. of the Indicator-Based Evolutionary Algorithm with cus- tomized reproduction operators and local search procedures The rest of the paper is organized as follows: in Section 2, to find a set of feasible, non-dominated solutions and anal- we briefly describe the production scheduling problem; in yse the required steps to achieve a digital twin. We show Section 3, we introduce required changes to create a digi- that with a multi-objective approach that is able to find tal twin; in Section 4, we present the main idea of Memetic high-quality solutions and flexibility of many “equal” solu- IBEA; in Section 5, we present the experimental environ- tions, the digital twin becomes a powerful tool for a decision ment and the results of the evaluation with the real-world maker. data; in Section 6, we present the usability study; and in Section 7, we draw conclusions and propose possible future Keywords work. multi objective, scheduling, optimisation, real world, digital twin 2. PRODUCTION SCHEDULING 1. INTRODUCTION PROBLEM Since production scheduling is important for smart special- The scheduling problem was introduced for a company that isation goals in factories of the future, we decided to take produces components for domestic appliances, including hot relevant results from [4], and apply them to see the impact plates, thermostats and heating elements. The fabrication of digital twins. A digital twin is a digital copy of physical process for components used in different types of plates is world (physical twin) in form of processes and systems. It similar, but due to clients’ demands the models differ in size provides both, the elements and dynamics of the real-word, (height, diameter), connector type, and power characteris- so one can simulate and predict the future events with an tics (wattage). For logistic reasons, the clients group dif- up-to-date model, which is relevant for a decision maker. ferent models of plates within the same order, implying the same due-dates for different products. As a consequence, In [4] we applied the multi-objective approach that uses spe- their production must be scheduled very carefully to fulfil cific local search procedures to the problem of production all the demands (quantities and due-dates), to maintain the scheduling. As the basic algorithm we used the Indicator- specified amounts of different models in stock, to optimally Based Evolutionary Algorithm (IBEA) [8]. We decided to occupy their workers, and to make efficient use of all the use the IBEA because it was shown that it can substantially production lines. The assignment of due-dates is usually outperform results generated by other multi-objective algo- performed separately and before the production scheduling, rithms, such as the improved Strength Pareto Evolutionary but since there are strong interactions between the two tasks, Algorithm [9] and NSGA-II [2], in terms of different per- using the proposed digital twin can allow for more accurate formance measures [8]. Due to the addition of local search arrangement of due-dates. For each order, the completion procedures, we called our approach the Memetic Indicator- time should be as close as possible to the due-date in order Based Evolutionary Algorithm (M-IBEA). As such it repre- to reduce the waiting time and costs [7]. Furthermore, not sents a synergy of the multi-objective evolutionary approach all the production lines are equal, since each of them can with separate, individual, learning or local improvement pro- produce only a few different models. A detailed formulation cedures (local searches). of the production scheduling problem is presented in [5]. 47 The required inputs to such a problem are: An indicator function assigns each Pareto-set approxima- tion a real value that reflects its quality. The optimisation goal becomes the identification of a Pareto-set approxima- • Production norms that specify which products are be- tion that minimizes an indicator function. The main advan- ing produced on each line and what is the changeover tage of the indicator concept is that no additional diversity- time from one product to another for each specific line. preservation mechanisms are required [1]. • Amount of stock for each product. The detailed description of the memetic IBEA can be found • Orders that need to be processed and their deadlines. in [4], but the main idea is presented as a pseudo code in Algorithm 1. In our implementation of the basic version, the • Number of planned shifts. IBEA is used to guide the local search procedures. Since we • Number of lines. are dealing with a combinatorial problem, we implemented problem-specific versions of the crossover and mutation op- erators. Additionally, we added different local search proce- Looking from the perspective of a simulation tool that is dures to enhance the efficiency of the algorithm. able to take into account all this inputs and evaluate the expected time of production for every order, it is a simple Algorithm 1 Memetic IBEA simulation tool. But such a tool alone lacks the dynamics of the real world, so it is not able to react “instantly” to the 1: SetInitialPopulation(P ) 2: Evaluate(P ) changes in the production environment. 3: while not EndingCondition() do 4: P 0 = MatingSelection(P ) 3. DIGITAL TWIN 5: Crossover(P 0, pc) For a simulation tool to become a digital twin, some capa- 6: Mutation(P 0, pm) bilities need to be added. Mainly, the interaction between 7: Evaluate(P 0) what is happening in the real world and the description of 8: LocalSearch(P 0) the problem instance. First of all, the relevant information, 9: P = CalculateFitness(P ∪ P 0) which defines the problem instance, can be gathered from 10: P = RemoveWorse(P ) the company’s information system. This allows receiving 11: end while up-to-date information about new orders, the current stock, and amount of products that were produced so far in the Compared to the basic version of the algorithm, the main day. With the way production companies are working, usu- difference is in the procedure LocalSearch(P 0). Here, not ally this needs to be done only once a day, since production only one but many problem specific local search procedures plans do not change for the current day (actually they are are applied [4]. fixed for up to several days in advance), due to the require- ments of having the required materials for producing orders Such a version of the algorithm is suitable for running a sim- at hand. The main reason for this is that an additional re- ulation based approach, but it lacks the required dynamicity quirement is also to have the stock of materials at the factory to actively adapt to changes in the production environment. as small as possible. We must be aware that any unneces- Two things need to change, first, the changes in the pro- sary stock is actually an expense that every company would duction environment should be transferred to the algorithm like to reduce or even remove. solution space, and second, the algorithm should be able to detect and adapt to such changes. Since the production is The simulation tool only takes into account the technical not a living system that changes every second and requires data provided by the company with regard to the above men- immediate changes (as mentioned above, the production is tioned required inputs. Though any changes in production fixed for several days in advance) this is not a crucial aspect, can be “detected” by the simulator through changes in in- since this changes could be applied to the algorithm on a puts (e.g., how many products were actually produced), this daily basis. But from the point of view of acquiring new or- does not provide a good baseline for predicting future pro- ders and providing potential deadlines to the customers, this duction with inclusion of predicting maintenance. For pre- is another matter. By providing a more dynamic system, a diction maintenance to be included in the digital twin a ma- product sales person could easily insert a new potential or- chine learning techniques should be used to estimate/model der and determine what would be the most efficient and safe any informalities that happen, but are not included in pro- deadline to be offered to the customer. And if a customer re- duction norms (e.g., failures on lines). All this is based on quires an earlier deadline, which could force other orders to previous experiences and requires to gather lots of data, so be put in jeopardy of missing the deadline, it allows a prod- the machine learning algorithm is able to be trained to de- uct sales person to better estimate the required higher price tect abnormal, correlated patterns in production, which will for covering the costs ocured from delays of other orders. lead to better predicting future production and provide in- The use of machine learning would also cover the irregular- sight into preventing maintenance, which will lead to further ities that happen in production. reducing of delays on production lines due to failures by ap- plying maintenance before a defect happens. 5. CASE STUDY 4. MEMETIC IBEA 5.1 Test cases The IBEA is a multi-objective version of a genetic algorithm, The algorithm was tested on two real order lists from the where the selection process is based on quality indicators. production company. Task 1 consisted of n = 470 orders 48 Table 1: Comparison of the BF (12 threads) and Table 2: Results of optimisation for Task 1. M-IBEA approach (1 thread). Statistics norders nworkers tchange ndays Evaluations Time Pareto Pareto min 18 631 353 127 n BF M-IBEA BF M-IBEA matching Pareto max 88 823 867 681 7 3.94 · 108 3.5 · 104 22 s 17 s 4/4 Single-objective 18 767 714 156 8 1.58 · 1010 5 · 105 15 min 33 s 5/5 9 7.09 · 1011 5 · 106 11 h 5 min 15/15 Table 3: Results of optimisation for Task 2. Statistics norders nworkers tchange ndays for 189 different products and Task 2 consisted of n = 393 Pareto min 16 538 355 59 orders for 175 different products. The number of orders n Pareto max 50 778 433 330 represents the problem dimension, with m = 5 representing Single-objective 15 702 443 155 the number of available production lines. To mimic the digital twin which is being updated with in- single-objective solution. Though we used the same num- formation once a day (after the end of the daily production) ber of evaluations, this single-objective solution does not we ran a task overnight and looked at the results. In this stand out with respect to any objective – quite the oppo- time, the algorithm made about 300 million evaluations, so site is the case. This can also be observed from Table 2, this was set as our stopping criterion for future tests. A where the single-objective solution returns an average qual- lexicographic evaluation [6] was used for presenting multi- ity solution on all the objectives except n objective solutions. In the simulation evaluation, the num- orders. The results are summarised in Tables 2 and 3, where the width of the ber of delayed orders (norders) was set as the most impor- Pareto approximation front is denoted with “Pareto min” tant objective, followed by the required number of workers and “Pareto max”. (nworkers), the sum of delayed days for all the delayed or- ders (ndays), and the sum of the change-over downtime in From the results we can conclude that using the Pareto-front minutes (tchange). This order was set according to the most approach gives us an expected greater versatility in choosing common objective hierarchy. a good solution, while at the same time we are not sacrificing one, likely the most important, objective. The only impor- 5.2 Evaluating the approach tant drawback is that multi-objective approaches need many To make sure that our proposed M-IBEA was working well, more evaluations than single-objective approaches. So, if we we ran a brute-force (BF) approach where all the possible do not have time to carry out enough evaluations, then the solutions were evaluated for n < 10 orders and the optimal single-objective approach is the only way. Pareto front was constructed for each of them. Table 1 shows a comparison of the number of problem evaluations, the ex- ecution time, and the matching of the Pareto front obtained 6. USABILITY OF MULTI-OBJECTIVE for n = 7, 8, 9. We did not include smaller n values, since SOLUTIONS in all cases a sub-one-second time was needed with perfect The multi-objective approach provides a set of feasible so- Pareto matching. From the obtained results it is clear that lutions, offering the possibility to select the final schedule with more than nine orders, the complexity increases well based on the specific decision maker needs. Since none of beyond an acceptable time (approximately two months) to the given solutions dominates the other solutions, all of them calculate all the solutions. Also, in all cases we were able are acceptable. Based on the current conditions, and accord- to acquire the same Pareto front using the BF and M-IBEA ing to the proposed set of solutions, a decision maker can approaches. When considering times, one must take into give more weight to some of the decision criteria. For this consideration that the BF was ran multithreaded with 12 an intuitive representation of the resulting solutions inside threads fully utilized, while the M-IBEA approach was sin- the GUI of the Planer application was provided, which is gle threaded. The perfect Pareto-front matching is unsur- presented in Figure 1. prising, since the IBEA already proved to be one of the best algorithms for solving multi-objective problems with more After the M-IBEA algorithm found the set of non-dominated than three objectives [3], which was also the main reason solutions, they are presented in the Planer application. In that we selected the IBEA as our base algorithm. the upper-right section there is a list of all the non-dominated solutions. In general, there might be up to several hundred 5.3 Results possible solutions. In [5], we optimized only according to the number of orders. To show that the multi-objective approach presented in [4] However, some of the criteria can be set tighter according is a better alternative, we compared the results with regard to the resulting range of each criterion, and according to the to the best result from the single-objective to the multi- current business conditions. In the specific example shown objective approach. The results showed that the single- in Figure 1, the initial set consisted of 518 solutions. The de- objective solution primarily concentrated on the number cision maker put the first objective into the range from 16 to of orders, while it neglected other objectives. But this is 17 out of 50, which in consequence moved the sliders of the not a surprise, as multi-objective solutions were able to find second objective from 697 to 738, the third objective from equally good solutions with regard to the number of orders 405 to 415, and the fourth objective from 60 to 111. So ir- and significantly better for other objectives, compared to regardless of which slider is moved, the ranges move accord- 49 ten not realistic. 8. ACKNOWLEDGMENTS The authors acknowledge the financial support from the Slovenian Research Agency (research core funding No. P2- 0098). This work has also received funding from the Eu- ropean Union’s Horizon 2020 research and innovation pro- gramme under grant agreements No 692286 (project SYN- ERGY) and No 722734 (under the Marie Sk lodowska-Curie action project UTOPIAE). 9. REFERENCES [1] M. Basseur and E. Burke. Indicator-based multi-objective local search. In Evolutionary Computation, 2007. CEC 2007. IEEE Congress on, pages 3100 –3107, sept. 2007. [2] K. Deb, S. Agrawal, A. Pratap, and T. Meyarivan. A fast elitist non-dominated sorting genetic algorithm for multi-objective optimization: Nsga-ii. In Figure 1: The graphical user interface. M. Schoenauer, K. Deb, G. Rudolph, X. Yao, E. Lutton, J. Merelo, and H.-P. Schwefel, editors, Parallel Problem Solving from Nature PPSN VI, ingly to the possible solutions of other objectives. Simulta- volume 1917 of Lecture Notes in Computer Science, neously, the list of possible solutions is updated to reflect the pages 849–858. Springer Berlin / Heidelberg, 2000. current setting of the objectives’ ranges. In the above ex- [3] H. Ishibuchi, N. Tsukamoto, and Y. Nojima. ample, the list narrowed down to 14 solutions. From them, Evolutionary many-objective optimization: A short the decision maker can select one solution which best fulfils review. In Evolutionary Computation, 2008. CEC 2008. the current demands. The visual representation consists of (IEEE World Congress on Computational Intelligence). all the relevant data, i.e., the production lines’ load, the or- IEEE Congress on, pages 2424–2431, june 2008. der types’ distribution, and change-over downtime lengths, [4] P. Korošec, U. Bole, and G. Papa. A multi-objective which are necessary to make the final decision. If the visual approach to the application of real-world production representation of the solution is accepted, it becomes the scheduling. Expert Systems with Applications, production schedule. By determining (using sliders), which 40(15):5839–5853, 2013. objective is the most important in the current situation and [5] G. Papa, V. Vukašinović, and P. Korošec. Guided to what extent, we automatically determine which part of restarting local search for production planning. Pareto front is important and at the same time disregard all Engineering Applications of Artificial Intelligence, the solutions from the Pareto front, which do not fulfil the 25(2):242–253, 2012. selected conditions. This way we are able to freely move the [6] M. Rentmeesters, W. Tsai, and K.-J. Lin. A theory of useful part of the Pareto front by moving sliders. lexicographic multi-criteria optimization. In Engineering of Complex Computer Systems, 1996. 7. CONCLUSION AND FUTURE WORK Proceedings., Second IEEE International Conference on, pages 76 –79, oct 1996. We presented what steps would be needed to make a memetic, multi-objective approach that used a simulation tool to asses [7] R. Zhang and C. Wu. A hybrid local search algorithm some real-world test cases of a production scheduling prob- for scheduling real-world job shops with batch-wise lem a more dynamic system by upgrading a simulation tool pending due dates. Engineering Applications of to a digital twin. From the perspective of the algorithm not Artificial Intelligence, 25(2):209 – 221, 2012. Special many changes would be required, since with a restart proce- Section: Local Search Algorithms for Real-World dure being already implemented any changes in the problem Scheduling and Planning. description could be “inserted” into the problem solving part. [8] E. Zitzler and S. Künzli. Indicator-based selection in multiobjective search. In in Proc. 8th International On the other hand, more substantial changes are required Conference on Parallel Problem Solving from Nature within the simulation tool. Primarily, how required inputs (PPSN VIII), pages 832–842. Springer, 2004. are being automatised (gathering data directly from the [9] E. Zitzler, M. Laumanns, and L. Thiele. SPEA2: company’s system). Additionally, an inclusion of some ma- Improving the strength pareto evolutionary algorithm chine learning algorithm, that would be able to detect and for multiobjective optimization. In K. C. Giannakoglou, predict failures on production lines, is foreseen for better D. T. Tsahalis, J. Périaux, K. D. Papailiou, and longterm estimation of production. T. Fogarty, editors, Evolutionary Methods for Design Optimization and Control with Applications to For future work, we are planning to implement the pro- Industrial Problems, pages 95–100, Athens, Greece, posed changes, which will enable for more real-life scenarios 2001. International Center for Numerical Methods in (including uncertainties-based worst-case scenarios), while Engineering. currently only “ideal” solutions are provided, which are of- 50 Indeks avtorjev / Author index Bartz-Beielstein Thomas .................................................................................................................................................. 11, 15, 43 Breiderhoff Beate ......................................................................................................................................................................... 15 Dovgan Erik ................................................................................................................................................................................. 27 Ellaia Rachid ................................................................................................................................................................................ 23 Filipič Bogdan .................................................................................................................................................................. 15, 27, 43 Fukumoto Hiroaki ........................................................................................................................................................................ 35 Korošec Peter ......................................................................................................................................................................... 39, 47 Monteiro Wellington Rodrigo ...................................................................................................................................................... 19 Naujoks Boris ............................................................................................................................................................................... 15 Oyama Akira ................................................................................................................................................................................ 35 Papa Gregor ............................................................................................................................................................................ 39, 47 Petelin Gašper .............................................................................................................................................................................. 39 Rahimi Mohammad ...................................................................................................................................................................... 31 Rehbach Frederik ......................................................................................................................................................................... 11 Reynoso-Meza Gilberto ............................................................................................................................................................... 19 Serrar Jihane ................................................................................................................................................................................. 23 Sodnik Jaka .................................................................................................................................................................................. 27 Stork Jörg ............................................................................................................................................................................... 11, 43 Talbi El-Ghazali ..................................................................................................................................................................... 23, 31 Tušar Tea ........................................................................................................................................................................................ 7 Vodopija Aljoša ........................................................................................................................................................................... 43 51 52 Konferenca / Conference Uredila / Edited by Mednarodna konferenca o visokozmogljivi optimizaciji v industriji, HPOI 2018 / International Conference on High-Performance Optimization in Industry, HPOI 2018 Bogdan Filipič, Thomas Bartz-Beielstein Document Outline 02 - Naslovnica - notranja - D 03 - Kolofon - D 04 - IS2018 - Predgovor 05 - IS2018 - Konferencni odbori 07 - Kazalo - D 08 - Naslovnica podkonference - D 09 - Predgovor podkonference - D 10 - Programski odbor podkonference - D 11 - Clanki - D 01_HPOI_2018_paper_8 02_HPOI_2018_paper_3 03_HPOI_2018_paper_7 04_HPOI_2018_paper_4 05_HPOI_2018_paper_11 06_HPOI_2018_paper_2 07_HPOI_2018_paper_1 08_HPOI_2018_paper_5 09_HPOI_2018_paper_6 10_HPOI_2018_paper_10 11_HPOI_2018_paper_9 12 - Index - D Blank Page Blank Page Blank Page Blank Page Blank Page Blank Page Blank Page Blank Page 04 - 05 - IS2018 - Predgovor in odbori.pdf 04 - IS2018 - Predgovor 05 - IS2018 - Konferencni odbori