COMPUTER-IMPLEMENTED METHOD FOR COMMUNITY ENERGY MANAGEMENT USING MIXED-INTEGER LINEAR PROGRAMMING, MODEL PREDICTIVE CONTROL, NON-INVASIVE LOAD MONITORING TECHNIQUES, ROBUST FORECASTING MODELS AND RELATED SYSTEM

Abstract

The present invention relates to a computer-implemented method and respective system (100) for community energy management in a community comprising energy consumers and energy producers. Such a community energy management method and system (100) allow optimal energy sharing within energy communities. The method of the present invention is based on mixed-Integer linear programming, operating under the receding horizon concept of Model Predictive Control.

Description

CROSS REFERENCE TO RELATED APPLICATIONS

This application claims the benefit of priority of Portuguese Patent Application No. 118847, filed Jul. 31, 2023, the contents of which are all incorporated herein by reference in their entirety.

TECHNICAL FIELD OF THE INVENTION

The present invention relates to the field of renewable energy communities and house energy management systems, specifically using prediction models based on artificial intelligence algorithms.

BACKGROUND OF THE INVENTION

In 2019, to face the challenges presented by climate change, the European Union (EU) presented the European Green Deal strategy, which according to the EU aims to transform EU the first climate neutral continent by 2050 and decouple economic growth from resource use. It is expected that the European Green Deal will transform the EU into a modern, resource-effective, and competitive economy. Within the scope of energy, simultaneously with promoting the use of renewable energy sources, the EU presented a set of initiatives focused on the energy end-users. These initiatives focused on the end energy end-users trend to accelerate the change of power system organization paradigm, previously essentially centralized. Among these initiatives there is the promotion of prosumers (i.e., who both consume and produce energy) and stated that the decentralized production of energy by prosumers must constitute a relevant component of EU's energy policy. The appearance of prosumers has enhanced the appearance of so-called energy communities. According to the EU, Citizen Energy Community “is a voluntarily legal entity established at the local level for the local production of electricity, distribution, supply, consumption, aggregation, energy efficiency services, peer-to-peer energy trading, electric vehicles, and energy storage, etc.” Although with similar definitions, a Renewable Energy Community considers only production from renewable sources, and includes different aspects such as electricity, thermal, gas and water. Other definitions have been proposed for energy communities, such as “a group of buildings with low-carbon technologies to supply, store and internally share/trade electricity and heating”, “a group of consumers and/or prosumers, that together share generation units and electricity storage”, and even “community energy” and “energy community” are found in literature. According to EU estimates, by 2030 energy communities can account for 17% of all installed wind capacity and 21% of all photovoltaic capacity installed in the EU. Thus, it is possible to anticipate the great impact of energy communities in the future. However, it is not only necessary to create energy communities, but more importantly, it is necessary to develop intelligent and sustainable communities, and for this it is necessary to develop efficient energy management systems for energy communities.

PRIOR ART

[1] W. Su, A. Huang, in “The energy internet: An open energy platform to transform legacy power systems into open innovation and global economic engines”, Woodhead Publishing, 2018 proposes an energy management system is developed for energy communities. In fact, two models are developed. In the first the model for house energy management systems is presented, i.e., when houses operate separately, and then a model where houses can exchange energy with each other. In both cases the objective is the minimization of operation costs. In the first case the objective is the minimization of house operation costs and the second the minimization of community costs. The results show that it is possible to reduce community costs when houses operate together, but some houses can see their cost increases in favour of a collective cost reduction goal.

[2] A. Manso-Burgos, D. Ribó-Pérez, T. Gómez-Navarro, M. Alcázar-Ortega, in “Local energy communities modelling and optimisation considering storage, demand configuration and sharing strategies: A case study in Valencia (Spain), Energy Reports.” 8 (2022) 10395-10408 proposes an optimization for local energy communities, with a case study in Spain. Different sharing strategies are considered among community participants through different allocation coefficients, such as static and variable coefficients. The authors consider that static coefficient fixed throughout the year and variable coefficient changing every hour. An analysis is made to assess the impact of increasing battery capacity on system operation.

Patent application CA2809011 A1, “Adaptative energy management system”, priority date Nov. 6, 2012, published May 6, 2014, MacMaster University, CA, discloses systems, methods and devices related to a microgrid system for providing power to a facility. A self-contained power system provides power to a facility using a combination of power storage elements and renewable energy sources. A connection to an external power grid may also be provided. The system optimizes power flow to the facility using power from the storage elements and the renewable energy sources and, if necessary, the external power grid. The optimization process predicts future power usage by the facility using power usage data from a predetermined time window. The optimization process can also take into account predicted energy generation amounts by the renewable energy sources. To optimize economic effects, the optimization process can also determine whether to buy and when to buy power from the external power grid.

Although there is a great interest in the subject of energy communities, the literature fails to present complete solutions, that are able to adapt to communities having houses with completely different situations in terms of consumption and production, or storage of energy. In addition, the discussion on the form of energy sharing within energy communities is in a very initial phase, namely when it comes to the allocation coefficients, which are addressed in this application.

SUMMARY OF THE INVENTION

The present invention refers to a computer-implemented method and respective system for community energy management in a community comprising energy consumers and energy producers. Such a community energy management method and system allow optimal energy sharing within energy communities, as it is a central system that makes the global management of the entire community. The method and system of the present invention is based on mixed-Integer linear programming (MILP), operating under the receding horizon concept of Model Predictive Control (MPC). A systematic classification of electric appliances used in the houses of the community, the use of external information such as weather information, as well as the use of intelligent forecasting techniques enables the present invention to achieve an excellent efficiency. It also allows for an easy installation of as well as a smooth scaling with an increasing number of houses. The community energy management method formulations include sharing of energy without restrictions, as well as employing different allocation coefficients strategies.

In a first aspect, the present invention refers to a computer-implemented method community energy management.

In a second aspect, the present invention relates to a system configured to implement the method of the first aspect.

In a third aspect, the invention refers to a computer program product comprising instructions that are able to carry out the method steps when executed in a computer.

In a fourth aspect, the invention refers to a computer-readable storage medium comprising instructions that are able to carry out the method steps when ran in a computer.

SOLUTION TO TECHNICAL PROBLEM

The present invention addresses the problem of energy management of a community of houses by proposing complete, highly efficient, and user-defined community technical solutions which implement sharing of electric energy produced and/or stored by its members. It is a complete solution since it details all hardware and software needed to implement community energy management, such hardware and software being described in the appended claims of the present application. It allows the community to decide the methodology of energy sharing among a set of pre-defined and implemented techniques, easily extendable to other approaches. Finally, for the chosen sharing technique, it provides a highly efficient solution as: a) MILP provides an optimal solution to the community economic cost, according to the data used, adapting naturally to changes of equipment usage; b) the forecasts necessary in the MILP formulation are among the more accurate available worldwide due to: i) the use of state-of-the-art robust forecasting models; ii) the division of electric appliances into four classes, which significantly reduces the uncertainty of power consumption in each house in the community.

ADVANTAGEOUS EFFECTS OF THE INVENTION

As advantageous effects of the present invention, the following can be listed:

• • A—The proposed invention enables to implement sharing of electric energy within a community, as well as the management of energy transfer with the grid, using several different criteria, whether pre-programmed or user-defined.
  • B—By classifying the types of Electrical Appliance between Monitored (M)/Not Monitored (NM), Controlled (C)/Not Controlled (NC), and Manually-Controlled (MC)/Computer-Controlled (CC), it enables, on one hand, a better forecasting of the power demand of each house within the community, and, on the other hand, allows the user to set suitable schedules for User-Controlled devices (for instance, thermo-accumulators), as well as settings and control algorithms for Computer-Controlled devices (for instance, Air-Conditioners or swimming-pool pumps).
  • C—The operation and power consumption of Monitored devices can be estimated and forecasted using Non-Invasive Load Monitoring (NILM) techniques.
  • D—Power consumption of Not Monitored devices, as well as power produced by the (photovoltaic apparatus) PV installations, must be forecasted within the Prediction Horizon (PH). Typically, with, for example, a time-step of 15 minutes, if the PH is 24 hours, this means that the forecasts must be available for 96 steps-ahead. To speed-up the calculations, this period is divided into two sub-periods, as an example: for the last 15 hours, naïve forecasting is used; for the first 9 hours, the forecasts are obtained using Robust Forecasting models, which are an ensemble of models, designed not only to minimize the forecasting errors within PH, but also ensuring Prediction Intervals.
  • E—C) and D) enable the power demand of each house, as well as its forecast within PH, to be divided between NM demand and M demand. As the larger uncertainty is related with NM demand, this decreases the total forecasting error.
  • F—The control of all CC devices, including inverters and/on batteries for all houses are obtained by solving, at each time-step, a MILP problem, formulated with the selected energy shared criterion, using the NM and M demand and PV forecasts, over PH, and energy selling/buying grid tariffs.
  • G—As at each time-step, a new Mixed-Integer Linear Programming (MILP) problem is solved, in a Receding Horizon Scheme, the energy management scheme adapts to new situations (for instance the initial operation of a dishwasher, which is not controlled, but monitored using NILM.
  • H—All this computation is done in an aggregator computational unit (a software executing in a computer), which can reside in one of the houses, or in a separate energy management company. Data is received/sent wirelessly between each house and the aggregator computational unit, using data acquisition systems in each house.

BRIEF DESCRIPTION OF THE FIGURES

Preferred embodiments of the present invention will be described with reference to the accompanying figures, which are to be construed as non-limiting the invention scope, which scope is defined by the appended claims, and wherein:

FIGS. 1A and 1B—illustrate a flowchart representing the method steps according to an embodiment of the present invention;

FIG. 2—illustrates a schematic diagram representing the system according to an embodiment of the present invention;

FIG. 3—illustrates schematically an aggregator computational unit according to an embodiment of the present invention;

FIG. 4—illustrates a schematic diagram representing an embodiment of a part of the system according to the present invention;

FIG. 5—illustrates a hybrid forecasting scheme used in embodiments of the present invention;

FIG. 6—illustrates a flowchart of a ApproxHull used in embodiments of the present invention;

FIG. 7—illustrates a chromosome and input space lookup table of a neural network used in embodiments of the present invention;

FIG. 8—illustrates a MOGA flow of operation of a neural network used in embodiments of the present invention;

FIG. 9—illustrates a model design cycle used in embodiments of the present invention;

FIG. 10—illustrates a chart plotting different strategies for identifying the best out of 10 training trials of a neural network used in embodiments of the present invention; and

FIG. 11—illustrates a chart Detail of objective space for a problem of a model used in embodiments of the present invention.

DETAILED DESCRIPTION OF THE INVENTION

The present invention relates to a computer-implemented method for community energy management and to a system configured to implement said method.

In the following detailed description, the method steps will be described, (method steps are indicated between brackets below) together with the system object of the present patent application. Also the figures annexed will be referred to, for better clarification.

1—Data d₀ acquisition (step a)), by an aggregator computational unit (110), during a minimum period of time of 15 days, typically during one month, wherein data do comprises at least one of the following data:

• • {tws}—Time basis for weather data;
  • {R}—Global Solar Radiation;
  • {T}—Atmospheric Air Temperature

For every house h:

• • {t_Cons}—Time basis for Consumption data;
  • {P_h}-Active power of house h
  • {R_h}—Reactive power of house h
  • {t_Inv}-Time basis for Inverter data;
  • {P_h^PV}—AC Power generated by the PV installation of house h
  • {O_h}—Daily Occupation of house h

For every house h that has Manually Controlled Appliances:

• • StP_h,j/EndP_h,j, jε{MC_h}—Start/end times of the operation of all Manually Controlled (MC) appliances of house h.
  • {P_h,j, J∈{MC_h}}—Active Power of all Manually Controlled appliances of house h.
  • {t_MCh}—Time basis for all Manually Controlled appliances of house h.

For every house h that has Computer Controlled Appliances (not Acs):

• • {U_h,j, j∈{CC_h}} Control Actions for all appliances of house h to be Computer Controlled (CC)
  • {P_h,j, j∈{CC}}—Active Power of all appliances of house h to be Computer Controlled
  • {t_cCh}—Time basis for all Computer Controlled appliances of house h.

For every house h that has Computer Controlled ACs:

• • {IT_h^room, room∈{AC_h}}—Air temperature for all ACs of house h to be Computer Controlled
  • {RHP_h^room, room∈{AC_h}} Relative humidity for all ACs of house h to be Computer Controlled
  • {WT_h^room, room∈{AC_h}} Wall temperature for all ACs of house h to be Computer Controlled
  • {P_h^room, room∈{AC_h}} Active Power for all ACs of house h to be Computer Controlled
  • {U_h^room, room∈{AC_h}} Control Actions for all ACs of house h to be Computer Controlled
  • {t_ACh}—Time basis for all ACs of house h to be Computer Controlled

For every house h that has Monitored but not Controlled (M) appliances:

• • StP_h,j/End_h,j, jε{M}—Start/end times of the operation of all Monitored (M) appliances of house h.
  • {P_h,j, j∈{M_h}}—Active Power of all Monitored appliances of house h.
  • {t_Mh}—Time basis for all Monitored appliances of house h.2—Model Design Procedure (step b))

The following procedure is employed for model design:

• • 1. Selection of the datasets from the available design data by using, for example, ApproxHull [3];
  • 2. Minimization of the different objectives, by formulating a Multi-Objective Genetic Algorithm (MOGA);
  • 3. Analysis of MOGA results, restricting some of the objectives and redefinition of the input space and/or the model topology;
  • 4. Usage of the preferential set of models obtained from the second MOGA execution;
  • 5. Determination of m_e^pref using (1);
  • 6. From m_e^pref obtain m_e^par using (2);
  • 7. If the model that should be designed is a Nonlinear Auto-Regressive model with Exogeneous (NARX) model, usage m_e^par for the exogenous variables;
  • 8. Obtaining the prediction intervals using (3)

The set m_e^pref is defined as:

$\begin{array}{cc}{m}_e^{\mathrm{pref}}=\left\{m_{\mathrm{pref}}:\mathrm{nw}\left(m_{\mathrm{pref}}\right)\le {\stackrel{⌣}{m}}_{\mathrm{nw}}^{\mathrm{pref}}\right\}\bigcap \left\{m_{\mathrm{pref}}:\mathrm{fore}\left(m_{\mathrm{pref}}\right)\le {\stackrel{⌣}{m}}_{\mathrm{fore}}^{\mathrm{pref}}\right\}& \left(1\right)\end{array}$

where m_pref is the set of preferred models obtained in the 2^nd MOGA execution, m̆_nw^pref is the median of the linear weight norm for each model in m_pref and m̆_fore^pref the median of the sum of the forecasts for the prediction time-series obtained by models m_pref.

The third subset of models, m_e^par, is obtained iteratively, by adding non-dominated models, considering nw(.) and fore(.) as the two criteria, until NS models are found. Initially, the set of models to be inspected, m_i is initialized to m_e^pref, and m_e^par to an empty set. In each iteration, both criteria are applied to m_i. Then, the non-dominated solutions found are added to m_e^par and removed from m_i.

$\begin{array}{cc}{m}_e^{\mathrm{par}}=\left\{\mathrm{nd}\left(\mathrm{fore}\left(m_e^{\mathrm{pref}}\right),\mathrm{nw}\left(m_e^{\mathrm{pref}}\right)\right):\#\mathrm{nd}=\mathrm{NS}\right\}& \left(2\right)\end{array}$

The prediction intervals are obtained as:

$\begin{array}{cc}{\underset{\_}{\stackrel{⌣}{y}}}_k={\stackrel{⌣}{y}}_k-\Delta {\stackrel{⌣}{y}}_k\le y_k\le {\stackrel{⌣}{y}}_k+\Delta {\stackrel{⌣}{y}}_k={\stackrel{\_}{\stackrel{⌣}{y}}}_k& \left(3\right)\end{array}$ $\Delta {\stackrel{⌣}{y}}_k=t_{1-\alpha /2,N-p}{\sigma }_{\mathrm{tot}}$

Where t_1-α2,N-p represents the α/2 quantile of a Student's t-distribution with N-p degrees of freedom. The total prediction variance for each prediction can be obtained as:

$\begin{array}{cc}{{\sigma }_{\mathrm{tot}}^2={\sigma }_{\epsilon }^2(1+{\phi }^T(x_k,v\text{*)}({\Gamma }^T(Z_{\mathrm{tr}},v\text{*)}\Gamma (Z_{\mathrm{tr}},v\text{*)})}^{-1}\phi (x_k,v\text{*)})& \left(4\right)\end{array}$

where v* denotes the optimal value of the nonlinear parameters and the data noise variance, σ_ε², can be estimated as:

$\begin{array}{cc}{\hat{\sigma }}_{\epsilon }^2=\frac{\sum _{k=1}^N{e}_k^2}{N-p}=\frac{\sum _{k=1}^N{\left(y_k-\stackrel{⌣}{y}\left(z_k\right)\right)}^2}{N-p},& \left(5\right)\end{array}$

Where N denotes the total number of samples and p the number of model parameters.

For classification models no forecasting criterion such as (62) is employed. In this case, criteria such as FP_tr, FN_tr, FP_te, FN_te are employed, where FP and FN denote the number of False Positives and False Negatives, respectively. Additionally, m_e^par will be the set with the smallest NS linear weight norm.

3—Determination of power profiles (step c)) for all monitored, manually controlled and computer-controlled appliances by computing the mean value of power consumption during all sampling intervals within a period of time, for example, within a day.4—Initializations (step d))

Before executing the system, the information described in the heading Parameters, below, must be supplied. For each battery bank, the initial SOC value must be supplied.

Additionally, as several data-based models are used, they need to be designed first. The ANN models described in the heading Models, below, need around 1 month of data to be designed, using the methods described before. This way, this amount of data needs to be acquired. Additionally, for each appliance that will be manually controlled, its consumption profile needs to be determined, which also requires data. All this data can be measured by Wi-fi Smart Plugs devices, available commercially, connected to the data acquisition of each house.

Parameters

For the community:

Energy sharing method, with possible coefficient values

For all Houses with PV:

• • P_h^Grid−/P_h^Grid+ Contracted bought(−)/sold(+) power from/to the grid by house hFor all houses with PV and Battery Bank:
  • P_h^B−/P_h^B+ Maximum discharged (charged) power of the battery bank of house h
  • SoC_h/SoC_h Minimum/Maximum value of the state of charge of the battery bank of house h
  • SoC_h,initial Initial state of charge of the battery bank of house h
  • RU_h/RD_h Maximum values of the ramp-up/ramp down rates for charging or discharging of the battery bank of house h
  • η_h^B−/η_h^B+ Charging/discharging efficiency of the battery bank of house h
  • Ē_h^B Rated capacity of battery bank of house hFor All Houses with Computer-Controlled Acs (Air Conditioning Systems):
  • t_s Start time when a room must be in thermal comfort
  • t_e End time when a room must be in thermal comfort
  • Δt Time before the start of the schedule when the AC might need to operate to ensure thermal comfort throughout the whole period
  • Θ_T PMV (Predicted Mean Vote) threshold, typically assuming the value of 0.5

Models

• • M^T Forecasting Temperature Model
  • M_h^P Forecasting Consumption Model for NMC appliances of house h
  • M^R Forecasting Radiation Model
  • M_h^PV Forecasting PV Power Model for house h
  • M_h,j^C Classifier model for monitored appliance j of house h5—Data acquisition (step e))

A data acquisition system was developed and installed in a house.

In the example of implementation of the present invention, the consumption of the house is measured by a Carlo Gavazzi (EM340) 3 phase imported energy meter. This meter is a certificated device, and electrical measurement is done using a 2-wires Modbus RTU connection. EM340 supplies 45 different electric variables, sampled at 1 Hz. Variables related with the energy produced by a photovoltaic apparatus (PV), stored in a battery and injected in the grid are available from a Kostal smart energy meter (SEM). In total, 21 variables are obtained by KI, at a sampling interval of 1 min. The data access is done using a cable IP network using the Modbus IP protocol. Smart Plugs (SP) are used for on/off control of some equipment. Additionally, they allow sockets belonging to the same CB to be measured individually. This type of devices connects to an existing wireless network using an initial access point and a manufacturer mobile app. They can be read/controlled using a cloud API or directly using an internal web service, which is the case here. A smart Plug (180) is represented in FIG. 4.

An Intelligent Weather Station (IWS) (120) represented in FIG. 2 is a device that measures the air temperature and relative humidity, and global solar radiation, and predicts their evolution within a self-defined prediction horizon. To enable these forecasts, a two stages strategy is employed. When little measured data are available, a nearest neighbor algorithm is employed. When enough data are available, neural networks predictive models are automatically off-line designed and uploaded to the IWS, for real-time use.

Finally, Self-Powered Wireless Sensors (SPWS) (150), represented in FIG. 4, are used for measuring room data, such as air temperature and relative humidity, status (open/close) of doors and windows, walls temperature, light and room movement. They are Ultra-Low-Power devices and communicate via ISM radio band working on 2.4 GHz or 868 MHz frequencies. They can be used to measure thermal comfort for usage in predictive control of air-conditioners in the house.

Gateways and a Technical Network are tasked with data transmission from/to the measurement devices. A technical IP-cabled and a wireless network have been created using a network router located inside an extension of the electric board case. This router separates the home network from the technical network. All devices, except the SPs and SPWS, are connected to that network. To perform the data acquisition from the existing devices EasyGateway devices are used. EasyGateway is a fault-tolerant IoT gateway that supports a variety of reception/acquisition protocols such as Modbus, SNMP, Easymodules and serial http, as well as a set of DataDelivered Connectors (DDC). The data transmission of the EasyGateways is always performed at a 1-min base, which means that the data acquired by the measurement devices related to each gateway are packed and transmitted at that rate. Control of electrical equipment is performed via Modbus.

The gateways can also send data to up to three different DDCs. Here, 2 DDCs and 2 IoT Platforms are used. One platform is used inside home, employing an EasyMqs DDC and another, using a Generic Ampq DDC, is in the cloud. Depending on where the aggregator resides, the former or the latter can be used. The same IoT platform is used inside home and in the cloud. It receives data from the configured message queue servers. The data arrived pass by a set of configured plug-ins for each type of entities configured on the platform. The application provides a web page where the end-user can configure a set of definitions related to data storage and management. It also allows data visualization using plots grouped by sensors category, and data download in 4 standard formats, csv, xlsx, mat and npz. For a more detailed description of the IoT platform, see [4].

6—Control Actions (step f) determined in the previous step, by the aggregator computational unit (110) sending to each data acquisition subsystem (140), which, on its turn, sends the respective control actions to the computer-controlled devices in each house of the community, via modbus;7—Data reception by an Intelligent Aggregator (steps g), h), i), i)

The concept of the intelligent aggregator for community energy management system is represented in FIGS. 2 and 3.

The aggregator computational unit (110) is a computational unit that receives information about:

• • weather parameters—Air Temperature (T) and Solar Radiation (R)—from possibly a weather station, such as described above or collected by an internet application that, on one hand, is valid for all houses in the community and, on the other hand, measures or interpolates these values at the sampling interval considered (15 minutes);
  • Tariffs employed in the energy prices, both from the point of view of purchasing as well as selling, for the current sample and for all future samples within a Prediction Horizon (PH), for example of 24 hours;
  • Daily occupation (number of occupants) for each house (current and next day), inserted by a user:
  • Electric information related with the electric consumption of each house, sent by the consumption energy meter:
    • P_h,k—Active power of house h at instant k;
    • R_h,k—Reactive power of house h at instant k
  • information regarding possible PV installation, sent by the inverter energy meter:
    • {circumflex over (P)}_h,k^PV—AC Power generated by the PV installation of house h at instant I
  • possible energy storage, sent by the inverter energy meter:
    • SoC_h,k—State of charge of the battery bank of house h at instant k
  • And possible climate information (Air temperature (IT_h,k^room), Relative humidity (RH_h,k^room), Wall temperature (WT_h,k^room) about each room in house h that is computer-controlled using Predicted Mean Vote (PMV), at instant k.

This information is sent wirelessly from the data acquisition of house h to the aggregator every sampling interval (for example, 15 minutes). Additionally, if there are appliances that should be monitored in house h, Active and Reactive Power sampled at, for example, 1 minute interval are also sent to the aggregator. If there are appliances that should be computer controlled in house h, then its data acquisition system receives from the aggregator, at each sampling interval, the corresponding actuation signals (in the case of battery control, the actions of charging the battery from the PV/discharging the battery; in the case a swimming pool pump, the on/off actuation signal; in the case of PMV control of air conditioners the on/off actuation signal and the reference temperature). The data acquisition system should supply these actuation signals to the corresponding device. Different acquisition systems can be used, provided they satisfy the above specifications. A data acquisition system according to the present invention is represented in FIG. 4.

8—Determination of the averaged active power consumption of each house (step k), by using:

$P_{h,k}=\frac{P_{h,k_1}}{N},N=\frac{\Delta k}{k_1}$

9—Forecasts (step I))

To implement the MILP-MPC approach, there is the need for forecasting both the consumption of each individual household and the PV generation of each PV system over the PH (obviously the community versions are just a sum of the values for the individual houses). Their forecasting, unless specified, is obtained as the joint usage of naïve forecasting and the use of multi-step NARX forecasting models, as FIG. 5 illustrates.

9.1—Electric consumption (substep I1) (EXAMPLE)

We have divided the appliances in four types, mainly because of the way that their forecasts are obtained. The forecast of the active Power of House h at instant k, within the Prediction Horizon, can be given as:

$\begin{array}{cc}{\hat{P}}_{h,\left\{k+1...k+\mathrm{PH}\right\}}={\hat{P}}_{h,\left\{k+1...k+\mathrm{PH}\right\}}^{\mathrm{NMC}}+{\hat{P}}_{h,\left\{k+1...k+\mathrm{PH}\right\}}^{\mathrm{NC}}+{\hat{P}}_{h,\left\{k+1...k+\mathrm{PH}\right\}}^{\mathrm{MC}}+{\hat{P}}_{h,\left\{k+1...k+\mathrm{PH}\right\}}^{\mathrm{CC}}& \left(6\right)\end{array}$

Forecasting of House h not Monitored Neither Controlled Consumption Within the Prediction Horizon

$\begin{array}{cc}{\hat{P}}_{h,\left\{k+1...k+\mathrm{PH}\right\}}^{\mathrm{NMC}}=\sum _{j=1}^{\#\left(\left\{{\mathrm{NM}}_i\right\}\bigcup \left\{{\mathrm{NC}}_i\right\}\right)}{\hat{P}}_{h,j,\left\{k+1...k+\mathrm{PH}\right\}}& \left(7\right)\end{array}$ $\begin{array}{cc}{\hat{P}}_{h,\left\{k+1...k+\mathrm{PH}\right\}}^{\mathrm{NMC}}={\hat{P}}_{h,\left\{k+1...k+36\right\}}^{\mathrm{NMC}}\bigcup {\hat{P}}_{h,\left\{k+37...k+\mathrm{PH}\right\}}^{\mathrm{NMC}}& \left(8\right)\end{array}$

Equation (7) states that the not monitored neither controlled consumption is the sum of the active power of corresponding appliances. The naïve forecasting is obtained as:

$\begin{array}{cc}{\hat{P}}_{h,\left\{k+37...k+\mathrm{PH}\right\}}^{\mathrm{NMC}}={\hat{P}}_{h,\left\{k-59...k\right\}}^{\mathrm{NMC}}& \left(9\right)\end{array}$

The first 9 hours of forecasting are obtained by a RBF NARX model:

$\begin{array}{cc}{\hat{P}}_{h,\left\{k+1...k+36\right\}}^{\mathrm{NMC}}=\left\{M_h^P\left({\stackrel{\_}{P}}_{h,l}^{\mathrm{NMC}},{\stackrel{\_}{T}}_l,O_{h,l},D_l\right),l=k+1\dots k+36\right\}& \left(10\right)\end{array}$

Equation (8) states that the forecasting of the first 36 steps for the NMC consumption of house h are obtained as a multi-step forecasting of model M_h^P, which uses past values (relative to I) of that consumption (P_h,l^NMC), past values of the atmospheric air temperature (T_l), the current occupation of the house, and a codification of the corresponding day (D_l). Notice that it is assumed that the community houses are near enough so that we can consider there is no spatial variation of the atmospheric temperature.

TABLE 1
Day of the week and holiday occurrence encoding values
	Day of the week	Regular day	Holiday	Special
	Monday	0.05	0.40	0.70
	Tuesday	0.10	0.80
	Wednesday	0.15	0.50
	Thursday	0.20	1.00
	Friday	0.25	0.60	0.90
	Saturday	0.30	0.30
	Sunday	0.35	0.35

The encoding, presented in table 1, distinguishes the day of the week and the occurrence and severity of holidays based on the day they occur. The regular day column shows the coding for the days of the week when these are not a holiday. The following column presents the values encoded when there is a holiday, and finally, the special column shows the values that substitute the regular day value in two special cases: for Mondays when Tuesday is a holiday; and, for Fridays when Thursday is a holiday.

As the forecasts of the NMC consumption depend on the forecasts of the atmospheric temperature, another model, in this case a NAR model must be used to obtain these forecasts.

$\begin{array}{cc}{T}_{\left\{k+1...k+36\right\}}=\left\{M^T\left({\stackrel{\_}{T}}_l\right),l=k+1\dots k+36\right\}& \left(11\right)\end{array}$

Forecasting of House h monitored but not controlled consumption within the Prediction Horizon

$\begin{array}{cc}{\hat{P}}_{h,\left\{k+1...k+\mathrm{PH}\right\}}^{\mathrm{NC}}=\sum _{j=1}^{\#\left(\left\{M_i\right\}\bigcap \left\{{\mathrm{NC}}_i\right\}\right)}{\hat{P}}_{h,j,\left\{k+1...k+\mathrm{PH}\right\}}& \left(12\right)\end{array}$

The forecasting can be obtained using Algorithm (Alg.) 1:

Alg. 1-Forecasting of House h monitored but not
controlled consumption of appliance j within the
Prediction Horizon
# Inputs
# Ph,{   ...k},,Rh,{   ...k},,OPh,{   ...k-1},,{circumflex over (P)}h,{   ...k-1}
# Outputs
# OPh,j,{k...k},{circumflex over (P)}h,j,{k−1...k},{circumflex over (P)}h,j,{k+1...k+PH}
StPh,j = 0
EndPh,j = 0
for i=1:15
k1 =(k−1)*15+i
[OPh,j,k1 ,{circumflex over (P)}h,j,k1] = NILM (h, j, Ph,{k1−10...k1−1}, Rh,{k1−10...k1−1},
{circumflex over (P)}h,j,{k1−10...k1−1})
if OPh,j,k1 =1
if OPh,j,k1−1 ≠ 1
StPh,j = k1
end
else
if OPh,j,k1−1 =1
EndPh,j = k1
end
end
end
{circumflex over (P)}h,j,{k+1...k+PH} = UpdateProfileMNC (i, j, StPh,j , EndPh,j , OPh,j,{k...k}
{circumflex over (P)}h,j,{k...k} , k)

In Alg. 1, OP_h,j,k denotes the Operation Status (0—Off; 1—On) of appliance j at house h at instant k. R_h,j,k is the reactive power of the main meter of house h at instant k. StP_h,j/EndP_h,j denote the start/end sample index of the operation of appliance j of house h.

It should be noted that NILM operates with data sampled at 1 minute, not with a 15 minutes sampling interval. The symbol k₁ denotes sample indexes with a sampling interval of 1 min. The relation between k₁ and k is k→k₁∈{(k−1)*15+1. . .k *15}. The first index of k₁ is denoted as k while the last one is k.

At each minute, the NILM operation for appliance j at house h is executed. If the appliance is operating (OP_h,j,k=1), then the corresponding power estimation is performed ({circumflex over (P)}_h,j,k). The algorithm then determines if the start of the operation and the end of the operation occur during this period. Then the profile of the appliance is updated.

Alg. 2-NILM Operation and Power Estimation of
appliance j in House h at instant k
# Inputs:
# h,j,Ph,{k1−10...k1−1},Rh,{k1−10...k1−1},{circumflex over (P)}h,j,{k1−10...k1−1}
# Outputs
# OP,{circumflex over (P)}
{circumflex over (P)} = 0
OP = Mh,jC(Ph,{k1−10...k1−1},Rh,{k1−10...k1−1})
if OP = 1
{circumflex over (P)} = Mh,jE(Ph,{k1−10...k1−1},Rh,{k1−10...k1−1},{circumflex over (P)}h,j,{k1−10...k1−1})
end

For each appliance there are two models that intervene in the NILM operation: first a classifier model (M_i,j^C) which outputs 1 if the operation is working, and 0 if not. If it is working, then another estimation model is executed, to estimate the active power consumed.

Alg. 3-Updates the Profile of the jth MNC appliance of house h
# Inputs
# i, j, StPh,j , EndPh,j , OPh,j,{1...15} {circumflex over (P)}h,j,{1...15} , k
# Outputs
# {circumflex over (P)}h,j,{k+1...k+PH}
k1 = (k−1)*15+1
[ StPbef , EndPbef , OP , P, P2]=LoadProfileMNC(h, j)
{circumflex over (P)}h,j,{k+1...k+PH} = 0
If ∃(OPh,j,{k...k}) ≠ 0
Beg = 1
BegP = StPh,j−k1+1
StP = StPh,j
If StPh,j = 0
Beg = k1 − StPbef
BegP=1;
StP = StPbef
end
LastP = 15
EndP = EndPbef
If EndPh,j ≠ 0
Last = EndPh,j −StPbef
LastP= EndPh,j −k1
EndP=EndPh,j
end
OP{Beg...Last} = OP{Beg...Last} + OPh,j,{BetP...LastP}
P{Beg...Last} = P{Beg...Last} + Ph,j,{BetP...LastP}
P2{Beg...Last} = P2{Beg...Last} + Ph,j,{BetP...LastP}·{circumflex over ( )}2
SaveProfileMNC(h, j, StP , EndP , OP , P, P2)
{circumflex over (P)}h,j,{k+1...k+m} ={mean(P{j...j}./ OP{j...j}),l=1...m}
end

For each appliance a profile is stored, with the start and end of the last (or current) operation of appliance j at house h, a vector (OP) describing the number of times where the appliance was on since its operation started, for each sample, a vector (P) storing the sum of the powers for each sample, and a vector (P2) storing the sum of the square of powers for each sample. All these variables assume a sampling interval of 1 sec.

Then, the mean power is simply obtained as P./O, where the symbol “.I” means division element-by-element. If necessary, the variance of the power can be obtained as P2./O-(P./O){circumflex over ( )}2.

The forecasting of the power over the prediction horizon, with a sampling interval of 15 minutes, is just obtained as the mean of the mean powers over each period of 15 minutes.

Forecasting of House h Manually Controlled Consumption Within the Prediction Horizon

${\hat{P}}_{h,\left\{k+1...k+\mathrm{PH}\right\}}^{\mathrm{MC}}=\sum _{m=1}^{\#\left(\left\{{\mathrm{MC}}_i\right\}\right)}{\hat{P}}_{h,m,\left\{k+1...k+\mathrm{PH}\right\}}$

Typically, for this kind of appliances, the power is obtained from a vector of daily mean power values, which are obtained offline through measurements. This way, for each appliance, we use:

${\hat{P}}_{h,m,\left\{k+1...k+\mathrm{PH}\right\}}=\mathrm{RetrieveProfileMC}\left(i,j,k\right)$

Alg. 4-Obtain the forecasting of the jth MC appliance
of house h
# Inputs:
# i, j, k
# Outputs
# {circumflex over (P)}h,j,{k+1...k+PH}
[P]=LoadProfileMC(h, j)
kbeg = dayindex(k)
{circumflex over (P)}h,j,{k+1 k+PH} = [Pkbeg+1...96 P1...kbeg]

Forecasting of House h Computer-Controlled Consumption Within the Prediction Horizon

${\hat{P}}_{h,\left\{k+1...k+\mathrm{PH}\right\}}^{\mathrm{CC}}=\sum _{n=1}^{\#\left(\left\{{\mathrm{CC}}_i\right\}\right)}{\hat{P}}_{h,n,\left\{k+1...k+\mathrm{PH}\right\}}$

For this kind of appliances, the power can be obtained as for MC, or using other formulations, which are appliance specific.

${\hat{P}}_{h,m,\left\{k+1...k+\mathrm{PH}\right\}}=\mathrm{RetrieveProfileCC}\left(h,m,k\right)$

9.2—Forecasting of PV Power generated (substep I2) (EXAMPLE)

Here we use also hybrid forecasting. The first 36 steps-ahead forecasting are obtained as:

$\begin{array}{cc}{\mathrm{PV}}_{h,\left\{k+1...k+36\right\}}=\left\{M_h^{\mathrm{PV}}\left(\stackrel{\_}{P}{\stackrel{\_}{V}}_{h,l},{\stackrel{\_}{T}}_l,{\stackrel{\_}{R}}_l\right),l=k+1\dots k+36\right\}& \left(13\right)\end{array}$

Where R denotes global solar radiation. This way, we need an additional forecasting model to determine the forecasts of solar radiation over PH.

$\begin{array}{cc}{R}_{\left\{k+1...k+36\right\}}=\left\{M^R\left({\stackrel{\_}{R}}_l\right),l=k+1\dots k+36\right\}& \left(14\right)\end{array}$

10—MILP-MPC formulation (step m)

After having described a data acquisition system (140) according to the invention, it will now be described the software part that resides in the aggregator computational unit (110).

At every sampling instant k, the intelligent aggregator solves a Mixed-Integer Linear Programming (MILP) problem, formulated as:

$\begin{array}{cc}\min {\hat{C}}^{\mathrm{Com}}=\min \sum _{h=1}^{\mathrm{Nh}}\sum _{s=1}^{\mathrm{PH}}\left({\hat{P}}_{h,l}^{\mathrm{Grid}-}+{\hat{P}}_{h,l}^{\mathrm{Com}-}\right){\pi }_l^{\mathrm{Grid}-}\Delta k-\left({\hat{P}}_{h,l}^{\mathrm{Grid}+}+{\hat{P}}_{h,l}^{\mathrm{Com}+}\right){\pi }_l^{\mathrm{Grid}+}\Delta k,l=k+s& \left(15\right)\end{array}$

In (15), Nh denotes the number of houses in the community, PH is the prediction horizon, and Δk is the sampling interval (in hours). Typical values for these two last parameters are 24 hours (96 steps) and ¼ hours, respectively. {circumflex over (P)}_h,l^Grid−/{circumflex over (P)}_h,l^Grid− are the estimates (in Kw) of the active power bought(−)/sold(+) from/to the grid at the sampling instant/by house h. In the same way, {circumflex over (P)}_h,l^Com−/{circumflex over (P)}_h,l^Com+ are the estimates (in Kw) of the active power bought(−)/sold(+) from/to the community at the sampling instant/by house h. Finally, π_l^Grid−/π_l^Grid+ are the price, per kWh of the energy bought(−)/sold(+) from/to the community or the grid at the sampling instant I. This way, the estimate of the cost of the energy bought from the grid by the community during PH is minimized, while the estimate of the profit of the energy sold to the grid by the community during PH is maximized.

Notice that the receding horizon principle is used here. The solution of (15) is a sequence of control actions (the binary values described later) over the whole steps-ahead within PH. But only the first (one-step ahead) of these actions are applied at that instant of time. In the next instant of time the same process takes place.

There are several constraints that must be satisfied:

$\begin{array}{cc}0\le {\hat{P}}_{h,l}^{\mathrm{Grid}+}\le {\stackrel{\_}{P}}_h^{\mathrm{Grid}+}\left(1-v_{h,l}^{\mathrm{Grid}}\right)& \left(16\right)\end{array}$ $\begin{array}{cc}0\le {\hat{P}}_{h,l}^{\mathrm{Grid}-}\le {\stackrel{\_}{P}}_h^{\mathrm{Grid}-}{v}_{h,l}^{\mathrm{Grid}}& \left(17\right)\end{array}$

In (16-17) it is defined that {circumflex over (P)}_h,l^Grid+/{circumflex over (P)}_h,l^Grid− are positive variables with maximum value equal to the contracted power from the grid {circumflex over (P)}_h^Grid+/{circumflex over (P)}_h^Grid+ by house h. v_h,l^Grid is the binary variable modelling the processes of selling/purchasing energy from house h to/from the grid, assuming value 1 if purchasing or 0 if selling.

$\begin{array}{cc}0\le {\hat{P}}_{h,l}^{\mathrm{Com}+}\le {\stackrel{\_}{P}}_h^{\mathrm{Com}+}\left(1-v_{h,l}^{\mathrm{Com}}\right)& \left(18\right)\end{array}$ $\begin{array}{cc}0\le {\hat{P}}_{h,l}^{\mathrm{Com}-}\le {\stackrel{\_}{P}}_h^{\mathrm{Com}-}{v}_{h,l}^{\mathrm{Com}}& \left(19\right)\end{array}$

In the same way, {circumflex over (P)}_h,l^Com+/{circumflex over (P)}_Com− are positive variables with maximum value equal to the maximum value exchangeable with the community P_h^Com+/P_h^Com− by house h. v_h,l^Com is the binary variable modelling the processes of selling/purchasing energy from house h to/from the community, assuming value 1 if selling or 0 if purchasing.

The energy balance can be given as:

$\begin{array}{cc}\sum _{h=1}^{\mathrm{Nh}}{\hat{P}}_{k,l}^{\mathrm{Com}-}=\sum _{h=1}^{\mathrm{Nh}}{\hat{P}}_{k,l}^{\mathrm{Com}+}& \left(20\right)\end{array}$

It is not allowed to buy energy from the grid and sell energy to the community at the same time:

$\begin{array}{cc}{v}_{h,l}^{\mathrm{Grid}}+v_{h,l}^{\mathrm{Com}}\le 1& \left(21\right)\end{array}$

10.1−Energy balance

Without Restrictions on the Transfer of Energy Within the Community

The energy balance of house h at instant k is:

$\begin{array}{cc}{\hat{P}}_{h,l}^{\mathrm{PV}}+{\hat{P}}_{h,l}^{B+}+{\hat{P}}_{h,l}^{\mathrm{Grid}-}+{\hat{P}}_{h,l}^{\mathrm{Com}-}={\hat{P}}_{h,l}^{B-}+{\hat{P}}_{h,l}^{\mathrm{Grid}+}+{\hat{P}}_{h,l}^{\mathrm{NMC}}+{\hat{P}}_{h,l}^{\mathrm{NC}}+{\hat{P}}_{h,l}^{\mathrm{MC}}+{\hat{P}}_{h,l}^{\mathrm{CC}}+{\hat{P}}_{h,l}^{\mathrm{Com}+}& \left(22\right)\end{array}$

In (22), {circumflex over (P)}_h,l^PV is the estimate of the AC Power generated by the PV installation of house h at instant I.

The active power of house h, at instant k is composed of four terms:

• • a) related to devices (NMC) that are not monitored (NM) neither controlled (NC).
  • b) related to devices that are monitored (M) but not controlled (NC).
  • c) related to devices that are manually controlled (MC).
  • d) related to devices that are computer controlled (CC).

$\begin{array}{cc}{P}_{h,k}=\sum _{j=1}^{\#\left(\left\{{\mathrm{NM}}_i\right\}\bigcup \left\{{\mathrm{NC}}_i\right\}\right)}{P}_{h,j,k}+\sum _{l=1}^{\#\left(\left\{M_i\right\}\bigcap \left\{{\mathrm{NC}}_i\right\}\right)}{P}_{h,l,k}+\sum _{m=1}^{\#\left(\left\{{\mathrm{MC}}_i\right\}\right)}{P}_{h,m,k}+\sum _{n=1}^{\#\left(\left\{{\mathrm{CC}}_i\right\}\right)}{P}_{h,n,k}=P_{h,k}^{\mathrm{NMC}}+P_{h,k}^{\mathrm{NC}}+P_{h,k}^{\mathrm{MC}}+P_{h,k}^{\mathrm{CC}}& \left(23\right)\end{array}$

Equation (22) is valid when it is possible to transfer energy between houses, without any restrictions. Several alternatives exist:

When There is No Transfer of Energy Between the Houses

$\begin{array}{cc}{\hat{P}}_{h,l}^{\mathrm{Com}-}={\hat{P}}_{h,l}^{\mathrm{Com}+}=0\Rightarrow & \left(24\right)\end{array}$ $\Rightarrow {\hat{P}}_{h,l}^{\mathrm{PV}}+{\hat{P}}_{h,l}^{B+}+{\hat{P}}_{h,l}^{\mathrm{Grid}-}={\hat{P}}_{h,l}^{B-}+{\hat{P}}_{h,l}^{\mathrm{Grid}+}+{\hat{P}}_{h,l}^{\mathrm{NMC}}+{\hat{P}}_{h,l}^{\mathrm{NC}}+{\hat{P}}_{h,l}^{\mathrm{MC}}+{\hat{P}}_{h,l}^{\mathrm{CC}}$

In this case it is still possible to use the MILP-MPC formulation for each house individually.

Energy is Shared Between Houses, Using Energy Sharing Coefficients

$\begin{array}{cc}{\hat{P}}_{h,l}^{\mathrm{Com}+}={\hat{P}}_{h,l}^{\mathrm{PV}}\Rightarrow & \left(25\right)\end{array}$ $\Rightarrow {\hat{P}}_{h,l}^{B+}+{\hat{P}}_{h,l}^{\mathrm{Grid}-}+{\hat{P}}_{h,l}^{\mathrm{Com}-}={\hat{P}}_{h,l}^{B-}+{\hat{P}}_{h,l}^{\mathrm{Grid}+}+{\hat{P}}_{h,l}^{\mathrm{NMC}}+{\hat{P}}_{h,l}^{\mathrm{NC}}+{\hat{P}}_{h,l}^{\mathrm{MC}}+{\hat{P}}_{h,l}^{\mathrm{CC}}$ $\begin{array}{cc}{\hat{P}}_{h,l}^{\mathrm{Com}-}={\beta }_{h,l}\sum _{h=1}^{\mathrm{Nh}}{\hat{P}}_{h,l}^{\mathrm{PV}},& \left(26\right)\end{array}$ $\sum _{h=1}^{\mathrm{Nh}}{\beta }_{h,l}=1,$ $l=1\dots \mathrm{PH}$

All PV power is shared by the community; its use for the individual houses is, however, governed by coefficients. In the case of fixed coefficients, the members of the community agree a certain fixed percentage of the total power of all PV systems to be allocated to each house. In the case of employing variable coefficients, the members agree that the optimal percentage of the total power of all PV systems to be allocated to house h is obtained through an optimization problem—here β_h,l are variables to be optimized. The use of consumption-based coefficients the members agree that the percentage to be allocated to each house is calculated through the ratio between the consumption of the house and the total power consumption of all houses of the community. Finally, the use of a hybrid scheme implies that the members agree that the percentage to be allocated to house h is given by a hybrid coefficient that combines a predefined fixed coefficient and the percentage of the consumption relative to all houses of the community.

$\begin{array}{cc}{\beta }_{h,l}=\{\begin{array}{c}{\mathrm{FC}}_{h,l}\in \left[0,1\right];\mathrm{Fixed}\\ {\mathrm{VC}}_{h,l}\in \left[0,1\right];\mathrm{Variable}\\ \frac{{\hat{P}}_{h,l}}{\sum _{l=1}^{\mathrm{Nh}}{\hat{P}}_{h,l}};\mathrm{Consumption}-\mathrm{based}\\ \frac{{\mathrm{FC}}_{h,l}{\hat{P}}_{h,l}}{\sum _{l=1}^{\mathrm{Nh}}{\mathrm{FC}}_{h,l}{\hat{P}}_{h,l}};\mathrm{Hybrid}\end{array}& \left(27\right)\end{array}$

The methodology of energy transfer between the members of community is defined by the community.

11—Examples of forecast of consumption for certain Computer Controlled devices: (step I4) calculation of forecasts of consumption)

Battery Bank

There is a series of restrictions related with the battery bank of house h at instant l:

$\begin{array}{cc}0\le {\hat{P}}_{k,l}^{B-}\le {\stackrel{\_}{P}}_h^{B-}{v}_{h,l}^B& \left(28\right)\end{array}$ $\begin{array}{cc}0\le {\hat{P}}_{k,l}^{B+}\le {\stackrel{\_}{P}}_h^{B+}\left(1-v_{h,l}^B\right)& \left(29\right)\end{array}$ $\begin{array}{cc}{\mathrm{SoC}}_{h,l}={\mathrm{SoC}}_{h,l-1}+\frac{{\eta }_h^{B-}{\hat{P}}_{h,l}^{B-}}{{\stackrel{\_}{E}}_h^B}\Delta k-\frac{{\hat{P}}_{h,l}^{B+}}{{\stackrel{\_}{E}}_h^B{\eta }_h^{B+}}\Delta k& \left(30\right)\end{array}$ $\begin{array}{cc}\underset{\_}{{\mathrm{SoC}}_h}\le {\mathrm{SoC}}_{h,l}\le \stackrel{\_}{{\mathrm{SoC}}_h}& \left(31\right)\end{array}$ $\begin{array}{cc}{\mathrm{SoC}}_{h,0}={\mathrm{SoC}}_{h,\mathrm{initial}}& \left(32\right)\end{array}$ $\begin{array}{cc}{v}_{h,l}^B+v_{h,l}^{\mathrm{Grid}}\le 1& \left(33\right)\end{array}$ $\begin{array}{cc}{\hat{P}}_{h,k}^{B+}-{\hat{P}}_{h,k-1}^{B+}\le {\mathrm{RU}}_h& \left(34\right)\end{array}$ $\begin{array}{cc}{\hat{P}}_{h,k}^{B-}-{\hat{P}}_{h,k-1}^{B-}\le {\mathrm{RU}}_h& \left(35\right)\end{array}$ $\begin{array}{cc}{\hat{P}}_{h,k-1}^{B+}-{\hat{P}}_{h,k}^{B+}\le {\mathrm{RD}}_h& \left(36\right)\end{array}$ $\begin{array}{cc}{\hat{P}}_{h,k-1}^{B-}-{\hat{P}}_{h,k}^{B-}\le {\mathrm{RD}}_h& \left(37\right)\end{array}$

In (28-29) it is defined that {circumflex over (P)}_k,l^B−/{circumflex over (P)}_k,l^B+ are positive variables with maximum value equal to the maximum charge/discharge power P_h^B−/P_h^B+ of the battery bank of house h. v_h,l^B is the binary variable modelling the processes of charge/discharge of battery bank of house h, assuming value 1 if charging or 0 if discharging. The minimum/maximum values of the state of charge of the battery bank, SoC_h/SoC_h are defined in (31). The initial state of charge for each, SoC_h,initial is defined in (32). In (33) it is stated that the binary variables v_h,l^B and. v_h,l^Grid cannot assume the value 1 simultaneously, meaning that the battery cannot be charged directly with energy from the grid. To decrease the battery wear, in (34-35) maximum values of the ramp-up/ramp down rates RU/RD for charging or discharging of the battery bank are defined.

The update of state of charge of battery bank of house h is defined in (30). η_h^B−/η_h^B+ are the charging/discharging efficiency of the battery bank. Ē_h^B is the rated capacity of battery bank of house h.

Swimming Pool Pump

This is an example of a shiftable load, that must obey the following assumptions:

• • the pump must operate during a user-specified time (multiple of T_s) within one day—N_Tspp
  • it must remain on at least for a specified period of time (multiple of Δk)—N_mspp
  • The power consumption is constant and equal to the rated power—P_i,spp

$\begin{array}{cc}\sum _{j=k}^{k+N_{\mathrm{mspp}}-1}{v}_{h,j}^{\mathrm{spp}-\mathrm{on}/\mathrm{off}}\ge N_{\mathrm{Tspp}}{v}_{h,k}^{\mathrm{spp}-\mathrm{start}}& \left(38\right)\end{array}$ $\begin{array}{cc}{P}_{i,\mathrm{spp},k}={\stackrel{\_}{P}}_{i,\mathrm{spp}}{v}_{h,k}^{\mathrm{spp}-\mathrm{on}/\mathrm{off}}& \left(39\right)\end{array}$ $\begin{array}{cc}\sum _{j=1}^K{P}_{i,\mathrm{spp},j}={\stackrel{\_}{P}}_{i,\mathrm{spp}}{N}_{\mathrm{Tspp}}& \left(40\right)\end{array}$ $\begin{array}{cc}{v}_{h,k}^{\mathrm{spp}-\mathrm{start}}-v_{h,k}^{\mathrm{spp}-\mathrm{shutdown}}=v_{h,k}^{\mathrm{spp}-\mathrm{on}/\mathrm{off}}-v_{h,k-1}^{\mathrm{spp}-\mathrm{on}/\mathrm{off}}& \left(41\right)\end{array}$ $\begin{array}{cc}{v}_{h,k}^{\mathrm{spp}-\mathrm{start}}+v_{h,k}^{\mathrm{spp}-\mathrm{shutdown}}\le 1& \left(42\right)\end{array}$

Restriction (38) ensures that the pump must remain on for at least N_mspp*T_s. Restriction (39) says that the power consumption is equal to the pump rated power consumption. Restriction (40) states that the pump must operate N_Tspp steps within a 24 hour period. Restriction (41) ensure the correct operation of the shiftable demand. Finally, (42) imposes that a start-up cannot occur at the same time of a shut-down. v_h,k^spp-on/off, v_h,k^spp-start and v_h,k^spp-shutdown are the binary variables for the state on/off, for the start-up and the shut-down binary variable of the swimming pool pump of house h.

Air Conditioner

It is assumed that the air-conditioner must ensure thermal comfort in a room between user-specified periods of time, or schedules. Assume that one of these periods is [t_s t_e], which corresponds to the sample range [k_s k_e]. The reference temperature (a discrete, integer value) to apply to AC is given as the solution of the following Model Predictive Control (MPC) problem:

$\begin{array}{cc}{\left(t<t_s\right)^\left(t>t_s-\Delta t\right)\underset{U\left[k\right]\in v_{\mathrm{PH}}}{\min }\left(\sum _{i=k+1}^{k_s}J\left[i\right]\right)❘}_{U\left[k\right]}+\underset{s.t.❘\hat{\Theta }\left[i\right]❘<{\Theta }_T,i=1:\mathrm{PH}}{{\underset{U\left[k\right]\in v_{\mathrm{PH}}}{\min }\left(\sum _{i=k_s+1}^{k+\mathrm{PH}}J\left[i\right]\right)❘}_{U\left[k\right]}}& \left(43\right)\end{array}$ $\begin{array}{cc}t+\mathrm{PH}>t_e& \underset{s.t.❘\hat{\Theta }\left[i\right]❘<{\Theta }_T,i=1:\mathrm{PH}}{{\underset{U\left[k\right]\in v_{\mathrm{PH}}}{\min }\left(\sum _{i=k+1}^{k_e}J\left[i\right]\right)❘}_{U\left[k\right]}}\end{array}$ $\begin{array}{cc}\mathrm{others}& \underset{s.t.❘\hat{\Theta }\left[i\right]❘<{\Theta }_T,i=1:\mathrm{PH}}{{\underset{U\left[k\right]\in v_{\mathrm{PH}}}{\min }\left(\sum _{i=k+1}^{k+\mathrm{PH}}J\left[i\right]\right)❘}_{U\left[k\right]}}\end{array}$

In (43), the solution is divided into 3 alternatives, depending on the instant of time considered: a) if it is before the start of the period, we only need to ensure that the room is in thermal comfort when, within the prediction horizon (PH), the time considered is after t_s; b) near the end of the period, we do not need to ensure thermal comfort when time is larger than after t_e; c) all the other cases we determine the optimal sequence of actions, between the admissible sequence of actions (v_PH) throughout PH. J[i] is the cost function at step i, whose profile must be supplied to the system. The thermal comfort (Θ) uses the Predicted Mean Vote (PMV) formulation, and Θ_T is a threshold, typically assuming the value of 0.5. To compute Θ_k for a room, the following information is needed: Air temperature (IT_h,k^room), Relative humidity (RH_h,k^room ) and Wall temperature (WT_h,k^room). Further details can be found in [5]. Δt, the time before the start of the schedule when the AC might need to operate to ensure thermal comfort throughout the whole period, must be specified.

With this explanation, the incorporation of the automatic control of an AC is given by the following equation:

$\begin{array}{cc}{P}_{h,\mathrm{AC},k}={\hat{J}}_{h,\mathrm{AC},k}{v}_{h,k}^{\mathrm{AC}-\mathrm{on}/\mathrm{off}}& \left(44\right)\end{array}$

In (44), Ĵ_H,AC,k denotes the estimated cost of controlling the air conditioner AC of house i at step k and v_h,k^AC-on/off is the binary variable controlling the operation of AC during [t_s-Δt,t_e].

Design of Classification or Forecasting Models (EXAMPLE)

Hereinafter it is described how good models can be obtained from a set of acquired or existent data.

Three sub-problems must be solved:

• • a) A pre-processing stage where, from the existing data, suitable sets (training, testing, validation, etc.) are obtained; this is known as a data selection problem;
  • b) Determine the “best” set of inputs (feature/delays selection) and network topology, given the above data sets;
  • c) Determine the “best” network parameters, given the data sets, inputs/delays and network topology

These sub-problems are solved by the application of a model design framework, composed of two tools. The first, denoted as ApproxHull, performs data selection, from the data available for design. Feature and topology search are solved by the evolutionary part of MOGA (Multi-Objective Genetic Algorithm), while parameter estimation is performed by the gradient part of MOGA.

ApproxHull

To design data driven models like RBFs (Radial Basis Functions), it is mandatory that the training set involves the samples that enclose the whole input-output range where the underlying process is supposed to operate. To determine such samples, called convex hull (CH) points, out of the whole dataset, convex hull algorithms can be applied.

The standard convex hull algorithms suffer from both exaggerated time and space complexity for high dimensions studies. To tackle these challenges in high dimensions, ApproxHull was proposed in [3] as a randomized approximation convex hull algorithm. To identify the convex hull points, ApproxHull employs two main computational geometry concepts: the hyperplane distance and the convex hull distance.

Given the point x=[x₁ . . . x_d]^T in a d-dimensional Euclidean space, and an hyperplane H, the hyperplane distance of x to His obtained by:

$\begin{array}{cc}\mathrm{ds}\left(x,H\right)=\frac{a_1{x}_1+\cdots +a_d{v}_d+b}{\sqrt{a_1^2+\cdots a_d^n}}& \left(45\right)\end{array}$

where n=[a₁, . . . a_d]^T and d are the normal vector and the offset of H, respectively.

Given a set X={x_i}_i=1ⁿ⊂ and a point x∈, the Euclidean distance between x and the convex hull of X, denoted by conv(X), can be computed by solving the following quadratic optimization problem:

$\begin{array}{cc}\underset{a}{\min }\left(\frac{a^T\mathrm{Qa}}{2}-c^{T}a\right)& \left(46\right)\end{array}$ $s.t.e^{T}a=1,$ $a\ge 0$

where e=[1, . . . 1]^T, Q=X^TX and c=X^Tx. Assuming that the optimal solution of (46) is a*, the distance of point x to conv(X) is given by:

$\begin{array}{cc}\mathrm{dc}\left(x,\mathrm{conv}\left(X\right)\right)=\sqrt{x^{T}x-2c^{T}a*+a{*}^T\mathrm{Qa}*}& \left(47\right)\end{array}$

ApproxHull consists of five main steps. In Step 1, each dimension of the input dataset is scaled to the range [−1, 1]. In Step 2, the maximum and minimum samples with respect to each dimension are identified and considered as the vertices of the initial convex hull. In Step 3, a population of k facets based on the current vertices of the convex hull is generated. In Step 4, the furthest points to each facet in the current population are identified using (45) and they are considered as the new vertices of the convex hull, if they have not been detected before. Finally, in Step 5, the current convex hull is updated by adding the newly found vertices into the current set of vertices. Step 3 to Step 5 are executed iteratively until no vertex found in Step 4 or the newly found vertices are very close to the current convex hull, thus not containing useful information. The closest points to the current convex hull are identified using the convex hull distance under an acceptable user-defined threshold.

The flowchart of AproxHull is shown in FIG. 6.

In a prior step before determining the CH points, ApproHull eliminates replicas and linear combinations of samples/features. After having identified the CH points, ApproxHull generates the training, test and validation sets to be used by MOGA, according to user-specifications, but incorporating the CH points in the training set.

MOGA

This framework is described in detail in [6], and it will be briefly discussed here.

MOGA designs static or dynamic ANN (Artificial Neural Network) models, for classification, approximation or forecasting purposes.

Parameter Separability

MOGA employs models that are linear-nonlinearly separable in the parameters. The output of this type of models, at time step k, is given as:

$\begin{array}{cc}\hat{y}\left(x_k,w\right)=u_0+\sum _{i=1}^n{u}_i{\phi }_i\left(x_k,v_i\right)=\phi \left(x_k,v\right)u& \left(48\right)\end{array}$

In (48), x_k is the ANN input at time step k, φis the basis functions vector, u is the (linear) output weights vector and v represents the nonlinear parameters. For simplicity, we shall assume here only one hidden layer, and v is composed of n vectors of parameters, each one for each neuron (v=[v₁. . . v_n]^T). This type of models comprises Multilayer Perceptrons, Radial Basis Function (RBF) networks, B-Spline and Asmod models, Wavelet networks, and Mamdani, Takagi and Takagi-Sugeno fuzzy models (satisfying certain assumptions).

This means that the model parameters can be divided into linear and nonlinear parameters:

$\begin{array}{cc}w=\left[\begin{array}{c}u\\ v\end{array}\right]& \left(49\right)\end{array}$

and that this separability can be exploited in the training algorithms. For a set of input patterns X, training the model means finding the values of w, such that the following criterion is minimized:

$\begin{array}{cc}\Omega \left(X,w\right)=\frac{{y-\hat{y}\left(X,w\right)}_2^2}{2}& \left(50\right)\end{array}$

where ∥·∥₂ denotes the Euclidean norm. Replacing (49) in (50), we have:

$\begin{array}{cc}\Omega \left(X,w\right)=\frac{{y-\Gamma \left(X,v\right)u}_2^2}{2}& \left(51\right)\end{array}$

where Γ(X, v)=[φ(x₁, v) . . . φ(x_m, v)]^T, m being the number of patterns in the training set. As (51) is a linear problem in u, its optimal solution is given as:

$\begin{array}{cc}u*={\Gamma }^{+}\left(X,v\right)y& \left(52\right)\end{array}$

Where the symbol ‘+’ denotes a pseudo-inverse operation. Replacing (52) in (51), we have a new criterion, that is only dependent on the nonlinear parameters:

$\begin{array}{cc}\Psi \left(X,v\right)=\frac{{y-\Gamma \left(X,v\right){\Gamma }^{+}\left(X,v\right)y}_2^2}{2}& \left(53\right)\end{array}$

The advantages of using (53) instead of (51) are threefold:

• • It lowers the problem dimensionality, as the number of model parameters to determine is reduced;
  • The initial value of (53) is much smaller than (51);
  • Typically, the rate of convergence of gradient algorithms using (53) is faster than using (51).

Training Algorithm

Any gradient algorithm can be used to minimize (53) or (51). First-order algorithms, error back-propagation (or steepest descent method), conjugate-gradient method, or their variants, or second-order methods, such as quasi-Newton, Gauss-Newton or Levenberg-Marquardt (LM) algorithms can be employed. For non-linear least-squares problems, the LM method is recognized to be the state-of-the-art method, as it exploits the sum-of-squares characteristics of the problem. The LM search direction is given as the solution of:

$\begin{array}{cc}\left(J_k^T{J}_k+{\lambda }_{k}I\right)p_k=-g_k& \left(54\right)\end{array}$

Where g_k is the gradient vector:

$\begin{array}{cc}{g}_k=\frac{{\mathrm{\delta \Omega }}_k}{\delta w_k}\mathrm{or}\frac{{\mathrm{\delta \Psi }}_k}{\delta v_k}& \left(55\right)\end{array}$

And J_k is the Jacobean matrix:

$\begin{array}{cc}{J}_k=\frac{\delta {\hat{y}}_k}{\delta w_k}\mathrm{or}\frac{\delta {\hat{y}}_k}{\delta v_k}& \left(56\right)\end{array}$

λis a regularization parameter, which enables the search direction to change between the steepest descent and the Gauss-Newton directions.

It has been proved that the gradient and Jacobean of criterion (53) can be obtained by: a) first computing (52); b) replacing the optimal values in the linear parameters; c) and determining the gradient and Jacobian in the usual way.

In MOGA, the training algorithm terminates after a predefined number of iterations or using an early-stopping technique.

Radial Basis Function Networks

In the case described here, RBF models will be employed. The basis functions employed by RBFs are typically Gaussian functions:

$\begin{array}{cc}{\phi }_i\left(x_k,v_i\right)=e^{-\frac{{x_k-c_i}_2^2}{2{\sigma }_i^2}}& \left(57\right)\end{array}$

Which means that the nonlinear parameters for each neuron are constituted by a centre vector, c_i, with as many components as the dimensions of the input, and a spread parameter, σ_i. The derivatives of the basis function with respect to the nonlinear parameters are:

$\begin{array}{cc}\frac{\partial {\phi }_i\left(x_k\right)}{\partial c_{i,j}}={\phi }_i\left(x_k\right)\frac{x_{k,j}-c_{i,j}}{{\sigma }_i^2}& \left(58\right)\end{array}$ $\begin{array}{cc}\frac{\partial {\phi }_i\left(x_k\right)}{\partial {\sigma }_i}={\phi }_i\left(x_k\right)\frac{{x_k-c_i}_2^2}{{\sigma }_i^3}& \left(59\right)\end{array}$

Associated with each model, there are always heuristics that can be used to obtain the initial values of the parameters. For RBBs, centre values can be obtained randomly from the input data, or from the range of the input data. Alternatively, clustering algorithms can be employed. Initial spreads can be chosen randomly, or using several heuristics, such as:

$\begin{array}{cc}\sigma =\frac{z_{\max }}{\sqrt{2b}}& \left(60\right)\end{array}$

where z_max represents the maximum distance between centres.

The Evolutionary Part

MOGA evolves ANN structures, whose parameters separate (in this case RBFs), each structure being trained by the algorithm described before. The text below is for the design of forecasting models; for classifiers we shall describe the changes in the end of the section.

As we shall be designing forecasting models, where we want to predict the evolution of a specific variable within a predefined PH (Prediction Horizon), the models should provide multi-step-ahead forecasting. This type of forecasts can be achieved in a direct mode, by having several one-step-ahead forecasting models, each providing the prediction of each-step ahead within a PH. An alternative, which is the one followed in this work, is to use a recursive version. In this case, only one model is used, but its inputs evolve with time. Consider the Nonlinear Auto-Regressive model with Exogeneous inputs (NARX), with just one input, for simplicity:

$\begin{array}{cc}{\hat{y}}_{k+1❘k}=f\left(z_k\right)=f\left(y_{k-d_{0_1}},\dots ,y_{k-d_{0_n}},x_{k-d_{i_1}},\dots ,x_{k-d_{i_m}}\right)& \left(61\right)\end{array}$

Where ŷ_k+1|k denotes the prediction for time-step k+1 given the measured data at time k, and d_ij the/h delay for variable i. This represents the 1-step-ahead prediction within a prediction horizon. As we iterate (61) over PH, some or all of the indices in the right hand-side will be larger than k, which means that the corresponding forecast must be employed. What has been said for NARX models is also valid for NAR models (with no exogeneous inputs).

The evolutionary part of MOGA evolves a population of ANN structures. Each topology comprises of the number of neurons in the single hidden layer (for an RBF model), and the model inputs or features. MOGA assumes that the number of neurons must be within user-specified bound, n∈[n_m,n_M]. Additionally, one needs to select the features to use for a specific model, i.e, must perform input selection. In MOGA we assume that, from a total number q of available features, denoted as F, each model must select the most representative d features within a user-specified interval, d∈[d_m,d_M], d_M≤q. For this reason, each ANN structure is codified as shown in FIG. 7. In that figure, C denotes a chromosome and F the input space.

The first component corresponds to the number of neurons. The next d_m represents the minimum number of features, while the last white ones are a variable number of inputs, up to the predefined maximum number. The λ_j values correspond to the indices of the features f_j in the columns of F.

The operation of MOGA is a typical evolutionary procedure, represented in FIG. 8.

We shall refer the reader to publication [6] regarding the genetic operators.

The model design cycle is illustrated in FIG. 9.

First, the search space should be defined. That includes the input variables to be considered, the lags to be considered for each variable, and the admissible range of neurons and inputs. The total input data, denoted as F, together with the target data, must then be partitioned into three different sets: training set, to estimate the model parameters, test set, to perform early-stopping, and validation set, to analyze the MOGA performance.

Secondly, the optimization objectives and goals need to be defined. Typical objectives are Root-Mean-Square Errors (RMSE) evaluated on the training set (ρ_tr), or on the test set (ρ_te), as well as the model complexity, #(v)—number of nonlinear parameters—or the norm of the linear parameters (∥u∥₂). For forecasting applications, as it is the case here, one criterion is also used to assess its performance. Assuming a time-series sim, a subset of the design data, with p data points. For each point, the model (61) is used to make predictions up to PH steps-ahead. Then an error matrix is built:

$\begin{array}{cc}E\left(\mathrm{sim},\mathrm{PH}\right)=\left[\begin{array}{cccc}e\left[1,1\right]& e\left[1,2\right]& \dots & e\left[1,\mathrm{PH}\right]\\ e\left[2,1\right]& e\left[2,2\right]& \dots & e\left[2,\mathrm{PH}\right]\\ ⋮& ⋮& \ddots & ⋮\\ e\left[p-\mathrm{ph},1\right]& e\left[p-\mathrm{ph},2\right]& \dots & e\left[p-\mathrm{ph},\mathrm{PH}\right]\end{array}\right]& \left(62\right)\end{array}$

where e[i,j] is the model forecasting error taken from instant i of sim, at step j within the prediction horizon. Denoting the RMS function operating over the i^th column of matrix E, by ρ_sim(.,i), the forecasting performance criterion is the sum of the RMS of the columns of E:

$\begin{array}{cc}{\rho }_{\mathrm{sim}}\left(\mathrm{PH}\right)=\sum _{i=1}^{\mathrm{PH}}\rho \left(E\left(\mathrm{sim},\mathrm{PH}\right),i\right)& \left(63\right)\end{array}$

One should notice that every performance criterion can be minimized, or set as a restriction, in the MOGA formulation.

After having formulated the optimization problem, and after setting other hyperparameters, such as the number of elements in the population (n_pop), number of iterations population (n_iter), and genetic algorithm parameters (proportion of random immigrants, selective pressure, crossover rate and survival rate), the hybrid evolutive-gradient method is executed.

Each element in the population corresponds to a certain RBF structure. As the model is nonlinear, a gradient algorithm such as the LM algorithm minimizing (53) is only guaranteed to converge to a local minimum. For this reason, the RBF model is trained a user-specified number of times, starting with different initial values for the nonlinear parameters. MOGA allows initial centers to be chosen from the heuristics, or using an adaptive clustering algorithm.

As the problem is multi-objective, there are several ways for identifying which training trial is the best one. One strategy is to select the training trial whose Euclidean distance from the origin is the smallest. Symbol “o” in FIG. 10 illustrates this situation for d=2. In the second strategy, the average of objective values for all training trails is calculated and then the trial whose value is the closest to the average value will be selected as the best one (i.e., symbol “x”—or nearest to average, in FIG. 10).

The other d strategies are to select the training trial which minimized the i^th objective (i.e., i=1,2, . . . , d) better than the other trials. As an example, the symbols ‘<’e‘>’ in FIG. 10 are the best training trials which minimized objective 1 and objective 2, respectively.

After having executed the specified number of iterations, we have performance values of n_pop*n_iter different models. As the problem is multi-objective, a subset of these models corresponds to non-dominated models (nd), or Pareto solutions. If one or more objectives is (are) set as restriction(s), a subset of nd, denoted as preferential solutions, pref, corresponds to the non-dominated solutions, which meet the goals. An example is shown in FIG. 11, wherein the small dots represent all models, the larger dots in black represent non-dominated solutions and the large dot in white represents the chosen solution. In the graph shown in FIG. 11, ε_t represents an objective defined as the root mean square error obtained in the training set and ε_g represents an objective defined as the root mean square error obtained in the generalization training data set.

The performance of MOGA models is assessed on the non-dominated models set, or in the preferential models set. If a single solution is sought, it will be chosen on the basis of the objective values of those model sets, performance criteria applied to the validation set, and possibly other criteria.

When the analysis of the solutions provided by the MOGA requires the process to be repeated, the problem definition steps should be revised. In this case, two major actions can be carried out: input space redefinition by removing or adding one or more features (variables and lagged input terms in the case of modelling problems) and improving the trade-off surface coverage by changing objectives or redefining goals. This process may be advantageous as usually the output of one run allows reducing the number of input terms (and possibly variables for modelling problems) by eliminating those not present in the resulting population. Also, it usually becomes possible to narrow the range for the number of neurons in face of the results obtained in one run. This results in a smaller search space in a subsequent run of the MOGA, possibly achieving a faster convergence and better approximation of the Pareto front.

Typically, for a specific problem, an initial MOGA execution is performed, minimizing all objectives. Then a second execution is run, where typically some objectives are set as restrictions.

Robust Model Ensemble

In the end of the 2^nd MOGA execution, we obtain a set of preferential models. A sub-set of those, denoted as outset, will be selected using the procedure described below, and the model output, for a forecasting application, will be given as:

$\begin{array}{cc}\stackrel{⌣}{y}(·)=\mathrm{median}\left(\left\{{\hat{y}}_i(·),i\in \mathrm{outset}\right\}\right)& \left(64\right)\end{array}$

For a classification application, the output is given as:

$\begin{array}{cc}\stackrel{⌣}{y}(·)=\mathrm{majority\_vote}\left(\left\{\mathrm{step}\left({\hat{y}}_i(·)\right),i\in \mathrm{outset}\right\}\right)& \left(65\right)\end{array}$

Where the step function is given as:

$\begin{array}{cc}\mathrm{step}\left(x\right)=\{\begin{array}{cc}1,& x\ge 0\\ -1,& x<0\end{array}& \left(66\right)\end{array}$

In certain embodiments of a first aspect, the invention relates to a computer implemented method for community energy management (1000), comprising the following steps:

• • a) acquisition, by an aggregator computational unit (110), of data d₀ that will be used for designing prediction models and appliance profiles, during a period of time, preferably during at least one month;
  • b) design, using a cloud-based server (115) of Radial Basis Function prediction models to use for predicting selected from the group consisting of: global solar radiation and atmospheric air temperature of the community, power consumption and power generation in each house h of the community, as well as operation and power consumption for all monitored appliances by Non-Invasive Load Monitoring models, wherein said design comprises the following substeps:
    • b1) selection of the data available for design from the data d₀ acquired in step a);
    • b2) minimization of the model objectives through a Multi-Objective Genetic Algorithm;
    • b3) analysis of the results of the first design, restricting some of the objectives minimized in step b2) and redefinition of the input space and/or the model topology;
    • b4) utilization of the preferential set of models obtained from a second Multi-Objective Genetic Algorithm execution;
    • b5) determination of a second set of the preferred models m_e^pref obtained in the previous substep through the following equation:

$m_e^{\mathrm{pref}}=\left\{m_{\mathrm{pref}}:\mathrm{nw}\left(m_{\mathrm{pref}}\right)\le {\stackrel{⌣}{m}}_{\mathrm{nw}}^{\mathrm{pref}}\right\}\bigcap \left\{m_{\mathrm{pref}}:\mathrm{fore}\left(m_{\mathrm{pref}}\right)\le {\stackrel{⌣}{m}}_{\mathrm{fore}}^{\mathrm{pref}}\right\}$

• • wherein m_pref is the set of preferred models obtained in the second Multi-Objective Genetic Algorithm execution, m̆_nw^pref is the median of the linear weight norm for each model in m_pref and m̆_fore^pref the median of the sum of the forecasts for the prediction time-series obtained by models m_pref;
    • b6) obtaining of a third set of models, m_e^par from the previous step by using the following equation:

$m_e^{\mathrm{par}}=\left\{\mathrm{nd}\left(\mathrm{fore}\left(m_e^{\mathrm{pref}}\right),\mathrm{nw}\left(m_e^{\mathrm{pref}}\right)\right):\#\mathrm{nd}=\mathrm{NS}\right\},$

• • wherein nd(.) stands for non-dominated solutions and NS is its number, preferably NS=25;
    • b7) if the model that should be designed is a Nonlinear Auto-Regressive model with Exogeneous inputs, usage, by the server of m_e^par for the exogenous variables;
  • c) determination of the power profiles for all monitored, manually controlled and computer-controlled appliances in the community, by the cloud-based server (115);
  • d) reception, by the aggregator computational unit (110), of all models designed in b), all power profiles determined in c), initial values v₀ and initial control actions, through a user-interface, as well as data history for all models supplied by a data acquisition subsystem (140) installed in each house h, one or more consumption energy meters (130) installed in at least one house h, one or more inverter energy meters (160) installed in at least one house h, one or more sensors (150) installed in at least one house h;
  • e) reception, by a data acquisition subsystem (140), of data related to house h, wherein said data is sent wirelessly and/or wired and wherein said reception includes the following substeps:
    • e1) reception of electric consumption data di sent by a consumption energy meter (130) installed at each house h, wherein data d₁ comprise active power being consumed in house h at all instants k₁ within the sampling interval Δk, P_h,k1, and reactive power being consumed in house h at all instants k₁ within the sampling interval Δk, R_h,k1;
    • e2) reception of data d₂ related to at least one energy producing photovoltaic subsystem installed in at least one house h, sent by an inverter energy meter (160) wherein data d₂ comprise the power generated by at least one energy producing photovoltaic subsystem of house h at an instant k, P_h,k^PV; the power imported from the grid P_h,k^Grid−, and the power exported to the grid P_h,k^Grid+, wherein PV refers to photovoltaic;
    • e3) reception of data d₃ related to at least one battery bank installed in at least one house h, sent by an inverter energy meter (160) wherein data d₃ comprises at least the state of charge of the battery;
  • f) sending of the control actions determined in the previous step, by the aggregator computational unit (110) to each data acquisition subsystem (140), which, on its turn, sends the respective control actions to the computer-controlled devices in each house of the community;
  • g) reception, by an aggregator computational unit (110), of the data d₁, d₂ and d₃ related to house h at an instant k, sent by the data acquisition subsystem (140);
  • h) reception, by the aggregator computational unit (110), of weather data d₄, related to the community, wherein said weather data d₄ comprises air temperature data and global solar radiation data;
  • i) reception of tariff-related data ds of each house h in the community, wherein said data d₅ is retrieved from a database (125) by the aggregator computational unit (110);
  • j) reception, by the aggregator computational unit (110), of occupation-related data d₆ of each house h, inserted in a database by each user via a user interface (145);
  • k) determination of the averaged active power consumption of each house:

$P_{n,k}=\frac{P_{h,k_1}}{N},N=\frac{\Delta k}{k_1};$

• • l) calculation of forecasts for a user-specified Prediction Horizon PH of both the energy consumption of each individual house h and the energy generation of each energy producing photovoltaic apparatus in the community over the prediction horizon, by using one or more of the methods of the group consisting of naive forecasting methods, multi-step Nonlinear Auto-Regressive Exogenous forecasting models, AutoRegressive Integrated Moving Average, ARIMA, Classification and Regression Trees, CART, Wavelet, Fuzzy or Neuro-Fuzzy models, and an hybrid combination thereof; wherein the present calculation step comprises the following substeps:
    • l1) calculation of forecasts of consumption for the devices in each house that are monitored and not controlled NMC, which is given as the sum of the active power of corresponding appliances in each house, given by:

$\begin{array}{cc}{\hat{P}}_{h,\left\{k+1\dots k+\mathrm{PH}\right\}}^{\mathrm{NMC}}=\sum _{j=1}^{\#\left(\left\{{\mathrm{NM}}_i\right\}\bigcup \left\{{\mathrm{NC}}_i\right\}\right)}{\hat{P}}_{h,j,\left\{k+1\dots k+\mathrm{PH}\right\}};& \left(I\right)\end{array}$

• • • for the devices in each house h that are monitored but not controlled NC by means of:

$\begin{array}{cc}{\hat{P}}_{h,\left\{k+1\dots k+\mathrm{PH}\right\}}^{\mathrm{NMC}}={\hat{P}}_{h,\left\{k+1\dots k+{\mathrm{PH}}_{\mathrm{NARX}}\right\}}^{\mathrm{NMC}}\bigcup {\hat{P}}_{h,\left\{k+{\mathrm{PH}}_{\mathrm{NARX}}+1\dots k+\mathrm{PH}\right\}}^{\mathrm{NMC}}& \left(\mathrm{II}\right)\end{array}$

• • • wherein the first term in the right-hand-side of (II) refers to the use of the Nonlinear Auto-Regressive model with Exogeneous NARX model and the second term to the naïve forecast technique;
    • wherein PH_NARX denotes the first steps within PH whose forecasts are obtained as the median of the outputs of the models m_e^pref:

$\begin{array}{cc}{\hat{P}}_{h,\left\{k+1\dots k+{\mathrm{PH}}_{\mathrm{NARX}}\right\}}^{\mathrm{NMC}}=\left\{M_h^P\left({\stackrel{\_}{P}}_{h,l}^{\mathrm{NMC}},{\stackrel{\_}{T}}_l,O_{h,l},D_l\right),l=k+1\dots k+{\mathrm{PH}}_{\mathrm{NARX}}\right\};& \left(\mathrm{III}\right)\end{array}$

• • • wherein the second term is given as:

$\begin{array}{cc}{\hat{P}}_{h,\left\{k+{\mathrm{PH}}_{\mathrm{NARX}}+1\dots k+\mathrm{PH}\right\}}^{\mathrm{NMC}}=P_{h,\left\{k+{\mathrm{PH}}_{\mathrm{NARX}}+1-\frac{24}{\Delta k}\dots k+\mathrm{PH}-\frac{24}{\Delta k}\right\},}^{\mathrm{NMC}}& \left(\mathrm{IV}\right)\end{array}$

• • • i.e., the value obtained one day ago;
    • wherein P_h,l^NMC denotes a set of past values relative to l of P_h,l^NMC, T_l denotes a set of past values relative to l of T_l O_h,l denotes the daily occupation relative to instant/of house h, and D_l corresponds to a codification of the day, within the week:
    • wherein as the multi-step iteration of (III) involves predictions of the atmospheric temperature, these forecasts will be given as a multi-step iteration of a NAR model:

$\begin{array}{cc}{\hat{T}}_{\left\{k+1\dots k+{\mathrm{PH}}_{\mathrm{NARX}}\right\}}=\left\{M^T\left({\stackrel{\_}{T}}_l\right),l=k+1\dots k+{\mathrm{PH}}_{\mathrm{NARX}}\right\};& \left(V\right)\end{array}$

• • l2) calculation of forecasts of consumption for the devices in each house that are monitored but not controlled NC, which is given as the sum of the active power of corresponding appliances in each house, given by:

$\begin{array}{cc}{\hat{P}}_{h,\left\{k+1\dots k+\mathrm{PH}\right\}}^{\mathrm{NC}}=\sum _{j=1}^{\#\left(\left\{M_i\right\}\bigcap \left\{{\mathrm{NC}}_i\right\}\right)}{\hat{P}}_{h,j,\left\{k+1\dots k+\mathrm{PH}\right\}};& \left(\mathrm{VI}\right)\end{array}$

• • • wherein the forecasts of consumption of each device, {circumflex over (P)}_{h,j,{k+1 . . .k+PH}}, are given by the joint application of Non-Invasive Load Monitoring NILM classification models configured to determine if the appliance is operating or is off for every instant in PH and, if the appliance is on, the estimate of consumption given by the NILM forecasting model, and the power profile of the appliance, determined offline and updated on-line;
    • l3) calculation of forecasts of consumption for the devices in each house that are manually controlled MC, which is given as the sum of the active power of corresponding appliances in each house, given by:

$\begin{array}{cc}{\hat{P}}_{h,\left\{k+1\dots k+\mathrm{PH}\right\}}^{\mathrm{MC}}=\sum _{m=1}^{\#\left(\left\{{\mathrm{MC}}_i\right\}\right)}{\hat{P}}_{h,j,\left\{k+1\dots k+\mathrm{PH}\right\}};& \left(\mathrm{VII}\right)\end{array}$

• • • wherein the forecasts of consumption of each device

$\sum _{m=1}^{\#\left(\left\{{\mathrm{MC}}_i\right\}\right)}{\hat{P}}_{h,m,\left\{k+1\dots k+\mathrm{PH}\right\}}$

• • • are obtained using the power profiles determined off-line;
    • l4) calculation of forecasts of consumption for the devices in each house that are computer controlled CC, which is given as the sum of the active power of corresponding appliances in each house, given by:

$\begin{array}{cc}{\hat{P}}_{h,\left\{k+1\dots k+\mathrm{PH}\right\}}^{\mathrm{CC}}=\sum _{n=1}^{\#\left(\left\{{\mathrm{CC}}_i\right\}\right)}{\hat{P}}_{h,n,\left\{k+1\dots k+\mathrm{PH}\right\}};& \left(\mathrm{VIII}\right)\end{array}$

• • • wherein the operation of each appliance is given as the solution of the Mixed-Integer Linear Programming MILP problem carried out in step m), and whose consumption estimation can be given using the power profiles determined off-line;
    • l5) calculation of forecasts of photovoltaic power generation by means of a combination of the naive forecasting method and multi-step Nonlinear Auto-Regressive Exogenous NARX forecasting models for each energy generation photovoltaic subsystem in each house h comprising an energy generation photovoltaic subsystem; the NARX forecasts are given by:

$\begin{array}{cc}{\mathrm{PV}}_{h,\left\{k+1\dots k+{\mathrm{PH}}_{\mathrm{NARX}}\right\}}=\left\{M_h^{\mathrm{PV}}\left(\stackrel{\_}{P}{\stackrel{\_}{V}}_{h,l},{\stackrel{\_}{T}}_l{\stackrel{\_}{,R}}_l\right),l=k+1\dots k+{\mathrm{PH}}_{\mathrm{NARX}}\right\};& \left(\mathrm{IX}\right)\end{array}$

wherein R denotes global solar radiation,wherein an additional forecasting model is needed for solar radiation over PH, to determine its evolution is used:

$\begin{array}{cc}{R}_{\left\{k+1\dots k+36\right\}}=M^R\left(\stackrel{\_}{R_l}\right),l=k+1\dots k+36\}& \left(X\right)\end{array}$

• • m) formulation, by the aggregator computational unit (110), of a Mixed-Integer Linear Programming problem as indicated in formula (XI) below

$\begin{array}{cc}\min {\hat{C}}^{\mathrm{Com}}=\min \sum _{h=1}^{\mathrm{Nh}}\sum _{s=1}^{\mathrm{PH}}\left({\hat{P}}_{h,l}^{\mathrm{Grid}-}+{\hat{P}}_{h,l}^{\mathrm{Com}-}\right){\pi }_l^{\mathrm{Grid}-}\Delta k-\left({\hat{P}}_{h,l}^{\mathrm{Grid}+}+{\hat{P}}_{h,l}^{\mathrm{Com}+}\right){\pi }_l^{\mathrm{Grid}+}\Delta k,l=k+s& \left(\mathrm{XI}\right)\end{array}$

• • • wherein
    • Nh denotes the number of houses in the community;
    • PH is a prediction horizon, and
    • Δk is the sampling interval;
    • {circumflex over (P)}_h,l^Grid−/{circumflex over (P)}_h,l^Grid+ are estimates of the active power P bought−/sold+ from/to an electrical energy grid at a sampling instant/by house h in the community;
    • {circumflex over (P)}_h,l^Com−/{circumflex over (P)}_h,l^Com+ are estimates of the active power P bought−/sold+ from/to the community at the sampling instant/by house h;
    • π_l^Grid−/π_l^{Grid +} are the price of the energy bought-/sold+from/to the community or the electrical energy grid at the sampling instant I;
    • wherein the following constraints (XII) and (XIII) related to the electrical energy grid are satisfied:

$\begin{array}{cc}0\le {\hat{P}}_{h,l}^{\mathrm{Grid}+}\le {\stackrel{\_}{P}}_h^{\mathrm{Grid}+}\left(1-v_{h,l}^{\mathrm{Grid}}\right)& \left(\mathrm{XII}\right)\end{array}$ $\begin{array}{cc}0\le {\hat{P}}_{h,l}^{\mathrm{Grid}-}\le {\stackrel{\_}{P}}_h^{\mathrm{Grid}-}{v}_{h,l}^{\mathrm{Grid}};& \left(\mathrm{XIII}\right)\end{array}$

• • • wherein
    • {circumflex over (P)}_h,l^Grid+/{circumflex over (P)}_h,l^Grid− are positive variables with a maximum value equal to the contracted power form the electrical energy grid, P_h^Grid+/P_h^Grid− by house h,
    • wherein v_h,l^Grid is the binary variable modelling the selling/purchasing energy from house h to/from the electrical energy grid, assuming 1 for purchasing and 0 for selling; and wherein the following constraints (XIV) and (XV) related to the community are satisfied:

$\begin{array}{cc}0\le {\hat{P}}_{h,l}^{\mathrm{Com}+}\le {\stackrel{\_}{P}}_h^{\mathrm{Com}+}\left(1-v_{h,l}^{\mathrm{Com}}\right)& \left(\mathrm{XIV}\right)\end{array}$ $\begin{array}{cc}0\le {\hat{P}}_{h,l}^{\mathrm{Com}-}\le {\stackrel{\_}{P}}_h^{\mathrm{Com}-}{v}_{h,l}^{\mathrm{Com}}& \left(\mathrm{XV}\right)\end{array}$

• • • wherein
    • P _h,l ^Com+/P_h,l^Com− are positive variables with a maximum value equal to the maximum value exchangeable with the community, P_h,l ^Com+/P_h^{Com −} by house h, wherein v_h,l^Com is the binary variable modelling the selling/purchasing energy from house h to/from the community, assuming 1 for selling and 0 for purchasing, and considering that it is not allowed to buy energy form the grid and sell energy to the community at the same time: v_h,l^Grid+v_h,l^Com≤1;
    • wherein the energy balance of the community is given as (XVI):

$\begin{array}{cc}\sum _{h=1}^{\mathrm{Nh}}{\hat{P}}_{k,l}^{\mathrm{Com}-}=\sum _{h=1}^{\mathrm{Nh}}{\hat{P}}_{k,l}^{\mathrm{Com}+};& \left(\mathrm{XVI}\right)\end{array}$

• • • wherein the energy balance of house h in the community can be obtained by several ways, selected by the community:
    • m1) if there are no restrictions on the transfer of energy between the community, the energy balance of house h is given by the equality restriction:

$\begin{array}{cc}{\hat{P}}_{h,l}^{\mathrm{PV}}+{\hat{P}}_{h,l}^{B+}+{\hat{P}}_{h,l}^{\mathrm{Grid}-}+{\hat{P}}_{h,l}^{\mathrm{Com}-}={\hat{P}}_{h,l}^{B-}+{\hat{P}}_{h,l}^{\mathrm{Grid}+}+{\hat{P}}_{h,l}^{\mathrm{NMC}}+{\hat{P}}_{h,l}^{\mathrm{NC}}+{\hat{P}}_{h,l}^{\mathrm{MC}}+{\hat{P}}_{h,l}^{\mathrm{CC}}+{\hat{P}}_{h,l}^{\mathrm{Com}+}& \left(\mathrm{XVII}\right)\end{array}$

• • • wherein {circumflex over (P)}_h,l^PV is the estimate of the alternate current power generated by the PV installation of house h at instant I; wherein {circumflex over (P)}_h,l^B+/{circumflex over (P)}_h,l^B− is the estimate of the discharged or charged power of the battery bank of house h at instant I;
    • wherein the {circumflex over (P)}_h/l^NMC term relates to the estimate of the active power of house h being consumed at instant I by devices that are not monitored neither controlled, the {circumflex over (P)}_h,l^NC term relates to the estimate of the active power of house h being consumed at instant/by devices that are monitored but not controlled, the {circumflex over (P)}_h/l^MC term relates to the estimate of the active power of house h being consumed at instant/by devices manually controlled, and the PCC term relates to the estimate of the active power of house h being consumed at instant/by devices that are computer controlled;
    • m2) energy balance of house h (VIII) in the community, if it is not allowed the transfer of energy between the community, is given by the equality restriction:

$\begin{array}{cc}\begin{array}{c}{\hat{P}}_{h,l}^{\mathrm{Com}-}=P_{h,l}^{\mathrm{Com}+}=0⟹\\ ⟹{\hat{P}}_{h,l}^{\mathrm{PV}}+{\hat{P}}_{h,l}^{B+}+{\hat{P}}_{h,l}^{\mathrm{Grid}-}={\hat{P}}_{h,l}^{B-}+{\hat{P}}_{h,l}^{\mathrm{Grid}+}+{\hat{P}}_{h,l}^{\mathrm{NMC}}+{\hat{P}}_{h,l}^{\mathrm{NC}}+{\hat{P}}_{h,l}^{\mathrm{MC}}+{\hat{P}}_{h,l}^{\mathrm{CC}}\end{array}& \left(\mathrm{XVIII}\right)\end{array}$

• • • m3) wherein the energy is shared between the houses using energy sharing coefficients Bh, as shown in the formulas (XIX) e (XX):

$\begin{array}{cc}\begin{array}{c}{\hat{P}}_{h,l}^{\mathrm{Com}+}={\hat{P}}_{h,l}^{\mathrm{PV}}⟹\\ ⟹{\hat{P}}_{h,l}^{B+}+{\hat{P}}_{h,l}^{\mathrm{Grid}-}+{\hat{P}}_{h,l}^{\mathrm{Com}-}={\hat{P}}_{h,l}^{B-}+{\hat{P}}_{h,l}^{\mathrm{Grid}+}+{\hat{P}}_{h,l}^{\mathrm{NMC}}+{\hat{P}}_{h,l}^{\mathrm{NC}}+{\hat{P}}_{h,l}^{\mathrm{MC}}+{\hat{P}}_{h,l}^{\mathrm{CC}}\end{array}& \left(\mathrm{XIX}\right)\end{array}$ $\begin{array}{cc}\begin{array}{cc}{\hat{P}}_{h,l}^{\mathrm{Com}-}={\beta }_{h,l}\sum _{h=1}^{\mathrm{Nh}}{\hat{P}}_{h,l}^{\mathrm{PV}},& \sum _{h=1}^{\mathrm{Nh}}{\beta }_{h,l}=1,l=1\dots \mathrm{PH},\end{array}& \left(\mathrm{XX}\right)\end{array}$

• • • wherein said coefficients follow one of the following different alternatives (XXI):

$\begin{array}{cc}{\beta }_{h,l}=\{\begin{array}{c}{\mathrm{FC}}_{h,l}\in \left[0,1\right];\mathrm{Fixed}\\ {\mathrm{VC}}_{h,l}\in \left[0,1\right];\mathrm{Variable}\\ \frac{{\hat{P}}_{h,l}}{\sum _{h=1}^{\mathrm{Nh}}{\hat{P}}_{h,l}};\mathrm{Consumption}-\mathrm{based}\\ \frac{{\mathrm{FC}}_{h,l}{\hat{P}}_{h,l}}{\sum _{l=1}^{\mathrm{Nh}}{\mathrm{FC}}_{h,l}{\hat{P}}_{h,l}};\mathrm{Hybrid}\end{array};& \left(\mathrm{XXI}\right)\end{array}$

• • • wherein FC_h,l represent fixed coefficients, and VC_h,l variable coefficients which are obtained through an optimization problem;
  • n) solution, by the aggregator computational unit (110), of the problem formulated in step m), using the receding horizon principle, wherein said solution includes the determination of all control actions, for the whole prediction horizon PH, but where only the ones calculated for the first step within PH will be sent to the corresponding data-acquisition systems (140), in step f);
  • o) update of the history for every model used;wherein steps e) through o) are repeated every sampling interval Δk.

In other embodiments of the first aspect of the present invention, data do acquired in step a) is selected from the group consisting of:

• • tws, time basis for weather data;
  • R, global Solar Radiation;
  • T, atmospheric air temperature; and, for every house h:
  • t_Cons, time basis for consumption data;
  • P_h, active power of house h;
  • R_h, reactive power of house h,
  • t_Inv, time basis for inverter data;
  • {P_h^PV}, alternate power generated by the PV installation of house h,
  • {O_h}, daily occupation of house h; and, for every house h that has manually controlled appliances:
  • StP_h,j/EndP_h,j, j∈{MC_h}, start/end times of the operation of all manually controlled appliances of house h;
  • {P_h,j, j∈{MC_h}}, active Power of all Manually Controlled appliances of house h;
  • t_MCh, time basis for all Manually Controlled appliances of house h; and for every house h that has Computer Controlled Appliances that are not air conditioning systems:
  • {U_h,j, j∈{CC_h}} control actions for all appliances of house h to be computer controlled;
  • {P_h,j, j∈{CC_h}}, active power of all appliances of house h to be Computer Controlled;
  • t_cCh, Time basis for all Computer Controlled appliances of house h; and for every house h that has Computer Controlled air conditioning systems:
  • {IT_h^room, room∈{AC_h}}, Air temperature for all air conditioning systems of house h to be Computer Controlled;
  • {RH_h^room, room∈{AC_h}}, Relative humidity for all air conditioning systems of house h to be Computer Controlled;
  • {WT_h^room, room∈{AC_h}} Wall temperature for all air conditioning systems of house h to be Computer Controlled;
  • {P_h^room, room∈{AC_h}}, Active Power for all air conditioning systems of house h to be Computer Controlled;
  • {U_h^room, room∈{AC_h}}, Control Actions for all air conditioning systems of house h to be Computer Controlled;
  • t_Ach, Time basis for all air conditioning systems of house h to be Computer Controlled; and for every house h that has Monitored but not Controlled appliances:
  • {StP_h,j/EndP_h,j, j∈{M_h}}, Start/end times of the operation of all Monitored appliances of house h;
  • {P_h,j, j∈{M_h}}, Active Power of all Monitored appliances of house h;
  • t_Mh—Time basis for all Monitored appliances of house h.

Yet in other embodiments of the first aspect of the present invention, the initial values v₀ of step b) are at least one selected from the group consisting of, for all houses with a photovoltaic energy generation subsystem:

• • P_h^Grid−/P_h^Grid+ contracted bought(−)/sold(+) power from/to an electrical energy grid by each house h; and for all houses with a photovoltaic energy generation subsystem and battery:
  • P_h^{B −}/P_h^B+ maximum discharged/charged power of a battery bank of at least one house h;
  • SOC_h,SoC_h minimum/maximum value of the state of charge of a battery bank of a house h;
  • SoC_h,initial initial state of charge of a battery bank of a house h;
  • RU_h/RD_h maximum values of ramp-up/ramp down rates for charging or discharging of a battery bank of a house h;
  • η_h^B−/η_h^B+ charging/discharging efficiency of a battery bank of a house h;
  • Ē_h^B rated capacity of a battery bank of a house h; and, for all houses with computer-controlled air conditioning systems:
  • t_s start time when a room must be in thermal comfort;
  • t_e end time when a room must be in thermal comfort;
  • Δt time interval before the start of the schedule when an air conditioning system needs to operate to ensure thermal comfort throughout a whole period;
  • Θ_T predicted mean value threshold.

Still in other embodiments of the first aspect of the present invention, the prediction horizon is between 1 hour and 30 days, preferably between 2 hours and 15 days, and most preferably the prediction horizon is 1 day.

Yet in other embodiments of the first aspect of the present invention, weather data d₄ is retrieved from a weather database and/or data d₄ is sent by a weather station (120).

In certain embodiments of the first aspect of the present invention, the computer implemented method for community energy management (1000) further comprises the following substep:

• • h4) reception of data d₇ related to climate information about at least one room in at least one house in the community, sent by at least one sensor (150), wherein the at least one sensor (150) is electrically powered; and wherein substep h4) is carried out within step h).

In other embodiments of the first aspect of the present invention, the sampling interval Δk belongs to the interval from 5 minutes to 1 hour, most preferably the sampling interval Δk is 15 minutes and that the sampling interval Δk₁ belongs to the interval between 10 seconds to 2 minutes, most preferably the sampling interval Δk₁ is 1 minute.

Yet in other embodiments of the first aspect of the present invention, the computer-controlled devices are selected from the group consisting of air conditioning systems (180), battery banks (190), swimming pool pumps (200), or any other electric device capable of being controlled remotely, namely using a smart plug (210).

In certain embodiments of a second aspect, the present invention relates to a community energy management system (100) comprising:

• • a source of weather-related data d₄ wherein such data d₄ is at least one of the group consisting of air temperature, relative humidity, global solar radiation;
  • an aggregator computational unit (110) for community energy management system comprising:
    • a communication unit (11), configured to communicate with one or more data acquisition systems (140), one or more consumption energy meter (130), one or more inverter energy meter (160), one or more sensors (150), one or more smart plugs (210);
    • a processing unit (12);
    • a data storage unit (13), configured to store data before, during or after processing by the processing unit (12);
  • at least one energy generation photovoltaic subsystem, comprised of photovoltaic panels and one inverter;
  • at least one battery bank (170), configured to store energy produced by the at least one energy generation photovoltaic subsystem;
  • at least one data acquisition subsystem (140) installed in each house h of the community;
  • at least one consumption energy meter (130) configured to measure the energy consumption at each house h;
  • at least one inverter energy meter (160) configured to measure the energy produced at each house h;
  • at least one smart plug (180), configured to switch on and off electrical equipment connected thereof;
  • at least one cloud-based server (115);wherein the aggregator computational unit (110) is configured to carry out the method according to the first aspect of the invention.

In certain embodiments of the second aspect, the community energy management system (100) according to the invention further comprises:

• • one or more indoor sensors (150) configured to measure at least one quantity from the group consisting of air temperature, relative humidity, wall temperature.

Yet in certain embodiments of the second aspect, the community energy management system (100) according to the invention further comprises:

• • the source of weather-related data d₄ is a weather station (120) configured to measure d₄ wherein such data d₄ is selected from the group consisting of air temperature, relative humidity, global solar radiation.

Yet in certain embodiments of the second aspect, the aggregator computational unit (110) is selected from the group consisting of a server, a desktop computer, a laptop computer, or a combination thereof.

In a third aspect, the invention relates to a computer program product comprising logic instructions which, when executed by a computer of the community energy management system (100) of the second aspect, cause the computer to carry out the method of the first aspect of the present invention.

In a fourth aspect, the invention relates to a computer-readable storage medium comprising instructions which, when executed by a computer of the community energy management system (100) of the second aspect, cause the computer to carry out the method of the first aspect of the present invention.

DEFINITIONS

Throughout this application, the terms “around” and “approximately” mean that the value to which they refer could also be within a range of plus and minus 10% of the mentioned value.

As used in this application, the term “or” is to be interpreted in an inclusive meaning instead of the exclusive meaning, unless otherwise clearly stated. That is, an expression like “X utilizes A or B” shall be interpreted as including all possible combinations, i.e., “X utilizes A”, “X utilizes B”, and “X utilizes A and B”.

As used in this application, the indefinite article “a”, “an”, shall be interpreted as including “one” or “one or more”, unless otherwise clearly stated.

Throughout this application, the examples provided shall be interpreted as having the purpose to illustrate one or more examples of embodiments of the present invention and shall not be interpreted as preferences, unless otherwise clearly stated.

As used throughout the present application, the terms “comprise/comprises”, “comprising”, “include/includes”, “including” specify the presence of the features, elements, components, steps, and related operations, and do not exclude whatsoever the presence of further features, elements, components, steps, and related operations.

The subject-matter above-described is provided as an illustration of the present invention and shall not be interpreted as limiting it. The terminology used with the purpose of describing specific embodiments according to the present invention, shall not be interpreted as a limitation of the invention.

REFERENCE SIGNS LIST

• • 11—communication unit
  • 12—processing unit
  • 13—data storage unit (memory, data base)
  • 100—system according to the invention
  • 110—aggregator computational unit
  • 115—cloud-based server
  • 120—weather station
  • 125—database (tariff-related data)
  • 130—consumption energy meter
  • 140—data acquisition subsystem
  • 145—user interface
  • 150—sensors
  • 160—inverter energy meter
  • 180—smart plug
  • 190—computer-controlled equipment
  • 1000—method according to the invention

CITATIONS LIST

In the following, the citations list is presented:

Patent Literature

CA2809011 A1, “Adaptative energy management system”, priority date Nov. 6, 2012, published May 6, 2014, MacMaster University, CA.

Non-Patent Literature

• • [1] W. Su, A. Huang, in “The energy internet: An open energy platform to transform legacy power systems into open innovation and global economic engines”, Woodhead Publishing, 2018.
  • [2] A. Manso-Burgos, D. Ribó-Pérez, T. Gómez-Navarro, M. Alcázar-Ortega, in “Local energy communities modelling and optimisation considering storage, demand configuration and sharing strategies: A case study in Valencia (Spain), Energy Reports.” 8 (2022) 10395-10408.
  • [3] Khosravani, H. R.; Ruano, A. E.; Ferreira, P. M. A convex hull-based data selection method for data driven models. Applied Soft Computing 2016, 47, 515-533, doi:10.1016/j.eswa.2016.06.028.
  • [4] Silva, S.; Ruano, A. The IMBPC HVAC system: Wireless Sensors and IoT Platform. IFAC-PapersOnLine 2018, 51 1-8, doi:https://doi.org/10.1016/j.ifacol.2018.06.227.
  • [5] Ruano, A. E.; Pesteh, S.; Silva, S.; Duarte, H.; Mestre, G.; Ferreira, P. M.; Khosravani, H. R.; Horta, R. The IMBPC HVAC system: A complete MBPC solution for existing HVAC systems. Energy Build. 2016, 120, 145-158, doi:10.1016/j.enbuild.2016.03.043.
  • [6] Ferreira, P.; Ruano, A. Evolutionary Multiobjective Neural Network Models Identification: Evolving Task-Optimised Models. In New Advances in Intelligent Signal Processing, Ruano, A., Várkonyi-Kóczy, A., Eds.; Studies in Computational Intelligence; Springer Berlin/Heidelberg: 2011; Volume 372, pp. 21-53.

Claims

1. A computer implemented method for community energy management (1000), characterized in that the method comprises the following steps: a) acquisition, by an aggregator computational unit (110), of data d0 that will be used for designing prediction models and appliance profiles, during a period of time, preferably during at least one month;b) design, using a cloud-based server (115) of Radial Basis Function prediction models to use for predicting selected from the group consisting of: global solar radiation and atmospheric air temperature of the community, power consumption and power generation in each house h of the community, as well as operation and power consumption for all monitored appliances by Non-Invasive Load Monitoring models, wherein said design comprises the following substeps: b1) selection of the data available for design from the data do acquired in step a);b2) minimization of the model objectives through a Multi-Objective Genetic Algorithm;b3) analysis of the results of the first design, restricting some of the objectives minimized in step b2) and redefinition of the input space and/or the model topology;b4) utilization of the preferential set of models obtained from a second Multi-Objective Genetic Algorithm execution;b5) determination of a second set of the preferred models mepref obtained in the previous substep through the following equation:

2. The computer implemented method for community energy management (1000) according to claim 1, characterized in that data d0 acquired in step a) is selected from the group consisting of: tws, time basis for weather data;R, global Solar Radiation;T, atmospheric air temperature; and, for every house h:tCons, time basis for consumption data;Ph, active power of house h;Rh, reactive power of house h,tInv, time basis for inverter data;{PhPV}, alternate power generated by the PV installation of house h,{Oh}, daily occupation of house h; and, for every house h that has manually controlled appliances:StPh,j/EndPh,j, j∈{MCh}, start/end times of the operation of all manually controlled appliances of house h;{Ph,j, j∈{MCh}}, active Power of all Manually Controlled appliances of house h;tMCh, time basis for all Manually Controlled appliances of house h; and for every house h that has Computer Controlled Appliances that are not air conditioning systems:{Uh,j, j∈{CCh}} control actions for all appliances of house h to be computer controlled;{Ph,j, j∈{CCh}}, active power of all appliances of house h to be Computer Controlled;tcCh, Time basis for all Computer Controlled appliances of house h; and for every house h that has Computer Controlled air conditioning systems:{IThroom, room∈{ACh}}, Air temperature for all air conditioning systems of house h to be Computer Controlled;{RHhroom, room∈{ACh}}, Relative humidity for all air conditioning systems of house h to be Computer Controlled;{WThroom, room∈{ACh}} Wall temperature for all air conditioning systems of house h to be Computer Controlled;{Phroom, room∈{ACh}}, Active Power for all air conditioning systems of house h to be Computer Controlled;{Uhroom, room∈{ACh}}, Control Actions for all air conditioning systems of house h to be Computer Controlled;tAch, Time basis for all air conditioning systems of house h to be Computer Controlled; and for every house h that has Monitored but not Controlled appliances:{StPh,j/EndPh,j, j∈{Mh}}, Start/end times of the operation of all Monitored appliances of house h;{Ph,j, j∈{Mh}}, Active Power of all Monitored appliances of house h;tMh—Time basis for all Monitored appliances of house h.

3. The computer implemented method for community energy management (1000) according to claim 1, characterized in that the initial values v0 of step b) are at least one selected from the group consisting of, for all houses with a photovoltaic energy generation subsystem: PhGrid−/PhGrid + contracted bought(−)/sold(+) power from/to an electrical energy grid by each house h; and for all houses with a photovoltaic energy generation subsystem and battery:PhB−/PhB+ maximum discharged/charged power of a battery bank of at least one house h;SoCh/SoCh minimum/maximum value of the state of charge of a battery bank of a house h;SOCh,initial initial state of charge of a battery bank of a house h;RUh/RDh maximum values of ramp-up/ramp down rates for charging or discharging of a battery bank of a house h;ηhB−/ηhB+ charging/discharging efficiency of a battery bank of a house h;ĒhB rated capacity of a battery bank of a house h; and, for all houses with computer-controlled air conditioning systems:ts start time when a room must be in thermal comfort;te end time when a room must be in thermal comfort;Δt time interval before the start of the schedule when an air conditioning system needs to operate to ensure thermal comfort throughout a whole period;ΘT predicted mean value threshold.

4. The computer implemented method for community energy management (1000) according to claim 1, characterized in that the prediction horizon is between 1 hour and 30 days, preferably between 2 hours and 15 days, and most preferably the prediction horizon is 1 day.

5. The computer implemented method for community energy management (1000) according to claim 1, characterized in that weather data d4 is retrieved from a weather database and/or data d4 is sent by a weather station (120).

6. The computer implemented method for community energy management (1000) according to claim 1, characterized in that the method further comprises the following substep: h4) reception of data d7 related to climate information about at least one room in at least one house in the community, sent by at least one sensor (150), wherein the at least one sensor (150) is electrically powered; andwherein substep h4) is carried out within step h).

7. The computer implemented method for community energy management (1000) according to claim 1, characterized in that the sampling interval Δk belongs to the interval from 5 minutes to 1 hour, most preferably the sampling interval Δk is 15 minutes and that the sampling interval Δk1 belongs to the interval between 10 seconds to 2 minutes, most preferably the sampling interval Δk1 is 1 minute.

8. The computer implemented method for community energy management (1000) according to claim 1, characterized in that the computer-controlled devices are selected from the group consisting of air conditioning systems (180), battery banks (190), swimming pool pumps (200), or any other electric device capable of being controlled remotely, namely using a smart plug (210).

9. A community energy management system (100) comprising: a source of weather-related data d4 wherein such data d4 is at least one of the group consisting of air temperature, relative humidity, global solar radiation;an aggregator computational unit (110) for community energy management system comprising: a communication unit (11), configured to communicate with one or more data acquisition systems (140), one or more consumption energy meter (130), one or more inverter energy meter (160), one or more sensors (150), one or more smart plugs (210);a processing unit (12);a data storage unit (13), configured to store data before, during or after processing by the processing unit (12);at least one energy generation photovoltaic subsystem, comprised of photovoltaic panels and one inverter;at least one battery bank (170), configured to store energy produced by the at least one energy generation photovoltaic subsystem;at least one data acquisition subsystem (140) installed in each house h of the community;at least one consumption energy meter (130) configured to measure the energy consumption at each house h;at least one inverter energy meter (160) configured to measure the energy produced at each house h;at least one smart plug (180), configured to switch on and off electrical equipment connected thereof;at least one cloud-based server (115);

10. The community energy management system (100) according to claim 9 characterized in that it further comprises: one or more indoor sensors (150) configured to measure at least one quantity from the group consisting of air temperature, relative humidity, wall temperature.

11. The community energy management system (100) according to claim 9, characterized in that it further comprises: the source of weather-related data d4 is a weather station (120) configured to measure d4 wherein such data d4 is selected from the group consisting of air temperature, relative humidity, global solar radiation.

12. The community energy management system (100) according to claim 9, characterized in that the aggregator computational unit (110) is selected from the group consisting of a server, a desktop computer, a laptop computer, or a combination thereof.

13. A computer program product characterized in that it comprises logic instructions which, when executed by a computer of a community energy management system (100), cause the computer to carry out the method of claim 1.

14. A computer-readable storage medium characterized in that it comprises instructions which, when executed by a computer of the community energy management system (100), cause the computer to carry out the method of claim 1.