Abstract
In the recent decade, renewable energy sources (RES) such as wind and photovoltaic (PV) power generations gained more attention. However, despite their proven role in the reduction of operating costs of the network, their integration has some serious challenges such as the inherent uncertainties and the energy balance at the power grid. Some options to face these challenges, as well as to enhance the network’s flexibility, are demand response (DR), distributed FACTS (D-FACTS) devices, and energy storage systems (ESS). This paper focuses on the solution of the unit commitment (UC) problem in such a way that in addition to increasing the transmission network’s flexibility, the operation cost of the network decreases. To do this, distributed series reactors (DSRs) are employed as promising D-FACTS devices to enhance the transmission network’s flexibility. The uncertainty of RES is handled by the scenario-based stochastic programming technique. A mixed-integer linear programming (MILP) approach is developed, with a guaranteed global optimal solution, by developing a convex model for power flow in the presence of DSRs. The proposed UC model is implemented on the IEEE RTS 24-bus system, and the obtained results in different cases show the ability of the abovementioned flexibilities to reduce the operation cost of the network.
1. Introduction
1.1. Background and Motivation
The environmental concerns about conventional fossil fuel-based generation units in power systems have been rising in recent years. Hence, many electricity networks have begun to increase the share of renewable energy sources (RES) in line with the national and international sustainability goals. Uncertainty in the generation of power by RES may upset the balance between power generation and demand and endanger system security and increase the network operation costs [1]. In addition, there is a need to more ramping capability by thermal units in the presence of high levels of RES, which ultimately reduces the proper utilization of the existing network. Accordingly, the major challenge for power system operators is how to operate the network to gain the benefits of RES, as well as to better utilize the existing transmission network. According to the existing studies, there is a need for a flexible power transmission grid to handle this challenge and determine the optimal hourly scheduling of both generation and transmission networks. Proposed solutions include the use of energy storage systems (ESS) [2], the use of demand-side programs such as demand response programs (DRPs) [3], and the use of flexible AC transmission system (FACTS) devices to add more flexibility to the transmission grid [4]. By proper coordination between these flexibilities via a unit commitment (UC) framework, transmission system congestion can be mitigated. Consequently, more RES can be injected to the network as well as the operational cost the entire system will be reduced.
1.2. Literature Review
Numerous studies have been conducted around the world, each addressing one aspect of this challenge facing the electricity industry. For instance, in [5], a stochastic computational model for D-FACTS devices and their study on power flow in power systems was presented. In [6], the application of distributed series reactors (DSRs) in reducing line congestion and consequently congestion costs in the electricity market has been introduced. In [7], by increasing the load on the lines in normal conditions and contingency of units, the effect of adding distributed FACTS (D-FACTS) devices on a specific part of the transmission network was discussed. Reference [8] deals with the optimal energy storage capacity in smart grid performance. In this paper, the capacity of energy storage utilization has been reduced via demand response (DR). In [9], the aim is to add wind farm generation that has the nature of uncertainty as one of the sources of the reservation to the problem. To match the nature of the uncertainty of wind products with the network, the DRP has been used. Reference [10] deals with how to combine demand dispatch with forecasting wind power generation, which helps to coordinate wind farm production with the power system in the operation of power systems. Reference [11] has performed the DRP with a direct load control method to reduce operating costs via mixed-integer programming. In [12], DRPs with plugged-in electric vehicles are considered as security constraint UC by a two-stage stochastic mixed-integer programming problem. The uncertain demand in the presence of modular FACTS devices was investigated in [13].
1.3. Contributions
In this paper, the DR, ESS, and DSR as D-FACTS are coordinated to improve the power system flexibility in the presence of uncertain RES. To handle the uncertainty of RES, scenario generation and reduction are used, and a scenario-based stochastic model is developed considering the impact of the flexibilities via a comprehensive network-constrained UC model, which aims at decreasing the operation cost and improving the network’s power transfer capability.
The aim is to develop a novel and computationally efficient model for UC by considering the flexible transmission network via DSR incorporation. By considering the DSR, power flow in transmission lines can be controlled in both directions, and consequently more flexibility will be available to obtain lower operational costs. However, handling power flow constraints in the UC model in the presence of DSR could be challenging, as these constraints will no longer be linear. With respect to the existing literature, this work proposes a novel linear and computationally efficient model for the network-constrained UC problem in the presence of D-FACTS (such as DSR). By using the developed model, it is possible to determine the optimal scheduling of D-FACTS in coordination with other flexibilities such as DR and ESS through a standard network-constrained UC mode. As the developed model is a mixed-integer linear programming (MILP) optimization model, the global optimal solution can be attained.
1.4. Paper Organizations
This paper is organized in five sections. In Section 2, the UC problem formulation is given by considering RES and ESS constraints. In Section 3, a convex model is developed for UC by including the DSRs in transmission network. Numerical results are investigated and analyzed in Section 4. Finally, Section 5 concludes the paper.
2. UC Problem Formulation by considering RES and ESS
In this section, a scenario-based stochastic model for the network-constrained UC problem is formulated, by considering the impact of RES, DR, and ESS.
2.1. Objective Functions
The objective function of the problem includes the total operation costs of thermal units, DRP incentive costs, and curtailment cost of renewables (i.e., wind and solar), which is expressed in (1)where 1, 2, and 3 are simply weighting factors to include the impact of each cost term in the overall objective function, OF. In (1), the first term consists of the fuel cost of thermal generation units that is described in a piecewise linear form in (2)–(10) [14] and also the shut down and startup costs of units, which are calculated from (11) to (12)
The second line in (1) is the overall incentive paid to the loads, which are participating in the DR program. This incentive function is expressed as follows [11]:
The economic model of DR is defined in (13). This model is used to evaluate the impact of consumer participation in DR on the load profile. It should be noted that the proposed model is based on an incentive DRP, which includes an emergency DRP, which has been more attended in recent years [11].
The third line in (1) is related to the curtailment of wind and solar power, where πs indicates the probability value of each scenario.
It is worth noting that as all terms in (1) are cost components, their simple summation makes the sense of overall operation cost of the network in the UC model. Considering different weights for different terms of cost in (1) means that only the cost coefficient (e.g., manipulating the cost coefficients for thermal power generation) changes. As the aim is not to solve a multi-objective optimization problem (in where there are some conflicting objective functions), considering different weightings for the cost components in Equation (1) does not give any further meaning, and it only means the cost coefficients of different components change by simply increasing or decreasing the energy cost or DRP service or RES curtailment cost. These weighting factors simply show the relative importance of each cost component for the system operator.
2.2. Problem Constraints
Ramp rate limits are given in equations (14)–(18)
Constraints (14) and (17) indicate the range of changes in power generation of each unit, and constraints (15), (16), and (10) indicate the limit on the amount of power change when the units are ON and OFF. The relationship is related to the ON/OFF status of the units, and the logical constraints related to it are expressed in (19)–(21)
If a thermal unit generates power (ON status), it must remain committed according to (22)–(25) to a specified time period
Also, if a unit is OFF, it must remain OFF according to (26)–(29) to a certain time horizon
DR constraints are given in equations (30)–(32)where:
Equations (30) and (31) show the limits on the amount of incentive paid to customers participating in the DRP, and (32) indicates the limit on the number of DR resources.
To model the uncertainty of wind and solar radiation, it is assumed that the behavior of these two models is based on the possible distribution functions of Kayal and Chanda [15]. The scenario is generated for each of the possible distributions using the Monte Carlo simulations, based on the standard deviation and the average of their predicted hourly values. Random wind and solar production quantities are generated using the Monte Carlo method, but due to the high number of scenarios and computational burden of the problem, by using scenario reduction methods, the large set of generated scenarios is reduced to a much smaller set and the problem will be solved faster. Here, the k-means data clustering method is used to reduce the number of scenarios [16]. For the wind farm or solar farm connected to bus b at time t and scenario s, the following constraints are considered:
Equations (33) and (34) are wind power curtailment and wind farm power generation capacity range in which indicates the available power of the wind farm connected to bus b at time interval t and scenario s. Equations (35) and (36) are also related to the solar power curtailment and the power output range of the PV system, where indicates the available solar power generation in the solar farm connected to bus b at time interval t and scenario s.
Also, the ESS charging and capacity constraints are given by (37)–(40) [17]
Equation (41) gives the nodal active power balance of the network, where the power flowing through the transmission lines is expressed in (42) as a linear function of nodal angles
In (41), the first line deals with the generating capacity of thermal units, wind and solar farms, and the power charge/discharge of ESS. The second line consists of two parts; the first part (i.e., ) is the loads that do not participate in the DRP and the remaining part is related to the loads participating in the DRP.
3. Network-Constrained UC Model in the Presence of DSR
The transmission system congestion increases the operating costs of the network. This led the system users to make decisions to prevent congestions in transmission lines or even in some cases from overload in system. D-FACTS devices are cheaper, smaller, and sampler than conventional FACTS devices and can typically be quickly installed on transmission lines. D-FACTS devices are basically divided into three types: distributed static series compensation (DSSC), DSR, and distributed series impedance (DSI) [18]. DSRs are capable of enhancing the available capacity of a transmission system by moving the flow from heavily loaded lines to the lines that have unused capacity. Hence, DSRs can add flexibility to the transmission network if their on/off operational status is determined via a coordinated manner.
In this paper, due to the presence of DSR on the lines, the line susceptance (i.e., ) changes according to the value of the DSR device installed in the line. This flexibility of the DSR can change both the magnitude and the sign of power flowing through the line. Therefore, (42) will be a non-linear equation, as the term is appeared. In the following, a novel formulation is proposed to linearize (42) and hence to keep the convexity of the overall model. First, by defining, the direction of power flowing through the line connecting buses b and n is written in a general form as follows:
The above equations imply that for a given scenario s and time interval t, if the direction of power is from bus b to bus n, then ; otherwise, , which means that the power is flowing in the opposite direction, i.e., from bus n to bus b.
Due to the presence of DSRs on the lines, the line susceptance changes with DSR operation as follows:
In (47), the binary variable is related to the DSRs ON/OFF status on the line connecting buses b and n. Similar to (43) and (44), as can be both positive (i.e., when the power flow direction is from bus b to bus n) or negative (i.e., when the power flow direction is from bus n to bus b), the following equations can describe the sign of :
Thus, from (42) and (48), the power flowing between buses b and n can be written as follows:
By defining , (47) can be rewritten as follows, by considering (50):
By considering (50), this means:
In the above inequality, the product of a binary variable (i.e., ) and a positive real variable (i.e., ) is appeared. Now, by defining , the above inequality can be linearized as follows:
In summary, the non-linear equation (42) is linearized by considering (43)–(46), (48), (49), and (53)–(57). It is worth noting that the proposed reformulation allows the optimization model to determine the optimal direction of power flowing through the branches, be optimal switching of DSRs.
Finally, the line power limit and the maximum number of DSRs that can be operated in the ON state are limited by the following constraint:
4. Simulation Results
4.1. Data
The studied network is the IEEE-RTS 24-bus network. According to Figure 1, this transmission system includes 12 power plants, 2 wind farms in buses 8 and 18, and 2 solar farms in buses 3 and 22, as well as 2 ESS in buses 15 and 19, 34 transmission lines, and 17 load buses [17]. The total network load in the peak hour is 2992.5 MW. A uniform hourly demand profile is considered for all load buses, and their hourly demand factor (as a percentage of the peak demand) is also given in Table 1. Moreover, the hourly price of electricity consumption is also given in this table.
The percentage of the nodal load participation in the DRP (i.e., ζ) is assumed to be 40%. It is also assumed that DSR is installed on all lines. The proposed UC model is implemented in the General Algebraic Modeling System (GAMS) software environment and solved by Gurobi solver [19], which is a high-quality solver to deal with MILP optimization models. DSR flexibility is assumed to be equal to [17].
4.2. Case Studies
In this paper, the aim is to reduce the operating costs defined in (1), so simulations are done in four cases as follows:(i)Case study 1 (UC): In the first case study, the conventional network-constrained UC problem is solved without considering RES, DR, and DSR. This case is the base case for impact evaluation of the aforementioned flexibilities.(ii)Case study 2 (UC-DR): In the second case, the DRP is added to the UC model and its effect on reducing operating costs and the amount of incentives for the reduction of nodal demand, as well as the peak loads, will be investigated.(iii)Case study 3 (UC-DR-DER): In the third study, the presence of distributed energy resources (DERs) including RES and ESS, is studied by considering the inherent uncertainty of RES via the above scenario-based stochastic programming model.(iv)Case study 4 (UC-DR-DER-DSR): In the last study, the proposed MILP model for the inclusion of DSRs in the network-constrained UC model is studied. In this study, all flexibilities such as DRP, RES, ESS, and DSR are considered.
It is noteworthy that in all cases the weighting factors in (1) are assumed to be unit (i.e., 1 = 2 = 3 = 1), which means the same importance of all cost components for the system operator.
4.2.1. Case Study 1(UC)
The total operating cost in this case is $ 507,859.42, which merely is the cost of power generation by the thermal units, as in this case there are no RES and DRP. Table 2 shows the on/off status of the units in different case studies. Units G1-G4, G6-G8, and G12 are committed in the entire operation horizon in all case studies, and there is no change in their commitment. Therefore, in Table 2, only the units that their on/off status has changed in different cases are given.
4.2.2. Case Study 2 (UC-DR)
In this case, the network-constrained UC model is solved by considering the DRP. Table 3 shows the hourly price sensitivity coefficients of the loads participating in DRP. These price elasticity coefficients are adopted from [11].
The optimal value of the total cost in this case is $ 487,808.90, which is 3.95% lower than the cost obtained in Case study 1, as a result of DRP utilization. The on/off status of thermal generation units in this case is given in Table 2. The entities shown in gray are the commitment states in this case, which are different from the base case (i.e., Case study 1). In this case, some units switched off compared to the base case; for example, unit G11 has been completely shut down, and the unit G9 switched off in the interval t16–t22. In this case, the binary variable φb,t is introduced in equations (30)–(32) to determine the buses participating in the DRP. In this case, all 17 load buses are participating in DRP in the entire horizon, which means φb,t = 1 for all 17 load buses. The incentive of customers who participated in the DRP in this case is $ 9991.96.
4.2.3. Case Study 3 (UC-DR-DER)
This study includes the presence of RES such as wind energy and solar energy, as well as ESS. The ESS and DRPs are used to handle the inherent uncertainty RES and their better integration. Two wind farms are connected to buses 8 and 18, each with a capacity of 200 MWm and two solar parks are connected to buses 3 and 22, each with the capacity of 50 MW. To reduce the RES power curtailment, two ESSs, each with a capacity of 200 MWh, are connected to buses 15 and 19. Both the wind curtailment cost () and photovoltaic curtailment cost () are assumed to be 50 $/MWh. ESS data are given in Table 4.
Probabilities related to solar radiation and wind blowing were selected based on the standard deviation and the average of their predicted hourly values in the summer season from Ref [15]. Random wind and solar production quantities are generated using the Monte Carlo method by generating the number of 1000 scenarios; by using scenario reduction methods, the large set of generated scenarios is reduced to 10 scenarios. Here, the k-means data clustering method is used to reduce the number of scenarios [16]. Besides, the probabilities of reduced scenarios are given in Table 5.
A similar wind profile is assumed for both WFs, and the same solar irradiance is considered for both PV farms. The available wind and PV power generations in different reduced scenarios are given in Figures 2 and 3, respectively.
The on/off status of thermal units in this case study is shown in Table 2. It can be seen that the units G10 and G11 have been completely switched off in this case, and the operating cost is $ 43,120.15, which shows a 15.09% decrease with respect to the base case (i.e., Case study 1). The use of DERs has a great impact on reducing the costs. Figure 4 shows the performance of the ESSs in different scenarios at buses 15 and 19. The total amount of the DRP incentive in the third case study is $ 10,167.91.
(a)
(b)
4.2.4. Case Study 4 (UC-DR-DER-DSR)
In this case, the network-constrained UC problem is solved by considering all flexibilities, including DRP, DERs, and DSRs. As aforementioned in Section 3, the scheduling binary variable of DSRs on the transmission lines is modeled by . The maximum number of lines where the DSRs can be ON at any specific interval, i.e., in (59), is assumed to be 29. The obtained optimal cost in this case is $ 427877.17, which shows a 15.74% reduction compared to Case study 1. The commitment status of the units in this case is also given in Table 2.
Table 6 shows the obtained on/off schedule of DSRs in the entire operation horizon. It can be seen that only the 12 lines’ DSRs are switched on/off in the entire horizon, and the remaining 17 lines’ DSRs are switched off entirely. Besides, Table 7 gives the optimal hourly dispatch of thermal units in this case, where the changes in optimal hourly schedules with respect to Case study 1 are highlighted.
Also, the modified hourly demand profile of the network for all cases is shown in Figure 5. According to this figure, it can be seen that the peak demand has been reduced by applying the various flexibilities and the role of DRP in the demand profile modification is considerable, as there are no significant changes in this profile in cases 3 and 4 with respect to Case study 2 (i.e., UC-DR).
Table 8 summarizes the total cost and incentives of DRP in different cases. According to this table, it can be observed that the operating costs of thermal units in the second, third, and fourth cases have decreased compared to the base case. In the third case, by utilization of DERs (i.e., RES and ESS), the operation cost of the network reduces significantly. It is also observed that the cost of encouraging the customers to participate in DRP in the third and fourth case studies has been decreased compared to that of the second case, which shows the impact of DERs and DSRs on the grid’s need to demand modification. This coordinated scheduling of thermal generation units and available flexibilities in the transmission network will result in the lowest cost solutions to the network operators.
4.3. Sensitivity Analysis
In this section, the impact of different weighting coefficients in Equation (1) is studied. Besides, the proposed UC model is solved for different number of reduced scenarios. All analyses are done for the model studied in Case study 4.
4.3.1. Impact of Different Weighting Factors
To check and analyze the results, different simulations are done considering different weight coefficients on the final model. In each simulation, the effect of each part of the OF on the final result can be seen and checked, and the results are given in the following.
The OF, which is the minimization of operating costs, consists of three main parts. The first part is related to the cost of switching on and off and the cost of operating the thermal units, the second part is related to the incentive cost that is given to the consumers who have participated in the DRP, and the third and last part is related to the operating costs and energy interruption of the solar and wind park.
According to Table 9, it can be seen that by considering different sets of the weighting coefficients for the OF in (1), significant changes can be observed in the costs. As the first term of OF, i.e., the cost of thermal generation units has the highest impact on the overall value of OF, it can be seen from Table 9 that in cases where 1 is higher than 2 and 3, the proposed UC problem aims to reduce the thermal generation cost, such that its minimum level is $393520.38 for 1 = 1, 2 = 3 = 0. For the other components of cost, i.e., DRP incentive and RES curtailment cost, the same behavior could be observed. For example, the minimum level of RES curtailment cost is $ 9236.04, which is obtained for 3 = 1, 1 = 2 = 0.
4.3.2. Impact of Different Sets of Reduced Scenarios
In general, the model is based on scenarios generated from solar radiation and wind, which are included in the simulation process using the scenario-based model and k-means scenario-reduction technique. Running the model without scenario reduction is not computationally feasible, despite the convexity of the proposed optimization model. Hence, as aforementioned, k-means clustering method is employed in this paper for scenario reduction. The simulations were done for various numbers of reduced scenarios (5, 10, 15, 20, and 25) for Case study 4 (which is the most comprehensive case by considering all DR, DER, and DSR). The overall cost (i.e., OF in (1)) obtained for different sets of scenarios is compared in Table 10. It is observed from this table that considering different sets of reduced scenarios will result in different values for OF, such that generally by increasing the number of scenarios, the obtained OF increases slightly.
5. Conclusions
In this paper, an MILP model is proposed for the network-constrained UC problem, by considering the impact of DSRs on adding more flexibility to the transmission grid. In addition to DSR, other flexibilities such as DRP, RES, and ESS are also considered to supply the demand, as well as to increase the system ability in congestion relief of transmission network. The conclusions can be outlined as follows:(i)According to the results, it was observed that DR is effective at the peak demand interval and reduced the peak load. Hence, the operational costs reduced by implementation of the incentive-based DR program.(ii)The presence of RES in the grid, which included PV and wind farms with ESS, resulted in a significant reduction in operating costs.(iii)Enhanced transmission network flexibility by DSRs through the proposed convex UC model can result in further cost reduction by the optimal scheduling of DSRs.(iv)By considering different weight coefficients, it is possible to achieve global optimal solutions different from the original solution, but the weight coefficients can be changed according to the opinion of the system operator.
Nomenclature
Set and indices| : | Set of the time period |
| : | Set of all generation units |
| : | Set of generation units connected to bus b |
| : | Set of all buses |
| : | Number of segments in linearizing |
| : | Hour index |
| : | Index of bus and node |
| : | Index of scenario |
| : | Index of block k in linearizing |
| : | Index for thermal generating units |
| : | On/off cost for thermal unit i ($/h) |
| : | Startup and shutdown limits for thermal unit i |
| : | Ramp-up and ramp-down limits in thermal unit i |
| : | Downtime and uptime limits in thermal unit i |
| : | Time periods that unit i has been on/off at the beginning of the scheduling horizon |
| : | Initial status of unit i |
| : | Number of thermal generating units |
| : | Maximum number of DSRs |
| : | Wind/PV curtailment cost ($/MWh) |
| : | Minimum possible variation from the initial suspension of the transmission line l |
| : | Minimum/maximum generation limit of unit i |
| : | Interval length k at the linear cost of unit i |
| : | Maximum/minimum generation limit of unit i at time t |
| : | The initial and final power values of section k in the linear cost of heating unit i |
| : | Initial and final cost values of section k, linear cost of unit i ($) |
| : | Cost slope in section k, related to the linear cost of unit i (MW/h) |
| : | i-th thermal unit’s fuel cost coefficients |
| : | Probability of scenario s |
| : | Maximum charging/discharging power of the ESS connected to bus b |
| : | Maximum power generation capacity of WTs/ PVs at bus b, time t, and scenario s |
| : | Wind/PV power curtailment at bus b, time t, and scenario s |
| : | The slope of the k segment for bus b at time t in the linearization of the DRP |
| : | The effective susceptance of the line connecting buses b and n at time t |
| : | Active power generated by the i-th thermal unit at time t |
| : | Power generated in the wind/solar farm, at bus b, time t, and scenario s |
| : | Auxiliary variables for modeling the DSR’s impact on power flowing through lines |
| : | A variable to categorize the penalty level |
| : | The demand of bus at time t |
| : | Minimum/maximum incentive to bus at time t ($/ MWh) |
| : | Incentive at the beginning/end of section k |
| : | Maximum/minimum power limit across the line connected to buses b and n |
| : | Price of electricity in bus b at time t ($/ MWh) |
| : | Percentage of load that participates in DRPs |
| : | Initial cost of electricity at bus b at time t ($/ MWh) |
| : | Maximum limit voltage angle across the line connecting buses b and n |
| : | Maximum number of demand response sources at time t |
| : | Initial susceptance of the lines connecting buses b and n |
| : | Maximum/minimum energy stored in the ESS connected to bus b |
| : | Maximum charge and discharge power of ESS in bus b |
| : | ESS charge and discharge efficiency (%) |
| : | Minimum/maximum generation limit of unit i at time t |
| : | Weighting factor for each cost term of the objective function |
| : | Self/cross-load elasticity factors |
| : | Start p and shutdown costs of the i-th thermal unit |
| : | State of charge for the ESS connected to bus b at time t and scenario s |
| : | The charging/discharging power of the ESS connected to bus b at time t and scenario s |
| : | Binary variable that shows the on/off status of unit i at time t |
| : | Binary variable that indicates the i-th unit’s startup status |
| : | Binary variable that indicates the i-th unit’s stop status |
| : | Voltage angle difference between buses b and n at time t at scenario s |
| : | Auxiliary variable defined for the voltage angle difference between buses b and n at time t at scenario s |
| : | Active power flowing the line that connects buses b and n at time t at scenario s |
| : | Operation planning of unit i in section k at time t |
| : | The demand of bus b at time t after the DR application |
| : | Optimal incentive to the buses ($/ MWh) |
| : | Optimal incentive to bus in section k ($/ MWh) |
| : | Total incentive for DRP |
| : | Binary variable shows the participation of bus b at time t in the DRP |
| : | The difference between consumer demand and demand participation in the DR program |
| : | Demand response time (h) |
| : | Auxiliary variable defined for the active power flow of the line connecting buses b and n at time t and scenario s |
| : | Auxiliary binary variable defined for modeling the active power in the line connecting buses b and n at time t and scenario s |
| : | Binary variable related to the DSRs on/off status on the line connecting buses b and n at time t |
| : | Fuel cost of thermal unit i at time t. |
Data Availability
Data are contained within the paper.
Disclosure
This is an academic research within the University of Zanjan.
Conflicts of Interest
The authors declare that there are no conflicts of interest regarding the publication of this paper.
Acknowledgments
The research was supported by the University of Zanjan as the primary affiliation of all authors.