Strojniški vestnik - Journal of Mechanical Engineering 62(2016)11, 685-694 © 2016 Journal of Mechanical Engineering. All rights reserved.
D0l:10.5545/sv-jme.2016.3994	Original Scientific Paper
Received for review: 2016-08-26 Received revised form: 2016-09-09 Accepted for publication: 2016-09-21
Comprehensive Energy Resource Management for Essential Reduction of the Total Cost
Liljana Ferbar Tratar* University of Ljubljana, Faculty of Economics, Slovenia
Comprehensive energy resource management presents a kind of global optimisation, which is much more cost-effective than separate treatment of energy-demand forecasting and energy-inventory data. In this paper, a two-stage supply chain is presented where the distributor of derivate energy products (e.g., hot water, steam, electricity) shares information with its supplier of primary energy sources (e.g., coal, crude oil, gas). Different methods of energy-demand forecasting can be used, although exponential smoothing methods are most often used in practice because they are simple, fast and inexpensive. In this study we analysed a modified Holt-Winters (HW) method and we describe a method for simultaneous optimisation of a forecasting method and a stock control policy. We compared the proposed joint optimisation of total cost with the optimisation, based solely on forecasting data. The total cost presents the joint cost of the distributor and supplier. The data consisted of 756 quarterly real series from M3-Competition and 300 simulated demand patterns. We have shown that the total cost can be reduced dramatically if we use joint optimisation instead of a separate treatment of a forecasting method and an inventory model. We obtained the best result with the modified HW method; the essential reduction of total cost was reached in the case of simulated data as well as in the case of real data.
Keywords: energy resource management, supply chain, forecasting, inventory, M3-competition, total cost Highlights
•	Comprehensive energy resource management is presented.
•	Two-stage centralised supply chain of distributor and supplier is observed.
•	The joint optimisation comprises energy-demand forecasting methods and energy-inventory data.
•	The result shows an essential reduction of the total cost as an objective function.
0 INTRODUCTION
Sustainable energy supply, reducing dependence on energy resources, efficient energy production, climate-friendly and less energy-consuming production are some of the most important challenges facing the European Union in present times [1]. The European commission encourages energy companies to follow suggested climate action, which comprises cost-efficient and less environmental damaging power generation [2] and [3]. Companies have to look constantly for new strategies and tools to improve processes, decrease cost and increase productivity and efficiency [4] and [5]. The appropriate demand forecasting methods, inventory data and sharing of information constitute comprehensive energy resource management, which may lead to essential reduction of the total cost, better planning performance with energy resources and much more efficient production.
Customer demand is probably one of the most important factors to be predicted through a forecasting system that is especially problematic when we are faced with a noisy demand. Selecting a forecasting method from the many available is a very complex task (see, for example, [6] and [7]). The main approach towards the selection and optimisation of
alternative methods relates to the minimisation of forecast error measures. Most of the stock control studies consider demand data as an input to the model without explicitly considering that they are the results of a demand forecasting system. Even though this weakness has been highlighted in the academic literature, little empirical work has been conducted to develop an understanding of the interaction between forecasting and stock control, except in [8] and [9], where the problem of the local optimisation of forecasting methods was elucidated. In general, separate evaluation of the forecasting method and the stock control policy may easily lead to poorer overall performances [10].
We present an example of a centralized supply chain with an order-up-to inventory policy and evaluate different forecasting methods together with the inventory model. In this joint model we optimise total inventory costs where holding and stock-out costs (for different penalties) are included. We demonstrate that smoothing and initial parameters of forecasting methods can be determined to minimise the total costs. The joint model can be implemented easily, even by using an Excel spreadsheet, and the proposed solution is easy-to-use for managers.
*Corr. Author's Address: University of Ljubljana, Faculty of Economics, Kardeljeva ploščad 17, 1000 Ljubljana, Slovenia, liljana.ferbar.tratar@ef.uni-lj.si
685
Strojniski vestnik - Journal of Mechanical Engineering 62(2016)11, 685-694
Exponential smoothing methods are a class of methods that produce forecasts with simple formulae, taking into account the trend and seasonal effects of the data [11]. Distinguished by their simplicity, their forecasts are comparable to those of more complex statistical time series models [12]. The main advantages of HW methods, a special subset of exponential methods, are low expenses, fast calculations and simplicity.
The aim of this paper is to present a modified HW method, which is computationally stable and can handle both additive and multiplicative seasonality, even when a time series contains a nonlinear trend and a large noise component. We also expose the problem of the local optimisation of forecasting methods. We show that in the case of a centralised supply chain the calculated forecasts of demand, determined by minimising the mean square error (MSE), are not optimal. We therefore propose a method for the simultaneous optimisation of demand forecasting and a stock control policy, and demonstrate that initial and smoothing parameters in the forecasting methods can be determined to minimise the total costs of the supply chain.
The idea of the joint optimization could be used for different objective functions such as the total energy efficiency and environmental impact concerns [13].
The remainder of the paper is organised as follows. In Section 1, the methodology and data are described. In Section 2 we present forecasting methods; the additive, multiplicative and modified HW method. In Section 3 we describe our model of the supply chain and calculate average costs for all different forecasts considered as an input to the model. In Section 4, we describe the proposed joint optimisation of total costs and present the results, which allows us to compare different forecasting methods and replenishment policy combinations. Finally, after the conclusions of our paper some further steps of research are suggested.
1 DATA AND METHODOLOGY
We used real seasonal time series from the M3-Competition and conducted a simulation study to evaluate the performance of the modified HW method. The analyses were carried out in the program R [14]. The function sbplx from the nonlinear optimization package "nloptr" [15] was used to estimate the smoothing parameters. The starting values in the minimization step were set to a0=^0=Yo = 0 5 and the maximum number of iterations was set to 25,000.
1.1 Real Time Series from the M3-Competition
The Makridakis Competitions, known in the literature as the M-Competitions, are empirical studies that have compared the performance of large number of major time series methods using recognized experts who provide forecasts for their method of expertise [12]. The first M-Competition (1982) used 1001 time series and 15 forecasting methods. The second M2-Competition (1993) used only 29 time series. The third M3-Competition (2000) was intended to both replicate and extend the features of the first two competitions. A total of 3003 time series were used. The data sets used refer mainly to business and economic time series, although the conclusions can be relevant to other disciplines as well. The original time series data can be found in the R package Mcomp [16].
In our study, we have analyzed quarterly series. They refer to five different disciplines, as shown in Table 1.
First we used the "ets" function from the R package forecast [17] and [18] to classify the series by the form of their trend, seasonality and noise. Table 1 also shows this classification, where 'A' stands for 'additive', 'M' for 'multiplicative', and 'N' for 'none'.
We applied different forecasting methods to each of the series independently of its discipline and ets classification.
Table 1. Classification of quarterly time series from the M3-Competition
Discipline		Noise	Trend	Season	
Micro	204	A	N	N	67
Industry	83	A	N	A	40
Macro	336	A	A	N	113
Finance	76	A	A	A	64
Demographic	57	M	N	N	90
Total	756	M	N	A	31
		M	N	M	37
		M	A	N	87
		M	A	A	49
		M	A	M	42
		M	M	N	76
		M	M	M	60
				Total	756
1.2 Simulated Time Series
We have simulated 300 time series corresponding to six different demand patterns using the following formula (based on [19]):
686
Ferbar Tratar, L.
Strojniski vestnik - Journal of Mechanical Engineering 62(2016)11, 685-694
Dt = base + slope x t2 x A + season x snd ^ x t x B + +noise x snormal() x t x C,	(1)
where Dt is the demand in time t, base is the average demand, slope captures trend in data (series with no trend are obtained by setting slope = 0), season is a seasonal factor, noise is the coefficient of demand variation, and snormal() is a standard normal random number generator. sind are seasonal indices, which are obtained for t = 1, 2, ..., 12 with snormal(). By varying the quantities A, B and C we obtain time series differing in the form of their trend, seasonality and noise. The possible demand patterns are the following. If we make A equal to 1, the equation provides a nonlinear (quadratic) trend; otherwise, if we make A equal to 1/t, we get a pattern with a linear trend. If we make B equal to 1/t, we get a demand pattern with additive seasonality; otherwise, if B is equal to 1, we get a pattern with multiplicative seasonality. The expression in the last term of the equation provides additive or multiplicative noise if we make C equal to 1/t or 1, respectively. Using the different combinations of parameters (base, slope, season, noise) we can generate the demand patterns.
Table 2. Descriptions and marks for different demand patterns
Description	Mark		noise (step)	slope	season
add. noise,	AAM	5 50	200-650 (50)	5	50
nonlinear trend,	AAM	5 100	200-650 (50)	5	100
mult. seasonality	AAM	8 100	200-650 (50)	8	100
mult. noise,	MAM	5 50	10-55 (5)	5	50
nonlinear trend,	MAM	5 100	10-55 (5)	5	100
mult. seasonality	MAM	_20_50	10-55 (5)	20	50
In this study we analysed in detail only demand patterns with additive and multiplicative noise, additive nonlinear trend and multiplicative seasonality. Base was selected to ensure that the average demand was approximately 2000 units; the parameters for six demand patterns are shown in Table 2. We used notation noise (step) to point out that for each selected parameter slope and season we chose 10 different noises (e. g., for the AAM_5_50, noise = 200, 250, ..., 650). The AAM produced a demand with additive noise, additive nonlinear trend and multiplicative seasonality (quantity C is equal to 1/t, A and B are equal to 1); the MAM produced a demand with multiplicative noise, nonlinear trend and multiplicative seasonality (quantities A, B and C are equal to 1). Each combination of the simulation study was replicated five times (for different snormalQ). Thus, the total number of simulation runs for the
experiment is equal to 6 (patterns) x 10 (different noises) x 5 (simulations) = 300.
1.3 Symmetric Relative Efficiency Measure
The efficiency of the modified HW (MoHW) method was measured in terms of the MSE of the in-sample one-step-ahead forecasts and compared to that of additive (AHW) and multiplicative (MHW) methods. Because the first two complete seasons were used to initialize the methods, these observations were excluded from the reported MSE:
MSE = 7-rZC+1 - Y+
T 2S t=2s
(2)
where Yt is the actual value at the time point t, Ft is the forecasted value at the time point t, s is the length of seasonality and T is the length of the observed time series. To compare the MoHW method with the other method, we first find their mean square errors, MSEMoHW and MSEmethod as defined above. We define the symmetric relative efficiency measure as
SREM
1 -
MoHW ' method
MSEM
MSE
MSE
MSEM
--1;
MSEMoHW < MSEmethod
MSEmohw ^ MSE
method
(3)
We prefer the SREM to the standard relative efficiency measure the REM (percentage increase or decrease of the MSE),
REM,
MSEMoHW MSEmethod
(4)
because it treats the methods symmetrically: e.g., when MSEMoHW = 10 and MSEAHW=20, the REM is -50 %, while when the MSEMoHW = 20 and the MSEAHW = 10, it is 100 %%. On the other hand the SREM is 50 %% in the first case and -50 %% in the second case, indicating that 'on average', none of the methods is preferable. The value of the SREM is bounded by the interval [-1, 1], which mitigates the possibility of an individual time series to substantially outweigh other series in the group. The interpretation does not depend on the number of series in the group, so the SREM can easily be applied to the M3-Competition data where different disciplines (or types) have different numbers of time series.
Also, we use the definition of the SREM to compare the MoHW method with other methods
Comprehensive Energy Resource Management for Essential Reduction of the Total Cost
687
Strojniski vestnik - Journal of Mechanical Engineering 62(2016)11, 685-694
regarding average costs AC (in these cases, the SREM1 measures the percentage increase or decrease of the average costs): SREMU™ =
1 -
method
ACMoH
AC
AC„
AC
--1;
ACMoHW < ACmethod
AC > AC
ACMoHW ACmethod
(5)
2 FORECASTING METHODS
The Holt-Winters method estimates three smoothing parameters associated with level, trend and seasonal factors. The seasonal variation can be of either an additive or multiplicative form. The multiplicative version is used more widely and on average works better than the additive one ([20]; of course, if a data series contains some values equal to zero, the multiplicative HW method cannot be used).
In the multiplicative seasonal form of the HW method (MHW) fundamental equations for level Lt, trend bt, seasonal factors St and forecast Ft are:
L =a(Y /^) + (l-a)(( + bM),	(6)
b, = P{L, -L,_i) + (l -p)b,_l,	(7)
= y(Yt /Lt) + (l -y)S,_,	(8)
F+m =(L + btm)St_s+m,	(9)
where a, p, y are smoothing parameters (which must lie in the interval [0, 1]), m is the number of the forecast ahead, 5 is the length of seasonality (e.g., the number of months or quarters in a year) and Yt is the observed data at the time point t. We have given values for 5, m and Yt, t = 1, 2, ..., T. To initialize the level, we set L5 = (Yj + Y2 + ... + Y5) / 5; to initialize the trend, we use b5 = (Y,+i - Yj +Y+- Y2 + ... + Y25 - Y5) / 52; and for the initial seasonal indices we calculate Sp=Yp / L5, p = 1, 2, ..., 5.
The additive seasonal form of the HW method (AHW) works with the following equations:
L = a(( -S,_s) + (l-a)^ + b-),	(10)
b, = P{L, -L,_i) + (l -P)b,_l,	(11)
=y(Yl - Lt ) + (l -y)S,_,	(12)
F+m = L, + btm + S,_s+m.	(13)
Eq. (11) is identical to (7). The only differences in the other equations are that the seasonal indices
are now added and subtracted instead of relying on products and ratios. The initial values for level and trend are identical to those for the multiplicative method. To initialize the seasonal indices we calculate Sp = Yp-L5, p = 1, 2, ..., 5..
The only difference between the AHW and the MoHW method is in the equation for seasonal factors:
=y{Y1 - L, ) + (l-ay)S, _
and forecast:
F = L
t+m t
+ btm + aSt _
(14)
(15)
The other equations for MoHW conform to the AHW format. Thus, when we minimize the MSE with respect to the smoothing parameters, the new effect is to smooth the seasonal factors by changing them less. The initial values for the level, trend and seasonal components are the same as in the case of the AHW method.
For each of the series we also used the "ets" function [17] and [18] to obtain the MSE, where we set opt.crit = 'mse', ic = 'aic', bounds = 'usual', so that the MSE was minimized to estimate the parameters of each model; the aic was used to select the best model, and the standard parameter restrictions were applied. We use the notation ETS method. It is a state space model that includes some transition equations that describe how the unobserved components or states (level, trend, seasonal) change over time. The classical decomposition method splits a time series into a trend and a seasonal component and projects them into the forecast horizon [21].
3 THE SUPPLY CHAIN COST MODEL
Consider a simple two-stage supply chain consisting of one distributor of derivate energy products (e.g., hot water, steam, electricity) and one supplier of primary energy sources (e.g., coal, crude oil, gas). The distributor holds an inventory in order to meet an external demand and places inventory replenishment orders to the supplier. At the time t, the last known value of the external demand is Dt-1. The distributor places an order Qt to the supplier. We assume that the order placed one period ago is received (delivery lead time is one period). After the order placement, the external demand Dt is observed and filled. At the end of each period, the inventory costs are evaluated. The unsatisfied demand is backlogged and causes penalty costs for the distributor. The supplier is able to supply any requested quantity.
686
Ferbar Tratar, L.
Strojniski vestnik - Journal of Mechanical Engineering 62(2016)11, 685-694
We assume that the distributor follows an order-up inventory policy. An order Qt placed by the distributor to the supplier can be expressed as Qt=Ft+1 - FSt, where Ft+1 is the forecasted demand for a time point t+1 and FSt is the final stock for the time point t (if FSt > 0 the distributor has on-hand inventory, if FSt < 0 the unsatisfied demand occurs). When it is Qt < 0, an order is not placed. The final stock is calculated as FSt=ISt- Dt, where the initial stock is obtained as ISt = Qt-1 - FSt-1.
As the supplier has information about the external demand (the centralized supply chain), it places the order, which is equal to the forecasted demand (less FSt , if FSt > 0). The missing amount of products supplied from the marketplace (assuming that a perfect substitute for the product exists) causes penalty costs for the supplier.
The costs of the supply chain are the sum of the holding costs and the penalty costs for all links in the supply chain. We assume the penalty costs to be higher than the holding costs, which is expressed by introducing a weight, penalty, that is greater than 1. In our analysis, for all calculations of total costs (average costs and minimised average costs) we assume that the penalty is equal to 3 or 5.
In other words, using the common notation X+ = max(X, 0), the supply chain costs Ct at the time
point t are expressed as (l is an individual link in the supply chain [in our case, l = 1, 2], n is the total number of links in the supply chain [n = 2]):
C = jjC =
1=1
= jj -D )+ + penalty ( - IS,)). (16)
Because the first two seasons were used to initialize the methods, the average costs (AC) are calculated as:
1 r
AC = -
T - 2s =
Z c.
(17)
3.1 Real Time Series from the M3-Competiton
For each observed method and series, the symmetric relative efficiency measures (the SREM and the SREM1) of the MoHW with respect to the AHW, MHW or ETS were computed. Also, the portion of time series for which the MoHW method outperforms each other method was recorded.
Table 3 shows averages of the SREM for quarterly time series. We can observe that with the MoHW method the MSE can be reduced on average by more
Table 3. Averages of the SREM and the SREM1 for time series from the M3-Competition; the series are grouped by disciplines or types
MSE ^ COST
SREM (MSE) [%]
SREM1 (AC) [%]
penalty = 3
penalty = 5
MOHW/ AHW
MOHW/ MHW
MOHW/ ETS
MOHW/ AHW
MOHW/ MHW
MOHW/ ETS
MOHW/ AHW
MOHW/ MHW
MOHW/ ETS
MICRO
MACRO
INDUSTRY
FINANCE
DEMOGRAPHIC
7.30
3.73
13.29
2.76
1.33
9.11
2.64
3.56
3.57
10.43
1.72
1.80
.47
1.77
1.94
-2.92
6.82
0.93
-0.34
3.99
1.07
0.42
0.24
10.38
0.57
0.74
3.74
0.74
-3.74
-3.24
1.07
-4.09
-3.37
-0.19
-4.55
1.30
-0.22
0.89
-3.86
9.69
.81
3.84
2.91
-0.72
ANN
ANA
AAN
AAA
MNN
MNA
MNM
MAN
MAA
MAM
MMN
MMM
0.21
1.27
14.65
-0.02
-0.33
10.54
-0.27
1.55
5.33
9.30
0.45
2.87
11.83
0.40
-1.96
-0.27
6.67
-2.07
-1.05
5.77
-2.39
10.63
14.39
20.63
3.82
4.86
11.47
3.99
1.33
-0.52
16.11
0.97
0.21
9.31
0.91
-0.56
-0.14
16.89
-0.36
0.23
6.91
-0.31
3.68
-4.83
.84
0.63
-1.21
11.17
0.58
-2.00
-2.27
4.18
-2.56
-2.04
1.57
-2.79
3.56
7.83
4.62
0.10
4.84
3.86
-0.34
14.07
1.60
5.84
7.32
0.99
3.80
7.52
-1.17
-2.91
4.49
0.02
-0.18
3.51
0.38
23.22
9.33
12.10
13.52
5.95
.74
14.27
-0.54
2.85
-1.32
4.96
0.20
0.26
-1.07
-2.31
4.59
1.31
0.19
6.66
10.99
13.18
6.23
11.19
9.27
5.73
12.18
1.17
3.90
3.72
3.74
9.24
Average
3.53
2.05
10.09
1.36
0.94
7.02
1.36
0.96
7.19
Portion [%]
67.59
50.26
78.44
55.56
50.53
72.22
54.76
50.79
70.63
Comprehensive Energy Resource Management for Essential Reduction of the Total Cost
687
Strojniski vestnik - Journal of Mechanical Engineering 62(2016)11, 685-694
than 3 % (2 %) in comparison with the AHW (MHW) method. Also, the MoHW method outperforms the ETS, on average by more than 10 %.
For all included disciplines the MoHW method performs better than the ETS method, indicating the universality of the MoHW method regarding the ETS which tries to select the most appropriate method for forecasting.
The MoHW method does not outperform the classical AHW and MHW methods for the classes with no seasonal component (xxN), except for the ANN type, but the MoHW substantially outperforms them for classes with additive/multiplicative trend and additive/multiplicative seasonality, irrespective of the noise. Surprisingly, the fit of the MoHW method is better even in the AA and AM classes, where the AHW and MHW methods are the theoretically correct methods.
Since demand data is usually considered as input to the model in stock control studies, the average costs (for the period t = T- s) for forecasts obtained with different forecasting methods were calculated. Table 3 also shows the averages of the SREM1 (the percentage of improvement of the average costs) of the MoHW with respect to the AHW, MHW and ETS.
Almost the same as we observe for the SREM holds for the SREM1. If the MoHW outperforms classical methods regarding the MSE, the MoHW outperforms them regarding the average costs as well as in this case the costs are calculated for forecasts considered as an input to the stock control model.
3.2 Simulated Time Series
We analyzed 300 simulated time series in the same way that we analyzed the real time series from the
M3-Competition. For each method and series, the symmetric relative efficiency measures (the SREM and the SREM1) of the MoHW with respect to the AHW, MHW or ETS were computed. From Table 4 we can observe that for the AAM patterns the MoHW method reduces the MSE on average by more than 71 % (55 %, 95 %) in comparison with the the AHW (the MHW, the ETS) method. For the MAM patterns the MoHW reduces the MSE on average by more than 35 % (9 %, 84 %) in comparison with the AHW (the MHW, the ETS) method.
From Table 4 we can also observe that for the AAM patterns at penalty = 3 the averages of the SREM1 are more than 53 %, 40 % and 80 % with respect to the AHW, MHW and ETS, respectively, and that for penalty = 5 the values of the SREM1 slightly increase. This holds also for the MAM patterns (averages of the SREM1 are now more than 27 %, 9 % and 61 % for penalty = 3). The MoHW method outperforms the AHW and ETS in all 150 cases of the MAM patterns and the MHW in 70 % of cases.
Applying all smoothing methods under consideration to the simulated data confirms the good performance of the MoHW method and yields another important insight. For the AAM patterns we can observe that the improvement of the MoHW with respect to other methods increases as seasonality increases and decreases as the slope increases. For the MAM patterns the improvement of the MoHW with respect to the AHW and ETS increases and with respect to the MHW decreases as seasonality increases and increases with respect to all methods as the slope increases.
Table 4. Averages of the SREM and the SREM1 for the simulated time series
		SREM (MSE) [%]					SREM1 (AC) [%]			
MSE ^ COST						penalty = 3			penalty = 5	
	MOHW/	MOHW/	MOHW/		MOHW/	MOHW/	MOHW/	MOHW/	MOHW/	MOHW/
	AHW	MHW	ETS		AHW	MHW	ETS	AHW	MHW	ETS
AAM 5 50	62.17	37.58	93.09		47.15	30.40	74.94	49.28	32.73	74.96
AAM 5 100	81.29	66.70	97.83		59.21	48.47	84.67	60.48	49.54	84.42
AAM 8 100	70.20	62.79	95.42		55.50	43.86	81.56	62.83	50.24	84.51
Average	71.22	55.69	95.45		53.95	40.91	80.39	57.53	44.17	81.30
Portion [%]	100.00	100.00	100.00		100.00	100.00	100.00	100.00	95.56	100.00
MAM 5 50	25.76	9.32	74.98		22.64	12.96	50.82	24.58	14.98	51.62
MAM 5 100	29.15	1.49	88.54		22.66	-1.63	65.96	25.29	-3.04	66.30
MAM 20 50	51.58	17.91	88.95		37.93	18.56	66.53	39.76	20.11	66.68
Average	35.50	9.58	84.16		27.74	9.96	61.10	29.87	10.68	61.53
Portion [%]	94.44	61.11	100.00		100.00	70.00	100.00	100.00	70.00	100.00
686
Ferbar Tratar, L.
Strojniski vestnik - Journal of Mechanical Engineering 62(2016)11, 685-694
Table 5. Averages of the SREM1 obtained with joint optimisation for quarterly time series from the M3-Competition
SREM1 (AC) [%]
	JOINT		penalty = 3			penalty = 5	
		JMOHW / JAHW	JMOHW / JMHW	JMOHW / ETS	JMOHW / JAHW	JMOHW / JMHW	JMOHW / ETS
	MICRO	6.25	4.43	33.28	9.49	6.31	46.89
e in	MACRO	9.76	9.42	30.51	14.50	14.23	44.13
ipl 'o	INDUSTRY	4.83	2.71	26.95	7.22	4.79	41.85
D	FINANCE	2.44	2.55	26.38	3.42	3.29	38.56
	DEMOGRAPHIC	2.57	3.42	18.70	6.59	7.85	30.39
	ANN	0.88	2.96	32.01	3.52	4.26	45.15
	ANA	2.99	4.66	35.09	3.36	3.59	49.26
	AAN	5.10	4.85	25.11	9.51	8.07	36.75
	AAA	12.24	15.83	36.42	15.26	19.72	49.52
	MNN	3.64	1.42	33.73	3.93	3.92	47.56
e p	MNA	1.70	0.97	30.15	0.80	2.31	44.37
ypT	MNM	5.63	1.30	37.39	9.21	4.66	51.21
	MAN	6.15	5.41	24.03	12.53	11.53	38.43
	MAA	5.98	10.21	27.93	10.53	13.23	42.17
	MAM	12.60	3.93	26.61	15.86	6.93	42.30
	MMN	7.86	8.23	23.09	13.12	12.51	34.30
	MMM	21.84	13.85	32.44	27.22	17.81	47.85
	Average	6.99	6.20	29.56	10.64	9.47	43.03
	Portion [%]	68.65	63.36	98.02	70.90	64.15	98.68
4 JOINT OPTIMIZATION
In this section we illustrate the advantage of joint optimisation over the optimisation based solely on forecasting data. Namely, the smoothing and initial parameters calculated by optimisation of the forecasting method (the minimisation of the MSE) are not optimal values for minimising the supply chain costs. Therefore, joint optimisation was used (notation J method) and the initial and smoothing parameters were optimised to minimise the average total costs of the supply chain (Eq. (17)).
4.1 Real Time Series from the M3-Competition
From the joint optimisation of supply chain (costs) model for the time series from the M3-Competition (see Table 5), we observe the following: joint optimization with the MoHW method (JMoHW) outperforms all other methods for all disciplines and it is particularly good for micro and macroeconomic time series. Also, the JMoHW method outperforms all other methods for all types and it is particularly good for the AAA and MMM patterns.
At the bottom of Table 5, we can see that the SREM1 and a portion of series with an improved SREM1 increase as the penalty increases.
Table 6. Averages of the SREM1 for the time series from the M3-Competition (a comparison between the joint model and the model where the forecasting/inventory problem is treated separately)
			SREM1 (AC) [%]		
	JOINT/notJoint	penalty	= 3	penalty	= 5
		JMOHW/ AHW	JMOHW/ MHW	JMOHW/ AHW	JMOHW/ MHW
	MICRO	28.38	27.15	42.60	41.61
e in	MACRO	25.27	25.44	39.70	39.86
ipl 'o	INDUSTRY	25.36	24.01	41.11	39.91
Di	FINANCE	24.27	24.55	37.84	38.11
	DEMOGRAPHIC	15.31	15.46	27.39	27.27
	ANN	23.76	23.31	38.01	37.49
	ANA	39.28	41.19	42.07	43.52
	AAN	43.62	45.00	30.35	31.22
	AAA	58.30	59.83	45.96	46.70
	MNN	41.08	41.51	43.33	42.73
e p	MNA	42.84	44.35	41.41	40.64
ypT	MNM	47.13	44.78	44.46	43.23
	MAN	41.94	44.19	35.66	36.04
	MAA	41.54	43.51	39.61	42.81
	MAM	42.97	41.84	45.61	40.85
	MMN	44.28	45.12	32.12	32.06
	MMM	36.14	30.54	50.86	46.41
	Average	25.27	24.90	39.52	39.21
Comprehensive Energy Resource Management for Essential Reduction of the Total Cost
687
Strojniski vestnik - Journal of Mechanical Engineering 62(2016)11, 685-694
Table 7. Averages of the SREM1 obtained with joint optimisation for the simulated time series
SREM1 (AC) [%]
JOINT	_penalty = 3_penalty = 5_
_JMOHW / JAHW JMOHW / JMHW JMOHW / ETS JMOHW / JAHW JMOHW / JMHW JMOHW / ETS
AAM_5_50_52.98_3671_80.46_65.61_58.67_85.74
AAM_5_100_6123_44.29_8680_6756_44.93_88.69
AAM_8_100_5813_53.95_85.14_7177_53.^_89.80
Average_57.45_44.98_8413_68:3_52.23_88.08
Portion [%]_100.00_94.44_100.00_100.00_in_100.00
MAM_5_50_29:80_37.00_65.41_45.03_39.08_75.36
MAM_5_100_16.85_9.23_75.66_14.44_8.52_79.43
MAM_20_50_5315_19.56_76.45_65.00_26.93_81.52
Average_3327_21.93_72.50_41.49_24.85_78.77
Portion [%]	96.67	95.56	100.00	98.89	91.11	100.00
Finally, if we use the JMoHW instead of the models where forecasts are calculated with the AHW or MHW method regarding minimising the MSE, we can observe the following (see Table 6): for penalty = 5 (penalty = 3) the JMoHW (on average) can reduce the average costs by at least 37.84 % (24.27 %) for all disciplines, except for a demographic series, where the JMoHW can reduce costs "only" by 27.27 % (15.31 %) in comparison with the AHW (the MHW) method.
For different types of series we can observe that for penalty = 3 the JMoHW can reduce the average costs for the ANN and MMM patterns "only" by 23 % to 36 % and for all other types by 39 % to 60 %. For penalty = 5, the JMoHW can reduce the costs for patterns with no season "only" by 30 % to 43 % and for other types by 40 % to 51 %.
4.2 Simulated Time Series
For the AAM patterns we can observe (see Table 7) that the improvement of the JMoHW with respect to all others methods increases as seasonality increases, except for the JMHW for penalty = 5. Also, for penalty = 5 the improvement of the JMoHW with respect to all others methods increases as the slope increases, and this also holds for the JMHW at penalty = 3.
For the MAM patterns the improvement of the JMoHW with respect to the ETS increases and with respect to the JAHW and the JMHW decreases as seasonality increases. Also, the improvement of the JMoHW with respect to the JAHW and the ETS increases and with respect to the JMHW decreases as the slope increases.
Finally, if we use JMoHW instead of the models where forecasts are calculated with the AHW or MHW method regarding minimising the MSE, we can observe (see Table 8) that the JMoHW can essentially reduce the average costs on average by more than 60 % (48 %) in comparison with the AHW (MHW) method. The improvement of the JMoHW with respect to the AHW and the MHW at the AAM patterns is greater than at the MAM patterns and increases as the penalty increases for both types of patterns.
Table 8. Averages of the SREM1 for a simulated time series (a comparison between the joint model and the model where the forecasting/inventory problem is treated separately)
		SREM1 (AC) [%]		
JOINT/	penalty	= 3	penalty	= 5
not_joint	JMOHW/	JMOHW/	JMOHW/	JMOHW/
	AHW	MHW	AHW	MHW
AAM 5 50	58.91	45.60	71.20	61.66
AAM 5 100	64.93	54.93	71.46	63.27
AAM 8 100	63.99	53.85	75.54	66.88
Average	62.61	51.46	72.73	63.94
MAM 5 50	44.78	36.45	61.70	54.89
MAM 5 100	44.50	23.94	55.24	34.51
MAM 20 50	55.50	40.61	66.30	54.10
Average	48.26	33.66	61.08	47.83
5 CONCLUSION AND FURTHER RESEARCH
The biggest challenges facing the European energy policy are sustainable energy supply, reducing dependence on energy resources and energy efficiency. Appropriate energy-demand forecasting methods offer many opportunities regarding process optimisation and appropriate strategic decisions.
686
Ferbar Tratar, L.
Strojniski vestnik - Journal of Mechanical Engineering 62(2016)11, 685-694
This paper exposes two problems: the problem of forecasting energy demand with large noise and the problem of the local optimisation of forecasting methods. We propose the modified HW method for a simultaneous optimisation of demand forecasting and total costs.
From our study of 756 quarterly real series from the M3-Competition and 300 simulated demand patterns we can conclude that the average costs can always be reduced if we use joint optimisation with the MoHW method. For real data they can be reduced on average by more than 24 % for penalty = 3 and by more than 39 % for penalty = 5 in comparison with the models where forecasts are calculated with the AHW or MHW methods regarding minimising the MSE and are treated separately from the cost model. When joint optimisation with the modified HW method for simulated data was used, the average costs were reduced on average by more than 55 % (42 %) for penalty = 3 and by more than 66 % (55 %) for penalty = 5 in comparison with the AHW (MHW).
Joint optimization ensures better planning performance with energy resources, better utility of power plant assests and much more efficient production of derivate energy products from primary energy sources. Finally, the result shows that sharing information and comprehensive energy resource management lead to essential reduction of the total cost.
Further research will include a detailed study of the impact of the modified HW method on the bullwhip effect if the proposed joint optimisation is used. Also, relaxation of the limitation for the lead time assumed to be only one period will be analysed for better evaluation of the joint model over longer horizons. It would be also interesting to investigate how the proposed joint optimization with the modified HW method would effect the systems described in [22].
6 REFERENCES
[1]	European Commission. 2050 Energy strategy, from http:// ec.europa.eu/energy/en/topics/energy-strategy/2050-energy-strategy, accessed on 2016-08-18.
[2]	European Commission. Global Europe 2050, from https:// ec.europa.eu/research/social-sciences/pdf/policy_reviews/ global-europe-2050-report_en.pdf, accessed on 2016-08-18.
[3]	European Commission. Energy roadmap 2050, from https:// ec. europa.eu/energy/sites/ener/files/documents/2012_ energy_roadmap_2050_en_0.pdf, accessed on 2016-08-18.
[4]	Gracanin, D., Lalic, B., Beker, I., Lalic, D., Buchmeister, B. (2013). Cost-time profile simulation for job shop scheduling
decisions. International Journal of Simulation Modelling, vol. 12, no. 4, p. 213-224, DOI:10.2507/IJSIMM12(4)1.237.
[5]	Kremljak, Z., Palcic, I., Kafol, C. (2014). Project evaluation using cost-time investment simulation. International Journal of Simulation Modelling, vol. 13(4), p. 447-457, DOI:10.2507/ IJSIMM13(4)5.279.
[6]	Ferbar Tratar, L., Strmčnik, E. (2016). The comparison of Holt-Winters method and Multiple regression method:a case study. Energy, vol. 109, p. 266-276, DOI:10.1016/j. energy.2016.04.115.
[7]	Potočnik, P., Strmčnik, E., Govekar, E. (2015). Linear and neural network-based models for short-term heat load forecasting. Strojniški vestnik - Journal of Mechanical Engineering, vol. 61, no. 9, p. 543-550, DOI:10.5545/sv-jme.2015.2548.
[8]	Strijbosch, L.W.G., Syntetos, A.A., Boylan, J.E., Janssen, E. (2011). On the interaction between forecasting and stock control: The case of non-stationary demand. International Journal of Production Economics, vol. 133, no. 1, p. 470-480, DOI:10.1016/j.ijpe.2009.10.032.
[9]	Ferbar Tratar, L. (2010). Joint optimisation of demand forecasting and stock control parameters. International Journal of Production Economics, vol. 127, no. 1, p. 173-179, DOI:10.1016/j.ijpe.2010.05.009.
[10]	Ma, Y., Wang, N., Che, A., Huang, Y., Xu, J. (2013). The bullwhip effect on product orders and inventory: a perspective of demand forecasting techniques. International Journal of Production Research, vol. 51, no. 1, p. 281-302, DOI:10.1080/ 00207543.2012.676682.
[11]	Gardner, Jr. E.S. (1985). Exponential smoothing: The state of the art. Journal of Forecasting, vol. 4, no. 1, p. 1-28, DOI:10.1002/for.3980040103.
[12]	Makridakis, S., Hibon, M. (2000). The M3-Competition: results, conclusions and implications. International Journal of Forecasting, vol. 16, no. 4, p. 451-476, DOI:10.1016/S0169-2070(00)00057-1.
[13]	Glavan, I., Prelec, Z., Pavkovic, B. (2015). Modelling, simulation and optimization of small-scale CCHP energy systems. International Journal of Simulation Modelling, vol. 14, no. 4, p. 683-696, DOI:10.1016/S0169-2070(00)00057-1.
[14]	R Core Team (2014). R: A language and environment for statistical computing. Version 3.1.0. R Foundation for Statistical Computing, Vienna.
[15]	Ypma, J., Borchers, H.W. (2014). nloptr: R interface to NLopt. https://CRAN. R-project. org/web/package=nloptr, accessed on 2016-10-03.
[16]	Hyndman, R.J., Akram, M., Bergmeir, C. (2013). Mcomp: Data from the M-competitions. R package version 2.05.
[17]	Hyndman, R.J., Athanasopoulos, G., Razbash, S., Schmidt, D., Zhou, Z., Khan, Y., Bergmeir, C., Wang, E. (2014). Forecast: Forecasting functions for time series and linear models. R package version 2.05.
[18]	Hyndman, R.J., Khandakar, Y. (2008). Automatic time series forecasting: The forecast package for R. Journal of Statistical Software, vol. 27, no. 3, DOI:10.18637/jss.v027.i03.
[19]	Zhao, X., Lee, T.S. (1993). Freezing the master production schedule in multilevel material requirements planning systems under demand uncertainty. Journal of Operations
Comprehensive Energy Resource Management for Essential Reduction of the Total Cost
687
Strojniski vestnik - Journal of Mechanical Engineering 62(2016)11, 685-694
Management, vol. 11, no. 2, p. 185-205, DOI:10.1016/0272-6963(93)90022-H.
[20]	Bermudez, J.D., Segura, J.V., Vercher, E. (2006). A decision support system methodology for forecasting of time series based on soft computing. Computational Statistics & Data Analysis, vol. 51, no. 1, p. 177 - 191, DOI:10.1016/j. csda.2006.02.010.
[21]	Escuin, D., Polo, L., Cipres, D. (2017). On the comparison of inventory replenishment policies with time-varying stochastic
demand for the paper industry. Journal of Computational and Applied Mathematics, vol. 309, p. 424-434, DOI:10.1016/j. cam.2016.03.027.
[22] Šljivac, D., Nakomčic-Smaragdakis, B., Vukobratovic, M., Tople, D., Čepic, Z. (2014). Cost-benefit comparison of on-grid photovoltaic systems in pannonian parts of Croatia and Serbia. Tehnički vjesnik - Technical Gazette, vol. 21, no. 5, p. 1149-1157.
686
Ferbar Tratar, L.