Electric Load Forecasting in the Presence of Active Demand
51st IEEE Conference on Decision and Control December 10-13, 2012. Maui, Hawaii, USA
Electric load forecasting in the presence of Active Demand Simone Paoletti, Member, IEEE, Andrea Garulli, Senior Member, IEEE, Antonio Vicino, Fellow, IEEE
Abstract— Active Demand (AD) is a new concept in smart grids developed within the EU project ADDRESS. It refers to the active participation of households and small commercial consumers in energy systems by means of the flexibility they can offer. Upon receiving real-time price/volume signals, consumers may find convenient to change their load profiles in return of a monetary reward. In this way, they can contribute to the provision of services to the different participants in the electricity system. Since AD causes modifications of the typical consumers’ behaviour, classical load forecasting tools not considering AD signals as inputs are expected to give inaccurate results when applied to load time series including AD effects. In this paper, we study this problem by comparing the prediction performances of several linear models of the load exploiting or not AD signals as inputs. The comparison shows that enhanced prediction results can be obtained by suitably combining the use of AD inputs and the extraction of seasonal characteristics. This is demonstrated by applying the considered approaches to simulated AD effects added to real measurements, representing the aggregated load of about 60 consumers from an Italian LV network. Index Terms— Load forecasting, active demand, smart grids.
I. I NTRODUCTION Small residential and commercial consumers may become “active” players of the electricity system not only through distributed generation, but also by offering their flexibility to the players of the energy markets. The latter concept is known as Active Demand (AD), and relies on the idea that consumers may be willing to change their consumption patterns in return of a monetary reward. In this way, it becomes possible to modify the load profiles at given nodes of the network, which represents a new opportunity for solving network constraints and supporting the development of renewable energy sources, potentially providing economic benefits to all the participants in the electricity system. Along with distributed generation and energy storage systems, it is expected that AD will play a crucial role in the smart grids of the future. The development of a comprehensive commercial and technical architecture to enable AD in the smart grids, has been one of the main targets of the EU Project ADDRESS [1]-[4]. This project, next to the concept of AD, has pursued the functional specification of a new player in the electricity system, called the aggregator. Indeed, since individual consumers do not have direct access to the energy markets, a new intermediary function is needed to bridge this gap. Moreover, S. Paoletti, A. Garulli and A. Vicino are with Dipartimento di Ingegneria dell’Informazione and Centro per lo Studio dei Sistemi Complessi, Universit`a degli Studi di Siena, via Roma 56, 53100 Siena, Italy (e-mail: [email protected]; [email protected]; [email protected]).
978-1-4673-2064-1/12/$31.00 ©2012 IEEE
to fully exploit the consumers’ flexibility, aggregation is a central concept. The role of the aggregator is therefore twice: i) gather the flexibilities provided by its pool of subscribers to form AD-based services; and, ii) offer the AD-based services to the electricity system participants through various markets. The availability of reliable tools for the forecasting of load and distributed generation is a key requirement in active distribution networks. These tools provide useful information to the Distribution System Operator (DSO) for the management of energy flows in the network in order to balance supply and demand (dispatching), thus ensuring continuity and reliability of service provision. This activity is becoming more and more important with the increasingly widespread diffusion of distributed generation, causing parts of the LV or MV network to become “active” when the generated power exceeds the loads. Besides the motivations described above, in the ADDRESS framework the need for load forecasts at distribution network level also arises due to the presence of AD. The DSO has the role to validate the AD-based services sent by the aggregators to the markets by checking if they are compatible with the network constraints and operation. To do this, the DSO needs new forecasting tools able to take into account the load modifications determined by AD. Load forecasting is a well consolidated subject of research [5]. In the wide range of methodologies offered by the literature, a quite general approach is to estimate a model of the load from past data, and then use this model to predict the future load. Linear models are based on the assumption that the current load is a linear function of past values of the load itself and, possibly, of exogenous inputs such as temperature, day type, etc. The most common linear models, borrowed from time series analysis, are the AutoRegressive Moving Average (ARMA) models [6], possibly with eXogenous inputs (ARMAX), or the more general class of transfer function models [7], [8]. These models are often coupled with suitable data pre-processing tools, in order to account for seasonal characteristics at different time-scales (day, week, year). In spite of their simplicity, linear models have been shown to perform satisfactorily in short-term load forecasting [6]. Alternative approaches are based on nonlinear models such as Nonlinear AutoRegressive models with eXogenus inputs (NARX), Neural Networks (NN) [9], and Support Vector Machines (SVM) [10]. Note, however, that none of the mentioned techniques considers AD signals as inputs to the load model. The main contribution of this paper is to analyse the implications of AD in load forecasting. Since AD causes modifications of the typical consumers’ behaviour, classical
2395
load forecasting tools not considering AD signals as inputs are expected to give inaccurate results when applied to load time series including AD effects. To show this, the prediction performances of several linear models of the load are compared in this paper. Since data sets including AD are unfortunately not yet available, the comparison is based on data sets obtained by adding simulated AD effects to real measurements representing the aggregated load of about 60 consumers from an Italian LV network. The compared models are the following: 1) a model with pre-processing, not using AD signals as inputs (this corresponds to treating a load time series as it would be free of AD effects); 2) a model without pre-processing, but using AD signals as inputs; 3) a model with pre-processing and using AD signals as inputs. Pre-processing consists in the preliminary extraction of the seasonal components via an exponential smoothing. The third model is a refinement of the model previously proposed in [11]. To the best of the authors’ knowledge, [11] is the first contribution where load forecasting in the presence of AD has been addressed. The presented comparison clearly shows that the third model provides the best prediction results. It also shows that, the more frequent the AD-based services are, the better the performance of the third model is with respect to the other two models. The paper is organized as follows. Section II describes the mechanism through which an AD-based service built by an aggregator turns out into an actual load modification. The problem of load forecasting in the presence of AD is then formulated and analysed in Section III, and three models of the load in the presence of AD are presented in Section IV. The three models are compared in Section V by using data sets with AD effects simulated as described in Section II. Finally, conclusions are drawn in Section VI.
5) The aggregator receives the validated AD-based service from the DSO (curtailment is possible). 6) The aggregator sends appropriate price/volume signals to the consumers. It is very important to distinguish the validated AD profile (step 5), which we denote by ad, and the actual AD profile provided by the consumers in response to the request received from the aggregator (step 6), which we denote by adtrue . Indeed, the aggregator builds the AD-based services by managing the flexibility of the consumers’ load profiles, but, in principle, there is no guarantee that the consumers will adjust their profiles accordingly. For this reason, adtrue will be different from ad, in general. Examples of AD profiles are shown in Fig. 1. The upper part of the figure shows a standard AD profile ad, composed of two parts. The first part is the actual AD service, whose profile is characterized by the volume Vad and the duration of the service Tdur . The second part is the energy payback effect, due to the fact that a demand modification may be followed by an opposite sign modification. The lower part of Fig. 1 shows what the corresponding actual AD profile adtrue could be. It is important to clarify that ad is the information available to the DSO that can be used for load forecasting. Therefore, the load models must include some kind of mechanism describing how ad is “transformed” into adtrue . Such a model can be based on realistic considerations about the consumers’ behaviour: •
•
An AD signal has a finite duration in time, and therefore it can be assumed that the duration of the corresponding consumers’ response is finite. The consumers may not comply exactly with the request. This results into a delayed and/or partial response with respect to the requested AD profile.
If the model of the consumers’ response is assumed to be linear and time-invariant, a simple FIR model is appropriate:
II. T HE MECHANISM OF ACTIVE D EMAND
adtrue (k) = B(q)ad(k) + v(k),
In this section, we briefly describe the mechanism through which AD-based services built by an aggregator turn out into actual load modifications. This is important in order to understand how AD should be modeled. Moreover, we provide a simple model of the consumers’ response to AD signals which will be used in Section V to generate artificial data sets including AD effects. The AD mechanism consists roughly in the following steps having the aggregator, the DSO and the consumers as actors: 1) The aggregator identifies the AD-based service (e.g load reduction/increase) to meet the needs of other electricity system participants. 2) Based on the consumers’ flexibility, the aggregator builds an AD-based service. 3) The aggregator offers the AD-based service on the market. 4) If the AD-based service is sold, the aggregator submits it to the DSO for technical validation.
where k is the discrete-time index, v(k) is a zero-mean white error process, modelling the random perturbations, and B(q) is a polynomial in the backward shift operator q −1 such that q −1 x(k) = x(k − 1): B(q) = b0 + b1 q −1 + . . . + bnb q −nb .
(1)
(2)
Note that the dc gain of (2) can be interpreted as the steadystate “responsiveness” of the consumers. If it is less than 1, it means that the consumers are typically not fully compliant with the AD signals. Remark 1: The main limitation of model (1) could be that the consumers’ response to AD signals is assumed to be time-invariant. Alternative models of the consumers’ behaviour could be of course considered, to take into account both nonlinear and time-varying effects. More reliable models of the consumers’ response will be estimated when real data sets including AD effects will become available.
2396
payback
V AD profile (kW)
“predictor” is then given by yˆ(k + h|k) = f (z(k); θ∗ ). Another approach, which is pursued in this paper, is based on the decomposition of the load into its components. Indeed, exploiting the structure of the problem is expected both to enhance the prediction accuracy and lead to more parsimonious models. In view of the discussion in Section II and of “good practices” in load forecasting, the following decompositions can be done: • AD profiles are expressed as variations with respect to the expected load profile yb if no AD request is sent. Therefore, the load can be primarily decomposed as follows: y(k) = yb (k) − adtrue (k). (4)
ad
0
Tdur
actual AD profile (kW)
time
0
time
Fig. 1. AD profiles: (top) standard AD profile composed by an ADbased service (with volume Vad and time duration Tdur ) and a consequent energy payback effect; (bottom) example of actual AD profile provided by the consumers.
•
III. P ROBLEM FORMULATION AND ANALYSIS In this section, the load forecasting problem in the presence of AD is formulated and analysed in view of the discussion in Section II. Basic notations are firstly recalled. The electric load in a given area of the distribution network is denoted by y. Ts denotes the sampling time in minutes. It is assumed that nh = 60/Ts is an integer number, representing the number of samples per hour. It follows that nd = 24 · nh is the number of samples per day, while nw = 7 · nd is the number of samples per week. Typical values in the ADDRESS framework are Ts = 15 min, nh = 4, nd = 96 and nw = 672. The sample of a variable x at time kTs from the time origin is denoted by x(k), where k = 0, 1, 2, . . . is the discrete-time index. The forecast of x(k + h), with h a positive integer, computed using only the information available up to time k, is denoted by x ˆ(k + h|k). We are now ready to state the considered load forecasting problem. Problem 1: For fixed lead time h > 0, predict the electric load at time k + h based on the following information available at time k: • load observations y up to time k; • AD signal ad up to time k + h. Note that future values of the AD signal ad(i) for i = k + 1, . . . , k + h are assumed to be known. This is consistent because AD schedules are known in advance to the DSO for validation purposes. Other input signals, such as the temperature, can be considered; they are omitted in this paper to simplify the presentation. Problem 1 could be in principle tackled in a completely black-box framework by choosing a parametric mapping f (·; θ) and select θ = θ∗ such that some norm of the errors e(k; θ) = y(k + h) − f (z(k); θ)
(3)
is minimized, where z(k) is a regression vector containing the information available at time k, including AD. The
Note that positive (negative) values for ad and adtrue mean a load decrease (increase) with respect to yb . Classical load time series (i.e. not including AD effects) show a strong seasonal behaviour. Preliminary extraction of known seasonal characteristics at different time scales (day, week, etc.) helps the model estimation procedure to better capture the stochastic component of the underlying data generation mechanism. Based on this, yb can in turn be decomposed as follows: yb (k) = b(k) + rb (k),
(5)
where b is the so-called base load (the seasonal component) and rb is the residual due to stochastic fluctuations. By substituting (5) into (4), one obtains: y(k) = b(k) + rb (k) − adtrue (k) = b(k) + r(k),
(6)
where the new residual r(k) = rb (k) − adtrue (k)
(7)
takes into account all perturbations to the base load b determined by both stochastic fluctuations and AD. IV. L OAD MODELS In this section, we describe three linear models of the electric load in the presence of AD. The first model does not exploit AD inputs. This corresponds to treating load time series with AD effects by means of currently available load forecasting tools, not taking AD into account. The second and third model both exploit AD inputs, but the third is also coupled with suitable pre-processing tools, in order to account for seasonal characteristics. A. Model with pre-processing, without AD inputs We first consider a model which does not exploit AD inputs. In this case, the actual load y in (4) is treated as yb only. The base load b in (6) is obtained by applying an exponential smoothing to y according to: b(k) = α y(k − n(k) ) + (1 − α) b(k − n(k) ),
(8)
where the delay n(k) depends on whether the time index k falls within workdays, Saturdays, or Sundays and holidays, and α ∈ (0, 1) is the smoothing parameter. Values of α close to one give greater weight to recent changes in the
2397
100
AD volume (kW)
AD volume (kW)
100 50 0 −50 −100
50 0 −50 −100
0 1 2 3 4 5 6 7 8 9 1011121314151617181920212223242526272829
500
400
400
300 200 100 0
0 1 2 3 4 5 6 7 8 9 1011121314151617181920212223242526272829
weeks
500
load (kW)
load (kW)
weeks
300 200 100 0
0 1 2 3 4 5 6 7 8 9 1011121314151617181920212223242526272829
weeks
weeks
(a) Data set #1: Tint ∈ [12, 72] hours. Fig. 2.
0 1 2 3 4 5 6 7 8 9 1011121314151617181920212223242526272829
(b) Data set #2: Tint ∈ [1, 24] hours.
Data sets used in the numerical example of Section V: (top) AD signal ad(k); (bottom) load signal y(k).
data, while values of α close to zero determine a stronger smoothing effect on the stochastic fluctuations. The residuals r in (6) are modeled through an ARMA model of the type: r(k) =
C(q) e(k), D(q)
(9)
C. Model with both pre-processing and AD inputs A load model which exploits both AD inputs and the extraction of seasonal characteristics can be defined as follows [11]. The idea is to roughly estimate yb as: yˆb (k) = y(k) + ad(k),
where e(k) is a zero-mean white stochastic process, and C(q) and D(q) are polynomials in q −1 : C(q) = 1 + c1 q −1 + . . . + cnc q −nc
(10a)
D(q) = 1 + d1 q −1 + . . . + dnd q −nd .
(10b)
i.e. by replacing adtrue with ad in (4). The base load b in (6) is then obtained by applying an exponential smoothing to yˆb , i.e. by replacing y with yˆb in (8). In view of (7), the residuals r in (6) are modeled through a transfer function model of the type:
Predicted values of the load are then computed as follows: yˆ(k + h|k) = b(k + h) + rˆ(k + h|k),
r(k) = −B(q) ad(k) +
(11)
where rˆ(k+h|k) is obtained by applying the linear minimum variance h-step ahead predictor derived from (9), and b(k+h) is calculated from (8). For simplicity of exposition, we assume h ≤ n(k) (for longer lead times, different strategies for computing b(k + h) can be adopted). The ARMA model (9) coupled with the exponential smoothing (8) will be referred to as ES-ARMA model. B. Model without pre-processing, with AD inputs
C(q) e(k), D(q)
(12)
where e(k) is a zero-mean white stochastic process, B(q) is defined in (2), and C(q) and D(q) are defined in (10). Hereafter we will denote model (12) as Transfer Function (TF) model. Note that model (12) describes the actual load y(k) without data pre-processing. This means that the ARMA C(q) e(k) in (12) models the component yb of y, which term D(q) includes both seasonal and stochastic terms. Predicted values yˆ(k + h|k) are computed by applying the linear minimum variance h-step ahead predictor derived from (12).
C(q) e(k), D(q)
(14)
where e(k) is a zero-mean white stochastic process, B(q) is defined in (2), and C(q) and D(q) are defined in (10). The resulting model will be hereafter referred to as ES-TF C(q) model. Note that the ARMA term D(q) e(k) in (14) models the residuals rb in (5), which mainly include stochastic contributions. Predicted values yˆ(k + h|k) are computed through (11), where, in particular, rˆ(k + h|k) is obtained by applying the linear minimum variance h-step ahead predictor derived from (14).
In view of (1) and (4), a linear model of the electric load which exploits AD inputs can be built as follows: y(k) = −B(q) ad(k) +
(13)
V. N UMERICAL TESTS The results presented in this section concern the comparison of the prediction performance of the models described in Section IV on two data sets with AD effects. The aim is to show that combining the use of AD inputs and the estimation of seasonal components, as the ES-TF model of Section IV-C does, is a promising direction for developing effective load forecasting tools in the presence of AD. Since the two data sets are generated by simulating AD effects with different frequencies of the AD-based services, the proposed comparison also shows that exploiting the AD inputs becomes more and more important as the AD-based services become more frequent.
2398
130 120 residual autocorrelation
0.8
MSE
110 100 90 80
0.4 0.2 0
70 60 0
0.6
0.2
0.4 0.6 smoothing parameter
0.8
−0.2 0
1
40
60
80
100
lag
(a) Choice of the smoothing parameter α. Fig. 3.
20
(b) Residual autocorrelation function.
Estimation of the ES-TF model using the data set #1.
The original data set contains real measurements representing the aggregated electric load of about 60 consumers from an Italian LV network. The data set covers 29 weeks from April to October 2008 with sampling time Ts = 15 min, for a total of 19488 samples. Temperature records are not available. AD effects are simulated using model (1) with nb = 2, b0 = 0.55, b1 = 0.30, b2 = 0.05, and σv2 = 25.00, where σv2 denotes the variance of the noise v(k). The input AD profiles ad(k) are generated by selecting randomly the AD volume Vad in the interval [30, 100] kW and the duration Tdur in the interval [1, 4] hours. The energy payback is proportional to Vad with a factor 0.3, and its duration is assumed to be half of Tdur . Two AD signals are generated, differing for the range of values of the random interval Tint between the start times of two consecutive AD-based services: • Tint ∈ [12, 72] hours in data set #1; • Tint ∈ [1, 24] hours in data set #2. AD effects are then added to the original load time series. Fig. 2 shows the resulting data sets, with the AD signals ad(k) in the upper part, and the load y(k) including AD effects in the lower part. A. Model estimation The data sets are partitioned in two parts. The first 14 weeks are used for model estimation, while the remaining 15 weeks are used for model validation. The estimation procedure is different for the TF model, and for the linear models coupled with the exponential smoothing, which require also the estimation of the smoothing parameter α. For fixed model orders nb , nc and nd , a TF model is estimated by minimising the Mean Square Error (MSE). In the case of the models ES-ARMA and ES-TF, a gridding of the interval [0, 1] is considered for the parameter α, and for each value of α an ARMA or TF model of the corresponding residuals r(k) = y(k) − b(k) is estimated by minimising the MSE. For each data set, the value of α
(and associated linear model) with smallest MSE is selected. Fig. 3(a) shows the MSE as a function of α and the resulting optimal choice α = 0.15, in the case of the ES-TF model with nb = 2 and nc = nd = 4, using the estimation data set #1. Fig. 3(b) shows the residual autocorrelation function for the selected ES-TF model. Quite nicely, all samples lie within or very close to the confidence regions. For all model types, order selection is carried out by estimating models for different combinations of the orders nb , nc and nd , and then selecting the one with the smallest MSE computed on validation data. B. Comparison The models of the three types selected in Section V-A are used to predict the load for different lead times h ranging from one sample (15 min ahead) to 96 samples (one day ahead). For fixed lead time h, the performance of each model is evaluated by computing a normalized form of the root MSE using validation data: √ M SE × 100%, (15) NRMSE = |y|
where y is the sample mean of y. Smaller values of NRMSE indicate a better agreement between measured and predicted values. The use of a non-dimensional form of the root MSE makes it possible to compare performances obtained on different data sets. The plots of the NRMSE versus the lead time are shown for all three model types in Fig. 4(a), when the NRMSE is computed using the validation data set #1, and Fig. 4(b), when the NRMSE is computed using the validation data set #2. It can be observed from both figures that the ES-TF model outperforms the other two models for all lead times, which confirms our claim at the beginning of Section V. It can be further observed that for very short lead times (less than one hour) the models of the three types show similar performance, and that, as expected, the performance is worse
2399
45
40
40
35
35
30
30
NRMSE (%)
NRMSE (%)
45
25 20
20
15
15
10
10
5
5
0 0
4
8
12 16 lead time (hours)
20
0 0
24
(a) Validation data set #1. Fig. 4.
25
4
8
12 16 lead time (hours)
20
24
(b) Validation data set #2.
Plots of the NRMSE versus the lead time for the three compared models: TF (solid), ES-ARMA (dash-dotted), ES-TF (dashed).
for all models when the lead time increases. However, the performance of both the TF model and the ES-ARMA model deteriorates more rapidly in the very short term than the performance of the ES+TF model, with the TF model which soon reaches unacceptable performance. It is worthwhile to note that the relative performance of the ES-TF model with respect to the ES-ARMA model is affected by the Tint parameter used to generate the ADbased services. In particular, it can be observed that when AD-based services are more frequent, as in data set #2, the enhancement of the prediction performance guaranteed by the ES-TF model with respect to the ES-ARMA model is much more significant. Similar effects can be observed if the duration and/or the volume of the AD services are increased. VI. C ONCLUSIONS Active Demand will be a key concept in the smart grids of the future. This calls for the development of new modeling, estimation and control tools, in order to properly take into account the AD effects on the whole electricity system. In this framework, this paper has presented a comparison of three simple models of the load in the presence of AD. The comparison has clearly shown that, by combining the use of AD inputs and the estimation of seasonal components, effective models for load forecasting in the presence of AD can be devised. The impact of AD on smart grids is expected to open up a number of new research and technological issues. Future work in the specific area of load forecasting will concern the use of nonlinear or time-varying models. In particular, models estimated by means of piecewise affine identification techniques [12] will be investigated. The choice of piecewise affine maps may allow one to model the seasonal and cyclical properties of loads and the influence of AD, through an appropriate partition of the model domain. Another fundamental step will concern the testing and validation of forecasting models on real-world data including AD effects. In this respect, first testbed will be provided by field tests on LV and MV networks in France and Italy, which are being
carried out within the ADDRESS project. As soon as data from the tests will be available (probably by April 2013), they will be used as a benchmark for the techniques proposed in this paper. ACKNOWLEDGEMENTS The research leading to these results has received funding from the European Community’s Seventh Framework Programme (FP7/2007-2013) under grant agreement no. 207643. R EFERENCES [1] “The ADDRESS project,” http://www.addressfp7.org/. [2] R. Belhomme, R. Cerero Real de Asua, G. Valtorta, A. Paice, F. Bouffard, R. Rooth, and A. Losi, “ADDRESS - Active Demand for the smart grids of the future,” in Proc. 2008 CIRED Seminar: Smart Grids for Distribution, 2008, paper no. 0080. [3] E. Peeters, R. Belhomme, C. Batlle, A. Paice, F. Bouffard, S. Karkkainen, D. Six, and M. Hommelberg, “ADDRESS - scenarios and architecture for Active Demand development in the smart grids of the future,” in Proc. 2009 CIRED - 20th Int. Conf. on Electricity Distribution, 2009, paper no. 0406. [4] R. Belhomme, R. Cerero Real de Asua, G. Valtorta, and P. Eyrolles, “The ADDRESS project: Developing Active Demand in smart power systems integrating renewables,” in Proc. IEEE/PES General Meeting, 2011, paper no. PESGM2011-000188. [5] H. Alfares and M. Nazeeruddin, “Electric load forecasting: literature survey and classification of methods,” Int. J. of Systems Science, vol. 33, no. 1, pp. 23–34, 2002. [6] S.-J. Huang and K.-R. Shih, “Short-term load forecasting via ARMA model identification including non-Gaussian process considerations,” IEEE Trans. on Power Systems, vol. 18, no. 2, pp. 673–679, 2003. [7] L. Ljung, System identification: theory for the user. Prentice-Hall, Upper Saddle River, NJ, 1999. [8] G. Box, G. Jenkins, and G. Reinsel, Time Series Analysis: Forecasting and Control, 4th ed. Wiley, 2008. [9] D. Park, M. El-Sharkawi, R. Marks II, L. Atlas, and M. Damborg, “Electric load forecasting using an artificial neural network,” IEEE Trans. on Power Systems, vol. 6, no. 2, pp. 442–449, 1991. [10] M. Espinoza, J. Suykens, R. Belmans, and B. De Moor, “Electric load forecasting using kernel-based modeling for nonlinear system identification,” IEEE Control Systems Magazine, vol. 27, no. 5, pp. 43–57, 2007. [11] S. Paoletti, M. Casini, A. Giannitrapani, A. Facchini, A. Garulli, and A. Vicino, “Load forecasting for active distribution networks,” in Proc. 2nd IEEE PES Innovative Smart Grid Technologies Europe, 2011. [12] A. Bemporad, A. Garulli, S. Paoletti, and A. Vicino, “A boundederror approach to piecewise affine system identification,” IEEE Trans. Autom. Control, vol. 50, no. 10, pp. 1567–1580, 2005.
2400
Electric load forecasting in the presence of Active Demand Simone Paoletti, Member, IEEE, Andrea Garulli, Senior Member, IEEE, Antonio Vicino, Fellow, IEEE
Abstract— Active Demand (AD) is a new concept in smart grids developed within the EU project ADDRESS. It refers to the active participation of households and small commercial consumers in energy systems by means of the flexibility they can offer. Upon receiving real-time price/volume signals, consumers may find convenient to change their load profiles in return of a monetary reward. In this way, they can contribute to the provision of services to the different participants in the electricity system. Since AD causes modifications of the typical consumers’ behaviour, classical load forecasting tools not considering AD signals as inputs are expected to give inaccurate results when applied to load time series including AD effects. In this paper, we study this problem by comparing the prediction performances of several linear models of the load exploiting or not AD signals as inputs. The comparison shows that enhanced prediction results can be obtained by suitably combining the use of AD inputs and the extraction of seasonal characteristics. This is demonstrated by applying the considered approaches to simulated AD effects added to real measurements, representing the aggregated load of about 60 consumers from an Italian LV network. Index Terms— Load forecasting, active demand, smart grids.
I. I NTRODUCTION Small residential and commercial consumers may become “active” players of the electricity system not only through distributed generation, but also by offering their flexibility to the players of the energy markets. The latter concept is known as Active Demand (AD), and relies on the idea that consumers may be willing to change their consumption patterns in return of a monetary reward. In this way, it becomes possible to modify the load profiles at given nodes of the network, which represents a new opportunity for solving network constraints and supporting the development of renewable energy sources, potentially providing economic benefits to all the participants in the electricity system. Along with distributed generation and energy storage systems, it is expected that AD will play a crucial role in the smart grids of the future. The development of a comprehensive commercial and technical architecture to enable AD in the smart grids, has been one of the main targets of the EU Project ADDRESS [1]-[4]. This project, next to the concept of AD, has pursued the functional specification of a new player in the electricity system, called the aggregator. Indeed, since individual consumers do not have direct access to the energy markets, a new intermediary function is needed to bridge this gap. Moreover, S. Paoletti, A. Garulli and A. Vicino are with Dipartimento di Ingegneria dell’Informazione and Centro per lo Studio dei Sistemi Complessi, Universit`a degli Studi di Siena, via Roma 56, 53100 Siena, Italy (e-mail: [email protected]; [email protected]; [email protected]).
978-1-4673-2064-1/12/$31.00 ©2012 IEEE
to fully exploit the consumers’ flexibility, aggregation is a central concept. The role of the aggregator is therefore twice: i) gather the flexibilities provided by its pool of subscribers to form AD-based services; and, ii) offer the AD-based services to the electricity system participants through various markets. The availability of reliable tools for the forecasting of load and distributed generation is a key requirement in active distribution networks. These tools provide useful information to the Distribution System Operator (DSO) for the management of energy flows in the network in order to balance supply and demand (dispatching), thus ensuring continuity and reliability of service provision. This activity is becoming more and more important with the increasingly widespread diffusion of distributed generation, causing parts of the LV or MV network to become “active” when the generated power exceeds the loads. Besides the motivations described above, in the ADDRESS framework the need for load forecasts at distribution network level also arises due to the presence of AD. The DSO has the role to validate the AD-based services sent by the aggregators to the markets by checking if they are compatible with the network constraints and operation. To do this, the DSO needs new forecasting tools able to take into account the load modifications determined by AD. Load forecasting is a well consolidated subject of research [5]. In the wide range of methodologies offered by the literature, a quite general approach is to estimate a model of the load from past data, and then use this model to predict the future load. Linear models are based on the assumption that the current load is a linear function of past values of the load itself and, possibly, of exogenous inputs such as temperature, day type, etc. The most common linear models, borrowed from time series analysis, are the AutoRegressive Moving Average (ARMA) models [6], possibly with eXogenous inputs (ARMAX), or the more general class of transfer function models [7], [8]. These models are often coupled with suitable data pre-processing tools, in order to account for seasonal characteristics at different time-scales (day, week, year). In spite of their simplicity, linear models have been shown to perform satisfactorily in short-term load forecasting [6]. Alternative approaches are based on nonlinear models such as Nonlinear AutoRegressive models with eXogenus inputs (NARX), Neural Networks (NN) [9], and Support Vector Machines (SVM) [10]. Note, however, that none of the mentioned techniques considers AD signals as inputs to the load model. The main contribution of this paper is to analyse the implications of AD in load forecasting. Since AD causes modifications of the typical consumers’ behaviour, classical
2395
load forecasting tools not considering AD signals as inputs are expected to give inaccurate results when applied to load time series including AD effects. To show this, the prediction performances of several linear models of the load are compared in this paper. Since data sets including AD are unfortunately not yet available, the comparison is based on data sets obtained by adding simulated AD effects to real measurements representing the aggregated load of about 60 consumers from an Italian LV network. The compared models are the following: 1) a model with pre-processing, not using AD signals as inputs (this corresponds to treating a load time series as it would be free of AD effects); 2) a model without pre-processing, but using AD signals as inputs; 3) a model with pre-processing and using AD signals as inputs. Pre-processing consists in the preliminary extraction of the seasonal components via an exponential smoothing. The third model is a refinement of the model previously proposed in [11]. To the best of the authors’ knowledge, [11] is the first contribution where load forecasting in the presence of AD has been addressed. The presented comparison clearly shows that the third model provides the best prediction results. It also shows that, the more frequent the AD-based services are, the better the performance of the third model is with respect to the other two models. The paper is organized as follows. Section II describes the mechanism through which an AD-based service built by an aggregator turns out into an actual load modification. The problem of load forecasting in the presence of AD is then formulated and analysed in Section III, and three models of the load in the presence of AD are presented in Section IV. The three models are compared in Section V by using data sets with AD effects simulated as described in Section II. Finally, conclusions are drawn in Section VI.
5) The aggregator receives the validated AD-based service from the DSO (curtailment is possible). 6) The aggregator sends appropriate price/volume signals to the consumers. It is very important to distinguish the validated AD profile (step 5), which we denote by ad, and the actual AD profile provided by the consumers in response to the request received from the aggregator (step 6), which we denote by adtrue . Indeed, the aggregator builds the AD-based services by managing the flexibility of the consumers’ load profiles, but, in principle, there is no guarantee that the consumers will adjust their profiles accordingly. For this reason, adtrue will be different from ad, in general. Examples of AD profiles are shown in Fig. 1. The upper part of the figure shows a standard AD profile ad, composed of two parts. The first part is the actual AD service, whose profile is characterized by the volume Vad and the duration of the service Tdur . The second part is the energy payback effect, due to the fact that a demand modification may be followed by an opposite sign modification. The lower part of Fig. 1 shows what the corresponding actual AD profile adtrue could be. It is important to clarify that ad is the information available to the DSO that can be used for load forecasting. Therefore, the load models must include some kind of mechanism describing how ad is “transformed” into adtrue . Such a model can be based on realistic considerations about the consumers’ behaviour: •
•
An AD signal has a finite duration in time, and therefore it can be assumed that the duration of the corresponding consumers’ response is finite. The consumers may not comply exactly with the request. This results into a delayed and/or partial response with respect to the requested AD profile.
If the model of the consumers’ response is assumed to be linear and time-invariant, a simple FIR model is appropriate:
II. T HE MECHANISM OF ACTIVE D EMAND
adtrue (k) = B(q)ad(k) + v(k),
In this section, we briefly describe the mechanism through which AD-based services built by an aggregator turn out into actual load modifications. This is important in order to understand how AD should be modeled. Moreover, we provide a simple model of the consumers’ response to AD signals which will be used in Section V to generate artificial data sets including AD effects. The AD mechanism consists roughly in the following steps having the aggregator, the DSO and the consumers as actors: 1) The aggregator identifies the AD-based service (e.g load reduction/increase) to meet the needs of other electricity system participants. 2) Based on the consumers’ flexibility, the aggregator builds an AD-based service. 3) The aggregator offers the AD-based service on the market. 4) If the AD-based service is sold, the aggregator submits it to the DSO for technical validation.
where k is the discrete-time index, v(k) is a zero-mean white error process, modelling the random perturbations, and B(q) is a polynomial in the backward shift operator q −1 such that q −1 x(k) = x(k − 1): B(q) = b0 + b1 q −1 + . . . + bnb q −nb .
(1)
(2)
Note that the dc gain of (2) can be interpreted as the steadystate “responsiveness” of the consumers. If it is less than 1, it means that the consumers are typically not fully compliant with the AD signals. Remark 1: The main limitation of model (1) could be that the consumers’ response to AD signals is assumed to be time-invariant. Alternative models of the consumers’ behaviour could be of course considered, to take into account both nonlinear and time-varying effects. More reliable models of the consumers’ response will be estimated when real data sets including AD effects will become available.
2396
payback
V AD profile (kW)
“predictor” is then given by yˆ(k + h|k) = f (z(k); θ∗ ). Another approach, which is pursued in this paper, is based on the decomposition of the load into its components. Indeed, exploiting the structure of the problem is expected both to enhance the prediction accuracy and lead to more parsimonious models. In view of the discussion in Section II and of “good practices” in load forecasting, the following decompositions can be done: • AD profiles are expressed as variations with respect to the expected load profile yb if no AD request is sent. Therefore, the load can be primarily decomposed as follows: y(k) = yb (k) − adtrue (k). (4)
ad
0
Tdur
actual AD profile (kW)
time
0
time
Fig. 1. AD profiles: (top) standard AD profile composed by an ADbased service (with volume Vad and time duration Tdur ) and a consequent energy payback effect; (bottom) example of actual AD profile provided by the consumers.
•
III. P ROBLEM FORMULATION AND ANALYSIS In this section, the load forecasting problem in the presence of AD is formulated and analysed in view of the discussion in Section II. Basic notations are firstly recalled. The electric load in a given area of the distribution network is denoted by y. Ts denotes the sampling time in minutes. It is assumed that nh = 60/Ts is an integer number, representing the number of samples per hour. It follows that nd = 24 · nh is the number of samples per day, while nw = 7 · nd is the number of samples per week. Typical values in the ADDRESS framework are Ts = 15 min, nh = 4, nd = 96 and nw = 672. The sample of a variable x at time kTs from the time origin is denoted by x(k), where k = 0, 1, 2, . . . is the discrete-time index. The forecast of x(k + h), with h a positive integer, computed using only the information available up to time k, is denoted by x ˆ(k + h|k). We are now ready to state the considered load forecasting problem. Problem 1: For fixed lead time h > 0, predict the electric load at time k + h based on the following information available at time k: • load observations y up to time k; • AD signal ad up to time k + h. Note that future values of the AD signal ad(i) for i = k + 1, . . . , k + h are assumed to be known. This is consistent because AD schedules are known in advance to the DSO for validation purposes. Other input signals, such as the temperature, can be considered; they are omitted in this paper to simplify the presentation. Problem 1 could be in principle tackled in a completely black-box framework by choosing a parametric mapping f (·; θ) and select θ = θ∗ such that some norm of the errors e(k; θ) = y(k + h) − f (z(k); θ)
(3)
is minimized, where z(k) is a regression vector containing the information available at time k, including AD. The
Note that positive (negative) values for ad and adtrue mean a load decrease (increase) with respect to yb . Classical load time series (i.e. not including AD effects) show a strong seasonal behaviour. Preliminary extraction of known seasonal characteristics at different time scales (day, week, etc.) helps the model estimation procedure to better capture the stochastic component of the underlying data generation mechanism. Based on this, yb can in turn be decomposed as follows: yb (k) = b(k) + rb (k),
(5)
where b is the so-called base load (the seasonal component) and rb is the residual due to stochastic fluctuations. By substituting (5) into (4), one obtains: y(k) = b(k) + rb (k) − adtrue (k) = b(k) + r(k),
(6)
where the new residual r(k) = rb (k) − adtrue (k)
(7)
takes into account all perturbations to the base load b determined by both stochastic fluctuations and AD. IV. L OAD MODELS In this section, we describe three linear models of the electric load in the presence of AD. The first model does not exploit AD inputs. This corresponds to treating load time series with AD effects by means of currently available load forecasting tools, not taking AD into account. The second and third model both exploit AD inputs, but the third is also coupled with suitable pre-processing tools, in order to account for seasonal characteristics. A. Model with pre-processing, without AD inputs We first consider a model which does not exploit AD inputs. In this case, the actual load y in (4) is treated as yb only. The base load b in (6) is obtained by applying an exponential smoothing to y according to: b(k) = α y(k − n(k) ) + (1 − α) b(k − n(k) ),
(8)
where the delay n(k) depends on whether the time index k falls within workdays, Saturdays, or Sundays and holidays, and α ∈ (0, 1) is the smoothing parameter. Values of α close to one give greater weight to recent changes in the
2397
100
AD volume (kW)
AD volume (kW)
100 50 0 −50 −100
50 0 −50 −100
0 1 2 3 4 5 6 7 8 9 1011121314151617181920212223242526272829
500
400
400
300 200 100 0
0 1 2 3 4 5 6 7 8 9 1011121314151617181920212223242526272829
weeks
500
load (kW)
load (kW)
weeks
300 200 100 0
0 1 2 3 4 5 6 7 8 9 1011121314151617181920212223242526272829
weeks
weeks
(a) Data set #1: Tint ∈ [12, 72] hours. Fig. 2.
0 1 2 3 4 5 6 7 8 9 1011121314151617181920212223242526272829
(b) Data set #2: Tint ∈ [1, 24] hours.
Data sets used in the numerical example of Section V: (top) AD signal ad(k); (bottom) load signal y(k).
data, while values of α close to zero determine a stronger smoothing effect on the stochastic fluctuations. The residuals r in (6) are modeled through an ARMA model of the type: r(k) =
C(q) e(k), D(q)
(9)
C. Model with both pre-processing and AD inputs A load model which exploits both AD inputs and the extraction of seasonal characteristics can be defined as follows [11]. The idea is to roughly estimate yb as: yˆb (k) = y(k) + ad(k),
where e(k) is a zero-mean white stochastic process, and C(q) and D(q) are polynomials in q −1 : C(q) = 1 + c1 q −1 + . . . + cnc q −nc
(10a)
D(q) = 1 + d1 q −1 + . . . + dnd q −nd .
(10b)
i.e. by replacing adtrue with ad in (4). The base load b in (6) is then obtained by applying an exponential smoothing to yˆb , i.e. by replacing y with yˆb in (8). In view of (7), the residuals r in (6) are modeled through a transfer function model of the type:
Predicted values of the load are then computed as follows: yˆ(k + h|k) = b(k + h) + rˆ(k + h|k),
r(k) = −B(q) ad(k) +
(11)
where rˆ(k+h|k) is obtained by applying the linear minimum variance h-step ahead predictor derived from (9), and b(k+h) is calculated from (8). For simplicity of exposition, we assume h ≤ n(k) (for longer lead times, different strategies for computing b(k + h) can be adopted). The ARMA model (9) coupled with the exponential smoothing (8) will be referred to as ES-ARMA model. B. Model without pre-processing, with AD inputs
C(q) e(k), D(q)
(12)
where e(k) is a zero-mean white stochastic process, B(q) is defined in (2), and C(q) and D(q) are defined in (10). Hereafter we will denote model (12) as Transfer Function (TF) model. Note that model (12) describes the actual load y(k) without data pre-processing. This means that the ARMA C(q) e(k) in (12) models the component yb of y, which term D(q) includes both seasonal and stochastic terms. Predicted values yˆ(k + h|k) are computed by applying the linear minimum variance h-step ahead predictor derived from (12).
C(q) e(k), D(q)
(14)
where e(k) is a zero-mean white stochastic process, B(q) is defined in (2), and C(q) and D(q) are defined in (10). The resulting model will be hereafter referred to as ES-TF C(q) model. Note that the ARMA term D(q) e(k) in (14) models the residuals rb in (5), which mainly include stochastic contributions. Predicted values yˆ(k + h|k) are computed through (11), where, in particular, rˆ(k + h|k) is obtained by applying the linear minimum variance h-step ahead predictor derived from (14).
In view of (1) and (4), a linear model of the electric load which exploits AD inputs can be built as follows: y(k) = −B(q) ad(k) +
(13)
V. N UMERICAL TESTS The results presented in this section concern the comparison of the prediction performance of the models described in Section IV on two data sets with AD effects. The aim is to show that combining the use of AD inputs and the estimation of seasonal components, as the ES-TF model of Section IV-C does, is a promising direction for developing effective load forecasting tools in the presence of AD. Since the two data sets are generated by simulating AD effects with different frequencies of the AD-based services, the proposed comparison also shows that exploiting the AD inputs becomes more and more important as the AD-based services become more frequent.
2398
130 120 residual autocorrelation
0.8
MSE
110 100 90 80
0.4 0.2 0
70 60 0
0.6
0.2
0.4 0.6 smoothing parameter
0.8
−0.2 0
1
40
60
80
100
lag
(a) Choice of the smoothing parameter α. Fig. 3.
20
(b) Residual autocorrelation function.
Estimation of the ES-TF model using the data set #1.
The original data set contains real measurements representing the aggregated electric load of about 60 consumers from an Italian LV network. The data set covers 29 weeks from April to October 2008 with sampling time Ts = 15 min, for a total of 19488 samples. Temperature records are not available. AD effects are simulated using model (1) with nb = 2, b0 = 0.55, b1 = 0.30, b2 = 0.05, and σv2 = 25.00, where σv2 denotes the variance of the noise v(k). The input AD profiles ad(k) are generated by selecting randomly the AD volume Vad in the interval [30, 100] kW and the duration Tdur in the interval [1, 4] hours. The energy payback is proportional to Vad with a factor 0.3, and its duration is assumed to be half of Tdur . Two AD signals are generated, differing for the range of values of the random interval Tint between the start times of two consecutive AD-based services: • Tint ∈ [12, 72] hours in data set #1; • Tint ∈ [1, 24] hours in data set #2. AD effects are then added to the original load time series. Fig. 2 shows the resulting data sets, with the AD signals ad(k) in the upper part, and the load y(k) including AD effects in the lower part. A. Model estimation The data sets are partitioned in two parts. The first 14 weeks are used for model estimation, while the remaining 15 weeks are used for model validation. The estimation procedure is different for the TF model, and for the linear models coupled with the exponential smoothing, which require also the estimation of the smoothing parameter α. For fixed model orders nb , nc and nd , a TF model is estimated by minimising the Mean Square Error (MSE). In the case of the models ES-ARMA and ES-TF, a gridding of the interval [0, 1] is considered for the parameter α, and for each value of α an ARMA or TF model of the corresponding residuals r(k) = y(k) − b(k) is estimated by minimising the MSE. For each data set, the value of α
(and associated linear model) with smallest MSE is selected. Fig. 3(a) shows the MSE as a function of α and the resulting optimal choice α = 0.15, in the case of the ES-TF model with nb = 2 and nc = nd = 4, using the estimation data set #1. Fig. 3(b) shows the residual autocorrelation function for the selected ES-TF model. Quite nicely, all samples lie within or very close to the confidence regions. For all model types, order selection is carried out by estimating models for different combinations of the orders nb , nc and nd , and then selecting the one with the smallest MSE computed on validation data. B. Comparison The models of the three types selected in Section V-A are used to predict the load for different lead times h ranging from one sample (15 min ahead) to 96 samples (one day ahead). For fixed lead time h, the performance of each model is evaluated by computing a normalized form of the root MSE using validation data: √ M SE × 100%, (15) NRMSE = |y|
where y is the sample mean of y. Smaller values of NRMSE indicate a better agreement between measured and predicted values. The use of a non-dimensional form of the root MSE makes it possible to compare performances obtained on different data sets. The plots of the NRMSE versus the lead time are shown for all three model types in Fig. 4(a), when the NRMSE is computed using the validation data set #1, and Fig. 4(b), when the NRMSE is computed using the validation data set #2. It can be observed from both figures that the ES-TF model outperforms the other two models for all lead times, which confirms our claim at the beginning of Section V. It can be further observed that for very short lead times (less than one hour) the models of the three types show similar performance, and that, as expected, the performance is worse
2399
45
40
40
35
35
30
30
NRMSE (%)
NRMSE (%)
45
25 20
20
15
15
10
10
5
5
0 0
4
8
12 16 lead time (hours)
20
0 0
24
(a) Validation data set #1. Fig. 4.
25
4
8
12 16 lead time (hours)
20
24
(b) Validation data set #2.
Plots of the NRMSE versus the lead time for the three compared models: TF (solid), ES-ARMA (dash-dotted), ES-TF (dashed).
for all models when the lead time increases. However, the performance of both the TF model and the ES-ARMA model deteriorates more rapidly in the very short term than the performance of the ES+TF model, with the TF model which soon reaches unacceptable performance. It is worthwhile to note that the relative performance of the ES-TF model with respect to the ES-ARMA model is affected by the Tint parameter used to generate the ADbased services. In particular, it can be observed that when AD-based services are more frequent, as in data set #2, the enhancement of the prediction performance guaranteed by the ES-TF model with respect to the ES-ARMA model is much more significant. Similar effects can be observed if the duration and/or the volume of the AD services are increased. VI. C ONCLUSIONS Active Demand will be a key concept in the smart grids of the future. This calls for the development of new modeling, estimation and control tools, in order to properly take into account the AD effects on the whole electricity system. In this framework, this paper has presented a comparison of three simple models of the load in the presence of AD. The comparison has clearly shown that, by combining the use of AD inputs and the estimation of seasonal components, effective models for load forecasting in the presence of AD can be devised. The impact of AD on smart grids is expected to open up a number of new research and technological issues. Future work in the specific area of load forecasting will concern the use of nonlinear or time-varying models. In particular, models estimated by means of piecewise affine identification techniques [12] will be investigated. The choice of piecewise affine maps may allow one to model the seasonal and cyclical properties of loads and the influence of AD, through an appropriate partition of the model domain. Another fundamental step will concern the testing and validation of forecasting models on real-world data including AD effects. In this respect, first testbed will be provided by field tests on LV and MV networks in France and Italy, which are being
carried out within the ADDRESS project. As soon as data from the tests will be available (probably by April 2013), they will be used as a benchmark for the techniques proposed in this paper. ACKNOWLEDGEMENTS The research leading to these results has received funding from the European Community’s Seventh Framework Programme (FP7/2007-2013) under grant agreement no. 207643. R EFERENCES [1] “The ADDRESS project,” http://www.addressfp7.org/. [2] R. Belhomme, R. Cerero Real de Asua, G. Valtorta, A. Paice, F. Bouffard, R. Rooth, and A. Losi, “ADDRESS - Active Demand for the smart grids of the future,” in Proc. 2008 CIRED Seminar: Smart Grids for Distribution, 2008, paper no. 0080. [3] E. Peeters, R. Belhomme, C. Batlle, A. Paice, F. Bouffard, S. Karkkainen, D. Six, and M. Hommelberg, “ADDRESS - scenarios and architecture for Active Demand development in the smart grids of the future,” in Proc. 2009 CIRED - 20th Int. Conf. on Electricity Distribution, 2009, paper no. 0406. [4] R. Belhomme, R. Cerero Real de Asua, G. Valtorta, and P. Eyrolles, “The ADDRESS project: Developing Active Demand in smart power systems integrating renewables,” in Proc. IEEE/PES General Meeting, 2011, paper no. PESGM2011-000188. [5] H. Alfares and M. Nazeeruddin, “Electric load forecasting: literature survey and classification of methods,” Int. J. of Systems Science, vol. 33, no. 1, pp. 23–34, 2002. [6] S.-J. Huang and K.-R. Shih, “Short-term load forecasting via ARMA model identification including non-Gaussian process considerations,” IEEE Trans. on Power Systems, vol. 18, no. 2, pp. 673–679, 2003. [7] L. Ljung, System identification: theory for the user. Prentice-Hall, Upper Saddle River, NJ, 1999. [8] G. Box, G. Jenkins, and G. Reinsel, Time Series Analysis: Forecasting and Control, 4th ed. Wiley, 2008. [9] D. Park, M. El-Sharkawi, R. Marks II, L. Atlas, and M. Damborg, “Electric load forecasting using an artificial neural network,” IEEE Trans. on Power Systems, vol. 6, no. 2, pp. 442–449, 1991. [10] M. Espinoza, J. Suykens, R. Belmans, and B. De Moor, “Electric load forecasting using kernel-based modeling for nonlinear system identification,” IEEE Control Systems Magazine, vol. 27, no. 5, pp. 43–57, 2007. [11] S. Paoletti, M. Casini, A. Giannitrapani, A. Facchini, A. Garulli, and A. Vicino, “Load forecasting for active distribution networks,” in Proc. 2nd IEEE PES Innovative Smart Grid Technologies Europe, 2011. [12] A. Bemporad, A. Garulli, S. Paoletti, and A. Vicino, “A boundederror approach to piecewise affine system identification,” IEEE Trans. Autom. Control, vol. 50, no. 10, pp. 1567–1580, 2005.
2400
