Energy Cost-Efficient Task Positioning in Manufacturing Systems

A problem to determine a production schedule which minimises the cost of energy used for manufacturing is studied. The scenario assumes that each production task has assigned constant power consumption, price of power from conventional electrical grid system is defined by time-of-use tariffs, and a component of free of charge renewable energy is available for the manufacturing system. The objective is to find the most cost-efficient production plan, subject to constraints involving predefined precedence relationships between the tasks and a bounded makespan. Two independent optimisation approaches have been developed, based on significantly different paradigms, namely mixed-integer linear programming and tabu search metaheuristic. Both of them have been verified and compared in extensive computational experiments. The tabu search-based approach has turned out to be generally more efficient in the sense of the obtained objective function values, but advantages of the use of linear programming have also been identified. The results confirm that it is possible to develop efficient computational methods to optimise energy cost under circumstances typical of manufacturing companies. The set of numerous benchmark instances and their solutions have been archived and it can be reused in further research.


Introduction
Energy cost-efficient manufacturing supports competitiveness and sustainability of contemporary industry. There are many areas in which energy-related cost can be optimised: energy sources, plant equipment, manufacturing technology, choice of raw materials and so on. In this work, energy cost-efficient production task scheduling is studied. This subject is comprehensively reviewed in [1,2].
There are numerous studies on optimisation of the cost of energy drawn for manufacturing from conventional electrical grid systems (EGS). They can be roughly split into two groups. In the first group, there are works concerning the optimisation based on schedule-dependent energy cost of the production [3][4][5][6][7]. In particular, the manufacturing with machines of variable speed [8,9] or having multiple working states [10][11][12][13] is considered in this group, where the speed or state determines instantaneous machine power consumption and the related cost. The second group includes problems related to manufacturing cost optimisation under time-of-use (TOU) price tariffs of EGS energy. This subject is quite well known and studied for a few decades, but it is continuously extended by new research in which some specific manufacturing environments and conditions are considered [10][11][12][14][15][16][17][18] or new approaches are proposed, e.g., optimal production planning for continuous power-intensive processes [19] or scheduling of machine unavailability related to preventive maintenance [20]. One can notice that a few of the mentioned works belong to the two groups simultaneously, namely [10][11][12].
Another research area comprises minimisation of energy cost in manufacturing systems by planning consumption of power from renewable energy sources (RES) on the basis of short-term forecast of its availability [21,22]. Intermittent character of the power from such sources is the main challenge in this subject. Therefore, production tasks have to be aligned to cost-efficient execution time windows.
The above-mentioned perspectives need to be merged to obtain practically useful solutions. A growing number of manufacturing environments simultaneously use energy from EGS charged on the basis of TOU tariffs and consume energy from RES. However, such a scenario seems to be an emerging research field and there are only a few related works [23][24][25]. Therefore, the respective knowledge is incomplete. In this work, two knowledge gaps are identified and addressed. First, there is no common useful problem formulation which could be an actual bridge between real-world planning software and dedicated optimisation algorithms. The algorithm proposed by Cui et al. [23] uses quite specialised parameters and it is dedicated for the Bernoulli production line; hence, it could be inadequate for modern production scheduling dealing with flexible and reconfigurable infrastructures. The solution introduced by Golari et al. [24] is designed for cost-efficient balancing of energy from RES and EGS on the level of production plants and longer time periods and it does not support detailed task-level scheduling. Finally, Küster et al. [25] do not specify precisely a scheduling problem model, but rather consider architecture of their original multi-agent scheduling software. The second gap is that none of the three related works analyse practical computational efficiency of the proposed algorithms. Exemplary use cases are only provided, based on individual and relatively small problem instances.
In this work, the first gap is filled by formulation of a simple and legible, but comprehensive optimisation model of energy cost-efficient task positioning problem (ECETPP) which combines the fundamental data regarding correlated elements: a manufacturing system, EGS and RES ( Figure 1). The TOU tariff price of EGS is described by a piecewise constant time function (Figure 1 plot A). The forecast available RES power is represented by a piecewise linear time function (Figure 1 plot B). For each production task, the model includes: processing time, consumed power and predecessors in the processing order ( Figure 1 charts C1-C3). A scheduling horizon H is also provided as a model parameter. The goal of the ECETPP optimisation is to minimise the cost of energy drawn from the EGS, maintaining completion time of the schedule not greater than the horizon H.
Two algorithms have been implemented for the ECETPP, based on the mixed-integer linear programming (MILP) and the tabu search (TS) metaheuristic. To address the second identified gap, extensive computational experiments have been performed involving 182 benchmark instances. The number of tasks, being the main determinant of an instance size, varies between 200 and 2000.
To consider the concept of the ECETPP closer, assume that the tasks presented in Figure 1 have the precedence constraints: (1, 2), (1, 5), (4, 2), (4, 5), (5, 6), (5,3), where (a, b) means that the task b starts not earlier than at the time when the task a finishes. The Gantt chart C1 in Figure 1 represents the left-shifted schedule, which reveals that the tasks 1, 5 and 6 belong to a critical path. The tasks 2, 3 and 4 are non-critical and have some positioning freedom preserving the schedule makespan C max , between the leftmost and rightmost positions illustrated in the charts C1 and C2, respectively. This freedom is exploited in the ECETPP optimisation to assign the best cost-efficient execution time window for each task. A user may want to relax the makespan in exchange for additional positioning freedom and generating a more cost-efficient schedule. It can be obtained by setting H > C max increasing the positioning freedom of each task by H − C max , compared to the case preserving makespan of the left-shifted schedule (Figure 1 chart C3).
One can observe that "non-overlap" constraints, specific to resource conflicts in shop-scheduling problems, are not incorporated into the ECETPP. It is a consequence of the assumption that such conflicts are to be resolved by a standard advanced planning and scheduling (APS) module, before the problem data are provided to the ECETPP optimisation. This assumption is reasonable, because APS modules are well known and refined applications for general schedule optimisation; thus, it is better to optimise energy cost-related features of a schedule after it is pre-scheduled by an APS. Such a hierarchical bi-objective scheduling is stated as one of the useful options in [2] and employed in a few algorithm implementations [14,15,18].
The TOU tariff price and RES power functions presented in Figure 1 are examples, but they are also the realistic propositions, which have been used in all computational experiments. The TOU price function of an EGS (Figure 1 plot A) is taken from a price list of the Polish Energetic Group PGE Distribution and it represents a 4-tariff plan signed as B24, according to Table 1. The RES power function (Figure 1 plot B) describes a hypothetical solar source. There is no power available in the night, and the peak is about midday. A quite sunny day 1 is forecast, which will provide about 80% of the maximum power in the peak, the day 2 will be cloudy with bimodal power profile, finally, the very sunny day 3 is forecast. The function of available RES power repeats periodically every 3 days to keep the data simple and compact.  The author of this work participated before in research aiming to develop efficient algorithms for structurally complex scheduling problems [26,27]. The past experience has been useful in this work in the preparation of the ECETPP model and implementation of the scheduling algorithms, especially the one based on the TS.
The work is organized as follows. The formal description of the optimisation problem is given in Section 2. The approach based on a discrete-time MILP model is described in Section 3. In the subsequent section, the TS algorithm is introduced. Computational experiments are presented in Section 5, in which a general overview of results for 182 benchmark instances are provided, as well as some detailed properties of one selected solution are analysed. Relationships between work assumptions, results and their consequences are commented in Section 6. Conclusion and some remarks on future work are given in the last section.

Problem Formulation
The ECETPP is defined as follows:

1.
A piecewise linear function of time R(t) defines the available RES power.

2.
The variable TOU price component of the energy from a conventional EGS is described by a piecewise constant function of time C(t).

4.
For each J i ∈ J, a respective processing time p i is defined.

5.
The task J i ∈ J needs a constant electrical power w i during its execution.

6.
A directed acyclic graph G = (J, A) is given, where A ⊂ J × J. It defines precedence constraints between tasks. If (J i , J k ) ∈ A, the task J k cannot start before the task J i is completed. 7.
The planning horizon H is defined. All tasks have to complete by the time H. 8.
The objective is to determine start times s i of all tasks J i ∈ J which minimise the total cost of electrical energy:

Discrete-Time Linear Programming Model
A discrete-time model has been prepared for the ECETPP. The planning horizon is divided into a finite number of time slots. The model makes it possible to determine duration of the tasks with accuracy of the slot sizes. It limits precision of the model but supports an easy implementation based on the linear programming approach. More precisely, one has to choose between greater accuracy, i.e., smaller slots, or better efficiency of optimisation, resulting from the lower number of the slots and related decision variables.
To obtain the discrete-time model, a value of time step (slot size) δ has to be provided. It is needed to derive the following parameters of the model: The values c t and r t are obtained from the functions C(t) and R(t), respectively. The parameter w i from the ECETPP is incorporated to the discrete-time model directly.
The following decision variables are used in the model: . . , D}-index of the first slot of the task J i , ω i ∈ {1, 2, . . . , D}-index of the last slot of the task J i , ε t ∈ [0, E ]-total consumption of the EGS energy in the slot t, β i,t ∈ {0, 1}-auxiliary variable equal to 1, if the task J i is executed in the slot t, γ i,t ∈ {0, 1}-auxiliary variable used in constraints related to β i,t , where i ∈ {1, 2, . . . , n}, t ∈ {1, 2, . . . , D}.
The objective of the optimisation is to minimise the function: subject to: The constraints Equation (6) assert proper relations between the start and completion slots of tasks, according to their processing times. The component −1 in Equation (6) is the result of the discrete time representation which implies that a task with a processing time p i starts and ends in slots α i and ω i , respectively, such that ω i − α i = p i − 1. The precedence relationships between tasks are imposed by the constraints Equation (7), which assure that after a task completed in the slot t, a next one can start at the earliest in the slot t + 1. The constraints Equations (8) and (9) force assignment β i,t = 1 if the task J i is processed in the slot t, i.e., α i ≤ t ≤ ω i . The auxiliary variables γ i,t connect the constraints Equations (8) and (9). The constraints Equation (10) provide that the value of ε t is not less than the consumption of the EGS energy in the slot t. The variables ε t are used in the Equation (5) representing total cost of the energy which value is minimised.
The model contains real decision variables ε t and integer decision variables α i , ω i , in particular the binary ones β i,t , γ i,t . The objective Equation (5), as well as the constraints Equations (6)-(10) are linear expressions. Therefore, the model has a form of the standard MILP problem and many popular solvers can deal with it.

Tabu Search Algorithm for ECETPP
The TS metaheuristic has been introduced by Glover [28]. It is successfully used for combinatorial optimisation, in particular, for many hard production scheduling problems, and it is continuously developed and extended [29]. In the following sections the main elements of the TS algorithm implementation dedicated for the ECETPP are presented, namely: neighbourhood, calculation of the objective function value and some aspects of the algorithm configuration.

Neighbourhood
Let S be the set of all feasible solutions of an ECETPP instance. Let J i ∈ J be an arbitrarily chosen task of the instance. The task J i has the start time s i in the solution S = (s 1 , s 2 , . . . , s n ) ∈ S. If J i has sufficient freedom in S, this start time can be modified without changing start times of other tasks to any value from the interval The interval R i will be referred to as a range of The neighbourhood of a solution S ∈ S in the TS algorithm is represented by the function N : S → 2 S , such that: In other words, the neighbourhood consists of all feasible modifications of a source solution, such that only one task is shifted at time.

Calculation of Objective Function Value
The value of the objective function can be easily calculated for a fixed solution, directly from the formula given in Equation (1). However, in the TS algorithm, minima of this function in the neighbourhood defined by Equation (14) need to be found, and this neighbourhood includes infinite subset of solutions, according to Equation (13). The objective function is obviously non-regular, i.e., it does not change monotonically with the task start times. Therefore, the inspection of the neighbourhood is a non-trivial optimisation problem itself. A procedure to search the neighbourhood efficiently and exhaustively and to find an exact minimum has been designed and implemented. Solution subspaces induced by shifting different tasks can be searched independently. The actual difficulty is the inspection of a given subspace S [ To do it, the procedure splits the range R i into parts on which local minima can be easily found and, then, it combines the local results to obtain the range-wide minimum. The idea is illustrated in Figure 2, in which a hypothetical range R i = s L i , s R i of a task J i is presented. The base power (orange colour) is the total power of all the tasks except for J i . It is represented by a fixed piecewise constant function, as the range R i is defined under the assumption that the tasks other than J i have fixed positions. The power consumed by J i itself (blue colour) is considered separately and it reflects the selected position of J i in R i . The available power of an RES (green colour) is a piecewise linear function, while the EGS TOU price function (red colour) has a piecewise constant form, according to the ECETPP specification. All the components presented in Figure 2, namely the base power, task power, available RES power and EGS energy price, have to be considered together by the algorithm. The computation of a minimum of the objective function value inside S [R i ] is based on events. An event is a point in time at which the value of the base power or the TOU price changes, or at which an interpolation node of the RES power function is located. These three types of events are marked in Figure 2 by circles of different colours. The algorithm keeps the events in a list sorted by time and this list is used for splitting a range into smaller intervals. If a few events are synchronous, just as the two events at s 4 i in Figure 2, then the related time point is processed once. According to the example given in Figure 2, the range R i is fragmented using the following nested three-level scheme:

1.
Range. It is the full interval of the allowable positions of a given task, providing that positions of other tasks are fixed. In the considered example, s L i is the start time of the leftmost position of J i in the range R i , because at least one predecessor of J i completes at s L i . Analogously, s R i is the start time of the rightmost position of J i in the range R i , because at least one successor of J i starts at s L i + p i . In particular, a range is bounded by the time 0 or H, according to Equations (11) and (12), if the related task has no predecessor or successor, respectively.

2.
Section. Let E be the set of all events. A pair are the largest intervals (in the sense of inclusion), such that Σ L ∩ E = ∅ and Σ R ∩ E = ∅, will be called a section. A task J i is inside a section Σ if and only if s i ∈ Σ L and s i + p i ∈ Σ R . In general, a range may include many sections Σ 1 , Σ 2 , . . . , Σ u , which can be ordered in the natural way, such that sup Σ L r−1 = inf Σ L r for each r ∈ {2, 3, . . . , u}. As an example, a section Σ r is presented in Figure 2. It is left and right bounded by the events at the moments c 1 i and s 4 i , respectively.

3.
Segment. One can notice ( Figure 2) that the function which connects positions of J i and its energy cost changes character in the intervals σ R q , σ R q+1 and σ R q+2 , even though they belong to the same section. Indeed, in σ R q , all required power is provided by a free of charge renewable source, in σ R q+1 the component of renewable power decreases linearly, and in σ R q+2 all energy consumed by the task is for pay. In general, to obtain the intervals on which the cost functions are defined homogeneously, one has to split sections into smaller parts at the points at which These parts are referred to as segments. By analogy to a section, the q-th segment of the range is represented by a pair of intervals σ q = σ L q , σ R q , such that a task J i is inside this segment if and only if s i ∈ σ L q and s i + p i ∈ σ R q .
base power renewable power energy price task power The procedure which inspects the solution subspace S [R i ] of N(S ) traverses the range R i segment by segment, starting from the original position of J i in S and going maximally to the right and then maximally to the left inside this range. In Figure 2, a traverse through the segment σ q is highlighted as an example. Observe that this single-segment traverse is equivalent to the transfer of a part of the task from , while the power consumption and its cost in the interval (s 2 i , c 1 i ) remain unchanged and do not need to be recalculated. According to the property of a segment, where inside a segment, i.e., for t ∈ σ q , ∈ {L, R}. The value P i (t) represents the paid component of instantaneous power of the task J i , i.e, the component not covered by an RES. The constants C , a and b depend on the current solution S and the selected segment, but additional annotations are omitted for legibility, assuming that these constants describe the segment σ q from Figure 2. Suppose that the cost of the J i execution amounts to C 1 i , providing that J i starts at s 1 i . Then, the cost c q,i (s i ) after shifting the start time to The formula c q,i (s i ) is well defined on the closed interval because the cost change is a continuous function of the shifting distance. According to Equation (16), c q,i (s i ) is a simple quadratic function and one gets immediately that it obtains its minimum on [s The procedure of neighbourhood inspection determines the values of required parameters and finds the local minimum inside the segment σ q using Equations (16) and (17). Additionally, the value c i (s 2 i ) is calculated, which represents the cost of the processing of J i if it starts at the beginning of the segment σ q+1 . Therefore, the calculation can be recursively repeated for the next segment, and so on, up to the end of the range. Analogously, the range is traversed to the left, starting from the original task position. The only difference is the exchange of the bounds, such that s 2 i and s 1 i become the initial and final positions, respectively, in a traversed segment σ q = s 1 i , s 2 i . In this way, each range of the neighbourhood is inspected. More precisely, for practical reasons, the actual algorithm implementation checks the cost value for a finite set of candidate positions, which includes all the positions connecting segments and the positions conforming the first variant of Equation (17). Notice that the global cost minimum of the neighbourhood has to be reached at some candidate position; thus, the neighbourhood search is well-formed.

Tabu Configuration
In addition to the neighbourhood itself, a tabu attribute, tabu list and aspiration criterion are other important elements configuring a TS procedure.
A tabu attribute is the property of a solution or a move which is inserted to the tabu list. In the considered implementation, the fact has been exploited that all events are ordered on the time axis, as illustrated in Figure 2. Let e i be the event of start of the task J i and let the triple ( * e i , e i , e * i ) denote that the event * e i directly precedes and the event e * i directly follows e i . The triple ( * e i , e i , e * i ) is used as the tabu attribute added to the tabu list after execution of the move in which the start time of the task J i is changed. As a result, further moves inserting start of J i directly between * e i and e * i are forbidden as long as the triple ( * e i , e i , e * i ) remains in the tabu list. Experiments have confirmed that such a tabu attribute is effective for the optimisation problem, i.e., it prevents blocking the algorithm in a cyclic sequence of solutions. The tabu list length is controlled dynamically to adapt it to different sizes of solved instances. It has been experimentally verified that the algorithm works fine if the length is set as proportional to the neighbourhood size. More precisely, since the k-th algorithm iteration in which the length is set to L, this length is used for L consecutive iterations, and the new length L * is calculated after these iterations, according to the formula where F i is the number of feasible solutions, i.e., the neighbourhood size, in the i-th algorithm iteration. Then, L is recursively replaced by L * . At the beginning of the algorithm execution, L = 200 is fixed arbitrarily. The standard aspiration criterion is used in the implementation, which accepts a new candidate solution if it improves the current global minimum of the cost, even if it is forbidden by the tabu list.

Organization of Experiments
Classic job shop problem (JSP) benchmark instances have been used as the basic data sources. Well known benchmark sets have been employed, namely: Abz07-09 [30], Dmu01-80 [31], Ta01-80 [32], Swv01-15 [33], Yn01-04 [34]. The set includes 182 instances having from 200 to 2000 tasks. The abstract time unit of JSP instances has been mapped to the interval of 10 min. As a result, the benchmarks instances represent schedules having the makespan approximately between 5 and 50 days. The ECETPP requires predefined precedence constraints between tasks, hence, solved instances of the JSP has been used, imported from the archive of Shylo [35]. The complete data representing an ECETPP benchmark instance have been composed of the following elements: • Set of tasks, processing times of the tasks and precedence constraints between them-acquired from the related solved JSP instance. • Power consumption of each task-randomly determined using the probability distribution U(600, 12, 000) [Watts]. • EGS TOU price function-determined by the PGE Distribution tariffs B24 (Table 1, Figure 1 plot A), the same for all instances. • Time profile of renewable power-based on the function presented in Figure 1 plot B, scaled according to the average power consumption.
For all the benchmark instances, eight experiments have been conducted, one for each configuration from the following set of 3-tuples The variants involving H equal to 1.0 or 1.1 of C max make it possible to compare the cases with minimum or increased freedom of tasks positioning. In addition to the complete form of the ECETPP, the simplified variant without RES power has been optimised as well for results comparison.
The IBM CPLEX has been used for MILP optimisation. A run time for each instance has been limited to 1 h. The total time limit for the TS has been set to 1 h per instance, as well, but the algorithm was started many times within this limit providing different results, due to stochastic elements of its implementation. Each execution was terminated according to the stop condition stop = n nimp > n limp ∧ (t > 60 s) (18) where: n limp -the number of TS main loop iterations from the beginning of execution to the last objective improvement, n nimp -the number of TS main loop iterations from the last objective improvement to the current iteration, t-execution time.

Overview of Results
Selected numerical results are presented in Tables 2-4, which relate to the MILP approach with H = C max , the MILP with H = 1.1C max and the TS, respectively. Each of the tables has two parts-part A for the case of available renewable energy and part B for the variant without it. Each part, in turn, contains the 12 best results first, followed by the 12 worst results, sorted by the value of the objective function. Other results are omitted in the tables for brevity, but they are appended to Supplementary Materials.
The column T [days] of Tables 2 and 3 where c orig represents the original cost of the consumed electrical energy, based on the input left-shifted schedule, and c opt represents the energy cost of the schedule optimised using the ECETPP. The subscripts A, B, C, D specify that the ratio is calculated only for the cost restricted to: (a) load valley, (b) other hours, (c) morning peak, (d) afternoon peak, respectively, as defined in Table 1.
The coefficients ∆ d and ∆ l are calculated for entire schedules according to Equation (19), but the first one is based on the discretised function of renewable power in its actual form used in the MILP model, while the second one is based on the original, more precise and more realistic piecewise linear interpolation of this function, used also in the TS. In the case without renewable energy, ∆ d = ∆ l = ∆ because the TOU tariff function (see Figure 1 plot A) does not change after discretisation with the sample time of 10 min. The coefficient R applies only to the case with renewable energy and it is defined according to the formula R = r opt − r orig r orig 100% (20) where r orig , r opt represent total consumption of the RES energy in the original and the optimised schedules, respectively. From Tables 2 and 3 (column T), one can observe that the best results are obtained for relatively short schedules and the worst ones are generated for the longest schedules, especially in the case H = C max . Satisfaction of the inequality gap ≤ 0.01% was considered as a practical proof of solution optimality, and further optimisation was stopped. Almost all from the 12 best results are optimal in the case H = C max (Table 2), except for Ta35 (both in the parts A and B). The cost is reduced by approximately 15-20% if the renewable energy is available and by approximately 9-12% in the opposite case. Even for the worst cases, the cost is reduced, but not more than by 4%. The values of R reveal that the optimisation process has effectively used RES energy, increasing its consumption by approximately 7-17% in the best cases. As a result of the optimisation, the cost has been relatively reduced for the hours of peaks, when the power is more expensive (∆ C , ∆ D ), and it has been relatively increased for the hours of load valley (∆ A ). The results for the case H = C max (Table 2) have been described so far. In the case H = 1.1 C max , the results are typically worse, especially if renewable energy is involved. Some solutions are even worse than the original reference schedules (see Table 3 part A, ∆ d > 0). Moreover, as a consequence of unfavourable computational complexity of the algorithm, no results have been obtained within one hour of optimisation for the instances Ta71 . . . Ta80 (H = 1.1 C max , no renewable energy). [%]  [%] The results related to the TS are given in Table 4. As the TS was executed many times for each instance, the basic statistical parameters of the objective function value are given, namely: minimum (best) value min(∆), mean value µ(∆), standard deviation σ(∆) and coefficient of variation where ∆ is calculated according to (19) for each solution. The value ∆ relates to the full range of a schedule, and the division into ∆ A , . . . , ∆ D is omitted in this table to be concise. The instances are sorted by min(∆). Table 4 additionally contains the number of the algorithm repetitions for a given instance (column n) and the total number of iterations of the algorithm main loop, summed over all the repetitions (column iter). In the case of H = 1.1 C max , in which the freedom of task positioning is greater, the number of repetitions is typically about two times less for a given instance. Although the TS slows down for H = 1.1 C max , it keeps good performance, much better than the MILP for the same case. In particular, the best results for the variant with renewable energy are much better (approximately 27-41% vs. 7-20% of the cost reduction for the TS and MILP, respectively), and the TS reduces the cost for all the instances, at least by a few percent. On the other hand, the MILP has provided a few better solutions for Part B (approximately 19-22% of the cost reduction), than obtained by the TS. Inst. The collective performance comparison is presented in Figures 3 and 4, for the case with available and unavailable renewable energy, respectively. The results for H = C max are quite balanced, they do not differ too much and the MILP is typically only a bit better for some instances, while the TS for others. However, the MILP is evidently worse for some groups of instances with the numbers about 40 and 170 in Figure 3. For the case without renewable energy and H = 1.1 C max , the MILP is almost always significantly worse. For the case with renewable energy available and H = 1.1 C max , the MILP usually does not optimise the cost at all, while the TS performs it effectively.

Selected Instance
One particular solution will be presented in more details, namely the solution for the instance Abz09 with renewable energy and H = C max , for which the best EGS energy cost reduction has been obtained (by the MILP). In Figure 5, the time functions of consumed electrical power are presented for the reference and optimised schedules. The plot shows very explicitly how the schedule has been improved. As a result of the tasks positioning, the power consumption increased on average in the intervals in which renewable energy is available (green areas) and, simultaneously, it decreased in the intervals in which the external paid power (EP) is very expensive (yellow strips). In other intervals, in which the EP is relatively cheap, the power consumption statistically increased as well, e.g., in the hours about: 25-30, 72-78, 95-102.  In Figure 6, Gantt charts of the Abz09 schedules are presented, namely: (1) the reference schedule based on the left-shifted makespan-optimised best known solution, (2) the schedule considered in this section, i.e., the best cost-efficient solution obtained for H = C max by the MILP, (3) the best cost-efficient schedule generated for H = 1.1 C max (obtained by the TS). Both the schedules (2) and (3) are optimised with an available RES. The schedule (3) is taken into account only in this figure, to make the comparison more complete. The colour map indicates how much a task is delayed (in hours) respective to its reference start time defined by the schedule (1). As one can expect, many tasks of the schedule (2) are not delayed at all, because they belong to critical paths. In particular, the critical paths end on the machines 3 and 11, hence, the final tasks on these machines are not shifted. One can observe that only a few tasks are delayed more than 6 h. In an alternative way, the same results are presented in Figure 7a,b. According to Figure 7a, the cost saving cumulates gradually over the schedule. The histogram in Figure 7b, which counts the delays with a resolution of 15 min, confirms that almost all task delays are not greater than 2 h and that there are only 4 tasks delayed between 6 and 10 h. The non-delayed tasks (121 of 300) are omitted on the histogram. In Figure 8, a hourly histogram of energy consumption is presented. It is calculated for the first four complete days of the Abz09 schedule, the incomplete fifth day is not included, so as not to distort the data balance. One can notice that the power consumption in the optimised schedule is remarkably reduced in the afternoon peak (19)(20)(21)(22). There is also noticeable reduction at 7 and 8, when the morning peak starts, but the level of the available RES power is low. The high power consumption is redirected, instead, to the hours inside the load valley (mainly 2-4) and to the time intervals on which the available RES power is highest (mainly 10-14).

Discussion
The experimental results indicate a substantial impact of the relation between H and C max on performance of the MILP approach. In the case of H = C max and available RES power, there are 94 (54%) solutions with proven optimality and the lower-upper bound gap is not greater than 21% for every instance, except for the largest ones (Ta70-80), for which the gap reaches up to 37%. If H = C max and RES power is unavailable, the MILP provides 88 (48%) optimal solutions and the gap does not exceed 16% in any case. In general, the results are very promising, and this approach under the considered configuration can be directly recommended for practical usage. In the case of H = 1.1 C max , no solution is proven to be optimal and the gap values are large, around 100%; moreover, the result are remarkably worse than those obtained by the TS or even worse than initial unprocessed schedules. Therefore, the MILP approach is useless in this case.
It is difficult to objectively and perceptively evaluate performance of the TS algorithm for H = 1.1 C max , as there is no reference optimal solutions. The TS is efficient at least in the relative sense, such that it uses the additional freedom of task positioning to find start times of the tasks which minimise energy cost better than for H = C max . It is clearly shown in Figures 3 and 4. Therefore, if a user wants to obtain less energy cost at the expense of longer makespan, the TS will produce the expected outcome. In any case, the TS seems to be the only practically useful approach for the ECETPP optimisation if H is greater than C max by about 10% or more.
In practical applications, a rolling scheduling horizon will be probably used. The trade-off between scheduling speed and quality of schedules can be additionally controlled by the length of this horizon, independently on the relation between H and C max . A similar trade-off is controlled by the slot size δ, but only in the case of the MILP approach. If δ is too large, the scheduling result is imprecise, but there are a less number of slots and a solver can provide a result faster. The value δ = 10 min used in this work has been chosen properly, because the discretisation error, represented by the differences ∆ d − ∆ l in Tables 2 and 3 (part A), is practically insignificant. The value of δ needs to be adjusted according to a form of target problem instances.
In this study, the same arbitrarily chosen but realistic functions of TOU tariffs and RES power availability have been used for all benchmark instances. One can generate and verify results for other functions in future research; however, a precise form of these functions has no major impact on the problem instance size and, as a consequence, on the optimisation performance. It is reasonable to expect that only the values of the cost reduction will change, according to new TOU-and RES-related functions, but the performance-related properties of the MILP and the TS will be similar, and these properties are crucial for real-world scheduling systems.

Conclusions
The research confirms that it is possible to develop a comprehensive model and related algorithms for production scheduling with TOU price tariffs and RES power minimising energy cost. The approach based on the MILP and configured with H = C max is ready for use in manufacturing systems conforming the ECETPP specification, as it assures efficient optimisation for real-world size problem instances of order of 1000 tasks.
Fixing H ≥ 1.1 C max changes the situation remarkably. The MILP is no longer the efficient and recommended approach, and only the TS algorithm is able to make use of the increase of H over C max . The TS can also be used in real-world applications, but one should keep in mind that its implementation is harder to modify and extend, and that the TS cannot prove optimality and typically provides suboptimal solutions.
From a theoretical point of view, the obtained results reveal that production scheduling regarding TOU tariffs and RES power is, in general, computationally hard and sensitive on values of some parameters, even for a relatively simple problem model. Therefore, further research are needed in this field. The results of this work provide as well a kind of warning that proposed scheduling approaches should be tested on numerous real-world size benchmark instances and a small simple demonstrative example, often used by authors, may be not sufficient.
Optimal solutions for H = 1.1 C max are not known and it remains a challenge for future work to construct algorithms obtaining the optimal or at least better results. Another research direction is to extend the ECETPP formulation by useful features of real systems and verify how it affects the optimisation. An onsite energy storage, a maintenance cost of an RES or a more complex task model, e.g., with variable processing speed, are the possible extensions. The results obtained in this study, especially the optimal ones, could be a good reference base for comparison after the model extension.
Supplementary Materials: The materials available online at http://www.mdpi.com/1996-1073/13/19/5034/s1 include benchmark instances used for computational experiments, the full set of the experiment results, and the source codes of the implemented programs.
Funding: This research received no external funding.

Conflicts of Interest:
The author declares no conflict of interest.

Abbreviations
The following abbreviations are used in this article:

APS
advanced planning and scheduling EGS electrical grid system ECETPP energy cost-efficient task positioning problem JSP job shop problem MILP mixed-integer linear programming RES renewable energy source TOU time-of-use TS tabu search