CZECH INSTITUTE OF INFORMATICS
ROBOTICS AND CYBERNETICS
INDUSTRIAL INFORMATICS DEPARTMENT
Energy Optimization of Robotic Cells
Libor Bukata, Pˇremysl ˇS˚ucha, Zdenˇek Hanz´alek, and Pavel Burget
DOI: https://doi.org/10.1109/TII.2016.2626472
Cite as: L. Bukata, P. ˇS˚ucha, Z. Hanz´alek, and P. Burget. Energy optimization
of robotic cells. IEEE Transactions on Industrial Informatics, 13(1):92–102,
Feb 2017
Source code: https://github.com/CTU-IIG/EnergyOptimizatorOfRoboticCells
c⃝2017 IEEE. Personal use of this material is permitted. Permission from IEEE must be
obtained for all other uses, in any current or future media, including reprinting/republishing
this material for advertising or promotional purposes, creating new collective works, for resale
or redistribution to servers or lists, or reuse of any copyrighted component of this work in other
works.
arXiv:1802.05925v1  [cs.RO]  16 Feb 2018
IEEE TRANSACTIONS ON INDUSTRIAL INFORMATICS
2
Energy Optimization of Robotic Cells
Libor Bukata, Pˇremysl ˇS˚ucha, Zdenˇek Hanz´alek, and Pavel Burget
Abstract—This study focuses on the energy optimization of in-
dustrial robotic cells, which is essential for sustainable production
in the long term. A holistic approach that considers a robotic cell
as a whole toward minimizing energy consumption is proposed.
The mathematical model, which takes into account various robot
speeds, positions, power-saving modes, and alternative orders
of operations, can be transformed into a mixed-integer linear
programming formulation that is, however, suitable only for
small instances. To optimize complex robotic cells, a hybrid
heuristic accelerated by using multi-core processors and the
Gurobi simplex method for piecewise linear convex functions is
implemented. The experimental results showed that the heuristic
solved 93 % of instances with a solution quality close to a proven
lower bound. Moreover, compared with the existing works, which
typically address problems with 3 to 4 robots, this study solved
real-size problem instances with up to 12 robots and considered
more optimization aspects. The proposed algorithms were also
applied on an existing robotic cell in ˇSkoda Auto. The outcomes,
based on simulations and measurements, indicate that, compared
with the previous state (at maximal robot speeds and without
deeper power-saving modes), the energy consumption can be
reduced by about 20 % merely by optimizing the robot speeds
and applying power-saving modes. All the software and generated
data sets used in this research are publicly available.
Index Terms—robotic cells, energy optimization, holistic ap-
proach, Digital Factory, Industry 4.0, parallel hybrid heuristic,
mixed-integer linear programming, simplex method for piecewise
linear convex functions
I. INTRODUCTION
R
OBOTIC CELLS, which are broadly used in heavy-duty
industries, are highly complex systems (see Figure 1)
that consist mainly of industrial robots, conveyors, (stationary)
welding/gluing guns, PLC controllers, and physical barriers
(e.g., fencing). Because the primary objective is to guarantee
the repeatability of production, the robots carry out the same
operations periodically; i.e., they follow fixed cyclic event-
based schedules, the execution of which can be conditioned
by some signals, such as an acknowledgment of a free work
space. Signaling is also used for inter-robot synchronization
to avoid collisions. Thus far, robotic cells are often optimized
with respect to their production rate (see, e.g., Dawande et
al. [1] for a survey); however, the pressure to reduce energy
Manuscript received October 08, 2015; revised March 07, 2016 and August
05, 2016; and accepted October 22, 2016
Libor Bukata and Zdenˇek Hanz´alek are with the Czech Institute of Informat-
ics, Robotics, and Cybernetics, Czech Technical University in Prague, Zikova
1903/4, 166 36, Prague 6, Czech Republic. (e-mail: bukatlib@fel.cvut.cz,
hanzalek@fel.cvut.cz)
Pˇremysl ˇS˚ucha and Pavel Burget are with the Faculty of Electrical En-
gineering, Department of Control Engineering, Czech Technical University
in Prague, Karlovo n´amˇest´ı 13, 121 35, Prague 2, Czech Republic. (e-mail:
suchap@fel.cvut.cz, burgetpa@fel.cvut.cz)
Copyright (c) 2016 IEEE. Personal use of this material is permitted.
However, permission to use this material for any other purposes must be
obtained from the IEEE by sending a request to pubs-permissions@ieee.org.
consumption in the industry is growing due to increasing
costs of electricity, the energy savings goal of the European
Union [2], and other initiatives and requirements of ecological
groups. As a result, there is a high demand for energy-efficient
solutions that enable sustainable development in the long term.
One of the areas in which a significant amount of energy
can be saved is the manufacturing process. For example,
industrial robots consume approximately 8 % of the energy
needed during the production of cars (Meike and Ribickis [3]).
Fig. 1. A typical robotic cell in the industry.
The present study addresses the energy optimization of
robotic cells not only in the automotive industry. A mathemat-
ical model of the robotic cell is presented, and optimization
algorithms are proposed. The primary goal of the optimization
is to minimize the energy consumption of a robotic cell
for a given production rate by changing the robot speeds,
positions (including the rotation of a gun), and order of
operations and by applying the robot power-saving modes.
These optimization aspects have been identified as important
for energy consumption based on measurements made on a
small industrial KUKA robot. Such optimization has signif-
icant potential for energy savings because robotic cells are
usually designed under time pressure with the main objective
of meeting a desired production rate. As a result, the speed of
robot movements is typically maximal, although these mostly
result in long idle times. Such performance leads to a power
profile that sharply fluctuates during production; for instance,
a regular 6-axis industrial robot with a 200 kg payload requires
between 0.5 and 20 kW of power (Meike and Ribickis [3]).
For stationary robots, the robot brakes can slightly alleviate
the energy burden if they are released early because the motors
do not consume energy in holding the robot in a certain
IEEE TRANSACTIONS ON INDUSTRIAL INFORMATICS
3
position, and the power needed to keep the brakes open does
not apply to the released brakes. On one hand, the frequent
use of brakes could result in significant energy savings; on the
other hand, caution must be exercised because the number of
brake cycles is limited to 5 million, and the expected lifetime
of robots is about 15 years (Meike et al. [4]). Some robots,
such as those from KUKA, also support the PROFIenergy
profile [5], which allows controlling the power consumption
by sending commands through a PROFINET network. This
enables deeper energy-saving modes for robots, such as bus
power-off or hibernate. Several experimental facilities have
been using PROFIenergy, such as the Mercedes-Benz plant
in Sindelfingen, in which various experiments have been
carried out [6]. Meike and Ribickis [3] provide more detailed
information on energy-saving techniques.
The current research on the energy optimization of robotic
cells often concentrates on the optimization of individual robot
trajectories (local optimization) with respect to the physical
limitations of robots and the obstacles to be avoided (e.g.,
[7], [8], and [9]). However, to reach the full energy-saving
potential, it is necessary to consider the robotic cell as a whole
(global optimization). Such energy optimization was proposed
in a pilot work of Wigstr¨om and Lennartson [10], in which a
nonlinear mathematical model was formulated and solved by
using various algorithms. Nevertheless, the authors neglected
the energy-saving modes and alternative positions of robots,
and their algorithms were benchmarked on an instance of job
shop scheduling problem.
In the work of Mashaei and Lennartson [11], an energy
model of a pallet-constrained flow shop problem was formu-
lated to find an optimal switching control strategy for achiev-
ing the desired throughput and minimal energy consumption.
Idle states of machines were also taken into account to reduce
the energy consumption when the machine was not working.
However, the model required a line with a particular structure,
i.e., a closed-loop pallet system; therefore, it was not generally
applicable to robotic cells.
A special structure, i.e., a rotationally arranged robotic cell,
was also considered by Foumani et al. [12], who proposed
a methodology for maximizing the production rate subject
to timing constraints, in which multifunction robots that per-
formed an extra operation during the movement were taken
into account to reduce the cycle time. For example, the
robot could measure the dimensions of a workpiece during
its transportation by a Grip-Gage-Go gripper. Sindiˇci´c et al.
[13] proposed a method that determined whether a flexible
manufacturing system was stable for given repeatable se-
quences. The proposed method used a machine-job incidence
matrix that was typically smaller than matrices describing the
structure of Petri nets, thus reducing the complexity of the
analysis. The authors applied their approach on a system with
three robots and four machines.
On the boundary of local and global optimization is the
work of Meike et al. [4], [14], in which the last movement
to the robot home position and the subsequent waiting period
were optimized to save energy. Compared with the work of
Wigstr¨om and Lennartson [10], in this study the robot brakes
were considered at the robot home position. Although only
a relatively small part of the robotic cell was considered, the
authors estimated that the energy savings for the robotic cell
reached 7.3 %.
To push the aforementioned research further, the present
work proposes a holistic mathematical model for the energy
optimization of robotic cells. In contrast to the existing works
[10], [15], [16], this model supports the energy-saving modes
and alternative positions of robots. Because the continuous
manufacturing process is considered and not the production-
free time, only robot brakes and bus power-off energy-saving
modes are applicable in practice due to the long transition time
of deeper energy-saving modes. In the proposed model, the
robot cycle time, which is known as the period from cyclic
scheduling, corresponds to the time interval (~1/throughput)
between two outgoing workpieces after the start-up phase.
The robot cycle time is usually determined according to the
desired production rate of a product. The duration of the start-
up phase, which is called the lead time in cyclic scheduling,
is the total time required for one (first) workpiece to be
processed by the robotic cell; or in other words, it is the
time difference between the entering and the leaving time of
a given workpiece. The cycle time is typically shorter than
the duration of the start-up phase; thus, there are usually
more unfinished workpieces in the robotic cell at a given
time. Such parallelism is analogous with the cyclic scheduling
of the processor pipeline (see, e.g., [17] and [18]). How-
ever, Wigstr¨om and Lennartson [10] presumed that the whole
robotic cell is dedicated to one workpiece during the entire
processing time, which the authors called the work cycle time;
hence, the overall throughput may be limited because there is
no parallelism. This different view of time distinguishes the
problem in the present work. Furthermore, the proposed model
is suitable for both existing and new robotic cells because the
number of optimized aspects, such as the selection of robot
positions, speeds, and energy-saving modes, can be adjusted.
Based on the model, the mixed-integer linear programming
(MILP) problem is formulated and subsequently solved by
applying the Gurobi solver. The experimental results, however,
show that only small problems are solvable close to optimality
and that bigger problems are intractable for the state-of-the-
art MILP solvers. Therefore, this research proposes a parallel
hybrid heuristic that solves some bigger instances in a short
time. The proposed heuristic is very efficient and fast, thanks
to the application of a special Gurobi simplex algorithm for
piecewise linear convex functions and the parallelization that
enables the algorithm to exploit the power of all processor
cores. The proposed algorithms are tested on a real robotic
cell from ˇSkoda Auto (see Section VI) and generated data
sets. The estimation, based on simulations and measurements
of power profiles in ˇSkoda Auto, indicate that up to 17.6 % of
energy can be saved merely by changing the speed of robot
movements and that even more savings can be achieved if
energy-saving modes are applied. These results are partially
supported by Paryanto et al. [19], who identified the speed of
robot movements as crucial for energy efficiency. In addition,
a generator of problem instances is implemented to enable
other researchers to evaluate their algorithms given the lack
of instances that are publicly available. Both the algorithms
IEEE TRANSACTIONS ON INDUSTRIAL INFORMATICS
4
and the generator are published as open-source software and
can be downloaded from https://github.com/CTU-IIG.
The main contributions of the present work are the extended
mathematical model, which considers the power-saving modes
and alternative positions of robots, and the design of the
parallel hybrid heuristic, which has been shown to be effective
in the experiments. The results of the case study from ˇSkoda
Auto indicate that the proposed approach may considerably
reduce the energy consumption of robotic cells. The rest of
the report is structured as follows. Section II sets out the
formal specification of the problem, Section III presents the
MILP model, and Section IV introduces the parallel heuristic
algorithm for solving the problem. The performance of the
proposed algorithms is tested on benchmark instances in
Section V, and the case study from ˇSkoda Auto is discussed
in Section VI. Finally, the conclusions are presented in Sec-
tion VII.
II. PROBLEM STATEMENT
The energy optimization problem of the robotic cell can be
defined as follows. There is a set of robots R = {r1, . . . , r|R|}
and their graphs denoted as Gr = (Vr, Er) where nodes Vr
are static activities (e.g., waiting or welding), and edges Er
are dynamic activities, which define the possible robot moves
between static activities. Let V = S
∀r∈R Vr, E = S
∀r∈R Er,
and Ar = Vr ∪Er.
Each static activity v ∈V has the assigned set Lv of pos-
sible robot positions, i.e., locations, in which the motionless
robot either works or waits for a signal. During this stationary
phase lasting from dv to dv, one of the robot power saving
modes m ∈Mr, including a dummy power-saving mode for
the robot held by motors, can be applied if the activity duration
dv ≥dm, where dm is the minimal time to switch the mode
on. The energy consumption of activity v is then Wv = pm
v,ldv,
where pm
v,l is the robot input power at location l ∈Lv for mode
m ∈Mr.
The dynamic activity, i.e., edge e ∈Er, consists of the set of
trajectories Te between v1 and v2 ∈Vr, where each trajectory
t ∈Te interconnects two locations, and its duration is from
dt
e to d
t
e. The energy consumption of activity e ∈Er is the
function of the selected trajectory t ∈Te, the duration of the
movement de, and a convex function f t
e mapping de to energy
consumption We. The form of the function (see Vergnano et
al. [20], p. 389, Eq. 12) is f t
e(de) = P5
i=1 Ct
e,id2−i
e
, where
Ct
e,i are constants.
Let A = V ∪E be the set of all activities; then, Sa = suc(a)
and Pa = pred(a) are the set of successors and predecessors
of activity a ∈A, respectively. If |Sv| > 1 for v ∈V , then
the dynamic activities e ∈Sv are optional because only one
trajectory t ∈Te will be selected as the leaving one. The
optional activities, denoted as EO ⊆E, can model alternative
orders of operations. Note that |Se| = 1 for ∀e ∈E.
The closed path of robot r ∈R, including the order
of operations, can be represented as a directed Hamiltonian
circuit, HCloc
r
= (l1, . . . , l|Vr|), through the activity locations,
where each v ∈Vr is visited only once, with the exception of a
closing activity. The order of operations (activities), regardless
of the selected locations, can be expressed in the form of a
directed Hamiltonian circuit, HCact
r
= (v1, . . . , v|Vr|), where
the last activity v|Vr|, denoted as vh
r , closes the cycle of robot
r. The last activity vh
r , regardless of the work of robot r,
is predetermined as home positions of robot r, where the
robot typically waits for its next cycle. Note that HCact
r
can be derived from HCloc
r ; however, the converse is not
true. The robot cycle time CT, i.e., the production cycle
time of the robotic cell, has to be fulfilled; in other words,
P
∀a∈A∗r da = CT for ∀r ∈R, where A∗
r ⊆Ar refers to all
activities on the related Hamiltonian circuit HCact
r .
To guarantee that the robots cooperate with each other at the
right time and location, there are global constraints for time
synchronization and spatial compatibility. Correct timing is
ensured by inter-robot time lags ETL, where each time lag
(a1, a2) ∈ETL has the fixed length la1,a2 ∈R and the
height ha1,a2 ∈Z. The length la1,a2 denotes a time shift
between a1 ∈Ari and a2 ∈Arj, and the height ha1,a2 enables
addressing the current and previous robot cycles. Both terms,
which are well-known in cyclic scheduling, define the time
relation as follows: sa2 ≥sa1 + la1,a2 −CTha1,a2, where sa1
and sa2 are the start times of activities a1 and a2, respectively.
To ensure that the workpiece is passed correctly from one
robot to another, i.e., that the workpiece is picked up from the
same place where it was put, it is necessary to define the spatial
compatibility as follows. Let v1 and v2 ∈V be two static activ-
ities in which an inter-robot handover is carried out; then, the
set of compatible location pairs is Qv1×v2 ⊆Lv1×Lv2, and set
Qli = {∀lj : (li, lj) ∈Qv1×v2 or (lj, li) ∈Qv1×v2} contains
all the locations (robot positions) compatible with location li.
Finally, because the concurrent work of robots may result in
collisions, time disjunctive combinations of trajectories and
locations are specified as quadruplets (a1, g1, a2, g2) ∈C,
where gi ∈Lai if ai ∈V ; otherwise, gi ∈Tai.
The goal of the optimization is to find the closed paths of
robots (i.e., HCloc
r ), determine the activity timing (i.e., the
assignment of sa and da for ∀a ∈A), and decide which
power-saving mode of the robot to apply for each static activity
v ∈V such that global constraints are satisfied, collisions are
prevented, and the energy consumption (i.e., P
∀a∈A Wa) is
minimized.
A. Example
Figure 2 shows a robotic cell consisting of two robots:
r1 and r2
∈
R. Robot r1 takes the weldment, carries
out two spot-welding operations, and subsequently puts the
weldment on a bench. Robot r2 takes the weldment from
the bench, carries out an additional spot-welding opera-
tion, and puts the weldment on a conveyor. The process
flow of each robot can be expressed as a graph Gr; such
as Gr1
= (Vr1, Er1), where Vr1
= {v1, v2, v3, v4} and
Er1 = {e1,2, e1,3, e2,3, e2,4, e3,2, e3,4, e4,1}. Nodes v2, v3, v6
correspond to spot-welding operations, whereas v1, v5 and
v4, v7 are take and put operations, respectively. The edges
define how to move from one operation to another, and the
edge style determines whether the move is mandatory (solid
line), or optional (dashed line). The dashed edges, i.e., the
IEEE TRANSACTIONS ON INDUSTRIAL INFORMATICS
5
Fig. 2. Robotic cell with two robots that perform welding operations.
optional dynamic activities EO, enable modeling alternative
orders of operations; in this case, there are two possible
orders: (v1, v2, v3, v4) and (v1, v3, v2, v4). The graphs shown
in Figure 2 have activity-level granularity because each edge is
a dynamic activity and each node is a static activity. However,
a lower-level of granularity is achievable because each node v
has locations Lv, and each edge e contains trajectories Te.
Because the weldment has to be passed at the right place
at the right time, spatial and time synchronization between
robots is necessary. For example, if Lv4 = {l4,1, l4,2} and
Lv5 = {l5,1, l5,2}, i.e., v4 and v5 have 2 locations each, then
the spatial compatibility can be defined as Ql4,1 = {l5,2} and
Ql4,2 = {l5,1}. To ensure correct timing, two time lags need to
be added: sv5 ≥se4,1 +le4,1,v5 and sv4 ≥se5,6 +le5,6,v4 −CT.
The first time lag ensures that robot r2 takes the weldment
from the bench after r1 has put it there. The second time
lag guarantees that robot r2 has left the bench before r1 puts
another weldment there. The lengths of the time lags are set
to the required time for a robot to safely leave the bench.
III. MIXED-INTEGER LINEAR PROGRAMMING MODEL
The problem is intrinsically nonlinear due to f t
e being
convex functions; nevertheless, a MILP model can be formu-
lated if the functions are piecewise linearized. The piecewise
linear functions ˆf t
e are propagated through constraints (3) to
the criterion (1), where W is an upper bound on energy
consumption, yt
e is a binary variable that determines whether
or not the trajectory t ∈Te is selected for activity e ∈E,
B is a set of indices, and kt
e,b, qt
e,b are coefficients of the b-
th linear function approximating f t
e. The number of segments
|B| of each ˆf t
e can be adjusted to meet the desired accuracy
of the approximation. There is no need to introduce binary
variables for each segment of the function because functions
f t
e are convex, and the criterion ensures that the right linear
function is active. The energy consumption of static activities
is calculated by constraints (2), where xl
v determines whether
or not location l is selected for v ∈V , and zm
v is set to true
if and only if mode m ∈Mr is applied for activity v ∈V .
Min
X
∀a∈A
Wa
(1)
s.t.
pm
v,ldv −W
2 −zm
v −xl
v


≤Wv
(2)
∀r ∈R, ∀v ∈Vr, ∀l ∈Lv, ∀m ∈Mr
kt
e,bde + qt
e,b −W
1 −yt
e


≤We
(3)
∀e ∈E, ∀t ∈Te, ∀b ∈B
Assignment constraints (4), (5), and (6) ensure the correct
selection of locations, power saving modes, and trajectories
for activities. Note that optional activities EO are omitted in
constraints (6) because they may or may not be carried out.
If an optional activity e ∈EO is not performed, then none of
its trajectories is selected. Therefore, We is zero because all
the constraints (3) for activity e are deactivated, and criterion
(1) pushes We down to zero due to the minimization.
X
∀l∈Lv
xl
v = 1
∀v ∈V
(4)
X
∀m∈Mr
zm
v = 1
∀r ∈R, ∀v ∈Vr
(5)
X
∀t∈Te
yt
e = 1
∀e ∈E \ EO
(6)
Flow constraints (7) and (8) state that if the robot moves
to location l then it also has to move away from the same
location.
X
∀e∈Pv
X
∀t=(lfrom,l)∈Te
yt
e = xl
v
∀v ∈V, ∀l ∈Lv
(7)
X
∀e∈Sv
X
∀t=(l,lto)∈Te
yt
e = xl
v
∀v ∈V, ∀l ∈Lv
(8)
Constraints (9) to (13) enforce the time ordering of activi-
ties. Some precedences of activities are mandatory (see con-
straints (9) and (10)) due to the structure of graph Gr(Vr, Er),
whereas others are optional (see constraints (11) to (13))
because some e ∈EO do not have to be carried out.
The binary variables we,v determine whether or not activity
e ∈EO is performed and which order of operations is selected.
sa2 = sa1 + da1
(9)
∀a1 ∈A \ EO, ∀a2 ∈Sa1, ∄vh
r = a1
se = svh
r + dvh
r −CT
(10)
∀vh
r ∈V, ∀e ∈Svh
r
X
∀t∈Te
yt
e = we,suc(e)
∀e ∈EO
(11)
sv + (1 −we,v) CT ≥se + de
(12)
∀e ∈EO, ∀v ∈Se
sv −(1 −we,v) CT ≤se + de
(13)
∀e ∈EO, ∀v ∈Se
Restrictions on the duration of static and dynamic activities
are reflected in constraints (14) and (15), respectively.
IEEE TRANSACTIONS ON INDUSTRIAL INFORMATICS
6
max(dv, dm)zm
v ≤dv ≤dv
(14)
∀r ∈R, ∀v ∈Vr, ∀m ∈Mr
dt
eyt
e ≤de ≤d
t
e + CT
1 −yt
e


(15)
∀e ∈E, ∀t ∈Te
The problem described by Equations (1) to (15) is robot
independent; i.e., the constraint matrix of the MILP formu-
lation is block diagonal; thus, each robot can be treated
separately. In the following constraints, however, the robots
are coupled together with regard to time synchronization and
spatial compatibility.
sa2 −sa1 ≥la1,a2 −CTha1,a2
∀(a1, a2) ∈ETL
(16)
xl1
v1 ≤
X
∀l2∈Lv2,l2∈Ql1
xl2
v2
(17)
∀v1 ∈V, ∀l1 ∈Lv1, |Ql1| > 0
Collisions between robots are resolved by constraints (18)
and (19) where ugi
ai is either xgi
ai or ygi
ai depending on whether
ai ∈V or ai ∈E, and the binary variables cn
o determine
the execution order of ai and aj for the multiples of cycle
time CT if the collision applies. As an example, consider xg1
a1
and yg2
a2 to be a colliding pair for a multiple of cycle time
n = 2. If location g1 and movement g2 are selected, i.e., if
the variables ug1
a1 and ug2
a2 are equal to one, then exactly one
of constraints (18) and (19) is active for this collision pair
and multiple n = 2; therefore, either sa2 + 2CT ≥sa1 + da1
or sa1 ≥sa2 + da2 + 2CT has to be satisfied. In case the
collision does not apply, i.e., if location g1 or movement g2 is
not selected, then these two equations remain inactive.
sa2 + nCT + 2 |R| CT(3 −cn
o −ug1
a1 −ug2
a2) ≥sa1 + da1
∀o = (a1, g1, a2, g2) ∈C, ∀n ∈{−|R| , . . . , |R|}
(18)
sa1 + 2 |R| CT(2 + cn
o −ug1
a1 −ug2
a2) ≥sa2 + da2 + nCT
∀o = (a1, g1, a2, g2) ∈C, ∀n ∈{−|R| , . . . , |R|}
(19)
All the variables of the model are either nonnegative floats
or binary variables, as expressed in Equation (20).
Wa, sa, da ∈R≥0
xp
v, zm
v , yt
e, cn
o , we,v ∈B
(20)
IV. PARALLEL HEURISTIC ALGORITHM
Motivated by the inability of MILP solvers to solve bigger
instances, this study proposes a parallel hybrid heuristic based
on a linear programming (LP) solver that iteratively solves
partially fixed problems, called tuples, for selected locations,
power-saving modes, and alternative orders. Although the
application of an LP solver to partially fixed problems has been
previously studied (see, e.g., [21], [22], and [23]), the present
work accelerates the heuristic by using a tailor-made Gurobi
simplex algorithm for piecewise linear convex functions and
the parallelization. The justification for the parallelization of
combinatorial problems can be found in [24], [25], and [26].
The flowchart of the heuristic, as shown in Figure 3, can be
divided into two parts: the control thread and worker threads.
The control thread generates various alternative orders in the
form of HCact
r , launches worker threads, and occasionally
combines elite solutions to create promising tuples. After the
specified time limit tmax, the control thread joins the worker
threads and prints the best solution if found. The process
flow of the worker thread consists of reading and generating
tuples and then optimizing them iteratively by alternating
between solving a reduced LP problem and carrying out
three sub-heuristics. In general, each sub-heuristic performs
a small modification of the tuple that may reduce the energy
consumption, and an LP solver evaluates a real impact of this
modification. The iterative optimization is stopped if there is
no significant improvement after Φmax iterations, in which
case the next tuple is read, or if the time limit tmax is exceeded,
in which case the thread is terminated.
The heuristic (see Figure 3) is accelerated by using multiple
worker threads that independently process tuples one by one
from a list. If the list becomes empty, then new tuples are
added to it, and the process continues. Because the tuples, al-
ternative orders HCact
r , and elite solutions are accessible from
all the threads, it is necessary to use synchronization primitives
to guarantee consistent reads and writes. Fortunately, the
introduced overhead is negligible (see the experiments in
Section V) because the access time is very short compared
with the time needed to solve the reduced LP problem. The key
parts of the heuristic are detailed in the following subsections.
A. Generation of Alternatives
The order of operations of robot r, i.e., an alternative, is
encoded as a Hamiltonian circuit HCact
r
in graph Gr(Vr, Er).
Because there can be more alternative orders of operations
(i.e., circuits), it is necessary to find some of them. For this pur-
pose, the finding of Hamiltonian circuits in graph Gr(Vr, Er)
is transformed into a search for Hamiltonian paths in graph
G′
r(V ′
r, E′
r) as follows. Nodes V ′
r = Vr\vh
r ∪{v→, v←}, where
v→, v←are the starting and ending nodes, respectively. The
edges leaving vh
r are modified such that they start from v→;
similarly the edges entering vh
r are changed such that they
end in v←. Both the nodes and edges of G′
r are weighted by
the minimum possible duration da of the related activity; i.e.,
min∀t∈Te dt
e for ∀e ∈E′
r and max(dv, min∀m∈Mr dm) for
∀v ∈V ′
r. The goal is then to find random Hamiltonian paths
from v→to v←such that P
∀a∈path da ≤CT. To ensure effec-
tive pruning during an iterative random tree search from node
v→, the obviously infeasible or completely searched sub-paths
from v→are detected and skipped. The sub-path (v→, . . . , ai)
is infeasible if there is an unvisited node au that cannot be
reached or |(v→, . . . , ai)| + U(ai, au) + U(au, v←) > CT,
where the first term is the sub-path length, and U(ai, aj)
is the length of the shortest path from ai to aj calculated
by the Floyd-Warshall algorithm. If a Hamiltonian path is
found, its fastest sequence through the activity locations can
be obtained by applying a shortest path algorithm for directed
acyclic graphs. If the duration of this sequence is longer than
the robot cycle time, then the related Hamiltonian path is
neglected because it cannot lead to a feasible solution for the
original problem. Note that the Hamiltonian path and its fastest
IEEE TRANSACTIONS ON INDUSTRIAL INFORMATICS
7
start
generate HCact
r
launch worker
threads
wait for timer
join worker
threads
print best
solution
end
combine
elite solutions
timeout
otherwise
control thread
start
read next tuple
generate tuples
end
change
locations
(de)select
power mode
change
path
solve reduced
LP problem
Φ = 0
Φ = Φ + 1
add
constraint
no tuple
available
timeout
otherwise
otherwise
otherwise
read next
sub-heuristics
worker thread
elites
Fig. 3. Flowchart of the parallel hybrid heuristic.
sequence through the locations can be easily transformed back
to HCact
r
and HCloc
r , respectively. Because the algorithm is
exact, i.e., it can prove the nonexistence of a Hamiltonian
circuit of a given length, it is possible to detect the infeasibility
of some instances. Moreover, the generation of alternatives is
robot independent; therefore, the work is distributed among
threads to accelerate the initialization of the heuristic.
B. Generation of Tuples
The generation of tuples, which are formally defined be-
low, is crucial for finding initial feasible solutions; therefore,
the emphasis is placed on feasibility rather than on energy
optimality during the generation. At the beginning of the gen-
eration, a random alternative HCact
r
and its fastest closed path
HCloc
r
through the locations (see Section IV-A) are assigned
to each robot r. Only the fastest power-saving mode of robot
r is considered. The closed paths are subsequently modified
to meet the spatial compatibility as follows. For each pair
of static activities with the violated spatial compatibility one
fixing pair q ∈Q is selected with respect to the prolongation
of the related closed paths; i.e., the paths with a minimal
duration close to the robot cycle time are penalized. After the
spatial compatibility is fixed, the tuple is created as a triple
T = (A , P, α : V →M), where A , P are the selected
alternatives and their closed paths, respectively, and function
α maps each activity v ∈V to its power mode m ∈M,
where M
= ∪∀r∈RMr. The process of tuple generation
is applicable to two blocks shown in Figure 3: generate
tuples and combine elite solutions. However, in the second
block, only the alternatives HCact
r
used in elite solutions are
considered to intensify the search process.
C. Reduced Linear Programming Problem
The timing of a partial problem, i.e. a tuple, is determined
by the reduced LP problem. If the resulting solution is feasible,
then it is feasible for the original problem and can be added
to the list of elite solutions if it ranks among the top solutions
in terms of energy consumption. Otherwise, the solution is
infeasible, and one of two cases arises: the related tuple is not
modified by any sub-heuristics, in which case another tuple is
read; or the previous tuple resulting in a feasible solution is
loaded, and the next sub-heuristic is selected.
Before formally describing the reduced LP problem, it is
necessary to define some sets extracted from tuple T . Set ET
contains all the dynamic activities in ∀HCact
r
∈A , where each
dynamic activity e ∈ET has to be carried out (is mandatory)
and has assigned its fixed trajectory t as a pair (e, t) ∈F1
T ;
similarly, each static activity v ∈V is linked to its location
l and power mode m as a triple (v, l, m) ∈F2
T . Then, the
reduced LP problem can be stated as follows. Criterion (21)
minimizes the sum of the piecewise linear convex functions
ˆf t
e passed to the Gurobi LP solver as a list of function points
or alternatively expressed similarly to constraints (2) and (3)
for other solvers. The equations marked with an asterisk are
the same as the original ones (without the asterisk); however
the sets from which the constraints are generated are differ-
ent. Constraints (9)*, (10)*, and (16)* establish precedences
between activities AT = V ∪ET , and constraints (14)* and
(15)* transform into the domain specification of da variables.
IEEE TRANSACTIONS ON INDUSTRIAL INFORMATICS
8
Min
X
∀(e,t)∈F1
T
ˆf t
e(de) +
X
∀(v,l,m)∈F2
T
pm
v,ldv
(21)
subject to: (9)*, (10)*, (14)*, (15)*, (16)*
sa2 + nCT ≥sa1 + da1
∀(a1, a2, n) ∈D≥
(22)
sa2 + da2 + nCT ≤sa1
∀(a1, a2, n) ∈D≤
(23)
All the aforementioned constraints are always present in
the above formulation; however, these alone may not ensure
feasibility because some collisions could occur. For this rea-
son, additional constraints (22) and/or (23) could be iteratively
added to heuristically resolve active collisions. Thus, if a
solution is feasible with respect to the current constraints of
the reduced LP problem but is not feasible for the original
problem due to some active collisions, then the worst collision
(formally defined later) is detected and resolved by adding
constraint (22) or (23). Afterward, the problem is resolved,
and the same procedure is repeated if collisions still occur
(see Figure 3).
To introduce the collision resolution more formally, con-
straints (22) and (23) have to be shown as specialized forms
of constraints (18) and (19), respectively. First, the variables
ugi
ai in constraints (18) and (19) are fixed for tuple T , i.e.,
ugi
ai = 1 if (ai, gi) ∈F1
T or (ai, gi, mi) ∈F2
T ; otherwise,
ugi
ai
= 0. Second, it is sufficient to consider set CT
=
{∀(ai, aj) : ugi
ai + ugj
aj = 2, (ai, gi, aj, gj) ∈C, ai, aj ∈AT }
instead of C because constraints (18) and (19) are neither
binding nor violated for unselected locations and movements
of tuple T , and each activity ai ∈AT has an assigned
movement or location. Note that if ugi
ai + ugj
aj = 2, then the
variable cn
o makes only one constraint, either (18) or (19),
active (“big M method”) depending on the decision whether
sa2 + nCT ≥sa1 + da1 or sa1 ≥sa2 + da2 + nCT, which
guarantees collision-free ordering. However, these decisions
correspond exactly to constraints (22) and (23), in which the
decision on ordering is given by adding a triple (ai, aj, n) to
the related set D≥or D≤, respectively; i.e., a constraint is
added to the problem. The constraint added is not removed
during the solving of the reduced LP problem; therefore, this
approach is heuristic. The key question is how to select which
constraint to add and how such constraint is related to the
worst collision. To find an answer, the maximal violation Γ
of constraints (18) and (19) for tuple T has to be defined as
follows.
υn
ai,aj = sai + dai −saj −nCT
(24)
µn
ai,aj = saj + daj + nCT −sai
(25)
Γ =
max
∀(ai,aj)∈CT
∀n∈{−|R|,··· ,|R|}
min(υn
ai,aj, µn
ai,aj)
(26)
Note that if Γ ≤0, then there are no active collisions,
and no extra constraints are needed; failing that, i.e., Γ > 0,
means that there is a pair of colliding activities that is not
time disjunctive. In that case, let (a∗
i , a∗
j, n∗) be an optimal
argument of the max function in Equation (26) (replace max
with argmax). This triple defines the worst collision, i.e., the
collision with the biggest time intersection, occurring in a
current solution. Whether this collision should be resolved by
adding a constraint of type (22) or (23) is determined by the
υ∗= υn∗
a∗
i ,a∗
j and µ∗= µn∗
a∗
i ,a∗
j values. If υ∗≤µ∗(µ∗≤υ∗),
then the worst collision is resolved by adding the triple to D≥
(D≤) because it may result in smaller changes in timing after
the problem is resolved with the added constraints.
D. Sub-heuristics
The aim of the sub-heuristics is to modify a given tuple
T and its timing calculated by the reduced LP problem such
that the energy consumption would be reduced in successive
LP calls. Because the performed modifications can result in a
violation of time lags or the occurrence of collisions, a penalty
based on the duration of breakage and average input power is
added to these modifications. If the modifications of an active
sub-heuristic do not lead to significant energy improvements,
then the next sub-heuristic is selected in round-robin order.
The goal of the (de)select power mode sub-heuristic, which
focuses on the application of the power-saving modes of
robots, is to find and apply an alternative power-saving mode
for some activity v ∈V . To select a suitable mode and activity,
the energy consumption is estimated for all applicable modes
of each v ∈V as follows. If a different mode m ∈Mr of robot
r ∈R is applied for activity v′, then some activities v ∈Vr are
uniformly shortened/prolonged to meet the cycle time. After
the timing is modified, both the energy consumption and the
above-mentioned penalty are determined. Finally, the choice
falls on mode m and activity v′ with the lowest sum of the
energy consumption and penalty. Note that the sub-heuristic
uses a Tabu list, i.e., a short-term memory with forbidden
modifications, to avoid cycling; therefore, some power-saving
modes may not be applicable for some activities.
The change locations sub-heuristic optimizes the closed
paths of robots, i.e., ∀HCloc
r
∈P, by modifying the go-
through locations as follows. Let ‘99K la
t1
−→lb
t2
−→lc 99K’ be
a part of HCloc
r , where t1 and t2 are the movements between
the locations; then, the question arises as to whether the inner
part, i.e., ‘
t1
−→lb
t2
−→’, can be replaced so as to achieve a
reduction in energy consumption. To find an answer, each
viable substitution ‘la
t⊢
−→l◦
t⊣
−→lc’ has to be evaluated in
terms of energy by solving the convex problem (27) to (28),
where t1, t⊢∈Te1, t2, t⊣∈Te2, e1, e2 ∈ET , lb, l◦∈Lvb,
vb ∈V , and K = dLP
e1 +dLP
e2 and dLP
vb are constants determined
from the LP solution. Because the problem is convex and one-
dimensional, the golden section search algorithm can be used
to find the optimal solution. Then, the most energy-friendly
substitution not breaking the spatial compatibility is applied,
and the process is repeated for other sub-paths.
minimize f t⊢
e1 (de1) + f t⊣
e2 (K −de1) + pm
vb,l◦dLP
vb
(27)
max(dt⊢
e1, K −d
t⊣
e2) ≤de1 ≤min(d
t⊢
e1, K −dt⊣
e2)
(28)
The last sub-heuristic, called change path, diversify the
search process to allow exploring some otherwise unreach-
able HCloc
r . The sub-heuristic selects one HCloc
r
∈P and
IEEE TRANSACTIONS ON INDUSTRIAL INFORMATICS
9
TABLE I
PERFORMANCE OF LP SOLVERS FOR THE HEURISTIC.
lp solve
Gurobi
Gurobi CF
S5
sequential
23.1
78.5
365
parallel
304
1111
4708
M8
sequential
7.9
44.2
170
parallel
93.4
488
2134
L12
sequential
3.3
30.8
97.9
parallel
42.6
371
1084
100
101
102
103
104
105
time [s]
75000
80000
85000
90000
95000
100000
energy [J]
sequential heuristic
parallel heuristic
parallel Gurobi ILP solver
Fig. 4. Progress of the heuristic and MILP solver on M8 8 instance.
randomly changes its go-through locations such that spatial
compatibility is achieved and the resulting closed path exists.
V. EXPERIMENTAL RESULTS
To evaluate the effectiveness and performance of the pro-
posed approaches, the algorithms were tested on benchmark
instances with 5, 8, and 12 robots. The experiments were
carried out on a Linux server with two Intel Xeon E5-2620 v2
2.10GHz processors (2 x 6 physical cores + hyper-threading)
and 64 GB of DDR3 memory, in which Gurobi 6.0.4, a
state-of-the-art MILP solver, and lp solve 5.5.2.0, an open-
source MILP solver, were installed to solve LP and MILP
problems. The optimization program, written in C++11 and
compiled by GCC 4.9.2, was benchmarked on three generated
data sets: S5 (small, 5 robots), M8 (medium, 8 robots), and
L12 (large, 12 robots), each of which contains 10 instances
of the problem. For the purposes of all experiments, each
energy function f t
e was approximated by 10 linear functions,
i.e., |B| = 10, and the minimal number of optimization
iterations per tuple Φmax was set to 100, 600, and 1000 for the
S5, M8, and L12 data sets, respectively. The data sets, their
configuration files, and the generator are publicly available on
https://github.com/CTU-IIGto enable comparisons with other
research.
The first experiment measures the effect of the paralleliza-
tion and Gurobi simplex method on the performance of the
hybrid heuristic. A good indicator of the performance is the
number of LP evaluations per second because more than 90 %
of the computational time is consumed by LP evaluations,
i.e., the building of the LP problem, its optimization, and the
extraction of a solution. Table I shows the average number
of LP evaluations per second for each data set. The figures
indicate that the parallel heuristic is about 12 times faster
than the sequential one, and the specialized Gurobi simplex
method, denoted as ‘Gurobi CF’, accelerates the heuristic
about 3 to 4 times compared with the regular one; therefore,
an overall speedup of about 36 to 48 can be expected for the
aforementioned configuration. For comparison, Table I also
presents the corresponding values for the lp solve open-source
solver. To show that a higher number of evaluated tuples
also has a positive impact on the quality of solutions, the
dependence of the criterion value, i.e., energy consumption,
on the time limit was plotted in Figure 4. The results on
an instance with 8 robots revealed that the parallel heuristic
with 24 threads (12 cores + hyper-threading) converged sig-
nificantly faster than the sequential version; therefore, a similar
solution quality was achievable in a fraction of the time. In
comparison to the Gurobi MILP solver, the heuristic seems
to be stronger in finding feasible solutions, compare 1.5 h
with 3.1 s (5.3 s) required by the parallel (sequential) heuristic,
and is more suitable for a short-term (re)optimization. On
the other hand, if the MILP solver is given enough time,
then better solutions may be found for the medium instances
(see Figure 4). The same experiment was repeated on L12 9
instance with 12 robots. The parallel (sequential) heuristic
found the first feasible solution in 10 m (4 h) and the best
criterion was 177276 J (186855 J) after the 10-h time limit.
The Gurobi MILP solver had run out of memory (installed
64 GB of memory) in less than 4 hours without having any
feasible solution.
The second experiment evaluates the quality of solutions in
terms of energy consumption (lower energy consumption is
better) for both the parallel heuristic and the parallel Gurobi
MILP solver. The quality of obtained solutions is compared
with a tight lower bound, calculated as the sum of the
lower bounds on energy consumption of individual robots in
which global constraints are neglected (see Equations (1) to
(15)), to obtain an upper estimate of the gap from optimal
solutions. The run time for the bound calculation was limited
to 5 hours. To ensure that the measured data are statistically
significant, the best, average, and worst qualities of solutions
are determined from 10 runs for data sets S5 and M8. Because
the biggest data set (L12) is too computationally demanding,
only 3 measurements were carried out for each instance, and
the criterion values were stated explicitly. Surprisingly, the
Gurobi MILP solver has a deterministic behavior and thus
provided the same quality of solutions for all runs. That is
why only the average is stated for the Gurobi MILP solver.
Tables II to IV show the qualities of solutions for the
S5, M8, and L12 data sets, respectively. Instance S5 2 is
omitted from the results because it was proved to be infeasible.
For small instances, the best solutions found by the heuristic
and the MILP solver were near the optimal ones because
the average gap from the lower bound was about 2 %. On
one hand, the Gurobi MILP solver could find almost optimal
solutions for some small instances, on the other hand, it had
difficulty in providing feasible solutions even for instances
IEEE TRANSACTIONS ON INDUSTRIAL INFORMATICS
10
TABLE II
QUALITY OF SOLUTIONS IN TERMS OF ENERGY CONSUMPTION FOR THE S5 DATA SET.
heuristic (tmax = 30 s)
MILP (tmax = 30 s)
heuristic (tmax = 1 h)
MILP (tmax = 1 h)
instance
best run
avg run
worst run
avg run
best run
avg run
worst run
avg run
lower bound
S5 0
40654.80
40788.45
40936.20
42521.9
40494.90
40637.85
40776.60
40272.8
39849.4
S5 1
31195.10
31335.18
31479.20
34582.6
31052.90
31111.81
31184.60
30623.3
30208.2
S5 3
47482.80
47665.48
47861.60
–
47434.10
47483.28
47674.90
–
46374.3
S5 4
47826.30
48039.38
48326.60
–
47541.10
47739.06
48067.60
–
45922.5
S5 5
41692.00
41902.74
42081.20
–
41654.70
41699.96
41741.00
41498.8
40349.3
S5 6
34920.30
35043.71
35128.20
36315.6
34760.00
34856.44
34918.90
34720.9
34398.7
S5 7
42511.60
42821.70
42987.60
45670.2
42384.40
42515.82
42627.60
42033.0
40801.4
S5 8
39181.00
39423.62
39602.40
41440.2
39174.00
39301.44
39464.10
38364.8
37772.9
S5 9
38864.90
39067.57
39235.90
39721.5
38649.80
38840.31
38924.80
38283.8
37923.5
TABLE III
QUALITY OF SOLUTIONS IN TERMS OF ENERGY CONSUMPTION FOR THE M8 DATA SET.
heuristic (tmax = 600 s)
MILP (tmax = 600 s)
heuristic (tmax = 1 h)
MILP (tmax = 1 h)
instance
best run
avg run
worst run
avg run
best run
avg run
worst run
avg run
lower bound
M8 0
86831.60
90580.66
95333.80
–
87365.40
90258.81
94670.50
–
78641.3
M8 1
88098.90
88491.50
88883.80
89297.4
87999.70
88236.49
88488.00
89182.7
83677.7
M8 2
90048.70
90667.98
91463.60
92991.0
89939.70
90600.33
91804.20
89744.0
82503.9
M8 3
83412.60
83948.53
84942.60
–
83224.70
83709.86
84741.30
83621.6
78564.3
M8 4
76914.70
77549.80
78191.60
82881.8
77058.90
77400.65
77723.00
76582.8
70838.1
M8 5
89239.40
90337.15
91486.00
–
89369.40
89873.88
90545.30
–
81856.3
M8 6
95314.50
96021.27
96629.40
97875.4
95638.50
95836.32
96134.60
94918.1
88023.1
M8 7
83777.20
85218.11
86951.50
–
83016.50
84881.71
86266.50
–
77038.8
M8 8
78514.40
79220.31
80530.70
–
78465.80
78859.41
79342.40
–
74447.2
M8 9
92875.60
93309.59
93736.40
93951.9
91960.60
92942.01
93654.10
92521.9
86900.0
TABLE IV
QUALITY OF SOLUTIONS FOR THE L12 DATA SET.
heuristic (tmax = 3 h)
instance
3 runs
lower bound
L12 0
–, –, 204464
167229
L12 1
–, 177731, –
154744
L12 2
211019, –, 211377
173457
L12 3
–, –, 188527
155121
L12 4
187163, –, 171860
148441
L12 5
–, –, –
–
L12 6
–, –, –
–
L12 7
190060, 202737, –
155598
L12 8
218844, –, –
176156
L12 9
173743, 176969, 176567
152531
with 5 robots. Compared with the MILP solver, the heuristic
obtained very good solutions in a short time, and its ability to
find feasible solutions seemed to be significantly better. The
results for medium instances indicated a similar trend. The
heuristic solved all the instances, whereas the MILP solver
found a solution for only 6 instances in the 1-h time limit.
Moreover, the quality of solutions achieved by the heuristic
in the 1-h time limit was comparable with that of the MILP
solver, and the average gap increased to about 7.5 % for the
best solutions found. The largest instances were intractable for
the MILP solver; therefore, only the results for the heuristic
are presented. The heuristic solved 8 of 10 instances in 3 runs
for the time limit of 3 hours. Two instances remained unsolved
because they were either too difficult to solve or infeasible.
The last experiment shows how the robot cycle time, which
was scaled by 1.0, 1.1, and 1.2 factors, respectively, influences
the performance of the heuristic and MILP solver. The sum-
mary of results for the S5 data set is in Table V, where each
figure is the average criterion value from 10 runs with tmax =
TABLE V
DEPENDENCE OF THE QUALITY OF SOLUTIONS ON THE CYCLE TIME.
CT
1.1 ∗CT
1.2 ∗CT
instance
heur.
MILP
heur.
MILP
heur.
MILP
S5 0
40717
40355
40435
40022
41381
40895
S5 1
31198
30812
29410
29269
29433
29239
S5 3
47579
–
44177
44338
43846
43284
S5 4
47749
–
43993
44233
43577
43400
S5 5
41730
–
38248
38227
37644
37566
S5 6
34880
34980
32359
32186
31744
31630
S5 7
42560
42619
42230
41867
42737
42571
S5 8
39360
38761
38407
38045
39051
38432
S5 9
38935
38412
38769
38233
39695
39302
600 s. The outcomes indicate that the heuristic outperforms the
MILP solver if the cycle time is tight, i.e., a feasible solution
is hard to find due to the limited time for the movements and
operations. However, if the cycle time is prolonged, then the
MILP solver gets ahead because its ability to find optimal
solutions by a systematic search becomes dominant for the
given time limit.
The approach proposed in this study cannot be directly
compared with existing works [16], [20], [15] because the
robot cycle time is considered instead of the work cycle
time. However, in general, much bigger instances were solved
compared with the aforementioned works, which considered
only one to four robots per robotic cell. Besides, the proposed
approach took into account additional optimization aspects,
such as the robot power-saving modes and positions. Finally,
the algorithms were tested on more problem instances than
those in others works. The data sets are available for others.
IEEE TRANSACTIONS ON INDUSTRIAL INFORMATICS
11
VI. CASE STUDY FROM ˇSKODA AUTO
This case study shows a potential impact of the optimiza-
tion on the energy consumption of existing robotic cells by
considering a long-operating robotic cell from ˇSkoda Auto
(see Figure 1 for a screenshot from the simulation), in
which a part of an automotive body is welded, glued, and
assembled by 6 industrial robots with a robot cycle time
of 56 seconds. The timing of individual robotic operations
was obtained from robotic programs. An alternative way is
to measure and subsequently identify the movements. The
optimized speed of movements can be easily entered into the
robotic programs. The energy function f t
e of a trajectory was
fitted from the points obtained from simulations in Siemens
Tecnomatix Process Simulate (Digital Factory - Industry 4.0),
in which the robot controller supported the calculation of the
energy consumption of movements according to the Realistic
Robot Simulation standard. It is also possible to carry out
the measurements with a physical robot to obtain the energy
functions; however, this is impractical in existing robotic
cells. Alternative approaches to the modeling and simulation
are found in [27], [28]. To ensure the repeatability of the
production process in terms of output quality, the welding,
gluing, and assembling operations remained the same, i.e.,
the fixed duration and energy consumption were extracted
from the measured power profiles. Only the robot speeds and
power saving-modes (at home position) were addressed in
the optimization because minimal intervention is desirable for
existing robotic cells. More information about the structure of
the robotic cell, the timing, and the synchronizations can be
obtained from a published instance file.
Based on the results, it was estimated that the original en-
ergy consumption of 500 kJ (maximal speeds, without power-
saving modes) could be decreased to 391 kJ (reduced speeds,
with power-saving modes) per cycle, resulting in about 20 %
energy savings. The power-saving modes of robots saved about
2.4 % of energy, whereas the remaining savings were attributed
to the optimization of speeds. If the cycle time is extended to
70 s and 80 s, then the application of energy-saving modes
would improve the consumption by about 6.8 % and 12 %,
respectively, compared with their nonuse. This finding may
be particularly useful during over production or production
cuts. Note that the inclusion of the bus power-off mode
is not always straightforward because it requires interaction
between the robot and a superior controller, which may not be
ready in existing cells. However, the effort to implement such
interaction is not high either. The outcomes of the optimization
indicate that a significant reduction in energy consumption can
be achieved for existing robotic cells and that even more can
be expected for planned robotic cells that will maximize the
potential of the optimization algorithm.
VII. CONCLUSION
Energy optimization is undoubtedly a current and important
problem for the industry because such optimization could lead
to a significant decrease in costs. Besides being able to save
money, an involved company could also improve its green
credentials and become more competitive.
This work presents a holistic approach to the energy op-
timization of robotic cells that considers many optimization
aspects. A universal mathematical model for describing robotic
cells is proposed, from which is derived an MILP formulation
that is directly solvable by MILP solvers. However, these
solvers can solve only small instances. Therefore, a hybrid
parallel heuristic is devised for instances with up to 12 robots.
The strength of the heuristic is that its performance scales
almost linearly up to 12 cores; thus, significant acceleration
is achievable on modern processors. Moreover, the heuristic
uses a specialized Gurobi simplex method for piecewise linear
convex functions that is about 3 or 4 times faster than
broadly used simplex methods. The merits of the heuristic
are verified by experiments, the results of which showed
that the heuristic found a solution to 27 of 29 problems
(+1 infeasible) and clearly outperformed the MILP solver
for the largest problems. Note that all the algorithms, the
generator, and the generated data sets are publicly available
from https://github.com/CTU-IIG; thus, other researchers can
experiment with the heuristic and explore the documented
open-source code.
Finally, a real robotic cell from ˇSkoda Auto is optimized.
The results indicate that energy savings of up to 20 % can
be achieved merely by changing the robot speeds and apply-
ing power-saving modes. In future studies, we plan to fully
optimize another robotic cell and integrate our approach into
industrial software.
ACKNOWLEDGMENT
This work was supported by the Grant Agency of the Czech
Republic under the Project FOREST GACR P103-16-23509S.
The authors are also very grateful to Antonin Novak from
the Czech Technical University in Prague for his advice of
using special Gurobi simplex method.
REFERENCES
[1] M. W. Dawande, H. N. Geismar, S. P. Sethi, and C. Sriskandarajah,
Throughput Optimization in Robotic Cells, 1st ed., ser. International
Series in Operations Research & Management Science.
Springer US,
2007, vol. 101.
[2] Energy
Efficiency
–
European
Commision.
[Online].
Available:
http://ec.europa.eu/energy/en/topics/energy-efficiency
[3] D. Meike and L. Ribickis, “Energy Efficient Use of Robotics in
the Automobile Industry,” in Advanced Robotics (ICAR), 2011 15th
International Conference on, June 2011, pp. 507–511.
[4] D. Meike, M. Pellicciari, and G. Berselli, “Energy Efficient Use of
Multirobot Production Lines in the Automotive Industry: Detailed Sys-
tem Modeling and Optimization,” Automation Science and Engineering,
IEEE Transactions on, vol. 11, no. 3, pp. 798–809, July 2014.
[5] PROFIenergy homepage. [Online]. Available: http://www.profibus.com/
technology/profienergy/
[6] U.
Schnabl.
Auch
Roboter
k¨onnen
jetzt
abschalten.
[On-
line].
Available:
http://www.szbz.de/nachrichten/artikel/detail/?tx
szbzallinone szbznews[news]=454255
[7] M. Pellicciari, G. Berselli, F. Leali, and A. Vergnano, “A Minimal
Touch Approach for Optimizing Energy Efficiency in Pick-and-Place
Manipulators,” in Advanced Robotics (ICAR), 2011 15th International
Conference on, June 2011, pp. 100–105.
[8] K. Paes, W. Dewulf, K. V. Elst, K. Kellens, and P. Slaets, “Energy
Efficient Trajectories for an Industrial ABB Robot,” Procedia CIRP,
vol. 15, no. 0, pp. 105 – 110, 2014, 21st CIRP Conference on Life
Cycle Engineering.
[9] P. Boˇzek, “Robot Path Optimization for Spot Welding Applications
in Automotive Industry.” Tehnicki vjesnik / Technical Gazette, vol. 20,
no. 5, pp. 913 – 917, 2013.
IEEE TRANSACTIONS ON INDUSTRIAL INFORMATICS
12
[10] O. Wigstrom and B. Lennartson, “Integrated OR/CP optimization for
Discrete Event Systems with nonlinear cost,” in 2013 IEEE 52nd Annual
Conference on Decision and Control (CDC), Dec 2013, pp. 7627–7633.
[11] M. Mashaei and B. Lennartson, “Energy Reduction in a Pallet-
Constrained Flow Shop Through On–Off Control of Idle Machines,”
IEEE Transactions on Automation Science and Engineering, vol. 10,
no. 1, pp. 45–56, Jan 2013.
[12] M. Foumani, I. Gunawan, K. Smith-Miles, and M. Y. Ibrahim, “Notes
on Feasibility and Optimality Conditions of Small-Scale Multifunction
Robotic Cell Scheduling Problems With Pickup Restrictions,” IEEE
Transactions on Industrial Informatics, vol. 11, no. 3, pp. 821–829,
June 2015.
[13] I. Sindicic, S. Bogdan, and T. Petrovic, “Resource Allocation in Free-
Choice Multiple Reentrant Manufacturing Systems Based on Machine-
Job Incidence Matrix,” IEEE Transactions on Industrial Informatics,
vol. 7, no. 1, pp. 105–114, Feb 2011.
[14] D. Meike, M. Pellicciari, G. Berselli, A. Vergnano, and L. Ribickis,
“Increasing the energy efficiency of multi-robot production lines in the
automotive industry,” in Automation Science and Engineering (CASE),
2012 IEEE International Conference on, Aug 2012, pp. 700–705.
[15] O. Wigstrom and B. Lennartson, “Sustainable Production Automation
- Energy Optimization of Robot Cells,” in Robotics and Automation
(ICRA), 2013 IEEE International Conference on, May 2013, pp. 252–
257.
[16] O. Wigstrom, B. Lennartson, A. Vergnano, and C. Breitholtz, “High-
Level Scheduling of Energy Optimal Trajectories,” Automation Science
and Engineering, IEEE Transactions on, vol. 10, no. 1, pp. 57–64, Jan
2013.
[17] M. Ayala, A. Benabid, C. Artigues, and C. Hanen, “The resource-
constrained modulo scheduling problem: an experimental study,” Com-
putational Optimization and Applications, vol. 54, no. 3, pp. 645–673,
2013.
[18] P. Sucha and Z. Hanzalek, “A cyclic scheduling problem with an
undetermined number of parallel identical processors,” Computational
Optimization and Applications, vol. 48, no. 1, pp. 71–90, 2011.
[19] Paryanto, M. Brossog, M. Bornschlegl, and J. Franke, “Reducing the
energy consumption of industrial robots in manufacturing systems,” The
International Journal of Advanced Manufacturing Technology, vol. 78,
no. 5-8, pp. 1315–1328, 2015.
[20] A. Vergnano, C. Thorstensson, B. Lennartson, P. Falkman, M. Pellicciari,
C. Yuan, S. Biller, and F. Leali, “Embedding detailed robot energy
optimization into high-level scheduling,” in Automation Science and
Engineering (CASE), 2010 IEEE Conference on, Aug 2010, pp. 386–
392.
[21] D. Kim and H. Shin, “A hybrid heuristic approach for production plan-
ning in supply chain networks,” The International Journal of Advanced
Manufacturing Technology, vol. 78, no. 1-4, pp. 395–406, 2015.
[22] B. Haddar, M. Khemakhem, S. Hanafi, and C. Wilbaut, “A hybrid
heuristic for the 0-1 Knapsack Sharing Problem,” Expert Systems with
Applications, vol. 42, no. 10, pp. 4653 – 4666, 2015.
[23] A. French and J. Wilson, “An LP-Based Hybrid Heuristic Procedure
for the Generalized Assignment Problem with Special Ordered Sets,”
in Hybrid Metaheuristics, ser. Lecture Notes in Computer Science,
M. Blesa, C. Blum, A. Roli, and M. Sampels, Eds.
Springer Berlin
Heidelberg, 2005, vol. 3636, pp. 12–20.
[24] W. Bo˙zejko, J. Pempera, and C. Smutnicki, “Parallel tabu search
algorithm for the hybrid flow shop problem,” Computers & Industrial
Engineering, vol. 65, no. 3, pp. 466 – 474, 2013.
[25] L. Bukata, P. ˇS˚ucha, and Z. Hanz´alek, “Solving the Resource Con-
strained Project Scheduling Problem using the parallel Tabu Search
designed for the CUDA platform,” Journal of Parallel and Distributed
Computing, vol. 77, no. 0, pp. 58 – 68, 2015.
[26] I. Chakroun and N. Melab, “Towards a heterogeneous and adaptive
parallel Branch-and-Bound algorithm,” Journal of Computer and System
Sciences, vol. 81, no. 1, pp. 72 – 84, 2015.
[27] L. Mikov´a, M. Kelemen, I. Virgala, and M. Michna, “Simulation
model of manipulator for model based design,” Applied Mechanics and
Materials, vol. 611, pp. 175–182, 2014.
[28] E. Prada, I. Virgala, G. Granosik, A. Gmiterko, and ˇS. Mrkva, “Sim-
ulation Analysis of Pneumatic Rubber Bellows for Segment of Hyper-
Redundant Robotic Mechanism,” Applied Mechanics and Materials, vol.
611, pp. 10–21, 2014.
Libor Bukata received his Master’s degree in Soft-
ware Engineering from the Czech Technical Univer-
sity in Prague (CTU) in 2012, where he currently
conducts the research on the energy optimization
of robotic cells as a post-graduate student. Up to
now, he has published 1 journal paper (Journal of
Parallel and Distributed Computing) and more than
3 contributions to international conferences (PDP
2013, SCOR 2014, ECCO 2014, MISTA 2015).
Moreover, he reviewed some papers from MISTA
2015 conference, where he also participated in or-
ganization. His research reflects his passion for scheduling, combinatorial
optimization, and parallel and distributed systems.
Pˇremysl ˇS˚ucha graduated in Technical Cybernet-
ics from the Czech Technical University in Prague
(CTU) in 2003 and received a PhD degree in Electri-
cal Engineering and Informatics from CTU in 2007.
From 2011 to 2012 he was in a postdoc position
at LAAS-CNRS, Toulouse, France. Since 2012 he
works as an Assistant Professor at CTU. His main
research activities are in the area of scheduling al-
gorithms and parallel programming. He participated
in several industrial projects (with e.g. Skoda, Air
Navigation Services of the Czech Republic) and
he was a national coordinator in the Design, Monitoring and Operation of
Adaptive Networked Embedded Systems project (DEMANES, ARTEMIS
Joint Undertaking). He published more than 40 papers in international journals
and conferences; he was a co-chair of the Multidisciplinary International
Scheduling Conference (MISTA 2015).
Zdenˇek Hanz´alek graduated in Electrical Engi-
neering at the Czech Technical University (CTU)
in Prague in 1990. He obtained his Ph.D. degree
in Industrial Informatics from the Universite Paul
Sabatier Toulouse and Ph.D. degree in Control En-
gineering from the CTU. He worked on optimization
of parallel algorithms at LAAS CNRS in Toulouse
and on discrete event dynamic systems at LAG INPG
in Grenoble. He founded and coordinated the Indus-
trial Informatics Group at CTU Prague focusing on
scheduling, combinatorial optimization algorithms,
real-time control systems and industrial communication protocols. From 2011
to 2014, he served as Senior Manager of newly established Mechatronics
group at Porsche Engineering Services in Prague. He has been involved in
many industrial contracts (e.g. Skoda, Porsche, Volkswagen, Rockwell, Air
Navigation Services, EATON). He is currently working on different aspects
of scheduling and optimization problems used in production and in time-
triggered communication protocols.
Pavel Burget received the M.Sc. degree in computer
engineering at Czech Technical University (CTU)
in Prague, Czech Republic in 1996. In 2008 he
received his Ph.D. degree in Control Engineering
and Robotics at the CTU in Prague. He is assistant
professor at the Department of Control Engineering
at CTU in Prague. He has been taking teaching and
research responsibilities in the areas of industrial
communications and control systems. He has had
several publications in journals (e.g. IEEE T on
Industrial Informatics, Control Engineering Practice)
and conferences focusing on industrial control systems (e.g. ETFA, CASE,
ECRTS, IFAC World Congress).