Advanced design and operation of Energy Hub for forest industry using reliability assessment

A large part of the refining heat production in the thermomechanical pulp mill can be recovered to supply the paper machine heat demand. This study introduces a novel approach for the heat integration of a thermo-mechanical pulp mill and paper machine using Energy Hub. An Energy Hub consisting of a steam generator heat pump and the electric boiler is integrated with the thermomechanical pulp mill to provide the heating demand of the paper machine. The advanced cost-efficient design and operation of the Energy Hub are investigated in this research by integrating thermo-economic analysis, reliability & availability assessment, and load profile prediction. The thermo-economic analysis combines economics and thermodynamics, which is necessary for energy system unit commitments. Reliability assessment will lead to more accurate modeling of real-life system operating conditions since system components ’ availability is considered in the design process. Load profile prediction estimates the Energy Hub load for the next hour, which helps with the optimal operation of the Energy Hub. Different state-of-the-art long-short-term memory (LSTM) neural network models have been developed to achieve the best time series model for refining heat prediction in the thermomechanical pulp mill. Results show that all the time series models are effective for refining heat prediction, while Bidirectional LSTM appears to perform better than others with the correlation coefficient and root mean square error of 0.9 and 0.15, respectively. In addition, the proposed Energy Hub design approach is compared with the conventional design method. The proposed design method offers a robust design that isn ’ t impacted by unsupplied demand penalty rates. Depending on the penalty rates, the total system cost could decrease by 14%-28% utilizing the proposed design method.


Introduction
Many industrial processes require heat in different parts of production.There are two major energy sources in the pulp & paper industry: electricity and heat.Energy savings can be achieved by heat recovery.For example, if the Thermo-mechanical pulp (TMP) mill and paper mill are close to each other, the heat produced in the refining process in the TMP mill can be recovered and used to supply the heat demand of the paper machine.As a significant amount of electricity is converted into heat by refiners in the TMP mill, it is necessary to recover the refininggenerated heat for other purposes [1].However, sustainable industrial heat production is vital for the utility system to avoid thermal load loss (unsupplied demand).For instance, if the paper mill's heat demand is not met, problems will arise during the drying procedure, resulting in financial and raw material losses.The heat produced by the refining process in the TMP mill depends on several factors mentioned in [2][3][4].A paper machine can benefit significantly from an Energy Hub (EH) concept and application for sustainable heat production and minimizing heat supply costs.Multi-carrier energy systems can be optimized with the help of an EH, which is a new and promising concept.Research work on the EH has begun over the last few years following its introduction in 2007 [5].EHs facilitate the move towards sustainable multi-energy systems and achieve the realization of energy system models [6].Several energy carriers are converted, stored, and supplied through the EH system for multiple energy types, meeting challenges such as energy consumption and resources, multiple energy infrastructure, and flexibility.
Industries are the biggest energy-consuming sectors.Waste heat recovery in industrial processes improves plant energy efficiency and reduces greenhouse emissions.High and medium-temperature industrial waste heat has a higher value than low-temperature industrial waste heat.On the other hand, several in-use technologies recover high and medium-temperature waste heat [7].In today's world, heat recovery from low-value thermal wastes is gaining attention.Heat pumps are promising technologies for recovering low-temperature industrial waste heat.Based on the analysis, today's heat pump technology is able to supply 16% of the industrial heat demand (up to 80℃) in Germany [8].However, these systems are complex and have high capital and operational costs.On the other hand, since they are made of many components, the system's availability is less than other electrical heatproducing technologies, such as electric boilers.Therefore, boiler technology is used in the heat distribution system as it offers high reliability and availability, ensuring stability in production and reducing initial and operating costs.
The optimum design and operation of EH are highly regarded nowadays.Distributed energy production units such as EH, which produce different types of energy simultaneously, enhance the efficiency and stability of the energy systems [9][10][11].The cost-effective and secure operation of EH depends on accurate load profile predictions.Operators can complete unit commitments and assess system stability with the help of load forecasting.Time series models are the best way to estimate the future value of a time-dependent process.In time series forecasting, future values are predicted based on the history of the previous values.Long-short-term memory (LSTM) networks are well suited to processing, prediction, and classification in analyzing time series data.LSTMs have recently received increased attention due to their ability to deal with highly complex, non-linear datasets.An LSTM algorithm is tested on a publicly available dataset of residential meters in ref [12].The results indicated that LSTM outperformed rival machine learning (ML) algorithms at forecasting residential loads in the short term.In ref [13], a dynamic model is trained using the LSTM algorithm to generate predictions for impulsive loads.LSTM is discussed in multiple configurations in ref [14], and the best architecture and optimum hyperparameters are proposed.In the other research, Taheri et al. [15] studied electricity demand time series forecasting based on a hybrid prediction model by integrating LSTM with empirical mode decomposition (EMD).Intrinsic mode functions (IMFs) are subseries of load time series data that are broken down by the EMD algorithm.Different LSTM models are developed for each IMF derived.For the final step of energy demand prediction, the outputs from each LSTM learner are fed into a meta-learner.In experiments on the California ISO dataset, their proposed hybrid model outperforms single LSTM, LR, and XGBoost algorithms.However, forest industry research is lacking when it comes to predicting electricity consumption and heat production.This study utilizes different LSTM structures to predict the TMP mill heat production.The analyzed LSTM structures are single LSTM with hidden dense layer, stacked LSTM with hidden dense layer, and Bidirectional LSTM.Results show that all the LSTM structures perform effectively for the TMP heat prediction.
It was common practice to plan and manage different energy systems independently in the past.Nevertheless, as technology advances, such as the development of efficient multi-generation systems, and integrated energy infrastructures such as electricity, natural gas, and district heating (DH) networks become increasingly advantageous, the movement towards multi-energy systems (MES) is accelerating.Such systems interact in a synergistic way between different energy carriers and systems.Nevertheless, the application of such a concept requires a management tool that integrates the components of the system.EHs, which can be described as places where different energy carriers are produced, converted, stored, and consumed, are promising for MES integration [16].Based on 200 articles published from 2007 to 2017, Maroufmashat et al. [17] conducted a systematic review of energy hubs.An overview of the articles is provided on the basis of modeling methodology, planning and operation, economic and environmental aspects, and applications of energy hubs.Techno-economic analysis of the heat pump system has been comprehensively studied in references [18,19].Reliability modeling has only been presented in a few reports for cogeneration systems.Haghifam and Manbachi [20] made a reliability and availability modeling of combined heat and power (CHP) systems.They classified combined heat and power subsystems and modeled the meantime-to-failure, availability, and reliability from the production unit to end-use.Their suggested CHP availability and reliability model is based on the state space and the continuous Markov method with heatgeneration, fuel distribution, and electricity-generation subsystems.As a result of this study, it appears that improving the gas distribution network, the water distribution network, and the water delivery network to the CHP, as well as optimizing the failure and repair rates of CHP systems, can enhance the reliability of a complete integrated system and play a major role in determining whether CHP systems are economically and technically feasible.A cooperative approach to energy management was presented by Tiwari and Singh [21] for multi-carrier energy hubs.The proposed scheme facilitates energy exchange between energy hubs, resulting in greater economic and environmental benefits.Results show that networked systems are more reliable and have lower operating costs and greenhouse gas emissions.Parisio et al. [22] proposed a robust solution to the EH operation scheduling problem.First, they made an integer dynamic model to establish a general EH modeling framework.Next, a short-term operating schedule is proposed in a deterministic environment and uncertain scenario.robust schedule of input power flows, which is less sensitive to uncertain converter efficiencies than the nominal schedule.In the other research, Talebjedi and Behbahaninia [23] conducted research to integrate thermo-economic analysis with reliability assessment for the optimal design of the EH consisting of a combined cooling, heating, and power (CCHP) system.Their proposed optimal design results in a 16% reduction in the system total cost over the EH life span.However, their research could have been improved for optimal operation by implementing EH load prediction.This research proposes an advanced cost-efficient design and operation of an EH in interaction with a TMP mill and paper machine.The proposed approach is based on the integration of thermo-economic analysis, reliability & availability analysis, and EH load prediction.The studied system is an EH consisting of a steam generator heat pump and electric boiler integrated with the TMP mill to provide the paper mill heating demand.A future trend in industrial processes, including heat production, is electrification.Hence, heat pumps and electric boilers would be viable options for energy converters.Paper mill heat demand is given to the Central Hub controller on an hourly basis.The paper mill's heat demand prediction is not necessary since the heat demand of the paper machine remains relatively constant and stable, averaging at 70 MWh with a standard deviation of 10 MWh.However, the TMP mill heat production depends on various factors and is predicted using the LSTM model.The Central Hub Controller receives information regarding the paper mill heat demand, TMP mill heat production, electricity market price, and the availability of each subsystem to perform the unit commitment.The major paper novelty and contribution are summarised as follows: • A step-by-step framework is developed based on the integration of thermo-economic analysis, reliability assessment, and load profile prediction for advanced cost-efficient design and operation of the EH (section 8).• Different LSTM structures are developed and tested for the best TMP mill heat prediction model.• An economic analysis is carried out to compare the conventional and suggested cost-efficient design approaches of an EH in the forest industry energy system.

Case study
The schematic of the Energy Hub, integrated with the paper machine, TMP mill, and electricity grid, is depicted in Fig. 1.The LSTM model predicts the 1-hour ahead TMP heat production and sends the information to the central hub controller (CHC).A central hub controller also receives the hourly heat demand from the paper machine.Heat demand rarely fluctuates at the paper mill, so developing a heat demand prediction model hasn't been considered in this research.At each time step, the CHC receives information about the state of the heat supply system illustrated in the state-space diagram (Fig. 3).The heat supply system consists of the EH and TMP mill.The EH's energy conversion components are the steam generator electric boiler and heat pump, where the heat pump includes a series of five subsystems: compressor, evaporator, condenser, expansion valve, and heat recovery steam generator.As a series system, the entire heat pump fails in case of failure of any subsystems.The electricity grid supplies the electric boiler and heat pump electricity demand at the market price.The EH unit commitment depends on various factors, such as the electricity market price, energy conversion components efficiency, the state of the system (Fig. 3), TMP heat production, and paper machine heat demand.

General modeling concept
Reliability represents the likelihood of an object being healthy under  a particular period.The reliability analysis techniques are typically categorized into quantitative and qualitative approaches.Failure modes and effects analysis are the most common qualitative approach.The major quantitative approaches are reliability block diagrams, fault tree analysis, and Markov analysis.The random behavior of repairable systems, such as electric boilers and heat pumps, can be analyzed by the Markov method [23].
The Markov chain is a stochastic random memoryless method.The probability distribution of the following state is only affected by the present state, and the history of the states is not influential.The Markov chain contains a set of state transitions.For the system reliability and availability analyses, it is necessary to characterize the state transitions for a single step or time interval in a matrix format.A stochastic matrix describes the probability distribution of these state transitions, where each row explains the different system operating conditions.In this study, the state-space method (SSM) is combined with the Markov chain to analyze each subsystem's random failure mode and evaluate the whole heat supply system's availability.The SSM is the best technique for quantitative reliability and availability assessment of repairable and complex systems [24].The SSM consists of three essential steps: 1-Defining the system boundaries and all the system's possible operating and failed states, 2-Building the state transition diagram, 3-Figuring the ultimate state probabilities during the system life span via different methods such as a frequency balance approach.
Fig. 2 shows the state-space diagram of a repairable device with a constant failure and repair rate consisting of two states.The system is referred to as 'up' when the system is healthy and operating as normal.In contrast, the state of a failed system is mentioned as a 'down' state.λ = failure rate.

Failure and the repair rate of series systems
The repair and failure rate of the system consisting of a series of components is discussed in this section.Each system component's failure and repair rate is essential to formulate the frequency balance equations.Ref [25] introduces the calculation method of the failure and repair rate of the series system with a constant transition rate.The calculation method of failure and repair rate of a system made up of Nseries components is provided in the equations below.
Where N is the number of series components, and j determines each subsystem component.The heat pump consists of a series of five components more likely to fail: compressor, evaporator, condenser, expansion valve, and heat recovery steam generator.The failure and repair rate of the whole heat pump can be calculated by substituting each subsystem's failure and repair rate in Eq. ( 1) and Eq. ( 2), respectively.

Frequency balance method
The frequency balance is a calculation method to compute the statespace diagram's limiting state or steady-state probabilities.The frequency balance concept is established on an ergodic system where each  state's frequency entering and leaving is the same.Frequency balance variables are limiting state probabilities where their sum is 1.The frequency balance implementation for a model consisting of two states (=up and down) is described in the following.The state-space diagram of the system is given in Fig. 2. The frequency balance equations for state 0 and state 1 are: The above equations are not independent, so the following equation should be added to solve the system of equations.The sum of steadystate probabilities is equal to 1: Solving the above system of equations gives the steady-state probabilities of:

Limiting state probabilities: Case study
The state transition graph of the case study is shown in Fig. 3. Applying the transition rates to the transition graph is the first step in forming the state-space diagram for the continuous and discrete Markov process to simplify the use of the frequency balance approach.Observing possible states is necessary for developing a state-space diagram.In general, no restrictions are imposed on the type or number of states and transitions.In reliability and availability analysis, preparing a proper state-space diagram is essential.In Fig. 3, we present the state-space diagram for the studied system, taking into account four central subsystems that are more likely to fail: steam generator heat pump (HP), electric boiler (EB), and refining lines in TMP mill (RL1 and RL2).The computational load can be reduced by omitting states with low limiting state probabilities from the state-space diagram.As an example, in the studied system, it is unlikely that three or more subsystems fail simultaneously; therefore, states with more than two failed components are ignored.There is a significant reduction in computational cost due to this modification.In the diagram, 11 states are included, of which states with more than two collapses at the same time have been removed due to a small limiting state probability.A state-space diagram consists of the most likely to fail components, such as an electric boiler, a heat pump, and refining lines in a TMP mill.As stated earlier, a heat pump is made of five series components where; the system's failure and repair rate calculation is discussed in section 3.2.A frequency balance analysis with constant failure and repair rates over time has been implemented to get the system's limiting state probabilities.The failure and repair rates in this study are based on the data we gathered from the pulp & paper mill and our previous research in reference [23].
The frequency balance equations for each state are given in Eq. ( 8)- (18).A technique like this is beneficial for the solution of relatively small systems that do not require a matrix formulation and subsequent solution by a digital computer.State probabilities (p i ) are the unknown parameters in these equations, which must add up to one.
As mentioned earlier, the above equations are not independent, and there is one dependent equation in the equation set.Since there are ten independent equations and 11 variables, The equation set needs to be expanded with another equation of the total probability theorem:

Steam generator heat pump (SGHP)
The industrial energy consumption share is 30% of the total energy production.Over half of the energy consumed by industries is used to generate steam at various temperatures.Therefore, various heat production technologies are developed in order to minimize the energy required to produce industrial steam.About half of the total industrial waste heat has a temperature range of 60-100 ℃.Nevertheless, discharged hot water should be cooled to below 40 ℃ to minimize the environmental impact.Therefore, industrial waste heat recovery development is essential to reduce primary industrial energy consumption [26].
A steam generator heat pump (SGHP) is widely used to generate high-temperature steam from low-temperature waste heat.SGHP is a viable technology to recover waste heat from diverse industries and generate steam for heating, drying, humidification, and other industrial applications that require heat.A vapor compression heat pump is a technology that uses electricity in a compressor to transfer heat from a lower to a higher temperature.Fig. 4 shows the schematic of an industrial steam generator heat pump system.In this research, the waste heat from the paper machine (typical temperature: 40-60 ℃) could be used as low-temperature waste heat for the heat pump evaporator to increase the heat pump efficiency by reducing the lifting temperature.However, since the information regarding the amount of the paper machine waste heat is not accessible, we do not consider this fact in the analysis, and we assume the evaporator temperature is as same as the outside temperature.
CO2, or R744, is an environmentally friendly refrigerant with perfect thermodynamic properties appropriate for numerous applications.Low critical temperature and high operating pressure make R744 different from other refrigerants.R744 is non-flammable with low toxicity, operating in high pressure.A special system design is required for the high operating pressure of R744, which makes the system more complex.The refrigerant cost is low, while the system is more costly than other conventional systems.The most prominent possible hazard of utilizing CO2 as the refrigerant is dry ice formation in case of reducing the pressure and temperature below the triple point (4.2 bar, −56 • C).

Long Short-Term memory (LSTM)
Classical neural networks cannot catch the dataset variables' temporal dependencies.Recurrent Neural Networks (RNNs) are developed to address this problem.As presented in Fig. 5, a recursive structure of the RNNs gives the ability to predict a series of potential values using observed historical inputs while keeping knowledge over a multitude of time horizons.
Recurrent Neural Networks (RNNs) are capable of using contextual information in the process of mapping the input and output sequence.However, the context range for the architecture of standard RNNs is limited.The effect of the input variables on the hidden and output layer deteriorates rapidly as it cycles around the network's recurrent links.This phenomenon is repeatedly mentioned in the literature as the vanishing gradient problem.Therefore, typical RNNs are unqualified for learning long-term temporal dependencies because of the vanishing gradients problem phenomenon.Ref [27] further explains and discusses this phenomenon.In 1997, Hochreiter and Schmidhuber [28] introduced the Long-Short-Term Memory (LSTM) architecture to address the problem of vanishing gradients for RNNs, shifting from normal to deep recurrent neural networks.The LSTM training approach is based on keeping the information regarding the previous states and can be qualified for state or memory recognition tasks.Utilizing the mentioned method, LSTM is capable of solving gradient bursting and vanishing problems.The LSTM cell architecture is presented in Fig. 6, where the memory cell mode block controls the signal through the forget gate, input gate, and output gate.Each LSTM gate has a specific function, and the following Equations show the procedure of computing each vector at time t.accuracy.In this regard, three main LSTM structures are studied and employed: 1-Single LSTM with hidden dense layer, 2-Stacked LSTM with hidden dense layer, and 3-Bidirectional LSTM.

Single LSTM with a hidden dense layer
The principal LSTM network is formed of a single hidden LSTM layer followed by a standard feedforward output layer.The LSTM structure and principal concept have been discussed in chapter 4. The structure of the employed single LSTM followed by a hidden dense layer is depicted in Fig. 8.

Stacked LSTM with a hidden dense layer
A stacked LSTM architecture refers to a model comprised of multiple layers of LSTM.The stacked LSTM is an extension of the principal LSTM model, which consists of a single hidden LSTM layer followed by a standard feedforward output layer.This model has multiple hidden LSTM layers, where each layer includes several memory cells.The stacked LSTM hidden layers make the model deeper by increasing the hidden layers and LSTM units.In contrast to single value output, an above LSTM layer provides a series of values instead of a single value to the below LSTM layer.They produce one output per input time step instead of having a single output for all input time steps.The objective of the different proposed LSTM network structures is to predict the next total heat production from the pulp mill for different time horizons.Fig. 9 illustrates the general structure of stacked LSTM with a hidden dense layer.

Bidirectional LSTM
The bidirectional LSTM is an augmentation of the standard LSTMs to enhance a model's prediction performance by utilizing the sequence of information in both forward and backward directions.In this regard, two LSTM networks are trained based on the original and a reversed replica of the data.This approach adds more meaning to the network to improve its performance.The bidirectional LSTM principal concept is simple.The observed data in its current format is assigned to the first recurrent layer input data.In contrast, the repeated replica of the observed data is considered an input to the copy of the first recurrent layer.Fig. 10 shows the principal concept of Bidirectional LSTM [29].
Based on Fig. 10, the backward hidden layer (h ), and the output (y) can be calculated as [30]: where w x h → and w xh ← are the feedforward and backward hidden weight matrices, respectively.b h and b h ← are the bias vectors, and H indicates the hidden layer.

Hyperparameter tuning
Model hyperparameter tuning is applied to the experiments to achieve a time-efficient and accurate prediction model.A good performance can be achieved with the LSTM network by choosing the optimum hyperparameters.Depending on hyperparameters, deep learning algorithms differ in performance, time cost, and memory usage.A hyperparameter selection can influence the performance of machine learning algorithms from state-of-the-art to mediocre [31].While there are some general guidelines regarding the appropriate values of these hyperparameters, hyperparameter tuning is still necessary because the optimal values depend on the data sets, the performance criteria, and several other variables [32].The hyperparameter tuning in this research has been made through the RandomizedSearchCV module from the Scikit-learn (Sklearn) library in Python.The parameter values are not all tested as in GridSearchCV, but instead, a fixed number of parameter settings are sampled from the specified distributions.
In addition to the learning rate, the batch size is an important hyperparameter.As the main hyperparameter of the training algorithm, the learning rate determines how quickly the weights are updated during each training cycle [33].Using a smaller learning rate reduces loss but results in slower convergence [34].Training speed can be increased with a larger batch size, but the convergence speed may be lower, and with an increased batch size, generalization ability may decrease.In the literature [35], an extensive discussion is made of the processes of LSTM hyperparameter tuning.The hyperparameter dictionary and optimal values for the RandomizedSearchCV are given in Fig. 11.

Economic parameters
Following is an overview of the economic parameters used in this research to design and operate a cost-effective EH system to supply forest industry heat demand.Over the system's lifetime, all annuities will be converted to 2020.
Net present value: Net present value (NPV) is the value of all future cash flows over the entire investment life span discounted to the present.The concept of NPV is employed in this research to convert all the system's annual costs to the reference year of the analysis where the investment has been made.The annuity factor (AF) can be utilized to calculate the net present value of all annuities, such as maintenance, operation, and load loss costs.A system's annuity factor can be expressed by Eq. ( 29), assuming k as the interest rate and m as the system life span.This value is 18.04 in our case study considering 0.01 as the interest rate, and 20 years for the system lifespan.
Energy carrier price: The electricity price is comprised of three components, namely the hourly Nord Pool spot price, the cost of distribution, and the tax on electricity.Table 1 displays the charges that the Caruna company in Southern Finland levies on industrial users for electricity distribution fees and taxation [36].
Energy Cost: Energy cost is the cost of electricity consumed by the EH to produce industrial heat through energy conversion units.Electricity is priced according to the market price in this study.The energy cost calculation for subsystem n, at time t, and state i in the state-space diagram is: where FC, EP, and MEP are the energy cost, electric power consumption, and market electricity price, respectively.Capital Cost: An initial capital cost is incurred when a hub designer decides to install a specific component.Every EH component has a capital cost determined by its design parameters.A design variable that influences capital costs the most is capacity.Linear regression can roughly represent the capital cost of each EH subsystem as a function of its capacity [37][38][39].
CC HP (€) = 400C HP (kW) Where C EB and C HP are the electric boiler and heat pump capacity in kW, respectively.

Load Loss Cost:
The EH is responsible for providing the heating load determined by the Central Hub Controller.The EH operator is charged for each MWh shortage in the heating supply.For each MWh of unsupplied heat, there is a charge (unsupplied demand penalty cost) that is significantly higher than the heating price.In this research, the unsupplied demand penalty cost (pc heat ) is considered to be ten times more than the heating price.However, a sensitivity analysis has been made in the result section to analyze the effect of different unsupplied demand penalty rates on the final total cost reduction.The heat penalty cost coefficient defines the extent of unsupplied demand penalty cost based on the selling price of the heat.The following equation shows the mathematical expression to calculate the load loss cost (LLC): LLC i,t = HLL i,t × pc heat t (33) where hll i,t (MWh) is the heat load loss at time t and state i in the statespace diagram (Fig. 3).pc heat ( € MWh ) is the unsupplied demand penalty cost for each MWh of unsupplied heat.
Operation costs: Energy cost and load loss cost (llc) make up the operation cost.System operation cost over the life span (OC Total ) is given by Eq. (34).p i defines the system operating condition in the state-space diagram and is the limiting state probability of the i th stage.The system operation cost is the summation of the system operation cost at each operating condition state multiplied by the limiting probability of the state.All the annuities have been converted to the refrence year by annuity factor (AF). T is the number of hours in a year (=8765).
System total costs: The system total cost in this article encompasses the investment, operation, and maintenance costs over the system life span.The system state and availability of the subsystems are considered in calculating the system's total cost.As mentioned earlier, system operation costs include unsupplied demand penalty costs and fuel (electricity) costs.System total cost is one of the measures to evaluate the cost-efficiency of the system.
In Eq. ( 35), we find the mathematical term that can be used to calculate the system's total cost (TC).The annual maintenance cost (MC) is considered 10% of the system capital cost.Annual maintenance cost is transferred to the refrence year by annuity factor (AF).The last two terms indicate the system's maintenance and capital cost where n indicates the energy conversion component in the EH.

LSTM performance assessment
Any type of forecasting requires models' performance assessment criteria to evaluate the prediction accuracy.Root mean square error (RMSE), and correlation coefficients (R) are the most well-known accuracy metrics to evaluate the prediction model performance.In this study, the prediction model dataset is separated by the first P portion of a percent to the training set and the following (1-P) percent as the testing set.Using a time series forecast model requires dividing training and test data systematically and preserving the data's sequence.Therefore, the division of training and test data cannot be done randomly, and the sequence of the data must be respected.The following equations show the main performance evaluation criteria for the time series prediction model.
Where ym j , yp j , ym, yp, and n are measured value, predicted value, mean of measured data, mean of predicted data, and the number of data, respectively.

Cost-efficient EH design
This section introduces an optimization model for the EH's costefficient design.A cost-efficient optimization model minimizes the system's total cost in the long run by determining the optimal capacity of energy conversion technologies inside the hub.The optimization model decision variables are the capacities of the electric boiler and heat pump.The mathematical model is established by defining linear constraints and the objective function.Therefore, the final optimization model is linear and is solved using the CPLEX algorithm.The model objective function and constraints for the proposed and conventional design approach are discussed in the following.

Optimization objective function:
In both conventional and proposed designing scenarios, the optimization model minimizes the EH total cost provided in Eq. ( 35) to achieve the optimum capacities of the EH components.However, The index i = 0 should be used in all equations pertaining to the conventional design approach (case 1) since state 0 represents the situation where all system components are in healthy condition.The heat load loss cost in the model objective function is zero in case 1 as a healthy system is supposed to be able to supply the thermal load.Re-writing the EH total cost for the state 0, Eq. ( 38) gives the optimization model objective function for the conventional design approach (case 1): MC, CC, andFC are the maintenance, capital, and energy costs, respec-tively.p denotes the state of the system.In the state-space diagram, P = 0 represents a zero state in Fig. 3 where all EH components are in healthy condition.n indicates the energy conversion component in the EH.
The second design approach (case 2) uses real operating conditions as the basis for system modeling.A reliability and availability assessment is conducted to observe what effect a subsystem's failure would have on the system's performance.In this modeling approach, the optimization objective function is the total cost by considering all the system states in the state-space diagram.By reorganizing Eq. ( 35) over all the system's possible states, the optimization objective function for the second optimization scenario (case 2) will be: Optimization constraints: In addition to the constraints relating to the EH components' capacities, optimization constraints also include hourly constraints on meeting the EH heat load.A coupling matrix describes how the hub transforms power from the input to the output; it can be derived from a hub's converter's efficiency and structure.Eq. ( 40) displays the EH input power vector, coupling matrix, and output power vector for the system depicted in Fig. 1.
Where EP is EH's total electric power consumption, and ϑ is the dispatch factor.HLL t is the heat load loss at time t in case of subsystem failure.The heat load loss is zero in the state = 0 since all the system components are in healthy condition.η EB , η HP , and η HEX are the efficiencies of the electric boiler, heat pump, and heat exchangers, respectively.TMP Heatproduction,t and PM Heatdemand,t are the produced heat by the thermomechanical pulp mill and paper machine heat at time t, respectively.The following equations provide additional constraints regarding the maximum capacity of system components.The maximum capacities for the electric boiler and heat pump are considered to be 120 MW.

Cost-efficient EH operation
EH load prediction is key to ensuring the cost-efficient operation of the system.The operation stage takes place after the system design, where the optimal capacity of all system components is determined in the design phase.A central energy hub controller can calculate the EH load based on the heat supply from the TMP mill and the heat demand from the paper machine.Due to the low fluctuation in the heat demand of the paper mill, the development of the heat demand prediction model has not been considered in this research to reduce the complexity.However, the TMP mill's heat production is affected by several time- dependent factors, including seasonal variations in wood chips quality, refining plate erosion, and wood chips moisture content.Different LSTM architectures studied in this research for TMP heat prediction are single LSTM with dense layer, stacked LSTM with dense layer, and Bidirectional LSTM.Fig. 12 depicts the step-by-step framework to achieve an accurate TMP heat prediction model.The best prediction model will be incorporated with the central hub controller to perform the unit commitment.The central hub controller performs the unit commitment for the next time step by maximizing the hourly profitability of the EH operator.The hourly heat pump and electric boiler power consumption are the optimization decision variables at each time step.The optimization objective function of the linear programming problem is provided in Eq. (43).Each time step must adhere to the same optimization constraint as the coupling matrix provided in Eq. (40).Using the CPLEX algorithm, the linear optimization problem can be solved at each time step.

TMP refining heat prediction: LSTM method
The LSTM network performance evaluation to predict TMP heat production is discussed in this section.LSTM network is a state-of-theart RNN that permits information to remain by sequential network construction.The LSTM learns how to map the sequence of previous observations as input to an output observation.Different LSTM models are analyzed in this research to achieve the best prediction accuracy.Table 2 shows the accuracy performance of different LSTM structures studied in this research.The employed and tested LSTM structures in this research to predict 1 h ahead TMP heat production are Single LSTM with a dense layer, Stacked LSTM with a dense layer, and Bidirectional LSTM.
All prediction models perform very closely, according to Table 2.However, the Bidirectional model seems to perform better than other models.The correlation coefficient and root mean square error for the best model, a Bidirectional LSTM structure, are 0.95 and 0.12 for the training dataset and 0.90 and 0.15 for the test samples, respectively.A higher prediction accuracy indicates an improved hour-ahead EH load estimate.This means that if a prediction model is able to accurately  forecast the EH load for the upcoming hour, then it is considered to be a better estimator of the EH load.The accuracy of the prediction is an important measure of the model's performance, and a higher accuracy implies that the model is more reliable and robust to heat load fluctuations.Therefore, a model with higher prediction accuracy is more likely to provide an accurate estimate of the EH load, which can help in making informed decisions about energy consumption, distribution, and management.LSTM model prediction on the test dataset is represented in red data.In machine learning, an algorithm is tested using a second (or tertiary) data set after being trained using an initial training data set.The smaller difference between the measured data and the predicted data indicates higher prediction accuracy.Results show that the studied time series models have high accuracy in estimating the heat produced, confirming the feasibility of using these models to optimize EH performance.The data used in this research are measured data from the TMP mill in a Nordic country.Due to the contract made with the industrial partner to maintain data security, the y-axis label has been removed from the graphs (Figs.[13][14][15] to keep the data secure.

Energy Hub optimal design
This section investigates the cost-efficient design of the EH to supply forest industry heat demand.With a regulatory mechanism, the EH  balances energy supply and demand by utilizing the energy conversion technologies inside the hub.The EH's cost-efficiency design requires thermo-economic analysis along with reliability & availability assessment.A reliability assessment factors in the availability of each system component in the system modeling.As each component of the system has a chance of failing in real operating conditions of a system; therefore, reliability assessment is necessary for making a more realistic model.The system optimum (cost-efficient) design approach is introduced in section 8.1.The results of the frequency balance technique to calculate the limiting state probabilities of the heat supply system in the state-space diagram (Fig. 3) are presented in Table 3.The heat supply system in the studied case consists of the EH and TMP mill.Limiting state probabilities are important values to define the system's limiting states and availability.In a Markov chain, the limiting state probabilities refer to the long-term probabilities of the system being in each state, as the chain continues indefinitely.These probabilities are reached when the system has reached a state of equilibrium and is no longer changing significantly.The limiting state probabilities provide valuable information about the behavior of the Markov chain over time.In the case study we conducted, the likelihood of a stable system being healthy in the long run is 94%.Table 3 gives the steady-state probabilities for the whole heat supply system in different system operating conditions.This paper investigates the advanced cost-efficient design of the EH based on the conventional (case 1) and proposed (case 2) design approach.The proposed design approach, which is based on the real system operating condition (case 2), is compared with the ideal hypothetical scenario where all the system components are healthy and 100% reliable (case 1).Table 4 presents the results of the optimal capacities for the conventional and proposed EH design approach for different unsupplied demand penalty costs (pc heat ).
For all different unsupplied demand penalty costs, the research demonstrates that the electric boiler capacity is higher in case 2. Nevertheless, there is a slight decrease in the heat pump capacity in case 2 when compared to case 1.As an example, when the penalty cost is set at 10( € MWh ), the electric boiler capacity in case 2 rises to 67 MW, which is a 346.66% rise compared to the conventional design where the capacity is only 15 MW.The optimization model has enabled this increase in capacity to reduce load loss costs.On the other hand, transitioning from case 1 to case 2 only results in a small change of 1.51% reduction in the heat pump capacity.The electric boiler has a higher availability leading to less load loss cost which explains the impressive capacity increase of the electric boiler moving from case 1 to case 2.
Table 5 illustrates the system's costs over the system life span for case 1 at various unsupplied demand penalty cost values.Case 1 represents the operating condition where there is no possibility of failure for the system components.Due to this assumption which is not supported by reality, the cost of load loss is significant, amounting to 102.63, 149.75, and 194.22 million € for penalties associated with unsupplied demand at rates of 10, 15, and 20 ( € MWh ), respectively.Consequently, if failure modes are not taken into account in the optimization model, the EH operator will face a significant penalty cost due to the load loss cost.The system's total cost is equal to 472.63 million € where energy cost and load loss cost contribute to the majority of the system's total cost.Due to the high load loss cost, the possibility of improving the system design to minimize the total cost by considering the real system operating condition is quite evident.Small changes in the system capital cost could result in higher system availability and less load loss cost.This issue has been investigated in the second approach of system design (case 2).
Table 6 presents a comprehensive analysis of the system costs for case 2 and a comparison with the system costs of case 1.The system's load loss cost experiences a substantial reduction of 98% across all unsupplied demand penalty rates compared to case 1.As an example, when the penalty cost is 10 (€/MWh), the load loss cost drops significantly from 102.63 million € in case 1 to 1.94 million € in case 2. In addition, reductions are also observed in the system total cost, ranging from 14% to 28%, depending on the unsupplied demand penalty rates.According to the results, the suggested method allows for a higher possibility of decreasing the system's total cost when the penalty rates are increased due.The reason for this phenomenon is that load loss increases rapidly as penalty rates are increased, and the proposed design approach aims to minimize the cost of load loss through the optimal design of the system.The system's total capital cost increases by a range   of 17% in case 2 for different unsupplied demand penalty rates.Most of the rise in the capital cost of the system is dedicated to raising the capacity of the electric boiler to provide heat in the event of a heat pump malfunction.This plays an important role in minimizing the system load loss cost which ultimately leads to a decrease in the overall system cost.
In addition, findings indicate that in the second modeling method (case 2), the overall cost of the system throughout its lifespan does not change significantly with variations in unsupplied demand penalty costs.As a result, the proposed approach offers a robust design strategy that is not significantly impacted by unsupplied demand penalty rates, which is a clear advantage of this approach.By subtracting the system total cost in case 1 from the system total cost in case 2, the proposed EH design approach would yield a total benefit of 68.38, 115.11, and 160 million € for penalty cost at the rate of rates of 10, 15, and 20 ( € MWh ), respectively.
For an unsupplied demand penalty rate of 10 (€/MWh), Fig. 16.displays the results of a sensitivity analysis conducted on the total cost of the system using various capacities of EH components.The graph indicates that smaller capacities, when deviating from the optimum point, lead to a greater increase in the overall cost of the system compared to larger capacities.The primary explanation for this occurrence is the rise in load loss cost due to inadequate heat supply.Therefore, if the optimization model for an energy system does not take into account the different ways in which components can fail, there could be significant penalties imposed on the energy system operator due to load loss.This is because the system may not be able to meet the energy demands if there is a failure is the system component, resulting in disruptions and outages, and ultimately causing a significant cost to the operator.As a result, it is important to consider failure modes when optimizing the energy system to ensure it operates reliably and efficiently, even in the event of component failures.
Sensitivity analysis was conducted to determine the influence of the spot electricity price in the reference year on the results of the proposed Fig. 16.Sensitivity analysis on the system total to examine how changes in the capacities of EH components affect the total cost of the system.

Table A1
Optimal Energy Hub subsystem capacities for the different design approaches (reference Year: 2020 & 2021).method, please see Appendix.Table A.1 shows the optimal capacities for the EH subsystems for the reference years 2020 and 2021.There would be minor changes in the optimal subsystems capacities while changing the reference year from 2020 to 2021.Table A.2 shows the system costs for the reference years 2020 and 2021.The average Nord Pool spot electricity prices for the years 2020 and 2021 were 28 €/MWh and 73 €/MWh, respectively.The system costs would be influenced by altering the reference year, as is apparent from the results.However, the efficiency of the proposed design method remains unaffected by altering the reference year electricity price, as Table A.2 demonstrates.The results of the proposed method indicate a reduction of nearly 98% in system load loss costs and 11.5% in total costs.

Conclusion
A study of the advanced cost-efficient design and operation of the Energy Hub as an industrial heat provider is presented in this paper.For this purpose, the concept of integrating thermo-economic analysis, reliability & availability assessment, and load profile prediction is employed.The thermo-economic analysis combines economics and thermodynamics, which is necessary for energy system resource allocation.Reliability assessment will lead to more accurate modeling of real-life system operating conditions since system components' availability is considered in the design process.Load profile prediction estimates the energy demand for the next hour, which helps with the optimal operation of the Energy Hub.
The studied Energy Hub consists of an electric boiler and heat pump that receives refining recovered heat from the thermo-mechanical pulp mill nearby to meet the paper mill heat demand.The Markov chain integrated with a state-space diagram is employed to model the reliability and availability of the system.The proposed Energy Hub design approach is compared with the conventional design to evaluate the superiority of the proposed method.The suggested design approach offers a strong and economical design that isn't influenced by unsupplied penalty rates.The results demonstrate that the suggested method allows for a bigger possibility of decreasing the total cost when the penalty rates are increased.Depending on the penalty rates, the system's overall cost could decrease by 14% − 28% if the suggested design approach is put into practice.As an example, the suggested design approach enables a decrease of 98.10% in the system load loss cost and a reduction of 14.46% in the system total cost for unsupplied demand penalty rates of 10 (€/MWh).
For the best time series prediction model, different long-short-term memory (LSTM) structures have been developed and tested for the best refining heat prediction model.Our experiments demonstrate that all studied LSTM structures are effective at predicting refining heat production.However, Bidirectional LSTM appears to perform better than other models, with correlation coefficient and root mean square error of 0.9 and 0.15, respectively.Due to an accurate load prediction model, the optimal operation of the EH becomes more robust when it comes to fluctuations in heat demand.

Fig. 2 .
Fig. 2. State-Space diagram of a repairable system consisting of one component.

Fig. 3 .
Fig. 3. State-space diagram of the heat supply system consisting of Energy Hub and TMP mill.

B
.Talebjedi et al.

Fig. 4 .
Fig. 4. Operating principle of industrial steam generator heat pump system.
) the ⊗ operator refers to the element-wise product.Where H is activation function, σ denotes the logistic sigmoid function, which is the recurrent step activation function.f t ,i t ,o t ,c t ,anda t , represent the forget gate, input gate, output gate, memory cell, and hidden vector, respectively.W l* = (W lf ,W li ,W la ,W lo ), and W m* = (W mf ,W mi ,W ma ,W mo ), represent trainable weights of the respective gates while b f , b i , b o , andb a are output biases.Fig. 7 shows the single LSTM layer, including the LSTM cells.Each layer in LSTM is composed of the same number of cells as the window size.The LSTM layer learns how time series data and sequence data are related over time.The layer enhances gradient flow over long sequences during training through additive interactions.Different LSTM structures have been examined in this study to achieve the best heat prediction

Fig. 8 .
Fig. 8. Overview of the Single LSTM with a hidden dense layer.

B
.Talebjedi et al.

Fig. 11 .
Fig. 11.The hyperparameter dictionary and optimal values for the RandomizedSearchCV: a)Single LSTM with a hidden dense layer, b)stacked LSTM with a hidden dense layer, c) Bidirectional LSTM.

Fig. 12 .
Fig. 12.The methodology flow chart to obtain the accurate TMP heat prediction model.

B
.Talebjedi et al.

Figs. 13 -
15 illustrate the prediction results of the different studied LSTM structures on the training and testing datasets.The blue data points represent the primary measurements, divided into two parts: training and testing dataset.Since the prediction model is a time series, test data are not chosen at random.The test data are placed in the continuation of the training data so that the training data includes the first seventy percent of the total data set and the rest of the data constructs the test samples.The green data are the predicted TMP heat production on the training data set.An initial set of data used to train a machine learning model is called training data (or training dataset).
They got a B.Talebjedi et al.

Table 1
Industrial taxation and electricity distribution fees charged by the Caruna company in Southern Finland.

Table 2
Performance evaluation of different studied LSTM structures.
Fig. 13.Single LSTM prediction results for train and test datasets.B.Talebjedi et al.

Table 3
Limiting state probabilities of the entire heat supply system.

Table 4
Optimal Energy Hub subsystem capacities for the different design approach.

Table 5
Energy Hub costs for the first design approach (case 1).

Table 6
Energy Hub costs for the second design approach (case 2).