Thermoeconomic multiobjective optimization of tobacco drying heat pump recovering waste heat from monocrystal silicon furnace based on SVR ANN model in Southwest China

Great quantities of energy are consumed by heat pump drying systems (HPDSs), such as tobacco drying. It can be reduced by increasing evaporating temperature using waste heat. In this paper, based on the case of two main industries in southwest China (monocrystalline silicon industry and tobacco industry), a novel integrated HPDS especially for the tobacco drying using waste heat from single crystal silicon furnaces is proposed. It is then analyzed from an energetic and economical perspective, investigating the effects of heat exchanger areas on system performance. To simultaneously minimize system CO2 emissions and total annual cost (TAC), a multiobjective optimization using nondominated sorting genetic algorithm is conducted based on the support vector regression model, which has been rarely reported for integrated HPDS improved by waste heat. The final optimum solution on the Pareto optimal front is selected using different decision‐making methods (LINMAP, TOPSIS, and fuzzy decision method), and are compared with the performance of the initial system. The results show that the TOPSIS method is most suitable in three decision‐making methods with a great reduction in m C O 2 ${m}_{C{O}_{2}}$ (−12.7%) and nearly unchanged in TAC (+0.5%) compared with the initial solution. And the decreases in m C O 2 ${m}_{C{O}_{2}}$ and TAC of the optimized system is as large as 50% and 1/3, compared with traditional vapor compression heat pump dryer without heat recovery of the monocrystal silicon furnace.

Heat pump drying system (HPDS) has attracted increasing attention due to its high efficiency, good drying quality, and strong control ability, compared with the fuel oil/gas heating dryers or electric heating dryers. 2,3 Duan et al. 4 proposed a novel enclosed cascade-like HPDS in the process of hawthorn cake drying, and compared it with a traditional hot air dryer (HAD). The comparison results showed that the overall performance of HPDS system reached 3.45 and the energy consumption decreased by 32.55%. Seyfi 5 did experimental analyzes on various agricultural products under different climatic conditions, which showed that HPDS achieved the best performance with no significant damage to dried products. Generally, the application of HPDS is helpful for promoting energy saving and emission reduction for realizing the goal of carbon neutrality. However, traditional heat pump performs poorly with low ambient temperature and high heating temperature, as they may lead to an increase of pressure ratio. Additionally, frost often unavoidably yields on the evaporator outdoor coils at low ambient temperatures, further degrading the heating performance.
Nowadays, more than 72% of the overall energy consumption is dissipated into the environment as waste heat, 6 considering that, it is a suitable way to strengthen system efficiency by means of recovering possible sources of waste heat to increase the evaporation temperature. Recovering waste heat energy is becoming an interesting domain of research to reduce the use of fossil fuel resources. Soylemez 7 carried out a thermoeconomic optimization analysis on HPDS with waste heat recovery, aiming to find the optimum operating temperatures and the optimum sizes of system components. Tan et al. 8 compared the energy performance of four types of heat pumps in low-temperature waste heat recovery, and investigated the impact of operating parameters on system coefficient of performance (COP), exergy efficiency and annual savings of standard coal. Ge et al. 9 designed a metal hydride heat pump for recovering heat from steel plants, and compared different alloys based on operating temperatures, system efficiencies, and thermodynamic equilibriums.
The hybrid of waste heat recovery can increase the system performance but will also increase the cost, as the system thermos-economic performance is largely determined by the heat exchangers. To have a deep understanding of the system comprehensively, it is necessary to conduct an optimization study. Nowadays, as rapidly emerging technology, the intelligence algorithm has attached much importance to system prediction and optimization in energy and process systems. Esen et al. 10 applied the support vector machine (SVM) method in the prediction of the COP of a ground source heat pump. They found that SVM had a faster computation speed and used fewer free parameters compared with other intelligent algorithms. Guo et al. 11 compared the performance of the support vector machine for regression (SVR), back propagation neural network (BPNN), multiple linear regression (MLR), and extreme learning machine (ELM) models in forecasting energy demand of building heating system, whose feature set is optimized through correlation analysis and least absolute shrinkage and selection operator (LASSO) approach. Jain et al. 12 developed a sensor-based forecasting model based on SVR to evaluate the influence of temporal and spatial granularity on the predictive power of our single-step model. Cao et al. 13 conducted a novel ejector enhanced heat pump system from the thermodynamic viewpoint, and optimized the system by nondominated sorting genetic algorithm (NSGA-II). Li et al. 14 established a simulation-based multiobjective optimization model for solar hybrid heat pump heating and cooling system based on NSGA-II, as the primary energy consumption and the annual average cost are selected as two objectives.
In southwest China, YunNan Province, the tobacco industry is a pillar industry, consuming a lot of energy in tobacco drying by HPDS with a high heating temperature around 80°C (Figure 1). Also in this area, there are many single crystal furnaces operated by PV manufacturers. Because they can be easily powered by several of the world's top hydroelectric plants here, with the great vertical drop provided by Qinghai-Tibet plateau. The large amounts of ultra-low grade waste heat can be recovered (<40°C) in single crystal silicon industry, that is carried by water cooling the steel shell of furnace.
Therefore, we proposed an integrated HPDS system for the drying application of tobacco, recovering waste heat from a typical single crystal silicon factory in southwest China. To the best of the author's knowledge, the optimization and parametric study on the waste heat F I G U R E 1 Tobacco and silicon industries in southwest China. recovery heat pump drying technology is still scarce. So, we first establish and experimentally validate the thermodynamic and economic model of the system, and then introduce SVM to predict the system performance, which can largely save the computation time. Finally, a multiobjective optimization study is conducted based on NSGA-II, in which the system's annual CO 2 emissions and total annual cost (TAC) are selected as two objective functions. The Pareto front is achieved and three different decision-making techniques are employed to find the suitable solutions. Furthermore, the optimum results are compared with the initial system (not optimized) and traditional HPDS, for having a deep understanding of the proposed system.

| Case information and system description
The single crystal silicon factory is divided into four workshops, each workshop holds 672 single crystal furnaces with a model number of KLS-160. The factory totally consumes about 160 MW electric power when it is working at full capacity. Each furnace can produce 4.5t single crystal silicon. It can continuously run 15 days and needs several days of preparation before next operation. Therefore, in the proposed integrated HPDS, every sour furnaces are in a group to keep the system operating stably and uninterruptedly, with one under preparation and three running. Average waste heat power for furnace group is about 165 kW. Three heat pump dryers are connected in parallel, recovering waste heat and providing high-temperature heat to three tobacco drying rooms simultaneously, as shown in Figure 2.
The heat pump cycle is a single-stage vapor compressor cycle, which includes a compressor, a condenser, an expansion valve, and an evaporator. Unlike traditional air-source heat pump, the evaporator employs waterwater heat exchange to recover the heat from a single crystal furnace. The water cycle consists of a water pump and a fan-coil unit, to deliver the heat from the condenser to the drying room for tobacco drying. It is noted that the heat pump in this paper is charged with refrigerant R1234zez, which can offer a high condenser temperature. The refrigerant is first compressed in the F I G U R E 2 Schematic diagram of the waste heat recovery integrated heat pump drying system. HOU ET AL.
| 2353 compressor to a superheated state and then transfer heat to the fan-coil unit via the condenser. Then, the refrigerant flows into the expansion device throttling to state 5. Finally, it absorbs the heat from the waste heat recovered from the single-crystal furnace in the evaporator to saturated states and returns to the compressor. The novel system recovers heat from wasted water, thus having high efficiency. The proposed hybrid HPDS system is expected to perform well in the drying application of tobacco.

| Thermodynamic model
Several main assumptions were followed in the system thermodynamic model: (1) The whole system was worked in steady state.
(2) Pressure drops of working substances are neglected in heat exchangers and pipelines. (3) The system was well insulated. (4) Throttling processes were isenthalpic. (5) The irreversibility of compressor was described by isentropic efficiency η s and the volumetric flowrate of compressor was described by volumetric efficiency η v , as given by Equations (1) and (2). 15 (6) The water supply and return temperatures in the condenser were both equal to those of the fan-coil unit, which means that the water temperature changes in the two heat exchangers were the same. (1) In Equations (8) and (9), Pr is compressor pressure ratio, which is determined by Equation (10) where P 1 and P 2 are pressures of inlet and outlet of compressor, respectively. With the isentropic efficiency of Equation (1), compressor power W com can be given by In Equation (4), m r is mass flowrate. With theoretical displacement q vs , and refrigerant density ρ 1 , it can be calculated by Another part that consists of the power of the whole system is pump power W pump , which is assumed to be 10% of the compressor power. 16 The heating capacity Q h of a single heat pump cycle in this paper is calculated by the below Equation The logarithmic mean temperature differences of condenser T Δ c , evaporator T Δ e , and fan-coil unit T Δ f are shown as ( ) where T is the temperature of state point, and subscripts of c, wo, wi, e, ri, ro, ai, and ao stand for condenser, water inlet, water outlet, evaporator, waste heat inlet, waste heat outlet, air inlet, and air outlet, respectively. The system annual electricity consumption E total and annual heating energy Q total are calculated as follows, respectively where W ele is electricity power the whole system, being the sum of W com and W pump . And t op denotes to the yearly operating time. System COP can then be expressed as Equation (12): Besides energy consumption, CO 2 emission (m CO 2 ) is another key indicator, which is one of the two optimization objectives. In this paper, m CO 2 consists of the emission from electricity consumption and the emission from the material of heat exchangers m Cu and m Al (copper and aluminum). As only the heat exchanger areas of evaporator, condenser, and fan coil unit are chosen as decision variable in the optimization, the emissions from the manufacture of other components need not be calculated. To calculate the carbon emission, emission conversion factor of electricity (λ ele ), copper (λ Cu ) and aluminum (λ Al ) are used, with their values set 0.968 kg/kW, 0.07 kg/kg, and 0.07 kg/kg, 17 respectively, as shown below

| Economic model
TAC is introduced as the other of the two optimization objectives, evaluating system feasibility. As shown in Equation (14), it is comprised by annual capital cost Z inv+main , annual operating cost Zȯ p , and carbon emission penalty Zṗ en .
The annual capital cost of the whole system can be obtained by evenly distributing the total capital cost into all operating years, according to the below equation where Z sum is the sum of investment costs of all components, which can converse to total capital cost by multiplying the maintenance factor φ (specified as 1.05). Capital recovery factor (CRF) can be determined with interest rate (i = 3%) and system life time (n = 10), according to literature. 18 Investment cost of condenser (Z gc ), evaporator (Z e ), fan-coil unit (Z fan ), compressor (Z com ), pump (Z pump ), water tank (Z wt ), and expansion valve (Z ev ) are determined by Equations (16) Additionally, the cost of accessories, such as check valves, and connecting pipes, is added to the total investment cost, which accounts for about 1% of the cost of the whole system. 22 The operational cost rate (Z op ) and penalty cost of CO 2 emissions (Z pen ) are determined by Equations (23) and (24), respectively.
where C ele is electricity price, and C CO 2 is the penalty tax of unit CO 2 emission (0.22$/kg). 23

| Model validation
To validate the established thermodynamic model, the calculated COP value is compared against the experimental COP value of a vapor compression heat pump, in the evaporating range of 15-35°C. Mean relative error (MRE) and the largest relative error (LRE) are shown in Equations (25) and (26) as two criteria, where K sim is the simulated value, and K exp is the experimental value. indicate that the established thermodynamic model is feasible.

| General optimization procedure
During multiobjective optimization process, thermodynamic and economic models will be called many times. As the models contain many nested iterations in themselves, particularly iterations for thermophysical properties, the calculation takes a very long time. Learning machines can quickly predict system performance after being trained by a data set, without iteration. Therefore, in this paper, first, heat pump performance data set are calculated by thermodynamic and economic models. SVR artificial neural network is then trained using the data set. Finally, the optimizing algorithm calls the trained SVR model during optimization, saving a lot of time.

| Support vector machine (SVM)
SVM is a powerful machine learning method based on statistical learning theory. The core principle of the SVM is presented to build an optimal hyperplane as a decision surface. Based on the principle of structural risk minimization, the nonlinear problem is transformed into a quadratic optimization problem. Because of its elegance in solving nonlinear problems with the "kernels" technique, good performance with small data set, SVM is used in this paper after pretest.
To employ SVM to solve regression problems, Vapnik extended it to estimate the function which is called SVR. SVR implies a nice resistance on overfitting and excellent generalization abilities. 24 As a result, it is considered as one of the best approaches in regression from little-size data sets.
The SVR method introduces the kernel function K x y ( , ) for addressing nonlinear regression problems. The introduction of kernel function can translate the nonlinear problem into a linear problem in highdimensional space, and decrease the computation time.
In this paper, we apply the radial basis function as the kernel function. It can be infinite polynomials by extending the RBF kernel function, which means that the original features can be mapped into infinite dimensions, therefore, RBF kernel function can perform better and use fewer parameters than polynomial kernel function. The equation is set as below: where γ represents the width of Gaussian kernel function In SVM, researchers need to find an ideal classification surface to separate two samples, however, the core principle of SVR is to search for the optimal hyperplane for minimizing error of all training samples from the classification surface. The function can be considered as a combination of empirical risk and structure risk: Subject to constraints: This function can be further transformed to a dual form which needs to maximize the following functions: Subject to constraints:

| Nondominated sort genetic algorithm
In previous studies, the system is often optimized based on the single-objective optimization method, which only aims at maximizing or minimizing a single objective. To simultaneously minimize the two key objectives, TAC and yearly CO 2 emissions (m CO 2 ), the NSGA-II is employed. The NSGA-II was first proposed by Deb et al. 25 in 2002, which was used to solve the disadvantages of NSGA, such as the large computation time and the difficulty in specifying the sharing parameter. NSGA-II is based on the Pareto optimal concept of Goldberg's thinking.
NSGA-II introduces a new selection mechanism based on nondominated sorting and crowding distance calculation. Generally, all individuals in the same population are ranked based on their nondominated level. It is defined that individual j dominates individual k when one or more objectives of j are better than those of k and the rest objectives of j equal those of k. Individuals in the same rank cannot dominate each other. For a series of individuals in the same nondominated rank, the individual with a higher crowding distance is easily selected to enhance the diversity of the population. The crowding distance of each individual in the same nondominated set is shown as Equation (36). Obviously, all crowding distance in a certain rank with no more than two individuals will be ∞.
To protect the population diversity and avoid the loss of favorable features, NSGA-II combines the parent solution and offspring solution to build a new population. Each individual in the new solution is compared with remaining individuals, to sort them based on the nondominated level and choose the best individual for the new parent population. The specific steps of NSGA-II in this paper are described in Figure 4, where SVR model is employed to replace numerical model to compute two objectives.
Different from a unique solution by the single-objective optimization, a set of all the optimal solutions is obtained by the multiobjective optimization, in which all the individuals are not dominated by any other individuals. This set is called Pareto front. However, not all feasible solutions in Pareto front are reasonable for practical applications, as decisionmaking processes are often employed to find the best compromise solution in the Pareto front. In this paper, the fuzzy decision method, TOPSIS and LINMAP are applied in parallel to find the optimal solution.
The fuzzy membership function is used to represent the optimization degree of the objective functions. 26 To describe the difference between multiple objective functions, it is essential to introduce objective weight depending on the membership degree of the objective function, and then calculate the corresponding solution. 27 Linear programming technique for multidimensional analysis of preference (LINMAP) is a method suitable for multiobjective decision-making, which can easily reflect the preference and experience of the decision makers. 14 As given by Equation (37), the individual with the minimum distance to the ideal point (d i+ ) is considered as the best solution.
TOPSIS method is one of the most used decision-making processes in the literature, which is based on an aggregating function representing the closest to the ideal solution point and farthest to the nonideal solution point (min d

| Optimization model with SVR
During the training process, 3289 data sets calculated from the thermodynamic and economic model are input into the SVR model to search for the optimal hyperplane. The datasets are divided into 2800 training sets and 489 testing sets proportionally. The input parameters are the aera of heat exchangers (evaporator area, condenser area, and fan-coil | 2357 unit area). The output parameters are TAC and CO 2 emissions. The optimized parameters of SVR model are presented in Table 1.
After applying the greedy search algorithm for RBF kernel, the obtained optimum c and g are 32 and 0.03125, respectively. As observed in Figure 5, it can be seen that the optimal SVR performance achieved a mean square error (MSE) of 0.0021738.
Finally, test datasets are evaluated in the established SVR model to examine the model accuracy and generalization. Figures 6 and 7 show the comparison between the predicted value from SVR and the true value calculated from the numeric model. It can be obviously observed that the predicted value has the same trend as the true value, which presents the high accuracy of the established SVR. R 2 of 0.9985 for TAC and 0.9998 for m CO 2 were obtained with the proposed SVR model. The results show that SVR is a good method for predicting the system performance with high accuracy and good generalization, and the established SVR model in this paper is applicable for the following optimization work in replacing the numerical model.
The calculation time of SVR is compared with that of thermodynamic model, to verify the benefit of using SVR. In each data set, there are 100 elements. Five different data sets are calculated by the two models, respectively. As shown in Table 2, the average time cost of SVR is 0.458 s, while the one of thermodynamic model is 129.55, which is more than 280 times larger than that of SVR. These results show SVR can greatly reduce calculating time.

| Parameter setting for NSGA-II
In this paper, the evaporator area, condenser area and fan-coil unit area have significant impacts on both T A B L E 1 Basic parameters of SVR.

Parameter Value
Kernel type RBF function Step c 312,606.1 Step g 369,357.7

Max c 10
Min c −10

Max g 10
Min c −10

Gamma Auto
Cross validation check 5 Step MSE 0.06

Shrinking True
F I G U R E 5 Variation of MSE with different parameters of c and g.
F I G U R E 6 Comparison results of TAC between true value and predicted value. TAC, total annual cost.
F I G U R E 7 Comparison results of m CO2 between true value and predicted value.
system energy efficiency and economic cost. They are therefore chosen as decision variables, which are presented in Table 3. In the initial traditional design of the drying system, their values are chosen according to the temperature difference between the hot stream and cold stream (for water-water heat exchanger, it is set to be 4°C, for air-water heat exchanger, it is set to be 8°C).
The multiobjective optimization in this study is carried out by MATLAB 2020a programming. The basic parameters of NSGA-II are set based on some pre-tests done by the authors to achieve a reliable optimal result with less computational time, as shown in Table 4. For operators, the selection method chooses the binary championship method, and the crossover method and mutation method employ the single-point crossover method and the single-point mutation method, respectively. The hypervolume diagram (HV) of NSGA-II in the pretest calculation is shown in Figure 8. In the 100th generation, the HV value reaches a stable area, indicating a maximum generation of 120 is reasonable.

| Parametric analysis
In this part, a parametric analysis is conducted to investigate the effects of the heat exchangers areas on the system performance, for having a full understanding of the proposed system. Table 5 shows basic parameters from a typical single furnace factory for simulation and optimization purposes. Figure 9 represents the effects of condenser area A c , fan-coil unit area A f , and evaporator area A e on COP, which indicates that COP increases with all the increases of them. It can be explained that there is a decrease in the temperature difference between refrigerant, water, and air with increasing heat exchanger areas, resulting in higher system efficiency. The trend of COP increases sharply at first and becomes slower gradually with the variation of heat exchanger areas because the decrease of temperature difference becomes slower when the heat exchanger area is large. Besides, the condenser area has the most significant impact on COP in given ranges, as  Figure 10 represents the effects of condenser area, fan-coil unit area, and evaporator area on TAC. Evidently, TAC rises quickly with the increase of A c and A e , since the growth of component costs is faster than the decrease of operational costs. The influence of A f on TAC is more complicated, as TAC decreases first and then increases, which indicates that there is an optimal fan-coil area for heat pump system to F I G U R E 9 The effects of the condenser area, fan-coil unit area, and evaporator area on coefficient of performance.
F I G U R E 10 The effects of condenser area, fan-coil unit area, and evaporator area on total annual cost (TAC). achieve a relatively low capital cost with high efficiency. Figure 11 represents the effects of the condenser area, fan-coil unit area, and evaporator area on life cycle CO 2 emissions. It can be observed obviously that m CO 2 decreases with the increase of heat exchange areas, which is related to the decrease in system efficiency. To be specific, when A c increases from 80 to 300 m 2 , A e increases from 80 to 200 m 2 , and A f increases from 300 to 700 m 2 , m CO 2 quickly decreases from 4.76 × 10 5 to 3.12 × 10 5 kg. Figure 12 depicts the Pareto optimal curve achieved by the NSGA-II algorithm, which clearly reveals the conflict of two objective functions: TAC and m CO 2 . As can be seen from the figure, the reduction in TAC is accompanied by the growth in m CO 2 in the Pareto front, with a range of 3.13 × 10 5 to 3.68 × 10 5 kg in m CO 2 and 6.04 × 10 4 to 7.67 × 10 4 $/year in TAC.

| Optimization results and analysis
To find the final optimum solution among the Pareto front, three decision-making processes introduced in Section 4.2 are utilized based on different needs and F I G U R E 11 The effects of the condenser area, fan-coil unit area, and evaporator area on m CO2 .

F I G U R E 12
The optimal Pareto front curve of multiobjective optimization.
interests. Since the two objective functions in this study have different dimensions, the Euclidian distance cannot be calculated directly. It is noted that the Euclidian technique is applied to convert the two objectives to nondimensional parameters, which is shown as Equation (39): (39) Figure 13 displays the nondimensional Pareto fronts, in which three different solutions were chosen for different decision-making methods. It can be seen that the lowest TAC of 63304.6 $/year can be achieved by LINMAP method while the lowest m CO 2 of 319,706.6 kg is achieved by fuzzy decision method.
The TOPSIS method can find a solution relatively more compromised, with a TAC of 64172.3/year and a m CO 2 of 326137.8 kg, in the respect of 140.71 m 2 A e , 118.82 m 2 A c , and 700 m 2 A r .
The scatter distribution is also achieved to find the variation of optimum design variables as to gain a deeper understanding of the system. From Figure 14A, the evaporator area A e floats in the given ranges (80-200 m 2 ), implying that two objectives have much conflict with the change of A e . Conversely, the optimum condenser area A c and radiator area A r tend to stay in the lower and higher bounds of the given ranges. It is therefore believed that decreasing A c and increasing A r will be helpful to enhance the system performance.

| Comparison of the optimized system against other systems
From Figure 15, it is obvious that proposed heat recovery system can reduce more than 50% carbon emission and 1/3 economic cost than traditional system, even without optimization. With the waste heat from cooling water of monocrystal silicon furnace, not only the carbon emission and operating cost are greatly reduced, but the capital cost of evaporator is also reduced. It is certain that the F I G U R E 13 The nondimensional optimal Pareto front curve of multiobjective optimization. integration of monocrystal silicon furnace and heat pump tobacco dryer can be recommended to achieve benefit in the economy and environment. Figure 15 also compares the solutions selected from pareto front, using the three different methods. The solution by fuzzy decision has slightly larger TAC (particularly a larger economic cost of condenser) and slightly smaller carbon emission. The solutions of TOPSIS and LINMAP seems similar. Generally, the difference between the three methods is not significant. But when these results are compared against the ones from the initially designed heat pump heat recovery system without optimization, a significant decrease in carbon emission can be found, as TAC almost keeps unchanged (about 0.5% variation). For example, carbon emission by TOPSIS is 3.22 × 10 5 kg, a 12.7% decrease is achieved, comparing with the initial system (3.69 × 10 5 kg). The results indicate the optimization work effectively improves the environmental performance of the heat recovery system, further increasing its application potential.

| CONCLUSIONS
A novel configuration of single-crystal furnace waste heat recovery HPDS is proposed, evaluated and optimized in this study. The analyzing results demonstrate that the proposed system is efficient and economical for the drying purpose of tobacco, which performs much better than traditional drying systems. The major conclusion of this paper is set as below: 1. A novel waste heat recovery hybrid HPDS is proposed and analyzed. The simulation results are comprehensively compared with the experimental results to ensure accuracy. The mean relative error and largest relative error between experimental and simulation results are 4.00% and 7.28%, respectively. Furthermore, the impacts of heat exchanger areas on the system performance are deeply investigated. 2. To mitigate the prolix iterations in the thermodynamic model, SVR is introduced for simplifying the calculation process. The best parameters of c and g are selected to search for the optimal MSE of 0.021738. The calculation time can also be largely reduced using SVR. 3. The multiobjective optimization is carried out using NSGA-Ⅱ for minimizing CO