A general MILP based optimization framework to design Energy Hubs

2.1 The framework

The proposed framework has four basic characteristics: object/component-oriented modelling, symbolic constraint definition, interchangeable optimization methods and GUI supported analysis.

The object-oriented modelling in general and the component-oriented modelling in particular permit reusability, interchangeability and collaborative work in the modelling process. The proposed framework applies two concepts of object orientation (see Fig. 1): abstraction and encapsulation (represented by the boxes) as well as inheritance of functionality (represented by the arrows). The basic element in the modelling process of an Energy Hub is the Model, which is a class containing Attributes (representing parameters, variables, results) and methods needed for initialization, preprocessing and postprocessing. Attributes provide metadata needed for the automatic processing and thus are the second important class. Connectors are interfaces to the environment of the model. Every model except the system model has at least one. Nodes represent the connection of at least two Connectors. All these classes are subclasses of TreeNode (arrows in Fig. 1 are pointing to it) enabling them to access other objects inside the resulting data structure.

During initialization of every Model, the used components (i. e., converters; Fig. 1: gas engine) and constraints are defined. Furthermore, Connectors and connections are specified. Inside the method’s preprocess and postprocess the necessary data processing is done before and after the optimization.

After the different models are initialized, the overall system model is built in the same way by adding and connecting instances of component models. The resulting data structure is a tree, in which the root node is the system object providing important data for all other nodes (e. g., number of time steps and time increment). Inner nodes of this tree are of type Model, Connector or Node and leaves are of type Attribute. Thus, every attribute of a model and the attributes of its connectors and nodes can be addressed by a unique and clear path.

Figure 1

UML class diagram representing the basic structure of the software framework.

Representing the Energy Hub structure with a tree has five major features: First, the depth of the tree is not limited as well as the size of the system model. Second, components can be clearly mapped, even if they are used more than once in a single system (e. g., several detailed building models in a district). Third, the unique and clear path prevents naming conflicts of the Nodes. Fourth, it allows to reduce the complexity of individual models (e. g., a piecewise linearization model can be included in different component models as a child). Fifth, the tree structure supports an automated building of the equation system.

Along with the symbolic formulation of the constraints the process of setting up the constraint matrices of the overall Energy Hub is simplified. Furthermore, a symbolic formulation is easier to read which increases the user-friendliness of the software framework. The design of the framework and how the system model is built allows to implement different optimization methods (e. g., Benders decomposition, time slice method). Hence, the proposed framework is not only able to optimize different configurations and system sizes but also allows a variability in the optimization method. This illustrates the flexibility of the proposed framework. The GUI allows a straightforward configuration and evaluation of the optimization results of the investigated system. Furthermore, it supports the user during debugging.

2.2 Objective function

The goal of the optimization is to determine the optimal design and operation of the components within the Energy Hub based on an objective function. In this paper the objective function maximizes the net present value which depends on the nominal capacities of the investigated components. The calculations are oriented towards the standard VDI 6025 [26]. The net present value is the sum of the present values of the revenues Rpv reduced by the sum of the present values of the costs Cpv for the number of investigated components ncomp.

The objective function of the system (root object) collects the necessary values from the components implemented in the Energy Hub. The nomenclature and composition of variables is indicated in Tab. 1 and Tab. 2.

2.3 Economic elements

During the preprocessing the coefficients of the present value of the costs cpv and revenues rpv are calculated on component level. They are used to scale investment costs and annual payments to the according present value. In the following the six economic elements are presented, which are included in the framework.

The present value of investments considers operation and maintenance Cpv,om (2), replacement Cpv,repl (4) and the remaining value RVpv (5) at the end of the period under review Trev. The present value is based on the interest rate ic, general rate of expected change in prices pg (i. e., inflation rate), the replacement factor frepl, the remaining value factor frv, the time of replacement krepl (3), the number of replacements nrepl and the depreciation period Tdp. The coefficient of the present value of investments cpv,inv can be determined in the preprocessing because the investment cost can be extracted from the following equations and is defined by (7).

The coefficient of the present value of annual payments includes the interest rate ic and a rate of price change p, which can be either a general rate of change in prices pg or change in energy cost pe. The costs from energy procurement C and the revenues R from energy sales are annualized (denoted by index a) to calculate the present value. In the following the coefficient of the costs is shown. Revenues are calculated accordingly.

A piecewise linear function is used to model the investment cost of the electrode boiler and the hot water storage (use case 1 and 2) as well as the gas grid fees. The method is based on special ordered sets of type 2 (SOS2), which is described in [27], [28]. The grid points represented by xgp,pwl and ygp,pwl are calculated in the preprocessing using an adaptive method. The variables xpwl and ypwl can be understood as the “input” and “output” of the model. The interpolation between two of the ngp grid points is done by vector λpwl (of size ngp) in (9)–(11). The following constraints force λpwl to be a SOS2 using the Boolean vector bpwl (of size ngp−1).

The power grid fees are calculated according to [29]. The annual grid fees consist of the power related (16) and the energy related costs (17). For each coefficient cp,pgf and ce,pgf two prices exist which depend on the ratio of annual energy to peak power. To avoid this nonlinear problem, the coefficients have to be chosen manually and the validity has to be checked during postprocessing. The peak load is defined by (18) and (19) as well as the annual exchanged energy by (20) in which Pout,pgf and Pin,pgf depict the power that is imported from or exported to higher voltage levels.

The gas grid fees are calculated according to [30]. In contrast to the power grid fees two piecewise linearization components define the power and energy related costs. The peak load is defined by (22) and (23) as well as the annual transmitted energy (24) in which H˙out,ggf and H˙in,ggf depict gas that is imported from and exported to upper gas networks. Because gas grid fees are based on the higher heating value and the chemical power flow of the connectors is based on the lower calorific value, fu,ggf is used in (25) and (26) as a conversion factor in the piecewise linearization.

The energy trading is modelled considering the current spot market price ceex in (28) and (29) as well as the energy related trading fees ctrad,eex in (30). Because trading of gas is based on the higher heating value [31] a conversion factor fu,eex is used to convert from lower to higher heating value in (28)–(30). Furthermore, a change in energy prices is only taken into account in cpv,eex. In the following equations X represents either P or H˙ depending whether power or gas trading is regarded. The overall cost and revenue of sold energy are determined with (31) and (32).

2.4 Components

In this chapter twelve available Energy Hub components and their respective constraints are described. The static energy demand and supply is simply described on component level by assigning the time series given by parameter Pparam to the connector variable Pconn.

The photovoltaic system and windfarm are modelled as scalable plants. The design variable is the number of subsystems nres (i. e., PV modules or wind turbines) that is optimized for a given electrical output vector P1 and a given nominal electrical power Pn,1 using the specific investment costs cinv. Using less or equal instead of equal in (34) enables curtailment, which is needed to limit the power feed-in (37). Eq. (37) can be added to the overall model if needed (e. g., PV system for houses). Here X represents the maximum allowed infeed:

The gas engine is modelled considering constant electric and thermal efficiencies ηel,ge and ηel,th in (38) and (39). Moreover, generation limits are implemented with a Boolean vector bge as an indicator of the operating state in (40) and (41). The constraints in (42) and (43) restrict the minimum uptime ton,min,ge and downtime toff,min,ge according to [32].

Figure 2

Schematic representation of the peak load boiler.

The peak load boiler ensures the required thermal power and supply temperature in a district heating system. During pre-processing, the specific enthalpies hplb of inflow (in), outflow (out) and return (r) pipes are calculated from its corresponding temperatures given as parameters (see Fig. 2). With these specific enthalpies it is possible to convert the thermal power flows of the connectors to mass flows inside the component ((44) and (45)) and finally to model a valid system behavior ((46) and (48)) using a bypass (bp). Eq. (47) represents the overall energy balance taking into account the efficiency of the peak load boiler ηplb.

The hot water storage is modelled by balancing the inflow and outflow of thermal energy for every time step (49) taking into account the initial state of charge SOC0,hws. The nominal capacity Qn,hws restricts the currently stored energy Qhws (50). The model assumes a stratified storage which has a constant head and bottom temperature. The resulting temperature lift defines the volumetric energy density qvol,hws. The investment costs are modelled by a piecewise linearization in (51).

The battery storage has two design variables: the nominal electrical power Pn,bat and the nominal electrical capacity En,bat. The former is modelled as a limitation of the charging and discharging power by (53) and (54). The latter is modelled by the energy balance in (55) and the minimum and maximum possible state of charge ((56) and (57)). The energy balance considers conversion losses by ηbat, which is the product of the charging/discharging as well as the inverter efficiency. The investment costs consist of power and energy related costs, which are represented by the coefficients cinv,p,bat and cinv,e,bat in (55).

The electrical boiler and electrode boiler are modelled as energy converters with an efficiency (60). Because large scale electrode boilers have a significant nonlinear cost function, as shown in [33], the piecewise linearization component is used in this model (62). The investment costs are represented by yeb in (63).

The heat pump model represents an air-water heat pump. The coefficient of performance COPhp is calculated during the pre-processing using a vector of the ambient temperature, the temperature difference between ambient and evaporation temperature, the temperature difference between condensation and supply temperature as well as a linear adjustment of the resulting Carnot coefficient to the real COP. The energy balance is given in (64). The nominal heating capacity is implicitly defined by limiting the current heating power in (65). Eq. (66) and (67) are semi-continuous due to the Boolean variable bhp, which states if the heat pump is installed or not. Hence, the nominal heating capacity can be zero or inside the specified range.

The electrolyzer converts power to hydrogen. In (69) the ratio of the lower calorific value of hydrogen Hi,h2 to the specific electrical energy demand eelz is used as the conversion efficiency. The nominal electrical power is implicitly defined by (70).

The fuel cell converts hydrogen to electrical power. This model does not consider waste heat of the fuel cell. The constraints are equal to those of the electrolyzer except the conversion efficiency, which is formed by the ratio of the specific hydrogen demand efc to the lower calorific value of hydrogen Hi,h2.

The compressor model allows to take the electrical energy demand into account, which is necessary because the hydrogen needs to be compressed before utilization. The specific energy demand ecmp for compression from 30 bar to 450 bar is about 180 to 200 Wh/Nm³ according to [34], [35].

The hydrogen storage is a basic storage model containing the energy balance for every time step (79) and the sizing constraint (80).

The hydrogen refueling station which depicts the hydrogen load is modelled based on an average daily hydrogen demand per car mhfs,car,d and a number of cars refueling at the station per day nhfs,car,d. With the aid of a relative daily demand profile m˙rel, the absolute demand profile can then be created. The chemical energy flow H˙hfs is finally calculated using the mass specific lower heating value of hydrogen Hi,h2, which is 33 kWh/kg.

2.5 Implementation of constraints

The implementation of the abovementioned constraints is shown in Fig. 3 with the example of (55). The constraint is implemented by the symbolic constraint definition. Here, the constraint contains inline functions (purple), parameters (blue), connectors (green) and variables (red). The inline-function “tril” is a lower triangular matrix of ones based on the number of time steps nt. Finally, the constraint matrix and vectors are generated by the framework.

Figure 3

Exemplary illustration of the implementation of one constraint (i. e., battery system) in the case of three time steps.

3.1 Use case 1 and 2 (city)

In use case 1 and 2 an Energy Hub is set up to meet the heat and hydrogen demand of an entire city, while the electric power demand is met by the grid. The used components are displayed in Fig. 4 (top). The Energy Hub buys the gas and electricity which is needed to fuel the components and sells the heat, hydrogen and electricity. The electric and district heating demand are taken from [36]. The hydrogen demand is determined by an hourly variation taken from [38] and 30 cars with a daily demand of 0.6 kg (around 40 km [8], [38]) each. The Energy Hub is connected to the power and gas grid. In both use cases time-varying energy prices [39], [40] are considered. Furthermore, in use case 1 grid fees [29], [30] are included, whereas use case 2 has no grid fees. It has to be noted that the gas engine and peak load boiler will not be designed. Only their operational mode will be optimized, since it is assumed that both components are already existent in the investigated energy system. The optimal size of the other components is shown in Tab. 4. Depending on the parameterization in GUROBI the computation time amounts to several hours.

The net present value of both Energy Hubs is positive (use case 1: 19.3 mio. €; use case 2: 8.4 mio. €) and thus both cases generate a greater benefit than the interest rate. Yet, use case 1 generates a greater benefit, which arises mainly from savings in the grid fees and energy cost. In use case 1 the electrode boiler is designed to reduce the maximum power which is fed into the power grid (see Fig. 5, middle left at time step 860). The black line indicates a reference case which contains only a gas engine and peak load boiler to meet the heating demand. Thus, the maximum power, which is exchanged with the grid, is limited to 20 MW resulting in a reduction of power related grid fees. Because of the electrode boiler the gas engine and peak load boiler run less time leading to a reduction in gas drawn from the grid compared to the reference case (see Fig. 5, middle right). In use case 2 no incentive is given to reduce the exchange with the power or gas grid. Hence, the electrode boiler is not used and heat is supplied by the gas engine and peak load boiler only. This in turn affects the size of the hot water storage, which is much smaller than in use case 1, because the generation adapts to the load. Moreover, using a battery is not reasonable because of the great investment cost. The fuel cell is not used either because of its low efficiency [41] and the overall inefficiency of producing hydrogen to later generate electricity.

It can be seen that including grid fees in the optimization has a great effect on the results because they provide a great incentive to optimize the operation of the Energy Hub.

Figure 5

Power and gas exchange with grids compared with reference case and electric power and gas consumption for use case 1 (middle) and use case 2 (bottom).

3.2 Use case 3 and 4 (single-family house)

In use case 3 and 4 an Energy Hub is set up depicting a single-family house. The used components are displayed in Fig. 4 (middle and bottom). The electric, hot water and space heating demand for a family (2 working people and 2 children) are taken from [42]. Both Energy Hubs are connected to the power grid. Furthermore, in use case 3 the Energy Hub is connected to the gas grid, while in use case 4 the heat demand is met via the electric power grid. In contrast to the use cases 1 and 2 (37) is added to the PV system to limit the PV grid infeed. The optimal size of the components is shown in Tab. 5. Depending on the parameterization in GUROBI the computation time amounts to several minutes.

The heat pump meets the space heating demand (temperature level at 45 °C) and the electric boiler supplies the hot water demand (temperature level at 70 °C).

The net present value of both Energy Hubs is negative (use case 3: −52470 €; use case 4: −73430 €). Hence, both systems are not economic compared to the interest rate. Yet, the absolute net present value of use case 3 is smaller, compared to use case 4. The poorer net present value in use case 4 arises because of a higher electric load (due to the heat pump and electric boiler). Hence, on the one hand more electricity needs to be supplied by the grid (8405 kWh to 4079 kWh, see Fig. 6 left, light blue bars) leading to a greater energy costs (2521 €/a to 2458 €/a) although no gas is needed. On the other hand, more PV generation is used directly. This leads to less PV generation which is fed into the grid in use case 4 (6760 kWh to 10570 kWh, see Fig. 6) resulting in a smaller annual remuneration (745 € to 1163 €). Consequently, the higher energy cost and a smaller income from the PV system and a greater investment cost lead to the poorer performance of use case 4. Though, the combination of a heat pump, electric boiler and a hot water storage adds more flexibility to the electric load resulting in a curtailment of 0 kWh/a compared to 400 kWh/a in use case 3 (see Fig. 6 right, dark red bars).

Apart from this, in both use cases the battery is not used because of its great specific cost and short lifetime. Increasing the lifetime to ten years would not change the presented results.

Figure 6

Electric load (positive values) and generation (negative values) for a winter (left) and summer day (right) for use case 3 (top) and use case 4 (bottom).

3.3 Discussion

The investigation of the four use cases demonstrates that the proposed framework is able to simplify the configuration and analysis of different Energy Hubs. Through its GUI, the configuration of the Energy Hub can be customized and with its object-oriented modelling and symbolic constraint definition the framework is capable to build the constraint matrices automatically and connect the components to their respective energy carriers.

This leads to great time savings in the analysis of Energy Hubs. The following statements are based on experiences during the development of the proposed framework and the implementation of many use cases. If a user wants to implement a use case in an LP or MILP it takes several days or up to weeks to code the constraint matrices, test the system and fix bugs. Applying the proposed framework, the user only needs a few minutes or hours to create a use case depending on the complexity of the investigated Energy Hub. Then the automatic building of the constraint matrices lies between seconds and several minutes depending on the size of the system and the number of time steps. The development and implementation of a new component model takes several hours depending on the degree of detail. Beside the time savings in the creation of use cases the possibility to apply small changes to the investigated system easily (i. e., adding constraints or components) has also a great benefit.

However, the investigation of use case 1 and 2 shows that the proposed framework is not able to solve particular system complexities in an appropriate computation time. Though, this does not necessarily arise from the investigated system size (see Tab. 6). Although, the number of constraints is greater in use case 4 the computation time was more than ten times faster than in use case 1 and 2. Hence, it is assumed that the complexity of the system has a greater impact than the number of constraints. Since more components are implemented the solution space expands and a greater number of possible configurations and component sizes arise. Furthermore, the computation time in use case 1 is more than two times faster than in use case 2. The difference between both use cases is the main incentive (use case 1: grid fees; use case 2: only time varying energy prices). Apparently, the grid fees present a much stronger incentive leading to a better behavior of the optimization problem. This leads to the presumption that besides the number of constraints the potency of the incentives has a great impact on the complexity and computation time of the optimization problem.

Therefore, to investigate complex systems with different incentives the complexity of the optimization problem needs to be reduced. Here, two options are seen as promising: a method using time slices and Benders decomposition. First, using time slices splits the optimization period into time slices (e. g., months) and performs firstly the design optimization and then the operational optimization. Second, the widely used Benders decomposition [43] can be applied splitting the original optimization problem into a master and different subproblems [44]. The master problem determines the optimal design and the subproblems define the optimal operation of the components. Hence, Energy Hubs consisting of a large number of components can be planned and scheduled on full annual simulations.