A Matheuristic Approach for Resource Scheduling and Design of a
Multi-energy System
A. Bartolini
1
, G. Comodi
1
, F. Marinelli
2
, A. Pizzuti
2
and R. Rosetti
2
1
DIISM - Dipartimento di Ingegneria Industriale e Scienze Matematiche, Universit Politecnica delle Marche,
Via Brecce Bianche, I-60131 Ancona, Italia
2
DII - Dipartimento di Ingegneria dell’Informazione, Universit Politecnica delle Marche,
rosetti@dii.univpm.it
Keywords:
Multi-energy System Management, District Design, MILP, Matheuristic.
Abstract:
Modern energy system are evolving due to the opportunities and challenges that new technologies pose in the
energy sector. These changes create the requirements of decision tools able to effectively sustain the processes
of design and retrofit of energy systems. In this paper a multi-energy system management problem is taken into
account and a mixed integer linear programming (MILP) formulation is proposed to model both the design and
the resource scheduling of energy districts. However, since the size of the formulation restricts its applicability
to small cases far from the application of interest, a matheuristic based on constraint relaxations and variable
fixing has been designed. Preliminary computational results show that the proposed solution strategy is able
to achieve good solutions (i.e., solutions with small optimality gaps) on restricted random instances, and to
solve in reasonable times instances derived from a real case study.
1 INTRODUCTION
The push towards a more rational and sustainable use
of primary energy resources is increasing worldwide
and many countries and organization take concrete ac-
tions and commitments in addressing such challenges,
e.g. (European Commission, 2012), (BNEF, 2017).
Among those actions, a common commitment is to
adopt higher shares of renewable energy generation
technologies into the energy mix to satisfy the users
demand for energy and/or other commodities. More-
over, several studies as (European Commision, 2017)
and (ISE, 2015) predict in the near future a strong
reduction of the purchase cost of such technologies
that would induces a dramatic change, with relative
challenges, in existing energy distribution infrastruc-
tures. Most of renewable energy generation technolo-
gies are in-fact characterized by an intermittent and
non-dispatchable electric generation capacity, making
the integration with current national and local distri-
bution grids challenging (California ISO, 2012).
A possible solution lies in regional large central-
ized energy storage systems, such as pumped hy-
dro, compressed air and large battery storage systems.
Such approach has been followed both in many stud-
ies and actual projects (J
¨
ulch, 2016). Another emerg-
ing approach to address such challenging integration
resides in decentralization with the aid of microgrids
(Speer et al., 2015). While microgrids are already a
well established solution to accommodate the needs
of rural off-grid communities, their deployment in al-
ready urbanized contexts of developed countries is
also gaining attention from researchers (Enea, 2017),
(Center for Climate and Energy Solutions, 2017).
The basic idea behind microgrids is strictly re-
lated to electric energy management, but their poten-
tial to handle and optimize the employment of re-
newables and primary energy resources greatly im-
proves in multi-energy systems (MES), that is sys-
tems in which multiple energy vectors (i.e., carriers
of energy as electricity) are concurrently considered
(Mancarella, 2014). To capture the potential opportu-
nities for efficient energy use (polygeneration, (Jana
et al., 2017)), MES models have to explicitly consider
the interactions across different energy vectors, in ac-
cordance with the real scenario under analysis, e.g. a
comprehensive assessment of users demanded com-
modities such as heat, cooling and electric energy.
In order to exploit all the potential energy savings
opportunities, such cross vector interactions have to
be taken into account since the design phase of both
new energy systems and the retrofit of already exist-
Bartolini, A., Comodi, G., Marinelli, F., Pizzuti, A. and Rosetti, R.
A Matheuristic Approach for Resource Scheduling and Design of a Multi-energy System.
DOI: 10.5220/0007574104510458
In Proceedings of the 8th International Conference on Operations Research and Enterprise Systems (ICORES 2019), pages 451-458
ISBN: 978-989-758-352-0
Copyright
c
2019 by SCITEPRESS Science and Technology Publications, Lda. All rights reserved
451
ing settings. Thus, MES models must preferably be
adapted to represent such opportunities and have a
correct assessment of what high integration of renew-
ables might imply.
In this paper a mixed integer linear programming
(MILP) formulation is proposed to model the design
of MES, where the optimal investment choices are
performed by evaluating also the required commodi-
ties from the exogenous energy networks. Specifi-
cally, the formulation aim at minimizing the set-up,
maintenance and procurement costs by taking into ac-
count all the operational constraints, i.e., the schedul-
ing of the operation performed by the available tech-
nologies in order to convert energy vectors. Due to
the inner complexity of the addressed problem and the
size of the mathematical programming formulation,
a matheuristic scheme based on fixing integer vari-
ables and relaxing constraints is designed to produce
good sub-optimal solutions in a reasonable amount of
time. To the best of our knowledge such scheme has
only been applied by (Triad
´
o-Aymerich et al., 2016)
for energy systems related problems, where the opti-
mal design for a rural electrical microgrid is evaluated
using different heuristic approaches, among which a
fix-and-relax strategy. In the last few years the opti-
mal design and operation management of distributed
MES has been addressed by means of several methods
(Singh and Sharma, 2017). The problem has been for-
malized with both single or multi-objective functions,
and solved by techniques that can be mainly clas-
sified into single monolithic and decomposition ap-
proaches. With respect to the former, (Mehleri et al.,
2012) and (Omu et al., 2013) exploit MILP formula-
tions to compute the optimal design of a small district
needing electricity and heat; the same is done in (Bis-
chi et al., 2014) where the demand of cooling power
is also considered. For the latter, decomposition ap-
proaches rely on keeping the planning phase, that is
the decision to purchase and deploy, separated by the
operational phase. In (Elsido et al., 2017) a mixed in-
teger non linear programming (MINLP) formulation
is presented and solved by a two phases approach,
where the first is handled by means of two evolu-
tionary algorithms and a discrete variable relaxation
technique, while the second is solved by exploiting
a linearized MILP. In (Li et al., 2017) the design of
an hydrogen based microgrid is performed by using
a genetic algorithm to determine the size of the sys-
tems, while the operational variables are set by solv-
ing a MILP. A similar technique is described in (Sachs
and Sawodny, 2016), where the same framework is
applied for a multi-objective variant of the problem.
The rest of the paper is structured as follows: in
Section 2 the approach to the energy hub and tech-
nologies modeling is described, in Section 3 the prob-
lem is formalized and the mathematical model is pre-
sented. Section 4 describes the matheuristic used; in
Section 5 preliminary computational results are pre-
sented, based on both randomly generated instances
and a real case study. Finally in Section 6 conclusions
and future perspectives are presented.
2 ENERGY HUB AND
TECHNOLOGY MODELLING
In this section the concept of energy hub and the spe-
cific setting related to the case study in §5.2 are de-
scribed. Furthermore, the technologies that are con-
sidered within the energy hub and their general func-
tion are depicted.
2.1 The Energy Hub
Among the suitable tools for the management of
MES, a relevant option can be found in energy hub
(Mohammadi et al., 2017). An energy hub can be
thought as an integrated unit within which the con-
version and storage of different energy carriers is un-
dertaken in order to satisfy different commodities de-
manded by an user. The energy hub can get some of
such energy vectors from outside of its boundaries,
as for example the withdrawal of electric energy or
natural gas from the respective national energy distri-
bution infrastructures. Then, within the boundaries
of the energy hub such vectors can be transformed
in order to meet the demand of users, this provided
that the technology needed to operate the conversion
is installed. Such technologies can either be already
available or be deployed within the energy hub by sus-
taining the relative purchase and maintenance cost.
Thus, the choices regarding which kind of technolo-
gies to deploy within an energy hub have to be as-
sessed by considering both the monetary investment
needed to purchase such technologies and how such
systems will work, providing the needs of the com-
munity of users tied to the energy hub.
Energy hub is then a general scheme that can be
adapted to various settings and the case study pre-
sented in §5.2 reflects the structure presented in Fig-
ure 1. It serves a community consisting in buildings of
the tertiary sector, which will need to properly func-
tion electric energy, space heating for the winter and
cooling for the summer. It has access to both the nat-
ural gas network and the national electric energy grid;
also has access to a solar resource to be used to power
the district.
ICORES 2019 - 8th International Conference on Operations Research and Enterprise Systems
452
Figure 1: Proposed energy hub model.
2.2 Technologies Modeling
The technologies deployed in the energy hub convert
energy among the different vectors, or store it for a
later use. The technologies considered in this study
include the majority of the most mature ones available
on the market for small urban districts. These can be
listed as:
Natural Gas Engine (CHP): operates in com-
bined heat and power (CHP) mode, thus produc-
ing electricity and simultaneously recovering the
waste heat, using natural gas as fuel;
Natural Gas Boiler (GB): production of heat
from natural gas;
Electric Chiller (EC): conversion of electric en-
ergy in cooling power;
Absorption Chiller (AC): conversion of heat in
cooling energy;
Heat Pumps (HP): allows the conversion of elec-
tric energy into (alternatively) cooling energy or
heat, based on the required output;
Photovoltaic Panels (PV ): production of electric
energy from solar radiation;
Batteries (EES): storage of electricity;
Hot Thermal Storage (HT ES): storage of heat;
Cold Thermal Storage (CT ES): storage of cool-
ing energy.
A small database of purchasable devices is con-
sidered for each of such technologies, distinguished
by size. Then for each particular model, parameters
representing costs and performance are specified. The
performance of conversion devices is described by the
energy conversion efficiency for each available size,
whereas the relevant parameters for the storage de-
vices are the efficiency of both charge and discharge
phases. Besides, conversion devices are enforced to
work above a fixed threshold value for a certain times-
pan, this to represent a lower bound of partialisation
(i.e., the minimum output power) for the technology.
The only exception concerns the photovoltaic panels,
whose produced electric energy is directly dependent
by the total panel surface and the solar exposition.
Hence, in absence of solar radiation the output of the
photovoltaic can be zero.
Regarding the cost, each model of technology is
characterised by its own purchase cost related to its
size and performance. Such cost has to be sustained
in order to make use of the technology for a represen-
tative lifetime in years, after which the purchase cost
has to be re-sustained. Moreover, additional costs re-
lated to the maintenance operations are considered.
In this way, all of the costs related to purchasing and
operating phases of a certain technology are taken
into account. Finally, further costs sustained within
the energy hub are given by the possible withdrawal
of energy from the two national infrastructures. The
withdrawal is priced with a fixed cost per kWh of pur-
chased energy, different for natural gas and electricity.
3 PROBLEM DEFINITION AND
MATHEMATICAL
FORMULATION
The problem can be summarized as follow: an en-
ergy district has to be designed and managed (by se-
lecting an equipment of technologies and by control-
ling the relative operations) in order to satisfy the de-
mand of different energy vectors (commodities). De-
cisions have to be performed such that the total costs
sustained by the energy district within a certain time
horizon are minimized, both in the planning and op-
erational phases. We assume that the time horizon is
cyclically repeated to describe the optimal dynamic of
a fully operational energy district, satisfying requests
that (hopefully) follow periodical patterns.
Let K be the set of commodities, with a fixed unit
purchasing and storing price, that is respectively c
k
and h
k
. Then, let us define Q as the set of technol-
ogy models that can be deployed in the district, e.g.
heat pump 170kW or electric chiller 120kW. More-
over, let S Q be the subset of all the storage systems
and S
k
S that of models that store the commodity
k. Parameters v
i
are defined for each i Q, stating
the costs for buying, installing and maintaining the i-
th appliance type. This cost is normalized along the
considered period: indeed, the approximated profile
of the demand of each commodity is reported for a
discrete time horizon T = {1,...,
˜
T } with intervals of
one hour; namely d
k
t
is the demand of the commodity
k K at time t T .
A Matheuristic Approach for Resource Scheduling and Design of a Multi-energy System
453
Therefore, the matrix v is normalized in this way:
given the purchasing and installing cost (IC), main-
tenance cost per year (MCY) and estimated lifetime
(expressed in year and hours respectively as ELY and
ELH), each element of v is defined by:
v
i
=
IC
i
+ MCY
i
× ELY
i
ELH
i
×
˜
T
Although fixed costs should be fully sustained in the
same moment, the definition of v
i
includes only the
portions of hourly costs related to the considered time
horizon T . Hence, the ratio between fixed and vari-
ables costs (e.g. c
k
) is contained, while the objective
function weights the total cost associated to the life
cycle of each device.
As previously mentioned, each device can require
and supply different types of commodity, according
to his features. For example a Gas Engine uses gas
and can supply heat or electricity, whereas an Elec-
tric Chiller needs electricity to work and supply cool.
A parameter ω
hk
i
is then defined only for conversion
technology i Q \ S, expressing the unitary conver-
sion multiplier to convert commodity h into commod-
ity k by means of appliance i. Some devices, such
as solar panels or wind farms, have a base production
that exploits renewable resources and are dependent
by weather conditions at a certain time. Thus, pa-
rameter b
k
it
states the base production for commodity
k K produced by the equipment i Q \ S at time t.
The amount of output k produced by the conversion
device i at each time t ranges within the maximum
rated power U
k
i
and the lower bound L
k
i
.
Regarding the storage technologies, each device i S
k
is able to collect commodity k up to a fixed capacity
C
k
i
, with a common efficiency φ
k
.
We define non-negative variables r and p to model
the flow of commodities converted by conversion de-
vices: at each time instant t, variable r
hk
it
R
+
ex-
presses the amount of commodity h converted into
commodity k by device i Q \ S; on the other hand,
variable p
k
it
R
+
models the output commodity k re-
sulting by the conversion operation related to r
hk
it
and
alternative base productions b
k
it
(e.g. solar energy).
Furthermore, for each k K and t T we define vari-
ables l
k
t
and f
k
t
, to describe for the time t respectively
the total amount of commodity k stored within the dis-
trict, and the quantity of k acquired by the exogenous
supplier. By assuming that each technology has an
estimated lifetime much larger than
˜
T , we can use bi-
nary variables y
i
i Q to model the selection of de-
vices, i.e. y
i
= 1 if and only if the district is provided
with equipment i since the beginning of T . Finally, for
conversion technologies i Q \ S and for each t T ,
variables z
it
{0, 1} are further employed to describe
the operation of device i at time t through z
it
= 1.
The proposed MILP formulation reads as follow:
min
tT
kK
(c
k
f
k
t
+ h
k
l
k
t
) +
iQ
v
i
y
i
(1)
d
k
t
+ l
k
t
+
iQ\S
hK:
h6=k
r
hk
it
=
φ
k
l
k
t1
+
iQ\S
p
k
it
+ f
k
t
k K,t T (2)
p
k
it
=
hK:
h6=k
ω
hk
i
r
hk
it
+ b
k
it
y
i
i Q \ S,k K, t T
(3)
L
k
i
z
it
p
k
it
U
k
i
z
it
i Q \ S,k K, t T
(4)
tT
z
it
|T |y
i
i Q \ S (5)
l
k
t
iS
k
y
i
C
k
i
k K,t T (6)
l
k
0
= l
k
˜
T
k K (7)
l
k
t
, f
k
t
0 k K,t T (8)
p
k
it
0 i Q \ S,k K, t T (9)
r
hk
it
0 i Q \ S,h, k K : h 6= k,t T
(10)
y
i
{0,1} ∀i Q (11)
z
it
{0,1} ∀i Q \ S,t T (12)
The objective function (1) aims to minimize the
total cost given by the cost related to deployed tech-
nologies, the cost of acquiring external commodi-
ties from public suppliers and the expense for storing
surpluses for further use. Constraints (2) are stock-
balancing constraints and state that the commodity k
demanded, stored and required by devices is equal to
the summation among the previously amount of com-
modity k stored at time t 1, the total one produced by
all deployed technologies and the quantity acquired
from external networks. The set of constraints (3)
defines the variables p, for each conversion technol-
ogy i Q \ S , as the summation between the required
commodity h necessary to produce commodity k and
a base production without consumption. This is due to
the feature of green technologies, such as photovoltaic
and wind farm, to produce commodities without con-
suming other resources. The bounds of variables p are
given by the set of constraints (4): if binary variable
z
it
is equal to 1, the corresponding variable p
k
it
is lim-
ited between values L
k
i
and U
k
i
. The link among vari-
ables z and y is coherently modelled by the set of con-
straints (5). Constraints (6) limit the stored quantity
of commodity k K at any time t T to the total sum
of capacities of the storage systems in S
k
installed. In
order to translate the cyclicality of the selected time
ICORES 2019 - 8th International Conference on Operations Research and Enterprise Systems
454
horizon, the quantities of commodities stored in the
first timestep (t = 1) are balanced in (2) with the quan-
tities stored in the last timestep (t =
˜
T ) by means of
equality (7). It is worth noting that whenever an ap-
pliance i is not able to convert a certain commodity h
into k, the corresponding parameter ω
hk
i
is null. In all
such cases, t T variables r
hk
it
are fixed to zero as a
preprocessing phase. Moreover, t T variables f
k
t
are set equal to zero for any commodity k K that is
not purchasable by external suppliers.
4 MATHEURISTIC APPROACH
Computational experiments showed that the mathe-
matical formulation proposed in §3 cannot be ex-
ploited to solve real-size instances in a reasonable
amount of time by means of standard solvers. This
can be (partially) related to the size of the formula-
tion, as the number of constraints is O(|Q||K||T |) and
the number of binary variable is O(|Q||T |), whereas
continuous variables are O(|Q||K|
2
|T |). Therefore,
we design a matheuristic algorithm to approach in-
stances that are meaningful in real context.
The matheuristic used to solve the proposed MILP
is based on relaxing integrality constraints and fix-
ing subsets of variables, similarly to the Fix & Re-
lax scheme formalized in (Escudero and Salmeron,
2005). Starting with the linear relaxation of the math-
ematical formulation, the idea is to restrict it by restor-
ing the integrality conditions only for a subset of the
original integer variables. Then, the restricted linear
relaxation is solved and the integer variables are fixed
to their values in the resulting solution; hence, the
number of variables is reduced and the space of re-
search shrunk. This process is iterated until a feasible
integer solution for the original problem is found.
A common challenge for this kind of approach
is to find an equilibrium between the algorithm effi-
ciency and quality of the solution. Indeed, while the
former usually requires to proceed by small subset of
variables, the latter could be compromised by such
scheme as the first steps can bring the current solution
far from the optimum. Moreover, a robust strategy is
one that can avoid the convergence to an unfeasible
solution, or is able to efficiently backtrack whenever
this occurs.
The matheuristic approach is actually a cus-
tomization of the Fix & Relax scheme for the energy
district design and management problem, obtained by
specifying the strategy of selection and fixing of the
subset of variables. In particular, the subsets of bi-
nary variables are obtained by separating the corre-
sponding technologies for their capacity of convert-
ing or storing commodities. Indeed, each technology
rarely produce (store) more than one commodity, thus
creating a natural separation between the relative vari-
ables. Moreover, as detailed below the fixing policy
is defined to allow a certain grade of flexibility to the
matheuristic and to reduce the possibility of reaching
unfeasible states.
The scheme can be summarized as follow. Let us
consider the mathematical formulation proposed in §3
as M . Let V
y
and V
z
be the sets of variables in M for
which the binary condition holds, related respectively
to y and z. In a similar way, let us define V
r
y
and V
r
z
as
the subsets of variables whose integrality is relaxed in
M . Finally, let K
r
be an auxiliary set initialized as K
and collecting all the commodities.
The algorithm starts by considering the linear re-
laxation of M , that is V
y
=
/
0 and V
z
=
/
0. Firstly, the
demands and costs for acquiring external resources
are used to define a priority among the commodities
k K. Specifically, at each iteration the commodity
s K
r
is selected accordingly to the following rule:
s = argmax
kK
r
,tT
{d
k
t
/c
k
} (13)
That is, the commodity chosen is that one with the
largest cost required to satisfy its own demand peak
through the external supply. Commodity s is then re-
moved from K
r
and the integrality is restored for the
variables y
i
V
r
y
modeling appliances able to con-
vert (store) the commodity s. In other words, for all
the conversion technologies i Q \ S such that exists
ω
hs
i
> 0, with h K : h 6= s, the binary condition for
y
i
is restored, that is y
i
is removed from V
r
y
and added
to V
y
; the same is done for the storage systems i S
with positive C
s
i
. Formulation M is then solved and
for each variable y
i
just added to V
y
, its value is fixed
to one if and only if y
i
= 1 in the resulting solution.
Furthermore, only for conversion devices i Q \ S,
the value of y
i
is set equal to zero if y
i
= 0 currently
holds. By excluding the variables related to the stor-
age devices, the flexibility of the current solution is
partially preserved. In fact, if in a later iteration re-
lated to h 6= s the requests r
sh
it
increase due to deploy-
ment of different appliances, then storage devices of
s are still valid options to avoid the increment of f
s
t
.
For each y
i
currently fixed to 1, all the corresponding
z
it
in V
r
z
are fixed to one if z
it
> 1 ε (with ε reason-
ably small), and moved from V
r
z
to V
z
. This process is
iterated until K
r
is empty.
Afterwards, the integrality for all the variables left
in V
r
y
and V
r
z
is restored and all the subsets are coher-
ently updated. In order to reduce the number of binary
variables in V
z
, a further fixing step is performed. For
each i Q and t T , let ¯z
it
be the values of variables
z
it
obtained in the last solution. Each ¯z
it
is compared
A Matheuristic Approach for Resource Scheduling and Design of a Multi-energy System
455
with a threshold θ [0, 1] and fixed to one if ¯z
it
θ,
then M is solved. If the current state is unfeasible,
θ is incremented by a small step α and all the z
it
for
which ¯z
it
< θ are unfixed; then M is solved again.
This is done until M is feasible or it results unfeasi-
ble for θ = 1, case in which the algorithm ends with
an unfeasible state. Although possible, during com-
putational experiments no unfeasibility state has been
reached and a feasible integer solution is achieved for
the starting value of θ. It is worth noting that for each
commodity k K that is not purchasable by an exter-
nal supplier, c
k
= is assumed to preserve the con-
sistency within (13). Moreover, selecting those com-
modities in later steps help to prevent unfeasibility, as
variables f
k
t
are enforced to be zero and cannot bal-
ance the possible increment of requests r
kh
it
, due to the
selection of device i in further steps of the algorithm.
5 PRELIMINARY EXPERIMENTS
The code was implemented in AMPL v.20180308
(MS VC++ 10.0, 64-bit) and experiments were per-
formed on a Intel
r
Core i7-7500U 2.90 GHz with
16Gb RAM. All the MILP were solved by IBM
r
CPLEX
r
12.8.0.0.
The matheuristic performance were tested on two
instances derived from a real case study and five ran-
domly generated instances, all based on the energy
hub scheme depicted in §2: 4 energy vectors uni-
formly measured in kWh are taken into account, and
a database of 37 conversion appliances and 18 storage
devices are considered. For conversion technologies,
the lower bound of partialisation is set equal to 30%
of the maximum rated power. Parameters were set to
ε = 0.1, θ = 0.7 and α = 0.1.
5.1 Random Instances
The five random instances were generated by consid-
ering a time horizon T of 48 time steps, simulating the
requests of commodities within a period of two days.
The values of demands were randomly chosen within
[0, 250] kWh for each time unit t T , whereas solar
radiation was derived by a real profile.
Table 1 shows the results obtained by the
matheuristic (H ) with respect to the optimal solution
provided by CPLEX after exactly solving the MILP
formulation (M ). In all the instances the percentage
optimality gap is smaller then 1%, with an average
value of 0.65%; then the CPU time required by H is
reasonable and exceeds two minutes only in one case,
while being widely smaller than directly solving M
on the mean.
Table 1: Solutions for random instances.
Instance gap% CPU
H
CPU
M
d
1
0.93 62.22 807.48
d
2
0.79 128.46 1592.91
d
3
0.36 83.45 1052.86
d
4
0.82 18.60 1481.12
d
5
0.37 65.42 2017.73
Average 0.65 71.63 1390.42
5.2 Case Study: Tertiary Urban District
In order to validate the proposed method, the data de-
scribing the needs of an existing urban district is used
to derive realistic instances. Simulated data based on
building energy consumption models is used to ob-
tain the demands of a district composed of tertiary use
buildings, specifically ten offices and a school. De-
mands are distinguished in three commodities: elec-
tricity, space heating and cooling. The user demands
are computed as if the buildings were situated in Italy,
specifically in Rome, for a timespan of one year with
a resolution of one hour. The building models are pro-
totypes developed by the U.S. Department of Energy
(U.S. Department of Energy, a) which are then sim-
ulated by means of the software Energy Plus (U.S.
Department of Energy, b). The solar radiation is also
given as an input by referring to an average year ob-
tained from historical data. While the withdrawal of
electric energy or natural gas from the respective na-
tional infrastructures comes with no limit in terms of
maximum capacity, this is not the same of the radia-
tion resource. The availability of this resource is then
represented by realistic data with the same hourly res-
olution as for the commodities demanded.
5.3 Realistic Instances
Given the scenario depicted as case study, two real-
istic instances are derived. The former instance r
1
counts 2160 hour as time horizon of simulated de-
mands, corresponding to three months of data that
cross the first two quarter of the year and have re-
quirements linked to both the winter and summer pe-
riod. The detail of the simulated demands is shown in
kWh in Figure 2; the availability of the solar resource
for a surface with an inclination of 30
in kW/m
2
is
shown in Figure 3. The latter instance r
2
is based on
the yearly demands aggregated within a daily resolu-
tion, so that the time horizon is translated into 365
days. This is done in order to keep manageable the
size of the MILP, while obtaining a reasonable coarse-
grained description of the urban district annual be-
haviour. To preserve the coherence, all technical pa-
rameters of systems are scaled within r
2
. While solv-
ICORES 2019 - 8th International Conference on Operations Research and Enterprise Systems
456
ing r
1
and r
2
, the relative optimality gap tolerance of
CPLEX was set to 2%.
(a) Electric demand
(b) Heating demand
(c) Cooling demand
Figure 2: Three-month energy demands of the district.
Table 2 summarizes the results obtained by the
matheuristic. For the two instances, the total cost
and the CPU time (in seconds) are reported; more-
over, column gap% reports the percentage gap ob-
tained by comparing the heuristic integer solution
with the lower bound resulting by the solution of the
MILP in the first step of H ; that is the linear relax-
ation of M restricted only by the integrality of the vari-
ables y
i
related to the technologies that convert (store)
the first selected commodity. In Table 3 the list of ap-
pliances deployed in each solutions is reported, with
subscripts specifying the maximum rated power of the
models (or rated capacity, in case of storage devices).
The matheuristic was able to achieve an integer solu-
Figure 3: Three-month solar radiation.
Table 2: Solutions for realistic instances.
Instance gap% CPU
H
r
1
22185 12.90 5804.09
r
2
74356 11.51 432.18
Table 3: Selected devices of realistic instance solutions.
Instance Devices
r
1
CHP
140
,EC
39
,EC
52
,EC
70
,EC
120
,EC
160
,HP
160
,HP
240
,PV
2000
HT ES
1000
,CT ES
500
,CT ES
1000
r
2
CHP
140
,EC
39
,EC
70
,AC
83
,HP
160
,HP
240
,PV
500
,PV
2000
HT ES
500
,CT ES
500
tion for the instance r
1
in a CPU time of 5804.09 sec-
onds. This solution has a total cost of 22185e, with
a gap% of 12.90%. For the instance r
2
the heuristic
solution has total cost of 74356efound in 432.18 sec-
onds, with a gap% of 11.51%. Given the dimension
of the instances and the size of the MILP on which the
matheuristic is based (e.g., only the binary variables
in instance r
1
are 79975), the time required to solve
the problem is reasonable and the values of gap% ob-
tained testify that the maximum distance from the true
optimal solution is not excessive.
6 CONCLUSION AND FUTURE
STEPS
In this paper a MILP formulation is described for the
optimal design of energy hubs, aimed at modeling the
energy systems of modern urban districts. The pro-
posed model wishes to represent the possible interac-
tions among different energy vectors, with the goal of
meeting the multiple commodity demands of users.
However, due to the size of the formulation, it is pos-
sible to directly exploit the MILP only to solve small
instances. In order to efficiently solve instances of
more relevant dimension, a matheuristic scheme re-
lying on the MILP has been designed. The strategy
behind the matheuristic is based on fixing the values
A Matheuristic Approach for Resource Scheduling and Design of a Multi-energy System
457
of subsets of binary variables while other integrality
constraints are relaxed, iteratively moving up to a fea-
sible integer solution. Preliminary experiments on a
restricted group of instances show that the matheuris-
tic is able to solve instances with limited time horizon
in restricted time, providing a solution with strict op-
timality gap. Moreover, two instances derived from
real case scenario were solved in reasonable time with
a limited gap from the best lower bound.
The matheuristic presented in this paper can be
seen as a first promising step in approaching a chal-
lenging problem as the management of MES, with a
strategy that can be customized in different ways. In-
deed, elements of a MILP can be divided into sub-
blocks following several criteria. Further work will
be dedicated to refine the proposed framework based
on variable separation and fixing, whereas alternative
methodology (e.g., the decomposition into linked sub-
blocks of constraints) will be explored. Concerning
the formulation, future efforts could aim to include
a more accurate modeling of the technology func-
tion, as for example switching from a performance
based on a single constant efficiency value to a piece-
wise linear efficiency-load curve. Finally, the prob-
lem could be further extended by considering the ex-
istence of multiple separated districts, thus evaluating
the deployment of technologies for connecting differ-
ent districts and modeling the energy exchange.
REFERENCES
Bischi, A., Taccari, L., Martelli, E., Amaldi, E., M. G.,
Silva, P., Campanari, S., and Macchi, E. (2014). A de-
tailed MILP optimization model for combined cool-
ing, heat and power system operation planning. En-
ergy, 74(C):12–26.
BNEF (2017). New Energy Outlook 2017. page 6.
California ISO (2012). What the duck curve tells us about
managing a green grid. Technical report.
Center for Climate and Energy Solutions (2017). Micro-
grids: What Every City Should Know. Technical re-
port.
Elsido, C., Bischi, A., Silva, P., and Martelli, E. (2017).
Two-stage MINLP algorithm for the optimal synthesis
and design of networks of CHP units. Energy, 121.
Enea (2017). Urban Microgrids. Technical Report January.
Escudero, L. and Salmeron, J. (2005). On a Fix-and-Relax
Framework for a Class of Project Scheduling Prob-
lems. J. Ann Oper Res, 140:163–188.
European Commision, J. S. H. (2017). Pv status report
2017. Technical report.
European Commission (2012). Roadmap 2050. Technical
Report April.
ISE, F. (2015). Current and future cost of photovoltaics,
long-term scenarios for market development, system
prices and lcoe of utility-scale pv systems. Technical
report.
Jana, K., Ray, A., Majoumerd, M. M., Assadi, M., and De,
S. (2017). Polygeneration as a future sustainable en-
ergy solution A comprehensive review. Applied En-
ergy, 202:88–111.
J
¨
ulch, V. (2016). Comparison of electricity storage options
using levelized cost of storage ( LCOS ) method. Ap-
plied Energy, 183:1594–1606.
Li, B., Roche, R., Paire, D., and Miraoui, A. (2017). Sizing
of a stand-alone microgrid considering electric power,
cooling/heating, hydrogen loads and hydrogen storage
degradation. Applied Energy, 205(April):1244–1259.
Mancarella, P. (2014). MES (multi-energy systems): An
overview of concepts and evaluation models. Energy,
65:1–17.
Mehleri, E. D., Sarimveis, H., Markatosa, N. C., and Papa-
georgiou, L. G. (2012). A mathematical programming
approach for optimal design of distributed energy sys-
tems at the neighbourhood level. Energy, 44(1):96–
104.
Mohammadi, M., Noorollahi, Y., Mohammadi-ivatloo, B.,
and Yousefi, H. (2017). Energy hub: From a model
to a concept A review. Renewable and Sustainable
Energy Reviews, 80(December 2016):1512–1527.
Omu, A., Choudhary, R., and Boies, A. (2013). Distributed
energy resource system optimisation using mixed in-
teger linear programming. Energy Policy, 61:249–
266.
Sachs, J. and Sawodny, O. (2016). Multi-objective three
stage design optimization for island microgrids. Ap-
plied Energy, 165:789–800.
Singh, B. and Sharma, J. (2017). A review on distributed
generation planning. Renewable and Sustainable En-
ergy Reviews, 76(March):529–544.
Speer, B., Miller, M., Renewable, N., States, U., Schaffer,
W., Gmbh, S. N., Gueran, L., Reuter, A., and Jang,
B. (2015). The Role of Smart Grids in Integrating
Renewable Energy. Technical report.
Triad
´
o-Aymerich, J., Ferrer-Mart
´
ı, L., Garc
´
ıa-Villoria, A.,
and Pastor, R. (2016). MILP-based heuristics for
the design of rural community electrification projects.
Computers and Operations Research, 71:90–99.
U.S. Department of Energy. Commercial prototype building
models.
U.S. Department of Energy. Energy plus v8.9.0.
ICORES 2019 - 8th International Conference on Operations Research and Enterprise Systems
458