Hybrid stochastic/robust flexible and reliable scheduling of secure networked microgrids with electric springs and electric vehicles

Electric spring (ES) as a novel concept in power electronics has been developed for the purpose of dealing with demand-side management. In this paper, to conquer the challenges imposed by intermittent nature of renewable energy sources (RESs) and other uncertainties for constructing a secure modern microgrid (MG), the hybrid distributed operation of ESs and electric vehicles (EVs) parking lot is suggested. The proposed approach is implemented in the context of a hybrid stochastic/robust optimization (HSRO) problem, where the stochastic programming based on unscented transformation (UT) method models the uncertainties associated with load, energy price, RESs, and availability of MG equipment. Also, the bounded uncertainty-based robust optimization (BURO) is employed to model the uncertain parameters of EVs parking lot to achieve the robust potentials of EVs in improving MG indices. In the subsequent stage, the proposed non-linear problem model is converted to linear approximated counterpart to obtain an optimal solution with low calculation time and error. Finally, the proposed power management strategy is analyzed on 32-bus test MG to investigate the hybrid cooperation of ESs and EVs parking lot capabilities in different cases. The numerical results corroborate the efficiency and feasibility of the proposed solution in modifying MG indices.


a. Aims and Motivation
In recent years, the development of renewable energy resources (RESs) and electric vehicles (EVs) into the power system especially traditional distribution networks to be more environmentally friendly has created serious challenges [1,2]. Due to uncertainties associated with RESs generation and the heterogeneity of EV users, there is power imbalance, voltage and frequency fluctuations and feeder overloading in this network [3]. Further, due to the existence of randomness in the market energy price and load demand, implementing a proper management in the smart active distribution network (SADN) becomes difficult. These high-impact-frequent oscillations originated from the mismanagement of sources and active loads (ALs) cannot be managed by traditional operational techniques, and infrastructural reinforcement measures [4,5]. The application of electric spring (ES) as a novel smart load that can be placed serially between non-critical and critical loads has recently been suggested for microgrid (MG) vulnerabilities ( Fig. 1) [6]. The ESs like EVs can appear in the role of smart ALs or flexibility resources (FRs) due to exchanging active and reactive powers with the MG. Now there is a great amount of concerns about the mismanagement of ES and other FRs such as EVs due to its impact on critical loads and operating conditions of MG [7]. To cope with the growing share of uncertainties and managing the complex interactions of FRs, this paper proposes a coordinated power management strategy (PMS), aiming to jointly optimizes the technical, i.e. operation, security, reliability and flexibility, and economic indices of the SADN.

b. Literature Review
Upgrading the structure and capabilities of the ES has been the subject of various studies exploring different generations of ES [7][8][9][10][11]. To provide active power control capability to the capacitor-based ES (ES-1), in the second generation (ES-2), by replacing battery with the capacitor, the high cost of the battery and the limited operating range have been imposed on the ES. In the next generation, to eliminate the limitations of active power control and low operating range of the two prototypes, the use of back-to-back converter (B2B) with isolation transformer is proposed for permanent supply of the ES from the grid (ES-B2B) [8]; While the investment cost of ES has increased significantly, the power quality problems caused by power converters have increased. In another photovoltaic based ES structure (ES-PV) [9], the utilization of PV energy was proposed to supply the required energy of the ES and reduce the energy received from the MG, in which due to the uncertainty in PV production there was still a need for high storage capacity [10]. In [11], to remove the DC/DC converter of ES-B2B and ES-PV, application of a bidirectional AC/DC converter with a battery was proposed. In another approach, instead of focusing on improving the structure of the ES to reduce the required capacity of the battery, the coordinated power control of the battery with non-critical load (NCL) and grid is proposed [12]. Generally speaking, prior studies have focused on the configuration and dynamic modeling of ES [13][14][15]. It should be noted that modeling of power system elements is done in static, quasi-static and dynamic forms. These types of categorizations are based on the time interval of studies. The shortest time period is the dynamical model, followed by quasi-static and static models. According to the literature, generally, the articles presented the dynamic model of ES in which time-varying ES and network variables were analyzed. However, given some of the important issues in the power system, such as optimal power flow (OPF) in the power system planning, operation, etc., due to the longer period of study, we need to develop the static model of ES. Additionally, in order to investigate the capability of ES on the technical and economic indices of a MG faced with various uncertain resources and loaded by EVs, it is necessary to define an OPF problem. This requires developing the static model of ES due to the operation of the network in a long study period (e.g. 24 h).
In the field of static model of ES as shown in Table 1, for the first time in [11], we presented a non-linear programming (NLP) formulation of the static model in the context of scenario-based stochastic programming (SBSP). To model the uncertain parameters the roulette wheel mechanism (RWM) for generating scenarios as well as the simultaneous backward method to decrease generated scenario samples have been employed. One of the major problems of SBSP is the need to generate a large number of scenarios to obtain a guaranteed response and accurate estimation of uncertain parameters' probability density function (PDF). The high number of generated scenarios, however, results in increasing computational complexity and reducing convergence speed. Therefore, to reduce scenario samples many researches have used scenario reduction techniques, which does not remarkably eliminate the mentioned shortcomings [16]. Therefore, to solve the mentioned issues of stochastic modeling, in this paper, the problem will be modeled as a hybrid stochastic-robust optimization (HSRO). In other words, for modeling uncertainty in the stochastic part, unscented transformation (UT) method is used which has a much lower number of generated scenarios than other techniques, without any need to lessen generated scenario samples [17]. Therefore, uncertain parameters such as load demand, energy price, RESs output and MG equipment availability are modeled by the UT method. In addition, for modeling EVs uncertain parameters consist of charging and discharging power rate, arrival and departure energy and charger capacity of EVs in the parking lot, bounded uncertainty-based robust optimization (BURO) will be used to investigate the EVs' robust potentials in improving MG indices technically and economically [18,19]. In more researches [20][21][22], EVs capabilities have been used to facilitate MG operating conditions, i.e. voltage and frequency stabilization, energy loss, lines overloading, etc. In this paper, however, EVs parking lot is coordinately utilized with ES to modify the network flexibility, reliability, security, operation and economic indices.
The mentioned purposes can be achieved if the MG uses a suitable energy management system (EMS) or PMS [23]. In the EMS method, the energy or active power control of different sources and ALs are done to improve network operation [24], flexibility [25], reliability [26], and economic [27]. However, in order to investigate network security index, it is necessary to control the reactive power [28], which is done in the PMS simultaneously with the active power control.
According to Table 1, no work has been performed to analyze the impacts of distributed ESs and EVs coordination on the network indices. Owing to the dependency of network indices on each other, simultaneous analyzing of them is of utmost importance. Hence, this paper implements PMS in the MG based on simultaneous modeling of the flexibility, reliability, operation, security and economic indices. As shown in Fig. 1, in the proposed PMS, operation of distributed renewable and flexible resources including ESs and EVs parking lots are centrally coordinated by MG operator (MGO), so that by implementing optimal scheduling of MG, maximum revenue is obtained.
Finally, the non-convex NLP format of the original PMS model which is normally solved by time-consuming numerical methods [29][30][31] or evolutionary algorithms [32,33] and leads to the locally optimal solution, will be converted to convex linearized format to reach a faster globally optimal solution. This process is done through conventional linearization approaches, as a type of convex relaxation methods [34][35][36][37].

c. Research Gaps and Contributions
In summary, according to the above-mentioned researches and Table 1, the main research gaps and the foremost contributions made in this paper are: -According to the literature, relatively extensive studies have been conducted in the field of control and structure of ES compared to its modeling. Although limited studies have been done on the dynamic modeling of ES [38], to the best of the authors' knowledge, the integrated static mathematical modeling of ES and EVs parking lot has not been presented yet, which is necessary for the technical and economic evaluation of the proposed plan capabilities. -In the field of uncertainty modeling, the stochastic modeling of uncertain parameters including energy price, load demand and RESs output have been expressed in [11]. While in this paper, in addition to the UT-based stochastic model of the aforementioned parameters and the MG equipment availability, EVs parking lot uncertain parameters including charging and discharging power of EVs fleet, arrival and departure EVs' energy and EVs' charger capacity are modeled using BURO method in the proposed PMS. Uncertainties as a vital issue in the flexible operation of MG have not been studied in such a way so far. -In most studies, the auxiliary services capabilities of ES including voltage and frequency compensation have been widely considered [39,40], whilst in this paper, the impact of cooperative performance of ESs and EVs parking lots on MG economic and various technical indices are simultaneously quantified. It can be considered as a novel solution in demand-side management to tackle uncertainties and contingencies. -For the first time, the linear-hybrid-coordinated mathematical framework of distributed ESs and EVs parking lot in the renewable MG as FRs and ALs are presented. The proposed framework speeds up finding the best feasible solution and reduces response error, by transforming the MINLP form of the optimization problem into the MILP form using the linear approximation-based convex relaxation methods. d. Paper Organization The remaining parts of the paper are outlined as follows: Section 2, presents the original and linearized modeling of the proposed optimal scheduling. Section 3 expresses the HSRO framework of the proposed solution in 2 separate parts of stochastic and BURO. In Section 4, the numerical results including dataset, system setup, and the study results are presented. Finally, Section 5 concludes the paper.

Original model
A) Objective function: The objective function of the proposed problem is stated in (1a) which consists of 4 parts. The first part refers to the expected energy cost received from the upstream network. Also, reliability cost (lost load) is mentioned in the second part of the objective function, which is proportional to the value of lost load (VOLL) and expected energy not supplied (EENS). In addition, the third and fourth parts are dealing with flexibility and security benefits, respectively. Note that in this paper, flexibility as an operational option means the ability of MG in response to small/large fluctuations, while following technical and economic requirements set by the MGO. Also, flexibility is considered from an economic perspective, which is the benefit of flexibility. This profit is considered as the product of flexibility incentive price (FIP) and flexibility energy (FE), and the FE depends on the flexibility capacity. Therefore, due to the purpose of the study, FE is considered as the upward and downward potential capacities of FRs. The security benefit is also calculated by multiplying the security incentive price (SIP) and security energy (SE). It is noteworthy that FE or SE is the amount of energy required to provide the demanded flexibility or security for a MG. Also, FIP and SIP are incentive mechanisms for encouraging FRs to participate in flexibility procurement.
(1a) B) Microgrid constraints: These constraints are including power flow equations, (2a)-(6a), and microgrid technical constraints, (7a)-(10a) [42]. The power flow equations also refer to the active and reactive power balance in each bus, (2a) and (3a), lines' active and reactive power, (4a) and (5a), and the amount of slack bus angle, (6a). In these equations, it is assumed that the MG at the slack bus (ref.) is connected to the upstream distribution station. Hence, the DS P and DS Q values in other buses are equal to zero. It should also be mentioned that ESs and EVs batteries are modeled as ALs in (2a) and (3a), where the positive/ negative value of their active power variables corresponding to the charge/discharge mode, and for their reactive power variables means inductive/capacitive mode.
Not that load models can be categorized into three types: static, dynamic, and composite models. The static characteristics of the load can be classified into constant impedance (Z), constant current (I) and constant power (P) load, depending on the power relation to the voltage. The static ZIP load model as a well-known model in the power system, represents the relationship between the active and reactive power as a function of the applied voltage. Therefore, the ZIP load model can be a suitable counterpart for the NCLs model, wherein it is assumed that the consumption of NCL serially connected to the ES output, can be managed by the ES voltage controlling as per equations (2a) and (3a). On the other hand, because the ES uses the IGBT bridge converter, the NCL connected to the bridge can be modeled as constant impedance/ current/power or a combination of these models [43][44][45]. Therefore, the ZIP model is used for the NCL as a comprehensive model that covers all types of NCL models [11,38]. Hence, the term B c is equal to one for buses with ES, owing to the ZIP model, otherwise, it is zero. Also, the technical limitations of the MG are expressed in (7a)-(10a), which indicate respectively the allowed range of lines capacity, distribution station capacity, buses voltage magnitude and the voltage limit of buses containing critical loads (CLs).
D) Electric vehicle parking constraints: Equations (17a)-(23a) state the mathematical model of the integration of EVs into the MG [46]. Constraints (17a)-(20a) represent respectively the power balance between the EVs' battery and its output, the EVs' charger losses, the energy stored in the EVs' battery and the EV's battery charge/discharge limit. Also, constraints (21a)-(23a) express EVs' battery energy at arrival and departure time, and EVs' charger capacity. It should be mentioned that the total charge/discharge rate of EVs parking lot at each moment is equal to the sum of the charge/discharge rate of each EV existed in the parking lot. Also, the charger capacity of the parking lot at time h is equal to the total capacity of the available EVs in the parking lot. In addition, EV E Arr,h is equal to

E) MG flexibility constraints:
This section is dedicated to introduce an index for calculating the flexibility amount of FRs. This index quantifies the flexibility of the FRs and network individually against uncertainty resources [47]. The flexibility index for FRs including ESs and EVs is characterized as upward and downward flexibility. If the output power of the FRs in scenario s is greater than the base scenario, there is upward flexibility, otherwise, it is downward flexibility. Therefore, the upward and downward power flexibilities for ES and NCL arising from ES battery and NCL control are calculated by (24a) and (25a), and for EVs, it is obtained based on (26a). Moreover, the flexibility energy (FE) of the MG is calculated by (27a). It should be stated that scenario "1 ′′ is pertained to the base case that uses forecasted (normal) value for uncertain parameters.
F) MG reliability constraints: In this paper, the EENS index for reliability assessment is used as presented in (28a). Also, the limitation of the unsupplied active load is expressed as (29a) [48].
G) Security constraints: Voltage security, as a vital index in optimal scheduling of MG, is defined by (30a)-(31a). It is determined by voltage stability margin (VSM) through different methods [42]. In this paper, the VSM is calculating via the voltage stability index, i.e., the SI-index method. In this method, the weakest voltage bus or worst SI (WSI) is firstly specified and then the WSI is determined by (30a). In constraint (30a), V wb-1 is the voltage at the sending end bus, L P wb−1,wb and L Q wb−1,wb are respectively active and reactive power in line between buses wb-1 and wb. The SI-index value changes between 0 (voltage collapse) and 1 per unit (p.u.) (no-load condition) in the hourly operation of MG. Also, it varies from 0 to 24 in daily operation. The voltage security is limited by (31a), and the required energy to improve MG security is formulated as (32a) in which RESs, ESs (their battery and NCL control part) and EVs form the security energy.

Linear model
In the original model, equations (2a)-(5a), (7a), (8a), (12a), (16a), (18a), (23a), (25a), (30a) and (32a) are non-linear and equations (4a) and (5a) are non-convex [49]. Therefore, the problem is a non-convex NLP that in the best condition can reach the locally optimal solution. It is generally calculated by numerical methods such as Newton Raphson or evolutionary algorithms, so it will be difficult and take considerable time to solve the problem. Besides, in large-scale and complex problems, achieving a feasible solution will be very time-consuming and hard. Therefore, the main problem is transformed into a linear counterpart by conventional methods. The details of this process are as follows: A) Linearization of power flow equations: Based on the linearization theorems expressed in the piecewise linearization method, and assuming that the voltage difference between two ending buses of a distribution line is always less than 6 degrees or 0.105 rad, the linear format of power flow equations, (2a)-(5a), and the voltage limit of the buses, (9a)-(10a), will be as (1b)-(6b). More details of the linearization process can be found in [49,50].
Where, ΔV U and ΔV CB U are equal to (V U − V D )/n u and / n u respectively, where n u is the total number of linear pieces. Equations (25a) and (32a) will also be transformed into linear equations, (7b)-(8b), based on the assumptions in this section.
B) Linearization of circular equations: Equations (7a), (8a), (16a) and (23a) are circular constraints in P and Q coordinates with the center of zero and the maximum radius of S. Based on [19], a circular plane can be approximated by a polygon plane with n k edge, which each edge corresponds to a line equation. This equation for the k-th edge of polygon is equal to the tangent line equation of edge k at k × 2π/n k angle. Therefore, the linear format of circular equations will be as follows:

C) Linearization of the absolute terms:
For the linear format of equations (12a) and (18a), which have the absolute term of active and reactive power of ESs and EVs, each variable is divided into two positive and negative components. The positive component of active and reactive power is introduced as discharge and capacitive powers, and the negative component of active and reactive power is considered as charge and inductive powers. It is noteworthy that these variables have positive value, so the linear model of equations (12a) and (18a) will be written as follows:

D) Linearization of voltage security equation:
In Equation)30a(, due to the small amount of MG lines' reactance against the resistance, the third term versus the second term has a small amount and can be ignored. Also, in the second part, V 2 is equal to (V L ) 2 + ∑ u∈L sl u ΔV u , but note that since the product of ΔV by the other variables leads to a very small value, so in the second part of (30a) only the term (V L ) 2 can be replaced by V 2 . Finally, according to the linear method, the WSI equation can be approximated as follows: Where, ml represents the slope of the line. Finally, the linear model of the proposed problem is defined as follows:

HSRO framework
In this section, a hybrid stochastic/robust model is intended for modeling uncertain parameters. In this model, parameters such as energy price, ρ DS , active and reactive power of load, D P and D Q , renewable energy sources power, RES, and equipment availability based on forced outage rates (FOR) are modeled in the stochastic form. This approach is due to the fact that to obtain a more accurate EENS index many scenarios need to be considered, and the EENS is dependent on the amount of these parameters. In addition, the EVs' uncertain parameters including the charging and discharging rate of EV, EV P B.U and EV P B.L , maximum EV charger capacity, EV S U , arrival and departure energy of EV, EV E Arr and EV E Dep , are analyzed in the robust framework. Because another purpose of this paper is to consider the robust ability of EVs in enhancing the flexibility, security, and reliability of the MG in the worst-case scenario of uncertainty.
A) Stochastic framework: This section is aimed to present the UT method as an efficient tool for uncertainty modeling in the stochastic framework. It should be stated that the UT method requires much fewer runs than other methods such as Monte Carlo simulation (MCS) and analytical methods for convergence; Therefore, the computational burden and calculation time are significantly reduced. Also, unlike analytical methods, the proposed approach does not require any mathematical assumptions for model simplification. Additionally, the UT method as one of the approximation-based methods which are suitable for modeling uncertainties due to its ability in non-linear correlated/ uncorrelated transitions [51,52], acceptable PDF estimation and simple coding is employed in this paper. In this case, the mentioned stochastic uncertain parameters are modeled by the UT technique in which n, as the dimension of input uncertain parameters vector (U), is equal to 5. In this method, the total number of scenarios created is 2n + 1, which it is equal to 11 scenarios for the proposed problem.
Therefore, due to the small number of scenarios produced in this method, there is no need to use scenario reduction methods to reduce the computational burden and calculation time. Further details on this method, including the formulation and implementation algorithm, are presented in detail in [51]. B) Uncertainty bounded robust optimization framework: In this robust method, unlike the SBSP method, there is no scenario sampling process, but instead, the space of uncertainty is limited to an uncertainty set. It reaches the best solution that is feasible for all realizations of uncertainties that lie in the uncertainty space under consideration. In other words, in this method, the uncertain parameters are modeled as the boundary range, and then according to the problem and position of the uncertain parameter in the original model, the amount of that parameter is obtained in the worst-case scenario. Finally, the main optimization process is performed for the worst-case scenario. It should be noted that in order to integrate the problem, determine the uncertain parameters and variables of the problem simultaneously, BURO constraints are added to the basic problem related to the worst-case scenario.
All uncertain EV's parameters in this method are in the range of (1 ± σ) × X in which σ⩾0 is considered as the uncertainty level of forecasted parameter (X). It should be noted that the problem has a robust solution if the following conditions are met [18,19]: Note that according to the stated robust theory, equations (21a)-(22a), (4c), (7d)-(10d) are respectively converting to equations (1g)-(9g).
It is noteworthy that according to Section 3-A, in equations (1g)-(8g) and (1h)-(2h) for the parameters corresponding to the forecasted scenario (normal), the subscript "1 ′′ is used for the EV's parameters in these equations. In addition, the proposed robust model is applied to the unequal constraints in the standard form, hence, the inequality (≥) is replaced with equality in (4g) and (5g), because, in a problem with minimizing the objective function, the equality constraint can be equivalent to inequality. Therefore, the final format of the proposed model can be expressed as follows: Subject to: Finally, the flowchart of the proposed model is expressed as Fig. 2.

Dataset and system setup
The power management strategy based on hybrid stochastic/robust optimization is implemented on 32-bus radial MG with the base values of 1 MW and 12.66 kV, as shown in Fig. 3 [53]. The specifications of lines and peak load hour (20:00) are stated in [11]. Also, the amount of load at other hours is equal to the peak load multiplied by the daily load factor curve, whose predicted load values at different hours are based on the load data of Rafsanjan city in Iran, which is depicted in Fig. 4(a). The MG has two types of RESs, solar and wind systems. The capacity of the solar system is 2 MW and the capacity of wind systems at buses 10 and 14 are 1.5 MW and 1.8 MW, respectively. Also, the daily energy price curve is taken from [11], wherein the energy price is 16 $/MWh during 1:00-8:00, 24 $/MWh at the periods of 9:00-17:00 and 23:00-24:00, and it is presumed 30 $/MWh for the hours from 18:00 to 22:00. The ESs' location and other specifications are assumed to be based on the performed optimal planning in [11], where the ESs are placed at buses 1, 2, 3, 10, 11, 12, 13, 14, 15, 18, 19 and 22, with the sizing in the range of 0-10 MWh. It is worth noting that the CLs are in all buses where the ESs exist, and critical buses consist of 6, 7, 8, 9, 20, 21, 25, 26, 27 and 28. The acceptable voltage limit in the MG on the basis of international Standards ANSI C84.1-2006 and IEEE std 1250-1995, for non-critical buses is specified from 0.9 to 1.05 p.u., while that of critical buses is intended in the range of 0.97-1 p.u [54]. It is also assumed that all buses in the MG have EVs parking lot, for which the number of EVs per bus is based on [19]. Other specifications of the EVs such as SOC, battery capacity, charging rate, charging capacity, etc. are mentioned in [32]. Also, the loss coefficients for EVs charger, λ and γ, are similar to ESs loss respectively. Finally, it should be noted that the value of VOLL, SIP and FIP are set at 100 $/MWh, 10 $/MWh and 10 $/MWh, respectively. Also, the FOR is considered to be 1% for each MG equipment such as distribution lines and distribution station.

Results and discussions
The proposed problem is simulated in the GAMS software environment [55]. It should be stated that the number of linearization segments of the piecewise method and circular constraints is 5 and 45, respectively [50].
A) Evaluation of different solvers: Table 2 enumerates the convergence results of the NLP and MILP models for the proposed problem in case of HSRO(σ = 0, δ = 0). For the non-linear model in this section, four solvers comprising CONOPT, IPOPT, SNOPT and MINOS are selected in the GAMS software. It can be seen from Table 2 that the NLP model solvers have different convergence iteration, computational time, objective function and model status, despite having the same number of equations and variables. Therefore, the NLP model solvers do not have a unique answer, it is also seen that the best solver for the nonlinear model of the proposed problem is IPOPT. Because it has been able to achieve the lowest objective function in the shortest possible time compared to other NLP solvers in the setting of locally optimal model status. Also, the selected solvers for the proposed MILP model in this section are CPLEX, CBC, and CONOPT. As it can be seen from Table 2 all the solvers have a unique response and globally optimal model status. However, the calculation time by the CPLEX solver is smaller than that

Table 2
The results of the different proposed formulations for the model convergence in HSRO(σ = 0, δ = 0). LO: Locally optimal, GO: Globally Optimal, I: Infeasible of CBC and CONOPT. Therefore, it can be stated that the MILP model has a more reliable response rather than the non-linear model and the best solver is the CPLEX which has the lowest calculation time. Table 3 shows the amounts of operational, flexibility, reliability, security and economic indices in the MG based on MILP model in HSRO (σ = 0, δ = 0) for the following case studies: The given table reveals that RESs in the MG (Case II) raises the MG energy losses compared to its absence (Case I), and it also reduces the voltage deviation in the MG, however, the overvoltage of the MG increases compared to case I due to the high active power injection of the RESs. The presence of EVs parking in Case III reduces the energy loss and maximum voltage deviation/overvoltage in the MG compared to Case II. The reason for this is that EVs due to energy transition and storage capabilities have dispatched RESs output power in a more effective way relative to Case II. Thus, the injected power into the MG from the upstream distribution station has been reduced. This case demonstrates the advantage of managing the EVs and local resources operation in the MG. Note that due to EVs demanded energy from the MG in Case III, the energy loss, in this case, is higher than that of Case I. €Finally, the presence of ESs in Case IV has able to obtain a better operating condition to the MG than any other cases, there is the least amount of energy loss and maximum voltage deviation in this case because of effective power management among the local components of the MG. Although the overvoltage in Case IV is slightly higher than Case I, it is notably reduced in comparison with the other two cases.

Model
As shown in Table 3, the highest amount of EENS occurs in Case I, which does not have any local resources such as RES, EVs and ESs, but with the presence of more local resources in the MG, the EENS rate is significantly reduced so that Case IV has the lowest amount of EENS (3.78 MWh) than any other cases. This analysis is also valid for the MG voltage security and flexibility indices, the increase of local and flexibility resources in the MG could increase the WSI index and provide notable SE and FE values for the MG. Hence, in Case IV, the maximum possible values of WSI, SE and FE are calculated due to the superior performance of ES in the MG. Also, the highest voltage deviation in Case I is related to bus 18, hence this bus is known as poor voltage bus (PVB) in Case I. But in other cases, considering the location of RESs, EVs and ESs, bus 33 is known as PVB.
Finally, economic indices such as the price of energy received from the upstream network (MG operating cost), reliability cost, security and flexibility benefits and total cost of the proposed scheme have been expressed in rows 9-13 of Table 3, respectively. As is presented, Case I due to the lack of local resources in the MG only contains operational and reliability costs. With the addition of RESs to the MG (Case II), in addition to reliability and operational costs that have been reduced compared to Case I, due to the injection of RESs into the MG, the proposed scheme has reached revenue from providing voltage security. Therefore, the total cost of Case II has reduced by approximately 92.86% with respect to Case I. In Case III, the operating cost, due to the increased demand of power from the MG by EVs, increases relative to Case II, but as EVs can reduce EENS, the cost of reliability in Case III is lower than Cases II and I. In contrast, in Case III because of the presence of RESs as a source of voltage security, and EVs as a source of flexibility and voltage security, the benefits of deploying security and flexibility in the MG decreases the total cost in this case with regard to Case II by 97.25%. Eventually, the entrance of ESs into the MG brings higher revenues and lower costs for Case IV than other cases, so that the revenues outweigh the costs and the proposed PMS reaches 615.08 $ profit from the optimal scheduling of local resources. Therefore, based on the results of Table 3, it can be stated that by applying appropriate PMS to the local MG resources, in addition to enhancing technical indices, can also acquire the financial gain for the MG in market environments.
Evaluation of indices in three HSRO(σ,δ) scenarios for Cases III and IV are listed in Table 4. The given table shows that in two cases, increasing the level of uncertainty (σ) in the HSRO model compared to the deterministic model, HSRO(0,0), reduces the maximum overvoltage and increases the energy loss and maximum voltage deviation in the MG. Concerning the reliability, voltage security and flexibility indices, increasing the level of uncertainty would increase EENS, and also decrease WSI, SE and FE rather than the HSRO(0,0) scenario. In relation to economic indices it is also observed that as the level of uncertainty increases, due to increasing energy and reliability costs and declining flexibility and security revenues, the total cost increases significantly. These results are arising from the reduction of RESs output, the increase of load and EVs demand, and reducing EVs active and reactive power discharging capacities in case of HSRO(0.1,0) corresponding to (1g)-(8g). On the other hand, increasing feasibility tolerance (δ) will increase the feasibility space of the problem, inasmuch as the results are the opposite of increasing σ. Finally, the presence of ESs in Case IV, scenario HSRO(0,0.1), brings the optimal values for all indices.

B) Analysis of MG indices:
In this section, the indices curves of security, flexibility, reliability and economic, based on the incentive prices (SIP and FIP) and VOLL penalty price corresponding to Figs. 5-7 are analyzed. Fig. 5(a-b) plots the SE of the various HSRO scenarios along with the operating cost and security benefit of the HSRO(0.1,0) scenario based on SIP changes. According to Fig. 5(a), the SE curve has three zones, the first zone is linear (SIP ∈ [0,5]) where the SE increases linearly with SIP. The second zone (SIP ∈ (5,35]) is a semi-saturated zone where SE has a non-linear relation with SIP. The third zone (SIP ∈ (35,50]) is a saturated zone, where the value of SE has a constant value for different values of SIP. In addition, increasing the level of uncertainty over the HSRO(0,0) scenario causes the SE curve shifts downward, contrary to the increase of feasibility tolerance which shifts upward. The reason for this is that by increasing σ/δ, the discharging capacity of EVs will decrease/increase compared to the HSRO(0,0) scenario, while their demand for energy will increase/decrease. Hence, the SE which is dependent on the charging and discharging power of EVs, according to (32a), decreases with increasing σ and increases with increasing δ. Also, with increasing σ, SE reaches the saturation point for lower SIP values than HSRO(0,0) scenario, but with increasing δ, it will reach the saturation point for higher SIP values. Furthermore, it can be clearly Table 3 The value of different indices in the MG based on the proposed MILP model in the HSRO(σ = 0, δ = 0).  Fig. 5(b) that as the security benefit is dependent on SIP and SE, (1a), it increases linearly when SIP increases even for cases in which SE is constant, SIP ∈ (35,50], while the operating cost increases nonlinearly. Fig. 6 also shows the trends of FE in different HSRO scenarios along with the operating cost and flexibility benefit in HSRO(0.1,0) based on FIP. It is seen that the fluctuations are similar to Fig. 5, it implies that FE has the strongest linear growth between 0 and 5 ($/MWh), then, it experiences a gradual nonlinear increase from 5 to 25 ($/MWh), after that, it flattens out at just under 40 for HSRO(0,0.1). Now, turning to the details increasing the level of uncertainty reduces the FE, and reaches the saturation point for lower FIP values than the HSRO(0,0) scenario. It is worth noticing that the increase in feasibility tolerance has a reverse effect rather than the uncertainty level. As well the amount of flexibility benefit in this figure increases when the FIP rising, even for constant FE.
Finally, regarding the reliability index, based on Fig. 7(a) it can be said that with the increase of VOLL penalty price EENS through the proposed PMS is sharply decreased. Also, by increasing the feasibility tolerance EENS reaches zero for lower VOLL values, in contrast to the uncertainty level rising., It is also seen from Fig. 7(b) that increasing VOLL reduces reliability cost and subsequently increases operation cost. Because with increasing VOLL, all local resources are trying to reduce EENS, while the operation cost due to receiving more active power from the upstream network. C) Evaluation of voltage of critical buses: Fig. 8 shows the daily mean voltage curve of critical buses in different HSRO scenarios. The given figure illustrates that the voltage variations over 24 h for different HSRO scenarios are about 0.002 p.u. or 38 V. It shows that the PMS ensures voltage stability of critical buses by transiting between the operating modes of FRs. Also, increasing the level of uncertainty in the HSRO(0.1,0) scenario compared to the HSRO(0,0) scenario causes the mean voltage curve of the critical buses shifts downwards, because if the uncertainty level increases, unlike the growth of EVs demand, its discharging capacity decreases. However, increasing feasibility tolerance has the opposite effect of increasing the level of uncertainty. D) Evaluation of flexible resources performance: In this section, the daily active and reactive power curves of FRs such as EVs parking and ESs are considered in different HSRO scenarios, which is shown in   Figs. 9-10, respectively. According to Fig. 9(a), EVs charging is performed at hours 1:00 to 7:00, which corresponds to low price energy hours to meet the needed energy of daily traveling. Also, they are recharged between 12:00 and 17:00, which are the average energy price hours to inject their stored energy over 18:00 to 22:00, which is intended as a high energy price period. This process lessens operating cost in accordance with the PMS while EVs perform charging and discharging operations respectively at off-peak and on-peak hours. The provided figure illustrates that the charging and discharging power of EVs will increase if the level of uncertainty (σ) increases (HSRO (0.1,0)), because the amount of energy consumed by EVs based on (1g)-(8g) increases under this condition relative to the HSRO (0,0) scenario. In addition, the reactive power curve of the sum of EVs is plotted in Fig. 9(b), it could be plainly viewed that when the injection power of RESs based on Fig. 4 is high, i.e., 9:00 to 13:00, EVs inject less reactive power into the MG than other operating hours. Because the MG overvoltage should be in the allowed range and at this time the level of active power injection of the RESs is high, so the EVs do not even inject the reactive power into the MG between 10:00 and 11:00. In contrast, from hours 1:00 to 8:00, the reactive power injection of EVs to compensate the voltage drop caused by the high energy demand of FRs charging in these hours is notable. Similarly, due to the high power demand of passive loads from the MG, the EVs highly inject reactive power to compensate the MG voltage drop in the range of 14:00 to 24:00. As a matter of fact, by increasing the level of uncertainty, EVs reactive power generation capacity decrease according to (1g)-(8g). Thus, the reactive power curve of EVs shifts upward compared to the results of the HSRO(0,0) model, and the rate of reactive power injection is reduced, contrary to the increase of feasibility tolerance that shifts downward. As shown in Fig. 10(a), ESs to minimize the cost of MG operation likewise EVs store the energy received from the upstream network and RESs at low and medium energy price hours, i.e., 1:00-17:00 and    23:00-24:00, and delivering to the MG at peak time and high energy prices, 18:00-22:00. As the uncertainty level increases, the ESs injection to support the increased energy demand of EVs rises in HSRO(0.1, 0) scenario, while increasing feasibility tolerance, due to the increase in solution space in accordance with (1g)-(8g), has the reverse impact. Furthermore, comparing Figs. 9(b) and 10(b) explicitly expresses the reactive power curves of ESs and EVs have the same trend, inasmuch as ESs inject lower reactive power at the period of 9:00 to 15:00 because of the notable generation of RESs and avoiding MG overvoltage. In other hours, ESs notably inject reactive power to compensate the voltage drop caused by the demands of ESs, EVs and passive loads from the upstream distribution network.
In conclusion, the presence of ES apart from demand-side management, given the results in this section, addresses the contradiction between the uncertainties and an increasing need of CL's power quality. It can quickly stabilize the voltage fluctuations and balance supplydemand by adaptively controlling the NCL. Also, it provides an evident impact on the efficiency of MG by providing auxiliary services during the period of network stressful operation.
Further, the proposed PMS enables MGO to effectively coordinate MG components in such a way the MG is enhanced to a highly reliable, flexible, and self-healing system. Moreover, optimal and coordinated management of FRs provides rapid, real-time, and appropriate actions to the uncertainties. Finally, the flexible character that ES brings to the MG operation has overcome the financial challenges.

Conclusion
In this paper, optimal scheduling of MG using ESs and EVs based on a flexi-reliable and secure power management strategy as a hybrid stochastic/robust optimization was presented. In the first step, the basic and non-linear model of the proposed scheme is formulated with the aim of minimizing the sum of the operating and reliability costs, and flexibility and security revenues. The proposed model constrained to the MG OPF and the corresponding constraints of flexibility, reliability and voltage security. Then, the proposed problem model was converted into MILP equivalent model using conventional linearization methods. In addition, scenario-based stochastic programming has been used to capture the uncertainty of load, market energy price, the maximum RES output, and MG equipment availability. Also, to achieve the robust capability of EVs in improving the flexibility and reliability of the secure-renewable MG, bounded uncertainty-based robust optimization was used to model EV's uncertain parameters. Finally, based on the numerical results, it was observed that the proposed MILP model is able to achieve a globally optimal point with a low computational time by using the CPLEX solver with lower computational error than the original model. Also, implementing optimal performance between local MG resources such as RESs, EVs and ESs based on the proposed strategy has led to MG even in the worst-case scenario remained stable with high reliability, voltage security and flexibility. Additionally, flexibility and security benefits in the MG outweigh the operating and reliability costs, and MGO has achieved a financial advantage for the MG. These results further advocate the valuable role of ES as a potential key component in the newly emerging smart grid. In future research, more types of FRs and technical flexibility indices, as well as FRs' storage degradation will be investigated to enhance this model to a more efficient and practical model.