*2.1. Optimization Method*

The novelty of the optimization method is that it can solve computationally expensive seasonal storage operation optimization problems with both integer and non-integer operation variables (variables that vary from time step to time step during operation). This means that the method can optimize energy storage systems containing multi-mode devices (e.g., on-off mode devices and reversible fuel cells) where integer variables are used to control operation modes. This type of optimization problem tends to become too computationally expensive for conventional optimization as the number of integer variables increases.

The basic idea of the method is to first perform a rough optimization of the annual operation, and then gradually fine-tune the solution by re-optimizing smaller time intervals. The first step of the method is, hence, to optimize the annual operation using a relatively small number of time steps, *N*, (*N* = 26 time steps are used in the case study). Thereafter, the energy storage load at the beginning, the middle and the end of the one-year period is fixed. The one-year period is then bisected, so that the first half of the year forms one new time period and the second half of the year forms another new time period.

In the second optimization step, the new time intervals are re-optimized separately, using the same number of time steps (*N*) as in the first optimization step, and the storage loads at their center points (not the starting points and the end points) are fixed. Then, the time intervals are bisected, and new intervals are formed. The second optimization step is repeated until the time steps reach the desired size (seven hours in the case study). The time step will reduce in size for every new interval bisection since the number of time steps remains the same, while the time interval size shrinks.

A graphical explanation of the first and second optimization steps are presented in Figure 1. In the figure, mi,j is the energy storage load at the central point of time interval j in optimization step i. Step 1 also shows the storage load at the beginning, a, and the end, b, of the year. To guarantee a net zero annual storage balance, the storage loads at point a and b must be equal to one another. Otherwise, the energy storage cannot operate in a sustainable way for several years.

After each optimization step, the new center points of the new time intervals will be added to the final solution of the annual system operation. This means that the storage load for all the points that are fixed will remain constant throughout the optimization process. Hence, the storage load for 2(n−1) (where n is the number of the optimization step) new points will be added to the final solution from each optimization step, except for the first optimization step, where three points (center point, starting point and end point) are added to the final solution.

Each time interval optimization is solved with a Matlab R2019b solver called intlinprog, which uses a branch-and-bound approach to tackle MILP problems [19]. The solver is also able to identify and eliminate futile subproblem candidates by adjusting the constraints of the optimization problem.

*Buildings* **2020**, *10*, x FOR PEER REVIEW 4 of 27

**Figure 1.** Graphical description of the first and second optimization step in the interval halving optimization method. **Figure 1.** Graphical description of the first and second optimization step in the interval halving optimization method.

#### *2.2. Case Study: Reversible Solid Oxide Cell and Hydrogen Storage System for Seasonal Storage of Solar 2.2. Case Study: Reversible Solid Oxide Cell and Hydrogen Storage System for Seasonal Storage of Solar Energy*

*Energy*  A reversible solid oxide cell (RSOC) is an electrochemical device that can operate either as a solid oxide electrolysis cell (SOEC) or as a solid oxide fuel cell (SOFC). In SOEC mode, the RSOC consumes electricity and water in order to produce hydrogen and oxygen through a redox reaction. In SOFC mode, the redox process is run in reverse, so that the RSOC consumes hydrogen and produces electricity, heat and water [20,21]. An RSOC together with a hydrogen compressor and a compressed hydrogen storage can be used as an energy storage system which is able to store energy in the form of hydrogen gas. Such a system is hereinafter referred to as an RSOC and hydrogen storage A reversible solid oxide cell (RSOC) is an electrochemical device that can operate either as a solid oxide electrolysis cell (SOEC) or as a solid oxide fuel cell (SOFC). In SOEC mode, the RSOC consumes electricity and water in order to produce hydrogen and oxygen through a redox reaction. In SOFC mode, the redox process is run in reverse, so that the RSOC consumes hydrogen and produces electricity, heat and water [20,21]. An RSOC together with a hydrogen compressor and a compressed hydrogen storage can be used as an energy storage system which is able to store energy in the form of hydrogen gas. Such a system is hereinafter referred to as an RSOC and hydrogen storage (RSOCHS) system. The power input capacity of the RSOC in SOEC mode is significantly higher than the power

(RSOCHS) system. The power input capacity of the RSOC in SOEC mode is significantly higher than

more time consuming than charging the storage.

output capacity in SOFC mode [22]. This means that discharging the RSOCHS storage is more time consuming than charging the storage. FutureHub, (Figure 2), located in Espoo, Finland. The building characteristics are presented in Table 1. The location of the building is 60°11'11.4" N 24°48'49.0" E and the Köppen climate classification of the site climate is Dfb [23].

large seasonal variations in solar PV generation for a 7874 m2 floor-area office building, VTT

*Buildings* **2020**, *10*, x FOR PEER REVIEW 5 of 27

the power output capacity in SOFC mode [22]. This means that discharging the RSOCHS storage is

In the context of the case study, the RSOCHS system is used in an energy system to balance out large seasonal variations in solar PV generation for a 7874 m<sup>2</sup> floor-area office building, VTT FutureHub, (Figure 2), located in Espoo, Finland. The building characteristics are presented in Table 1. The location of the building is 60◦11011.4" N 24◦48049.0" E and the Köppen climate classification of the site climate is Dfb [23]. The assumed energy system of the office building (depicted in Figure 3) includes solar PV panels, an RSOCHS system, a water tank, a boiler with a hydrogen combustor for heat generation as well as connections to the electricity grid and the district heating network (for both import and export). The annual heating and electricity demand profiles of the office building are produced by IDA ICE simulations [24]. These profiles are presented in Figure 4.

**Figure 2.** The office building, VTT FutureHub, in Espoo, Finland. **Figure 2.** The office building, VTT FutureHub, in Espoo, Finland.


**Table 1.** VTT FutureHub building characteristics. **Table 1.** VTT FutureHub building characteristics.

The assumed energy system of the office building (depicted in Figure 3) includes solar PV panels, an RSOCHS system, a water tank, a boiler with a hydrogen combustor for heat generation as well as connections to the electricity grid and the district heating network (for both import and export). The annual heating and electricity demand profiles of the office building are produced by IDA ICE simulations [24]. These profiles are presented in Figure 4.

The optimization method presented in this paper is used to optimize the system operation, where the objective function is to minimize the annual operating expense (OPEX). Due to the differences in the two RSOC operation modes, they must be modelled as separate functions in the optimization problem. A binary decision variable, which determines the operation mode, must also be introduced at each time step in the optimization problem. The other variables in the optimization problem consist of energy transfer rates, which can be described as non-integer variables. The optimization problem is considered as a computationally expensive MILP problem due to its high number of integer and non-integer variables. The problem characteristics are, hence, ideal for the interval halving optimization method presented in this paper. Table 2 presents the dimensions of the case study optimization problem.

*Buildings* **2020**, *10*, x FOR PEER REVIEW 6 of 27

*Buildings* **2020**, *10*, x FOR PEER REVIEW 6 of 27

**Figure 3.** The reversible solid oxide cell and hydrogen storage system (RSOCHS). **Figure 3.** The reversible solid oxide cell and hydrogen storage system (RSOCHS). **Figure 3.** The reversible solid oxide cell and hydrogen storage system (RSOCHS).

**Figure 4.** Electricity and heating demand profiles of the office building. **Figure 4.** Electricity and heating demand profiles of the office building.

**Figure 4.** Electricity and heating demand profiles of the office building. The optimization method presented in this paper is used to optimize the system operation, where the objective function is to minimize the annual operating expense (OPEX). Due to the **Table 2.** Statistics of the optimization problem as well as the central processing unit (CPU) time required to solve the problem.


optimization problem is considered as a computationally expensive MILP problem due to its high

#### problem consist of energy transfer rates, which can be described as non-integer variables. The optimization problem is considered as a computationally expensive MILP problem due to its high number of integer and non-integer variables. The problem characteristics are, hence, ideal for the interval halving optimization method presented in this paper. Table 2 presents the dimensions of the 2.2.1. Calculated Cases

required to solve the problem.

number of integer and non-integer variables. The problem characteristics are, hence, ideal for the interval halving optimization method presented in this paper. Table 2 presents the dimensions of the case study optimization problem. case study optimization problem. **Table 2.** Statistics of the optimization problem as well as the central processing unit (CPU) time required to solve the problem. The optimal operation of three RSOCHS systems with different RSOC sizes is generated for different solar PV areas and hydrogen storage capacities. The RSOC sizes used for the RSOCHS systems are 20/80 kW (20 kW maximum power output in SOFC mode and 80 kW maximum power input in SOEC mode), 50/200 kW and 100/400 kW. All tested case parameters are presented in Table 3.

**Table 2.** Statistics of the optimization problem as well as the central processing unit (CPU) time

**Number of Binary Variables 1248**  Number of non-integer variables 9985 Number of constraints 6240 CPU time 15 s

**Number of Binary Variables 1248**  Number of non-integer variables 9985 Number of constraints 6240 CPU time 15 s


**Table 3.** Tested case parameters.

The RSOCHS systems are compared with a reference system without any energy storage in order to identify the strengths and weaknesses of the RSOCHS system. By comparing the energy systems, it is also possible to examine whether the RSOCHS system is a suitable solution for seasonal storage or not.

#### 2.2.2. Objective Function

The objective of the optimization problem is to minimize the annual OPEX of the RSOCHS system by using the operating mode and the different energy transfer rates as variables. The OPEX of the system in this context is defined as the sum of the electricity and district heating costs minus the income generated by energy export to the grid and the district heating network. The objective function is, thus, expressed as follows:

$$\min \sum\_{i=1}^{N} (\dot{Q}\_{im,i} \mathbb{C}\_{Qim,i} + P\_{im,i} \mathbb{C}\_{Pim,i} - \dot{Q}\_{ex,i} \mathbb{C}\_{Qex,i} - P\_{ex,i} \mathbb{C}\_{Pex,i}) t\_{\prime} \tag{1}$$

where *CQim*,*<sup>i</sup>* is the district heating import price at time step *i*, and *CPim*,*<sup>i</sup>* is the electricity import price at time step *i*, which is sum of the electricity spot price, the electricity distribution tariff and the grid tax. *CQex*,*<sup>i</sup>* and *CPex*,*<sup>i</sup>* are the export prices for heat and electricity at time step *i*. The variables . *Qim*,*<sup>i</sup>* and . *Qex*,*<sup>i</sup>* are the imported and exported heat rate at time step *i*, while the variables *Pim*,*<sup>i</sup>* and *Pex*,*<sup>i</sup>* are the imported and exported electrical power at step *i*, respectively. The time step size, *t*, and the number of time steps, *N*, are selected so that the sum of all time steps is equal to one year.

#### 2.2.3. RSOC Functions

In SOEC mode, the RSOC is modelled with one single function that describes the hydrogen output as a linear function of the electrical power input, *PE*,*in*,*<sup>i</sup>* , and the binary decision variable, δ*<sup>i</sup>* , for each time step *i*. The SOEC function is expressed as:

$$f\_{\mathbb{E}}(P\_{\mathbb{E},\text{in},\text{in}}\delta\_{i}) = a\_{\mathbb{E}}\delta\_{i} + k\_{\mathbb{E}}P\_{\mathbb{E},\text{in},\text{in}} \tag{2}$$

where the coefficients *a<sup>E</sup>* and *k<sup>E</sup>* are dependent on the size and operating range of the RSOC device.

In SOFC mode, two different linear functions are used; one for the hydrogen input, *fF*, and one function for the thermal energy output, *fG*. Both of these functions are dependent on the electrical power output of the RSOC. These functions are expressed as follows:

$$f\_{\mathcal{F}}(P\_{\mathcal{F},\text{out},i\prime}\delta\_{i}) = a\_{\mathcal{F}}\delta\_{i} + k\_{\mathcal{F}}P\_{\mathcal{F},\text{out},i\prime} \tag{3}$$

$$f\_{\mathbb{G}}(P\_{\mathbb{F}\rho\text{ut},i\prime}\delta\_{\mathbb{i}}) = a\_{\mathbb{G}}\delta\_{\mathbb{i}} + k\_{\mathbb{G}}P\_{\mathbb{F}\rho\text{ut},\mathbb{i}}.\tag{4}$$

The coefficients *aF*, *aG*, *k<sup>F</sup>* and *k<sup>G</sup>* are here also dependent on the size and operating range of the RSOC device.

#### 2.2.4. Constraints

There are five constraints for each time step in the optimization problem. These constraints consist of two energy balances, one for electrical energy and one for thermal energy, two constraints for the capacity of the RSOC and one constraint for the hydrogen storage.

The electrical energy balance is the sum of all electrical energy transfer rates and is expressed as:

$$P\_{\rm F,out,i} - P\_{\rm E,in,i} + P\_{\rm im,i} - P\_{\rm ex,i} - P\_{d,i} + P\_{PV,i} - \frac{h\_{\rm c,out} - h\_{\rm c,in}}{\eta\_{\rm c} \rm LHV} f\_{\rm E}(\delta\_{i\nu} P\_{\rm E,in}) = 0 \,\,\forall i \tag{5}$$

where the parameters *Pd*,*<sup>i</sup>* and *PPV*,*<sup>i</sup>* are the electrical power demand of the building and the power generated by the solar PV panels at time step *i*. The last part of the equation describes the power used by the hydrogen compressor at time step *i*, where η*<sup>c</sup>* is the total efficiency of the hydrogen compressor, LHV is the lower heating value of hydrogen gas. *hc*,*in* and *hc*,*out* are the specific enthalpies of the hydrogen gas at the inlet and outlet of the hydrogen compressor.

The thermal energy balance is expressed as follows:

$$f\_{\rm G}(\delta\_{i\prime}P\_{\rm F,out,i}) + \dot{Q}\_{b,i} + \dot{Q}\_{\rm im,i} - \dot{Q}\_{\rm ex,i} - \dot{Q}\_{d,i} = 0 \,\,\forall i \tag{6}$$

where . *Qb*,*<sup>i</sup>* is the heat generated through combustion of hydrogen, and . *Qd*,*<sup>i</sup>* is the heating demand of the building at time step *i*.

The input, *PE*,*in*,*<sup>i</sup>* , and output, *PF*,*out*,*<sup>i</sup>* , power ranges of the RSOC device are dependent on the size of the RSOC device as well as the operating mode. The RSOC power input and output constraints are thus controlled by the decision variables for each time step. These constraints are expressed as:

$$P\_{\rm E,in,max} \delta\_{\rm i} \ge P\_{\rm E,in,i} \ge P\_{\rm E,in,min} \delta\_{\rm i} \,\forall i \tag{7}$$

$$P\_{\mathcal{F},out,max}(1-\delta\_{\mathrm{i}}) \ge P\_{\mathcal{F},out,\mathrm{i}} \ge P\_{\mathcal{F},out,min}(1-\delta\_{\mathrm{i}}) \,\,\forall \,\mathrm{i} \tag{8}$$

The hydrogen constraint is dependent on the storage load at the beginning of each time step, which is the cumulative sum of the produced hydrogen plus the initial stored hydrogen, *EH*,0, minus the cumulative sum of consumed hydrogen. The hydrogen storage constraint is thus dependent on the system operation of earlier time steps, which is the reason why each time step cannot be optimized individually. Hydrogen storage constraints are expressed as:

$$\mathbb{E}\_{\rm H,up} \ge \mathbb{E}\_{\rm H,0} + t \sum\_{s=1}^{i} \left( f\_{\mathbb{E}}(\mathbb{P}\_{\rm E,in,s}, \delta\_{i}) - f\_{\mathbb{P}}(\mathbb{P}\_{\rm F,out,s}, \delta\_{i}) - \mathbb{Q}\_{\mathbb{B},\delta} \right) \ge 0 \,\,\forall \, i \in \{1 \ldots (N-1)\} \tag{9}$$

$$E\_{\rm H,amp} \ge E\_{\rm H,0} + t\sum\_{s=1}^{N} \left( f\_{\rm E}(P\_{\rm E,in,s}, \delta\_{i}) - f\_{\rm F}(P\_{\rm F,out,s}, \delta\_{i}) - Q\_{\rm b,s} \right) \ge E\_{\rm H,0} \tag{10}$$

where the parameter *EH*,*cap* is the maximum capacity of the hydrogen storage. The constraint for the last time step in Equation (10) is different from the constraint for the other time steps in Equation (9) since it only accepts a final storage load above the initial stored energy.

#### 2.2.5. Assumptions

The parameters in the case study are comprised of the annual energy demand, the solar PV generation and the energy price profiles as well as some technical data regarding the RSOC, the hydrogen compressor and the solar PV panels. The RSOC performance data are based on confidential data provided by the RSOC development team at the Technical Research Centre of Finland and are therefore not presented in this paper.

According to the US Department of Energy, the isentropic efficiency, η*is*, of hydrogen compressors used for small hydrogen gas terminals is about 65% [25]. The electrical motor efficiency, η*m*, of the compressor is assumed to be 95%. The total efficiency of the hydrogen compressor is thus assumed to be: *Buildings* **2020**, *10*, x FOR PEER REVIEW 9 of 27

$$
\eta\_{\rm c} = \eta\_{\rm is} \eta\_{\rm m} \approx 62\% \tag{11}
$$

The solar PV production capacity is assumed to be 0.17 kWp/m<sup>2</sup> and the PV generation profile is based on simulated data for a system where 50% of the solar panels are facing east and 50% are facing west [26]. The solar PV production capacity is assumed to be 0.17 kWp/m2 and the PV generation profile is based on simulated data for a system where 50% of the solar panels are facing east and 50% are facing west [26].

As a power and heat price scenario base, we use current prices and tariffs for Espoo, Finland. As mentioned earlier, the electricity import price is the sum of the electricity spot price, the distribution tariff and the grid tax. The electricity spot price used in this study is based on the Nord Pool spot price data from 2017 (Figure 5) [27] and the distribution tariff is selected according to the pricing of Caruna Oy, which is the electricity distributor in Espoo [28]. The grid tax in Finland is 2.253 c/kWh for non-industrial consumers [29]. The electricity export price is based on the price the Finnish energy company Fortum Oy pays for excess electricity generated by households [3]. This price is the hourly spot price minus a 0.24 c/kWh transfer fee. The average export and import electricity prices are summarized in Table 4. As a power and heat price scenario base, we use current prices and tariffs for Espoo, Finland. As mentioned earlier, the electricity import price is the sum of the electricity spot price, the distribution tariff and the grid tax. The electricity spot price used in this study is based on the Nord Pool spot price data from 2017 (Figure 5) [27] and the distribution tariff is selected according to the pricing of Caruna Oy, which is the electricity distributor in Espoo [28]. The grid tax in Finland is 2.253 c/kWh for non-industrial consumers [29]. The electricity export price is based on the price the Finnish energy company Fortum Oy pays for excess electricity generated by households [3]. This price is the hourly spot price minus a 0.24 c/kWh transfer fee. The average export and import electricity prices are summarized in Table 4.

**Figure 5.** Nord Pool electricity spot prices 2017 [27]. **Figure 5.** Nord Pool electricity spot prices 2017 [27].



**Total average price 3.06 9.12**  <sup>1</sup> Electricity distribution tariff for a 400 V connection; <sup>2</sup> Daytime distribution tariff, winter: Mon–Sat 7am–10pm, Nov–Mar.

1 Electricity distribution tariff for a 400 V connection

2 Daytime distribution tariff, winter: Mon–Sat 7am–10pm, Nov–Mar. The district heating import and export prices presented in Figure 6 are also based on pricing by Fortum Oy, who operates the district heating network in Espoo [4]. These prices vary from month to month, and are usually much higher during the winter months, when the heating demand is higher. The district heating import and export prices presented in Figure 6 are also based on pricing by Fortum Oy, who operates the district heating network in Espoo [4]. These prices vary from month to month, and are usually much higher during the winter months, when the heating demand is higher. The export prices are dependent on the temperature of the provided heat, but since temperatures are

export price for each month.

The export prices are dependent on the temperature of the provided heat, but since temperatures are not considered in the optimization model, the export prices in Figure 6 are based on the average

not considered in the optimization model, the export prices in Figure 6 are based on the average export *Buildings*  price for each month. **2020**, *10*, x FOR PEER REVIEW 10 of 27

**Figure 6.** District heating import and export prices in Espoo, Finland [4,5]. **Figure 6.** District heating import and export prices in Espoo, Finland [4,5].

#### **3. Results 3. Results**

installations.
