Abstract
Considering the increasing penetration of renewable energy sources (RESs) into power grids, adopting efficient energy management strategies is vital to mitigate the uncertainty issues resulting from the intermittent nature of their output power. This paper aims to resolve the energy management problem by presenting an adaptive robust optimization (ARO) model in which uncertainties associated with solar and wind powers, consumer demand, and electricity prices are considered. The proposed model comprises three stages, one master and two subproblems, using both linearized and convexified power flow equations. Meanwhile, a new optimization-based bound tightening (OBBT) method is also presented to strengthen the relaxation. The proposed solution method solves the microgrid (MG) operation management problem with high accuracy due to considering the convexification method in a reasonable computational time and reducing operational costs compared to conventional models. The numerical results indicate the benefits of the proposed ARO and solution approach over traditional methods.
1. Introduction
During recent years, owing to many grounds such as environmental and economic concerns, microgrids (MGs) have been receiving the attention of both academia and industry. MGs usually consist of several distributed generations, including renewable and nonrenewable sources and energy storage systems such as batteries and DC/AC loads [1]. The reliable and economic operation of MGs is of paramount importance to system operators, and the existing uncertainties are a major challenge that should be addressed [2, 3]. The intermittent nature of renewable energy sources (RESs) combined with unforeseen variations in consumer demand and electricity price, which might be resulted from forecast errors, are the common sources of uncertainties in the MG operation. These uncertainties complicate the operation of the MG and impose serious problems on the system operators.
1.1. Research Gap
MG operation optimization problem can be solved using metaheuristic and mathematical methods. Metaheuristic optimization methods involve using optimization algorithms inspired by natural processes, such as genetic algorithms, grey wolf optimization, and cuckoo search [4, 5]. This includes optimizing DG dispatch, such as PV, WT, and energy storage systems, to minimize costs, reduce carbon emissions, and ensure a stable power supply. The metaheuristic method is a type of soft computing that is a branch of artificial intelligence that deals with developing intelligent systems that can reason, learn, and adapt in uncertain and imprecise environments [6, 7]. Metaheuristic optimization algorithms are heuristic search methods that do not guarantee to find the optimal solution. On the other hand, mathematical modeling can provide a more precise and accurate solution, especially for problems that can be accurately represented mathematically and for which the solution space is well-defined and tractable.
In the presence of uncertainties, deterministic techniques are unable to ensure a reliable operation meeting the consumers’ demand [3, 8]. Some techniques have been recently used to model the discussed uncertainties, such as stochastic programming (SP) [9–12]. This method requires using probability distribution functions (PDFs) to generate numerous scenarios. It is worth mentioning that although SP approaches can study uncertainties well through generated scenarios, it remains a challenge to find the PDF of some uncertain parameters primarily because there is not enough historical information about uncertain parameters. Besides, a vast range of scenarios is required to achieve an efficient model, which in turn enormously increases the computational time.
Compared with the SP, which requires the PDFs of all uncertain parameters, the RO approach focuses on the largest extent of variability among uncertain parameters. Hence, it is particularly attractive for the continuous uncertain parameters such as PV and WT energy, consumer demand, and price of electricity [13]. This approach does not necessitate generating a significant number of scenarios [14, 15]. The authors in [16] suggest an energy management system (EMS) based on robust optimization (RO) that aims to reduce the impact of uncertainties associated with both RES and load demands in MGs. The uncertainties of RESs, demand, and price of electricity are addressed in [17] with the RO approach, where the expected profit of the MG is maximized. In [18], a new p-robust optimization technique is developed for the energy management of a combined power and heat system. Their proposed technique strengthens the robustness of the MG power management problem, and interestingly enough, no substantial fall in the MG expected revenue occurs. In [19], the authors succeeded in enhancing the robustness of power scheduling for the MG while maintaining expected revenue levels with only minor reductions. Besides, the authors in [20] conduct a survey on recent developments in RO and their application for the optimal planning and operation of power systems. DG placement in MGs is studied in [21] using a two-stage RO framework. The authors of [21] prove that the RO approach is more practical than conventional MG planning approaches and can well model the probabilistic nature of the load, PV, and WT powers. The uncertainty of loads is of the essence, particularly when demand response programs are applied. In [22], a RO scheme is presented to model the uncertainty of WT generation and elasticity of demands. Uncertainties about electric vehicles (EVs) are addressed with the aid of robust adaptive optimization (ARO) in [23], and sensitivity analysis of the number of EVs and uncertainty level are considered. An ARO application for security-constrained unit commitment problems is discussed in [24]. This study demonstrates that using an ARO approach lowers the expected operation cost compared to the deterministic models.
1.2. Motivation
MG operation can be challenging due to uncertainty in factors such as load demand, RES generation, and weather conditions. This uncertainty can cause imbalances between energy supply and demand, leading to issues such as power outages or overloads. Parametric uncertainties in load, solar, and wind generation can arise due to various factors and can affect the performance of renewable energy systems. ARO is a mathematical optimization technique that can help address these uncertainties by designing control strategies that are robust to the worst-case values of uncertain parameters. This approach can ensure that the system operates effectively and efficiently even in the presence of uncertainties, by balancing power generation and consumption, utilizing energy storage systems, and smoothing out fluctuations in supply and demand.
ARO is an enhanced version of RO that utilizes a two-stage formulation. By solving the second subproblem, ARO identifies the most unfavorable conditions for uncertain parameters, where the objective function value is optimized under these conditions. In this paper, the ARO technique is used to solve the operation management of a MG, considering the uncertainty of RES power, consumer demand, and electricity price. The goal is to minimize operational costs under worst-case uncertainty, considering different uncertainty levels (ULs). A min-max problem is formulated to achieve this goal, and then, Karush–Kuhn–Tucker (KKT) conditions are utilized to transform it into a single Max problem. It is worth mentioning that the ARO problem is complicated, particularly when modeled as a nonlinear problem. Also, the electricity price uncertainty makes the ARO problem nonlinear, like the nonlinear ARO models introduced in [25, 26]. In this paper, linearization techniques are used to linearize the ARO problem and make it tractable to solve it in a very reasonable solution time.
Furthermore, the linearized model does not accurately represent MG constraints, including power losses and limitations on AC power flow. It is important to note that the nonconvex nature of the AC power flow equations can result in problematic and uncertain optimal solutions. Therefore, in this energy management problem, the MG constraints are precisely modeled using convexified AC power flow equations. Moreover, an optimization-based bound tightening (OBBT) technique is introduced to strengthen the relaxation, thereby gaining an accurate model for the MG scheduling problem.
1.3. Contributions
The previous sections discussed a literature review of various approaches to optimizing microgrids, including metaheuristic, deterministic, and stochastic methods. Despite a thorough review of existing literature, it was found that no similar paper considers all uncertain parameters, such as wind and solar generation, load demand, and electricity prices, while also considering convexified AC power flow equations with an OBBT technique. Accordingly, the key contributions of this paper are outlined as follows:(1)Development of an effective three-stage ARO algorithm to manage the day-ahead operations of a MG, with a focus on minimizing total operational costs while accounting for uncertainties related to load, PV and WT power, and electricity prices(2)Utilization of a linear formulation for the ARO problem in MGs that have the ability to exchange power with the upstream network(3)Precise modeling of AC power flow constraints through the adoption of convexified AC power flow equations and the introduction of new duality cuts based on the convexified subproblem(4)Introduction of a novel OBBT technique to enhance the relaxation of the convexified subproblem
1.4. Paper Organization
The remaining sections of this paper are structured as follows. Section 2 presents the optimization problem for MG operation, while Section 3 outlines the proposed solution algorithm. The simulation results, obtained by applying the proposed method to a modified IEEE 33-bus test study, are presented and analyzed in Section 4. Lastly, the conclusions of this study are briefly discussed in the Section 5.
2. Modelling and Formulation
This paper aims to introduce an ARO model that can effectively handle uncertainties in MG scheduling. The formulation of this model can be written as follows:where is the total operation cost which equals sum of the first and third stage costs, i.e., and . Sets , , and are variable sets of first-, second-, and third-stage problems; represents startup and shutdown constraints, indicates the uncertainty set, and ensures the feasibility of variables in set . includes the third-stage variables as well as uncertain parameters which are variables in the ARO subproblem.
A hierarchical three-stage decomposition algorithm is introduced as (1) that cannot directly be solved; a hierarchical three-stage decomposition algorithm is introduced. The first stage minimizes the start-up costs to find the best unit commitment status (modeling with binary variables), which are considered as constant parameters in other stages. The second stage is a maximization problem determining the worst state of uncertain parameters considered as fixed values in the third stage. Despite the previous works, this paper introduces a linear model for the ARO problem of a MG that can buy/sell power from/to the upstream network depending on the real-time electricity price. As mentioned, the uncertainties in the MG operation stem from variations in PV and WT powers, electricity prices, and demand. An uncertainty bound is defined for each uncertain parameter, and the ARO subproblem aims to resolve the worst-case situation. The third stage contains a nonlinear and convexified economic dispatch problem. The third level considers the convexified AC power flows to accurately model a MG. The solution strategy for solving the resulting optimization problem will be presented in Section 3. The proposed model is formulated as a three-stage min-max-min optimization problem as follows.
2.1. MG Scheduling Formulation
In the initial phase of the optimization problem, the determination of the unit commitment status is accomplished through the application of an objective function (OF) and constraints as follows:and subject to
Equation (2) shows the OF of the first stage, which minimizes the start-up and constant costs of the DGs. The constraints defined by are related to the start-up and shutdown constraints of the DG units [27].
The robust optimization approach takes into account the most adverse scenario of uncertainty. The conventional RO model assumes that uncertain parameters vary from the expected values over the scheduling horizon. This realization does rarely happen due to its low probability. The ARO approach considers different ULs to adjust the robustness of the solution. To apply ARO, the uncertainties of PV and WT powers, electricity demand, and price are considered as polyhedral uncertainty sets for each uncertain parameter, i.e., set presented in [27].
Following the determination of the unit commitment status in the initial stage and the identification of the worst-case scenario in the second stage, the third stage problem is introduced. The OF of the third stage is depicted in equation (4), and is utilized to specify the feasibility of other optimization variables in set . This OF and feasibility sets are described as follows using AC power flow equations:and subject to:where (5) describes the OF of the third stage. It is worth mentioning that OF of the second stage is the same as the third stage by this difference that uncertainty variables are not fixed in the second stage. The output power of DGs is restricted to their bounds with respect to their commitment status in (6)–(9) utilized to restrict the active and reactive powers of the grid to their respective limits. Battery charging and discharging constraints are described in (10) and (11). The energy storage limits for batteries are restricted by (12). The battery is charged/discharged with respect to charging/discharging efficiencies by (13). Constraint (14) is used to avoid depletion of the batteries over the scheduling horizon and forces the final state of charge (SOC) to be equal or greater than the initial SOC. Spillage power of PV and WT systems is formulated in (15) and (16). Load curtailments are limited to the related load values in (17). The constraints (18) and (19) express the equations for balancing active and reactive powers, while (20) and (21) specify the power flow equations. Line power flow limitation is described by (22). Voltage magnitude bounds are also presented in (23). Moreover, (24) represents the upper and lower bounds of voltage angle difference values, and by (25), voltage angle of the reference bus is set to 0.
2.2. Linearized and Convexified AC Power Flow Equations
In order to formulate the ARO subproblem in a linear form, an approximation of the AC power flow equations is required. As such, (20) and (21) are substituted with a linear approximation of the power flow equations in the second stage, as described by the following equation [27]:
In this representation, maximum line power flows are neglected and voltage magnitudes are considered as important variables for a MG. Contrary to DC power flow equations widely used in power system studies, this approximation allows us to impose a limitation on the voltage magnitude. However, linear estimations are not the exact equivalent of AC power flow equations and fail to calculate power losses in the MG. Full AC power flow constraints are highly nonconvex constraints, making the problem complex. Therefore, convexified equations of AC power flow are used in this paper. Convexification techniques include convex relaxations, which involve replacing nonconvex constraints with convex ones and convex envelopes, which involve replacing nonconvex constraints with convex ones that preserve the original function’s properties. Hence, an effective convexification method is utilized for trigonometric and bilinear terms of power flow equations. This is carried out by a relaxation technique relying on the convex envelopes around nonconvex functions, i.e., the following terms in (20) and (21):
For that purpose, the following constraints are used for the relaxation of AC power flow equations instead of (20) and (21):
For additional information about these convexification envelopes, see Appendix [28, 29].
3. Solution Algorithm
3.1. Three-Stage Decomposition Algorithm
As solving problem (1) directly is challenging, a three-stage decomposition algorithm is employed to achieve an optimal solution for the optimization problem. The process flow of the suggested algorithm is depicted in Figure 1.

The master problem represents the initial stage problem that identifies unit commitment and involves a mixed-integer linear programming (MILP) approach that incorporates both ARO and AC cuts. The primary objective of this stage is to determine the optimal status of DGs. The commitment results are used as fixed parameters in the second and third stages. In the second stage of the decomposition model, known as the ARO subproblem, a max-min problem is formulated based on linearized power flow constraints to obtain the worst-case realization of the uncertain parameters. It should be noted that the KKT conditions replace the minimum component of the problem, effectively transforming the second stage into a single maximization problem. Since the ARO and AC power flow constraints impact the MG scheduling, additional cuts need to be incorporated into the master problem. This is achieved by iteratively applying the ARO and AC cutting planes. The procedure is outlined as follows:subject towhere (43) and (44) specify the inequality and equality constraints for the third-stage problem based on linearized power flow equations. Constraints (39) and (40) correspond to the ARO and AC cutting planes, respectively, which are introduced to the master problem. Accordingly, auxiliary variable is applied in constraints (41) and (42) for adding new cuts. Stage one is solved with respect to both ARO and AC cuts, and after that, ARO and AC lower bounds are calculated as follows:
The unit commitment status is determined after solving the first-stage problem. The second stage is the ARO subproblem, a maximization problem formulated as [27]. The iteration between the first-stage problem and the ARO subproblem will be in process until the satisfaction of the convergence criterion. To do so, at each iteration, the ARO upper bound is calculated as follows:
After achieving the convergence criterion between the initial two stages, the third stage problem employs predetermined values for robust commitment statuses and the most unfavorable uncertainty scenario. This stage, called the AC subproblem, performs an economic dispatch with convexified AC power flow equations for finding optimal MG operating conditions. The AC subproblem formulation is introduced by the following constraints as well as set:subject to
By solving the third-stage problem, the AC upper bound is calculated as follows:
The convergence criterion for the first and third stages is assessed, and if it is met, the solving algorithm concludes and provides the optimal solution for the MG energy management problem. If the convergence criterion is not met, a supplementary duality cut for the AC subproblem is created and sent to the master problem.
3.2. Bound Tightening Relaxation Method
The relaxation method’s characteristics and effectiveness are closely linked to how tightly the voltage magnitudes and angle differences are constrained. In this sense, an OBBT method is employed to firstly reduce the domains of the variables and secondly strengthen the relaxation. This tightening technique is an iterative algorithm in which the best variable bounds are determined by means of minimization/maximization of lower/upper bounds of variables considering the optimal power flow (OPF) problem in the third stage. The proposed optimization-based bound tightening relaxation algorithm is illustrated in Figure 2. This procedure is carried out as follows: Step 1. Solve the nonlinear problem: First, the nonlinear OPF problem is solved by nonlinear programming (NLP). Since the solution to this nonlinear problem may not be the global optimum, it is treated as an upper bound for the objective function of the third stage in the iterations described in (51). Step 2. Solve max(min) problem: assume is an optimization variable whose tight bounds should be determined. This variable is limited between its lower bound and upper bound . In order to find the best tight bounds, the following optimization problem is solved. In the abovementioned minimization/maximization problem, the variable takes its lower/upper value; and this value gives way to the lower/upper bound of the variable in the next iteration. Constraint (54) defines the upper bound of OF of the third stage , equal to the nonlinear solution in the first iteration. In the next iterations, takes OF value in the previous iterations. Step 3. Check stop criterion: After solving the tight-bound optimization problems, the lower bound and upper bound of variables are determined. The stopping criterion is checked by evaluating the value OF and its previous value. There are two possible conditions:(1)Provided that OF of the third stage is lower than , the parameter is updated by the new , calculated in the second step. Then, we go back to the second step and continue the OBBT algorithm.(2)Provided that OF of the third stage is equal to , the OBBT algorithm is terminated, and the achieved lower/upper bound is used for the solution algorithm.

4. Numerical Study
4.1. Case Study System
With the objective of exhibiting the efficacy of the proposed model and solution algorithm, a 33-bus test system is evaluated. This benchmark system is identical to the ref [27] which includes 4 dispatchable units, 3 PV stations, 2 WTs, and 1 battery energy storage system. This system was chosen in order to compare the results with the reference paper. Predicted load, PV, and WT generation and electricity price for the scheduled timeframe are illustrated in Figure 3 [27, 30]. Five different ULs are considered, including 0, 25, 50, 75, and 100 percent. Once UL equals zero, all the uncertain parameters are considered as forecasted values (deterministic case), and UL 100% means a conventional robust model. All simulations are done on 64-bit Windows with a core i7 processor and 8 GB RAM in GAMS software. Cplex, Gurobi, and Conopt solvers solve the linear, convex, and nonlinear models.

4.2. Results
Two different cases are addressed as follows: Case I: In the two-stage adaptive robust optimization problem that accounts for linear power flow equations, the MG energy management problem is presented as a linear ARO problem that neglects AC power flow equations. The master problem determines the optimal unit commitment, while the ARO subproblem identifies the worst-case scenario of uncertain parameters. The AC subproblem is not incorporated in the optimization problem in this scenario. Case II: In the three-stage adaptive robust optimization problem, the convexified power flow equations are considered, including the AC power flow constraints, and all three stages of the optimization problem are solved. Dual AC subproblem cuts are transmitted to the master problem using (55). To achieve the optimal result of the proposed optimization model, the AC subproblem is addressed concurrently with the ARO and sends dual cuts to the master problem to modify the unit commitment status and, consequently, the worst-case realization.
Table 1 reports the solution and run time of case I and case II. The ARO subproblem has a linear model, while the third stage considers the AC power flow equations using convexification envelopes. The proposed method finds the optimal solution in a reasonable computational time. As shown in Table 1, the more UL we consider, the more total cost we obtain. That is why the robust model considers the worst case of uncertainty. Operation costs are more in the robust approach compared to the deterministic approach because higher ULs cause broader ranges of uncertainty, and the most extreme uncertainty scenarios result in higher costs. Since generated duality cut by the third stage updates the unit commitment status and worst realization of uncertainty parameters, the third stage of the proposed method influences the optimal solution and eventually results in lower operation costs.
Figure 4 illustrates the commitment results of DGs with UL = 25% in Case I and Case II. As can be seen, DG1 and DG3 are committed most of the time over the time horizon because they are more economical compared to DG2 and DG4. Case II consists of convexified AC power flow equations and describes the MG constraints well, while case I uses a linear approximation of power flow constraints and does not model power losses and maximum line power flow constraints. As it is seen, the third stage affects the unit commitment results. The reason is that the dual AC cuts update unit commitment status in the master problem and result in lower operation costs.

Figure 5 shows the realization of electricity price at the connection bus between the MG and the power grid with UL = 25%. Moreover, Figure 6 exhibits the exchanged power between the MG and the power grid. Negative exchanged power shows that there is an excess of power in the MG that is sold to the main grid. In contrast, positive exchanged power refers to the power shortage in the MG, thereby purchasing power from the grid. When the MG sells power to the grid, the uncertain price parameter reaches its lower bound, while it ends up with an upper bound when the power is purchased from the upstream network. In Case I, the worst uncertainty realization appears at hours 5, 6, 19, 20, 23, and 24, reaching its upper bound. This worst uncertainty realization in Case II happens at hours 12, 13, and 14, hitting its lower bound and hitting its upper bound at hours 6, 20, and 23. The worst realization of the electricity price depends on two factors, the electricity price and the exchanged power between the MG and the upstream grid. The ARO aims to find the worst realization of uncertain parameters. Accordingly, it is expected that the worst realization of the electricity price occurs at hours in which power transaction is greater and also the electricity price is high. In other words, the ARO chooses the hours in which the sensitivity of OF to price is maximum.


Figure 7 illustrates the battery SOC in both cases. Charging/discharging the battery is one of the decisions that is made in the third stage. The battery is discharged when the electricity price is high, and the MG sells power to the grid, whereas the battery is charged when the electricity price is low. As it is seen in Figure 7, the initial and final SOC are equal.

According to Table 1, as UL increases, so does the associated operational cost. In order to demonstrate the privilege of the ARO approach, the expected cost through numerous scenarios is studied. Therefore, a Monte-Carlo simulation (MCS) is carried out, and many scenarios are generated. These realizations of uncertain parameters are used to evaluate the expected cost for each UL shown in Figure 8. According to this figure, all robust solutions have lower expected costs than the deterministic model (i.e., UL = 0). The best solution for Case I occurs with UL = 25%, but according to the proposed three-stage adaptive robust model, the lowest anticipated cost can be attained with a UL setting of 50%. Since uncertainty is not considered in the deterministic model, once the demand load exceeds the anticipated value or renewable generation power is less than the predicted value, we might lose some loads, which cause a heavy penalty for operation costs. In contrast, sufficient reserved powers will be considered in the ARO approach in case the uncertain parameters violate their predicted values. Therefore, the ARO model results in less load curtailment and lower expected cost. The superior solution of the proposed ARO model was attained for UL = 50%. When the UL equals 100%, we have the conventional RO model that is too conservative. This means that the ARO model has an optimal UL, and there is no need to increase the UL by 100%, since the expected cost does not improve for higher ULs.

Table 2 provides a comparison of the expected cost for the deterministic solution, conventional RO, and ARO of both cases. The AC-constrained model results in lower expected costs compared to the linearized model, as the latter neglects power losses and AC constraints. The proposed convexification method accurately models AC power flow equations and results in better cost performance than the linearized model. Among the deterministic, RO, and ARO models, the deterministic model performs the worst as it does not account for uncertainties. The conventional RO approach is overly pessimistic about uncertain parameters, though it has a lower expected cost than the deterministic model. However, it does not outweigh the ARO’s best solution. Therefore, the ARO model results in the best expected cost performance among all models due to its consideration of different uncertainty levels and optimization of solution robustness.
5. Conclusion
A three-stage adaptive robust optimization model for microgrids operation, considering the uncertainties of PV and WT generation, consumer demand, and price of electric power, was presented in this paper. To solve the model, a decomposition algorithm was proposed, which includes a master problem as well as two subproblems. The first subproblem solves a linear adaptive robust model to specify the worst case of uncertain parameters, whereas the second subproblem considers convexified AC power flow equations. It was verified that although higher ULs result in higher operation costs, they lead to a more robust solution for operation management problems. This robust solution provides support power in the presence of uncertainties and results in a lower expected cost. A MCS analysis was also performed to evaluate the superiority of the proposed model and the proposed solution algorithm. It was shown that the proposed ARO model outweighs deterministic, conventional RO, and linearized models in terms of expected cost since deterministic models neglect uncertainties, RO is too pessimistic, and the linearized model neglects power losses and AC power flow constraints.
The studied MG includes dispatchable DGs, WTs, PV stations, and batteries. This paper did not discuss the presence of EVs in the MG optimization problem. The growing presence of EVs in the power system has increased uncertainties in MG operations. This poses a challenge for managing and optimizing MG performance, as the behavior of EVs can be difficult to predict and control. Hence, for the future scope of this research study, it is highly recommended to consider EVs and their inherent uncertainties with the ARO approach.
Appendix
The convex envelopes for the square and bilinear functions (e.g., for variables and ), as well as trigonometric terms, can be expressed as follows [31]:where
Nomenclature
Indices: | Battery |
: | Demand |
: | Dispatchable generator (DG) |
: | Bus |
: | Iteration |
: | Line |
: | Network |
: | Photovoltaic (PV) system |
: | Time |
: | Wind turbine (WT) system |
: | Susceptance/conductance of line; |
: | Cost of DG, PV, and WT system |
: | Resistance and inductance of line |
: | Maximum allowable power of the line |
: | The end of scheduling time |
: | Start-up/constant cost of DG unit |
: | Upper bound of objective function in bound tightening relaxation problem |
: | Load curtailment cost |
: | Charge/discharge efficiencies of battery |
: | Stored energy of battery |
: | Active and reactive curtailed load |
: | Objective function of main problem |
: | The objective function of each stage |
: | Charge and discharge rates of battery |
: | Active and reactive power of line |
: | Active and reactive power of network |
: | Power of DG, PV, and WT system |
: | Spillage of PV/WT system |
: | Start-up cost of DG unit |
: | Voltage magnitude |
: | Commitment status of DG |
: | Dual variable |
: | Voltage phase angle |
: | Active and reactive power price |
: | Auxiliary variables for cuts sent to master |
: | First-stage decision variables |
: | Third-stage decision variables |
: | Uncertain parameters |
: | Constraints of each stage. |
Data Availability
The numerical data used to support the findings of this study are included within the article.
Additional Points
In whole paper, minimum and maximum of a variable are denoted by under bar and over bar , respectively; fixed values are shown by superscript and is used for convex envelopes.
Conflicts of Interest
The authors declare that there are no conflicts of interest regarding the publication of this paper.