Energy Optimization for Distributed Energy Resources Scheduling with Enhancements in Voltage Stability Margin

. The need for developing new methodologies in order to improve power system stability has increased due to the recent growth of distributed energy resources. In this paper, the inclusion of a voltage stability index in distributed energy resources scheduling is proposed. Two techniques were used to evaluate the resulting multiobjective optimization problem: the sum-weighted Pareto front and an adapted goal programming methodology. With this new methodology, the system operators can consider both the costs and voltage stability. Priority can be assigned to one objective function according to the operating scenario. Additionally, it is possible to evaluate the impact of the distributed generation and the electric vehicles in the management of voltage stability in the future electric networks. One detailed case study considering a distribution network with high penetration of distributed energy resources is presented to analyse the proposed methodology. Additionally, the methodology is tested in a real distribution network.


Introduction
The growing use of distributed generation (DG) in different voltage levels has been changing the power systems operation concept. To support the network operation and also to take advantage of the distribution energy resources, it is important to develop new operation and management methodologies.
One key aspect to guarantee adequate service levels is the power system stability margin that is ensured by adequate ancillary services. These services are traditionally provided by centralized power plants with high power capacity and coordinated by the system operators in the transmission level. However, in a near future and considering the growing penetration of distributed energy resources in medium and low voltage distribution network, the system stability should also be ensured by the system operator in the distribution level, as well as by the aggregators (e.g., virtual power plants) that manage the distributed energy resources [1,2].
Behind the DG units, the consumers, storage systems, and electric vehicles (EVs) are also important to support the power system stability [3]. At the distribution level, the distribution system operator and the aggregators can participate in several ancillary services [1,4] such as primary, secondary, tertiary frequency, and voltage control; fault-ride-through capability; the congestion management; the power losses minimization in distribution networks; the monitoring the waveform quality; and the islanded operation of networks.
Extensive reviews on voltage stability indexes can be found in the literature [5][6][7][8][9], with special focus on online assessment methods for voltage stability. Nevertheless, using this index to enhance scheduling, reconfiguration, and dispatch solutions has shown its potential to improve the solutions regarding the voltage stability limitation [10][11][12]. These approaches make the -index a possible option to enhance the solutions, which was initially proposed in [13] based on the power flow equations; however, recent attention to its use, application, and possible improvements has been reported in [14]. A thorough comparison with other indexes can be found in [8], as well as a discussion on -index limitations [9].

Mathematical Problems in Engineering
The line stability indexes like the [15] or the VCPI [16] are more accurate than the -index to predict the voltage collapse proximity in a real time operation. However, in the day-ahead scheduling optimization, the main goal of the voltage stability index is not based on determining the proximity to the collapse, but to influence the distributed resources scheduling for contributing to the system stability. As stated in [16], the evolution of the -index is similar as the other suggested indexes (when the -index increases the and VCPI also increase), meaning if we optimize the -index we are also improving the and VCPI indexes. The -index is used in the present study because it is easier to integrate in the optimization problem than the other suggested indexes. With the -index, the objective function does not depend on the use of and consequently not on due to the load consumption in each bus (also depends on the resources scheduling such as distributed generation and electric vehicles).
The paper proposes an energy resource scheduling problem with a multiobjective function incorporating the operation cost and the voltage stability. This multiobjective optimization problem will be applied for a scenario in a distribution network with a high penetration of distributed energy resources, mainly an intensive EVs penetration. In the operation cost different technologies of DG and the use of EVs with griddable capability are considered, also known as vehicle-to-grid (V2G). The use of -index is proposed to deal with the voltage stability in the joint optimization problem. The -index was initially proposed in [13] based on the power flow equations. Two optimization techniques are proposed in this paper for solving the proposed multiobjective energy resource scheduling problem. These techniques will determine the nondominated solutions of the multiobjective optimization problem, namely, the weighted-sum method and an adapted goal programming methodology. Thus, the nondominated solutions represent the Pareto front that was proposed in [17], yet the application in engineering and science fields only began in the end of the seventies [18]. Furthermore, the goal programming methodology can be very useful in real application due to the complex characteristics of the objective functions. In a regular power system operation, the system operators can establish a predefined range in the operation cost objective function, and, in the critical situations (operation near from boundaries), the system operators can limit the objective function concerning the voltage stability index.
To demonstrate the effectiveness of the proposed methodologies, concerning voltage stability, two studies were included: in the first one, the sensitivity analysis is performed considering variations in the power demand, in the voltage angle, and in the voltage magnitude on the slack bus (from the distribution network's point of view the reference bus is the connection in an upstream network). In the second analysis, the loadability limit is determined for an hour considering three different scheduling objective functions (operation cost, -index, and multiobjective), allowing the determination of the maximum load that can be supplied (voltage stability boundary) considering the voltage control constraints. This approach is equivalent to the bifurcations determined with continuation power flow algorithms that allow to calculate the loadability limit for the power system [19][20][21]. Both analyses show the improvements in the energy resource scheduling problem through the incorporation of -index as another objective function. In addition, both methodologies are tested in a distribution network with high penetration of distributed energy resources, considering the use of electric vehicles allowing the -index and the operation cost evaluation. The weighted-sum method is also applied to a real distribution network to evaluate its performance in a large network.
After the Introduction, Section 2 presents an overview concerning the energy resource scheduling problem. Section 3 focuses on the mathematical formulation and on the implementation of the proposed methodologies. Section 4 shows the case study considering a 33-bus distribution network, and finally the most important conclusions are presented in Section 5.

Energy Resource Scheduling Overview and Contributions
The development of energy resources scheduling methods considering the distributed resources in different voltage levels of power systems is an important research topic. Typically, the energy resource scheduling consists in an optimization problem to determine the best scheduling to minimize the operation cost of the available resources [22]. However, in a smart grid context it is also important to take into account other aspects than just the economic one, such as power quality, voltage stability, environmental aspects, or the load diagram profile. Therefore, all these aspects can be included in the energy resource scheduling providing different solutions to help the system operators in the decision making process. Several authors have proposed different methodologies to deal with the energy resource scheduling considering distributed energy resources, such as DG and active consumers with demand response programs and the network operation. In [23], it is described a framework for aggregators to determine the energy resource scheduling based on the concept of quality-of-service in power system. A more complex negotiation perspective is presented in [24] considering multilevel negotiation layers between aggregators and electricity market participation. For a microgrid level perspective, [25] proposes a multiagent base platform allowing the scheduling of the distributed energy resources.
Other works deal with the energy resource scheduling to integrate the electric vehicles with V2G capability. A comprehensive and exhaustive review is presented in [26] concerning the impact of EVs in the distribution network. In [27], the authors proved that EVs can improve the management of intermittent renewable resources such as wind farms, and in [28] it is shown that EVs can be used to level the daily load diagram. Wu et al. [29] claim that the charging control in EVs is required for a well accommodation in the power system. To handle the large number of electric vehicles, several artificial intelligence algorithms have been proposed [30][31][32] to provide the scheduling of charge and discharge energy from EVs batteries. Another innovative perspective is proposed in [33] considering a hierarchical model to coordinate the energy resource scheduling in smart grid with electric vehicles. The integration of plug-in hybrid electric vehicles in microgrids resource scheduling is proposed [34].
The use of multiobjective functions in the energy resources scheduling problems is an important challenge to improve the quality of the obtained solutions. Some approaches are proposed considering the environment aspects [30,35] or to levelling the load diagram [28,36] in the energy resource scheduling problem. However, as is possible to see in [26], few work was developed considering the contribution of the distributed energy resources and mainly the electric vehicles to the ancillary services like the voltage stability. The inclusion of a voltage stability index in the energy resource scheduling problem turns into a multiobjective function, because it is a competing objective with the operation cost. The main contributions of this work are as follows: (1) To propose a multiobjective model to deal with the operation cost and voltage stability in the energy resource scheduling problem.
(2) To use distributed energy resources, namely, distributed generation and electric vehicles, for contributing to the power system voltage stability.
(3) To apply the weighted-sum methodology and to adapt the goal programming methodology to determine the Pareto front of the proposed multiobjective energy resource scheduling problem.
(4) Test the proposed multiobjective approach in a real distribution network.

Energy Resource Scheduling in Distribution Network
The energy resource scheduling is an important task in the present and the future power systems operation. The growing penetration of distributed generation and other energy resources, such as the electric vehicles, increases significantly the problem complexity [37]. The energy resource scheduling can consider several objective functions, most of them based on the energy costs or on the entities profits. However, technical aspects, such as the system stability, are becoming more important in new operation paradigm of the future distribution networks. In this paper, it is proposed a multiobjective energy resource scheduling for the distributed energy resources, considering two objective functions, namely, the operation cost and the voltage stability, using two different methodologies. The first methodology, called weighted-sum, is one of the most popular methods to solve multiobjectives problems. The second implemented approach is the modified weighted goal programming which is also used in several real applications. These methodologies can be used by an aggregator with the responsibility to control different distributed resources as well as part of the distribution network. The goal programming approach can result in non-Pareto optimal solutions [38], and the execution time for each simulation should be higher due to the increased number of constraints (one of the objective function is formulated as constraint) [39]. On the other hand, it is possible to obtain an approach of Pareto front with few simulations. Therefore, the use of goal programming approach was selected considering the characteristics of the objective functions. In fact, when the system is operating normally, the system operators can establish a predefined range in the operation cost objective function, and in critical situations (operation near to boundaries) the system operators can define the objective function concerning the voltage stability index. (1)

Operation Cost
For the DG units, a quadratic function is used, which is commonly employed for fossil fuel units [40]. In DG units based on renewable sources (e.g., wind or solar), the linear term ( ) of the quadratic function is the only one considered. The cost with energy acquisition to external suppliers is also considered ( SP ) that allows the balance between the DG, EVs, and demand in the distribution network. In this formulation the cost with EVs discharge ( Dch(EV, ) ) is considered and also the benefit to the aggregator from charging EVs ( Ch(EV, ) ). In addition, the battery degradation cost ( Deg(EV) ) is considered during the EVs discharging process [41,42]. Finally, two penalization costs are considered. The first one, NSD( , ) , penalizes the aggregator when nonsupplied demand situations occur. The second one, GCP(DG, ) , refers to "takeor-pay" contracts violation. These contracts are considered mainly for wind and solar units, and the penalization occurs when generation curtailment is necessary. The penalization terms are important to make a robust mathematical formulation in order to handle with critical situations from high consumer demands or high power generation from DG units. ( 2 ). In the proposed mathematical formulation, the voltage stability is achieved considering the load index ( -index) minimization. In [13] the following expression is proposed that determines theindex ( ), considering bus as a generation bus and bus as the load bus:

Voltage Stability Objective Function
In [43] and most recently in [44], a new expression is proposed to determine the -index using measurements of voltage phasors at the bus and is defined as The -index value is between 0 and 1, and the optimal value is close to 0. If the maximum -index in the system is less than 1, the system is stable in terms of voltage level. The system is unstable if the -index value is above 1 [13]. From the optimization point of view, the goal is to minimize the maximum value of -index in all buses. Basically, the -index minimization involves taking into account the bus far from the stressed condition boundaries.
The evaluation of -index implies the use of expression (3) in all consumption buses. However, in the future distribution networks there will be generation connected in several buses, changing the consumption buses to the generation buses in some periods of the day depending on the distributed energy resources installed and on the generation and load forecast in each one. Therefore, all buses are compared with the bus connected to the high voltage level in order to determine the -index, where the minimization of the function 2 , which is the maximum -index in each period , is formulated: Function 2 is a nonconvex function, which requires more time to find the optimal solution. The epigraph variable Aux is used to turn the 2 function into a convex one: where the epigraph variable Aux removes the nonconvexity of function 2 (the maximum -index) turning the optimization problem simpler to be solved. The use of the epigraph variables is detailed, explained, and illustrated in [45], turning a nonlinear optimization problem into a linear optimization problem.

Multiobjective Function: Weighted-Sum Approach ( ).
The weighted-sum method [46] transforms the multiobjective function into a single one by summing all functions ( 1 and 2 ), where each function is multiplied by a different weight ( and ), as it is formulated: where the weight factors are between 0 and 1 for giving more or less relevance to each objective function. Additionally, the sum of the two weight factors must be equal to 1. To uniform the objective functions the voltage stability price factor (SF) is included. The voltage stability can be quantified as a price signal meaning that the multiobjective function can be treated as a single objective function to optimize the cost.
In the present paper, the value of SF is equal to the energy cost of the most expensive distributed resource, as given by where the Res Sche contains the prices of all resources scheduled (DG, external suppliers, and EVs) solving the optimization problem with just the operation cost function 1 . For the DG units that use a quadratic function, it is considered an average price determined by the multiplication of the DG maximum generation power and the coefficients of the quadratic function and then divided by the same maximum generation power. Typically, the price selected will be the most expensive resource scheduled in the peak periods, because in those periods it has the highest consumption power. However different expression can be also used depending on the aggregator's strategies and on the normal network operation conditions. The weighted-sum method is the most traditional and popular method that parametrically changes the weights among objective functions to obtain the Pareto front [47].

Multiobjective Function: Goal Programming
Considering the Utopia Point Approach ( ). The goal programming was firstly proposed in [48,49] and it is used in a large range of problems in different areas [50]. Several variations of the original method have been proposed, such as the reference goal programming [51], or the Archimedean goal programming (also known as weighted goal programming) [52]. The goal programming consists in the definition of a goal for the objective function, converting the original objective function into a constraint, as it is described subjected to: In order to cope with variations in the initial goal, positive ( + ) and negative ( − ) deviation variables for each objective function should be added to the new constraint. The objective in (8) is to minimize the positive and negative deviation variables [39]. Additionally, a weight factor can be multiplied in each deviation variable turning the method into a weighted goal programming. The Pareto front can be also obtained by this method through changing the weights of the positive and negative deviation in each simulation [39,53].
The proposed methodology is based on the goal programming method with additional changes in order to  adapt this approach to the envisaged problem. The proposed methodology uses some principles of the normal-boundary intersection proposed in [54]. The main idea is to establish only one objective function as constraint, trying to optimize the other objective function. The main advantage is the ability to execute simulations only in the predefined ranges for each parameter. Figure 1 presents the flowchart of the proposed goal programming methodology.
In the first step, the optimal solution for each objective function is determined using the expressions (1) and (5), respectively. The utopia point (Up) is determined using the obtained individual points ( 1 ) and ( 2 ). The utopia point, also known as ideal point, can be described as the point with the best result for both objective functions of the multiobjective optimization problem [55]. Considering the utopia point (Up) and the results obtained in the single objective functions, the users can reduce function limits in order to guarantee more adequate solution that results in the points 1 and 2 . These values can also be determined by functions according to the operation contexts.
Afterwards, the fixed values for functions 1 and 2 are obtained using functions steps ( 1 step ) and ( 2 step ) that are calculated as where num steps is the total number of steps that is defined by the users and they can be different for each objective function. In this paper, a total number of 19 steps are considered for each objective function. Finally, the fixed values for each objective function are used to determine new solutions ( 1 ) and ( 2 ): In (11) the objective is to minimize the function 1 obtained with (1) plus the positive ( 2 + ) and negative ( 2 − ) deviations in objective function 2 . The negative deviation is multiplied by 0.1 which means an incentive to the -index reduction. The positive deviation is multiplied by 10 to penalize the increase of the -index, because the positive deviation is only used to guarantee the problem feasibility. Additionally, it is necessary to transform the objective function 2 (5) into a constraint: The same approach is used to obtain a new solution ( 2 ) optimizing the objective function 2 (12), considering the expression (5) for function 2 plus the positive ( 1 + ) and negative ( 1 − ) deviations in objective function 1 . Then, the objective function (1) is converted into a constraint: 3.5. Problem Constraints. The energy resource scheduling should use an accurate model of the network to achieve scheduling results that are feasible in the electric network (avoiding lines congestion and bus voltage violations). For this purpose, an AC power flow is included as constraint in the energy resource scheduling problem. The active power balance equation establishes that the active power injected in each bus is equal to the active power generation minus the active power demand in the same bus: The active power injected in bus is defined as the sum of all power flow through the lines that connect to this bus; the power generation is the sum of the power generation from DG and external suppliers and EVs discharge; the active power demand is equal to the sum of the EVs charge plus the forecast consumers demand in each bus less than the nonsupplied energy. A more detailed loads model should be used for dynamic evaluation of voltage stability. However, in the scheduling process, the use of forecast values allows determining the impact of the power demand in the voltage stability indexes [5][6][7][8][9].
In addition, the reactive power injected in bus is also equal to the reactive power generation minus the active power demand: where ( ) is the voltage angle difference between bus and . contains the set of all lines that are connected to the bus . and represent the real and imaginary part of the admittance matrix corresponding to the row and column, respectively. Engineering   7 For the AC power flow model, it is also to establish the maximum and minimum limits for the voltage magnitude and angle, respectively:

Mathematical Problems in
Before solving the energy resource scheduling problem, a slack bus is necessarily selected in the distribution network. For this slack bus, a fixed value for the voltage magnitude and angle is specified. Finally, the line thermal limit (upper limit) is established for the power flow from bus to bus and vice versa, as is defined: The distribution network can be connected to upstream networks by transformers that change the voltage level from high voltage (HV) to medium voltage (MV). These HV/MV transformers have an upper limit (maximum capacity): Similarly, the buses in the distribution network can have transformers from MV to low voltage (LV) connecting small distributed energy resources, such as photovoltaic units and EVs, to the MV side of bus : The TFR MV LV and TFR MV LV are determined by Regarding the DG units, the minimum and maximum limits for active, reactive, and apparent power generation are considered: where DG(DG, ) is used to control if the DG unit will be turned on or off considering the optimized solution.
Regarding DG units with "take-or-pay" contracts with the system operator, mainly renewable sources, the system operator is mandatory to dispatch all the forecasted power that is given by In the case of external supplier, a maximum limit for the active and reactive generation, respectively, is defined as where Trip(EV, ) corresponds to the typical daily travel profile for reducing the energy stored in the battery when the EV is in travel. Thus, the system operator must ensure the energy required for the EV user to travel in the time horizon of the energy resource scheduling problem. A trip forecast can be considered in the Trip(EV, ) according to the history consumption profile for each EV [56]. The charge and discharge efficiency( (EV) and (EV) ) are included in the batteries balance equation (25).
The energy resource scheduling problem also considers the maximum and minimum energy stored in the EVs batteries: where BatMax(EV, ) is related to the battery's capacity and BatMin(EV, ) corresponds to the minimum amount of energy stored in the battery and in which period that energy must be guaranteed. EV users and system operator must use an adequate communication system to exchange information about several parameters, such as the BatMin(EV, ) [57].
The charge/discharge rates have their own upper limits that are formulated as where the maximum charge/discharge limits depend on the connection points (normal charge or fast charge points). For instance, EVs can be connected in a single phase (e.g., at home); therefore the charge/discharge limit is lower than when EVs are connected in three-phase mode (e.g., a parking lot at the work).
Finally, the sum of the two binary variables for charge and discharge must be lower or equal to 1 to avoid a simultaneous charge and discharge in the same period : The on-load tap changer (OLTC) in the HV/MV power transformer is modeled (29) considering different steps (STP) and a correspondent binary variable ( TFR(STP, ) ) for each one. Equation (30) assures that only one step is used and the power voltage in the power transformer ( 0( ) ) is defined in (31). Consider the following: 3.6. Software and Solvers Used. Both methodologies have been implemented in MATLAB software interconnected with the general algebraic modeling system (GAMS) [58]. MAT-LAB is used to process all the data regarding the resources characteristics and contracts and afterwards to organize all the results. GAMS is used to run the optimization algorithms. GAMS offers a large set of solvers in the same platform. The proposed energy resource scheduling is classified as a mixedinteger nonlinear programming (MINLP) problem. In GAMS software, the DICOPT solver was used [59], because it solves the MINLP problems by splitting them into mixed-integer programming (MIP) and nonlinear programming (NLP) subproblems. The coordination between these two subproblems is important to obtain the optimal solution of a MINLP problem, and DICOPT coordinates MIP and NLP solutions through "Outer approximation," "Equality relaxation," and "augmented penalty." These coordination mechanisms will create and handle relaxed problems to be solved by the MIP and NLP solvers, afterwards the obtained solutions are penalized, then the relaxed problems are decreased until the stopping criteria of DICOPT is reached.
The two solvers used to solve the two subproblems are CPLEX for the MIP subproblems and the CONOPT for the NLP subproblems. DICOPT uses an iterative process that stops when the MIP and NLP subproblems return solutions with a difference less than a predefined error that has been fixed at 0.01%. The local optima solutions are the main  obstacle to overcome by DICOPT due to the nonconvexities that characterize a MINLP problem. Therefore, the solver does not guarantee the global optimum even incorporating algorithms to handle this kind of obstacle.

Case Study
The present section shows the main results of the proposed methodologies and it is divided into five subsections. Section 4.1 presents the information and input data used in the first case study of this paper. Section 4.2 is related to the Pareto front results of the two proposed methodologies to solve the multiobjective energy resource scheduling problem. Section 4.3 presents an evaluation concerning the voltage stability margin of different Pareto front solutions determined by the adapted goal programming methodology. In Section 4.4 the behaviour of the adapted goal programming methodology in different operation scenarios is evaluated, such as voltage angle variation, voltage magnitude variation, and load consumption variation. Section 4.5 presents the Pareto front results of the proposed weighted-sum methodology for a real scenario.

33-Bus Distribution Network Description.
In the first case study, a 33-bus distribution network in [60] is used. The network supplies 218 consumers, including domestic, commercial, and industrial consumers. The network has 66 DG units spread over the buses: 32 photovoltaic (PV), 15 combined heat and power (CHP), 8 fuel cell, 5 wind, 3 biomass, 2 small hydro, and 1 waste-to-energy (WTE) unit. The network is connected to a HV upstream network through bus 0. The distribution system operator or the aggregator can negotiate energy with external suppliers in bilateral negotiations and/or electricity markets. The negotiated energy flows to distribution network through bus 0. Figure 2 shows the 33-bus distribution network. In Table 1, the power capacity and the energy cost of each generation technology, of external suppliers and of EVs, are presented. Regarding the EVs, the management of 2000 EVs is considered that can come and go from the network. Table 2 presents the seven EVs models that are used in this case study [61][62][63][64][65][66][67]. A simulation tool [56] is used to generate the daily traveling profiles for the 2000 EVs. This simulator obtains the bus location that each EV will have to connect in the distribution network. In terms of the batteries' degradation cost, a cost of 0.03 m.u./kWh considering the work proposed in [41,42] has been defined. In addition, the EVs discharge cost was established at 0.04 m.u./kWh. This value is based on the profits of the EV's owner; however, an extra incentive established in the contracts should be considered to stimulate the participation in these events.

33-Bus Distribution Network Results.
Considering the described scenario, two proposed methodologies were performed, namely, the weighted-sum and the adapted goal programming methodologies. In the weighted-sum methodology, 500 runs have been made and in each run the weighted factors changed from 0 to 1 in steps of 0.002 (considering that in each run the sum of the two weights must be equal to 1). For the adapted goal programming methodology 40 points were tested. The simulations were performed on a computer with two processors Intel5 Xeon5 E5645 2.40 GHz, each one with two cores, 24 GB of random-access-memory. Figure 3 shows the Pareto front obtained by the weightedsum methodology.
The weighted-sum method has been able to find 360 nondominated solutions, most of them with very close values in both objective functions. Considering these results, it is possible to conclude that 500 weights is an excessive number that influences the execution time of this methodology. The operation cost changed between 6933 and 9054 m.u, and the -index changed between 0.0305 and 0.1613. In Figure 4     Considering the excessive number of weights tested in the first simulation (Figure 3), 40 weights were tested for the weighted-sum method. Using the same number of simulations in both methods allows getting a better comparison. As shown in Figure 4, the curves are overlapped and it is very hard to see differences between the two methodologies. The overlapping of Pareto front curves allows to conclude that both proposed methodologies are suitable to solve the multiobjective energy resource scheduling considering the operation cost and the voltage stability. In the goal programming methodology, the reduction of the solutions ranges ( 1 = 1 and 2 = 2 ) was not considered. However, in a real application it is possible to reduce the range of the operation costs, due to the low variation presented by -index. In this paper, the normalized distance to the utopia point was used. The normalized distance of each obtained solution to the utopia point is presented in Figure 5.
The Pareto front gives useful multiple choices to the decision maker (i.e., distribution system operator or resource aggregator). However, the decision maker must choose the "best" nondominated solution that satisfies its requirements. The choice should be made according to a specific strategy depending of each decision maker. Some methods have been proposed to select the best choice of the Pareto front [68]. In this method, the "best" solution corresponds to an operation   cost of 7236 m.u. and to an -index of 0.0642, as is shown by Figure 5. In Figures 6 and 7 the comparisons of the operation cost and the -index for the "best" solution and for the solutions considering each objective function are presented.
Regarding to the operation cost, a high variation occurs in successive periods when the minimization of the -index objective function is only considered. However, the -index is constant throughout the day with a very small value of 0.0305. The -index value, represented in Figure 7, corresponds to the maximum value in each hour, and it can be obtained in different buses in each period. Results of the "best" solution also present a constant behaviour of -index, and the operation cost is according to the expected values for this parameter.
Looking more carefully at Figure 6, it is possible to see that in some periods the operation cost obtained in the "best" solution is higher than the cost obtained by the optimal solution in the minimization of -index. This happens because the optimization process considers the 24 periods and not an optimization for each period. In fact, one of the most difficult aspects in the optimization is to schedule the charge and discharge periods of EVs. For instance, the operation cost in period 16 is higher in the "best" solution than in  the optimal solution for the individual -index optimization due to the higher amount of scheduled EVs charging in the "best" solution, which increase the energy stored in the EVs batteries. This energy will be used in the remaining periods avoiding the use of more expensive resources. The voltage magnitude in each bus is presented in Figure 8.
As expected, the use of -index function results in a more stable voltage profile in each bus. Figure 9 presents the energy resource scheduling by each technology for the "best" solution. Figure 9 shows the high impact of the PV generation during the hours with a high solar radiation (around midday). Another important aspect is the use of electric vehicles discharge at the end of the day (i.e., peak periods). The use of EVs discharge in these peak periods can be restricted by the use of these vehicles, because approximately 96% of the time cars are parked, and in only 4% of the time cars are used for transportation [69]. In the peak periods, many vehicles are traveling. However, many of them are parked with the possibility of being connected in the electric network  Figure 10 shows the EVs charge and discharge and energy stored in the battery for the "best" solution.
In Figure 10 one can see that EVs charge their battery during the night (i.e., off-peak periods), but also in hours with high PV generation, because this case study considered that PV has "take-or-pay" contracts, so the aggregator or distribution system operator must fully dispatch the energy generated by PV. Another important aspect is the use of charge and discharge processes in the same periods. In fact, each EV only charges or discharges in each period. However, the optimization schedules the charge of some EVs and the discharge of other EVs in order to guarantee lower -index values. This means that EVs can be used as a resource to improve the voltage stability margin in the future distribution networks. Some detailed information concerning the performance and execution time of the proposed methodologies is shown in Table 3.
The high number of variables, including the discrete ones, leads to a very complex problem with high execution time. The execution time for one run (a simulation with specific weights) is, in average, of around 25 minutes. However, the execution time of a simulation can change between 13 and 58 minutes. The total execution time is higher than 27 hours for the weighted-sum method with 500 different weights (or simulations). When 40 different weights are used, the time becomes more acceptable to a little bit more than 2 hours. These execution times are only possible due to the use of parallel processing (8 cores in this case). In the adapted goal programming method, the time is higher (3.3 hours), because it is necessary to determine the results for each objective function regardless of defining the simulation steps and continuing with the rest of the method's process (see Table 3: Performance and execution time.

Performance indicators
Weighted-sum method Proposed goal programming 500 weights 40 weights Average execution time for one run (in minutes) 24 Figure 1). In this case, the parallel processing only can be used after determining the simulation steps and fixed values. The memory space is not a critical aspect in the optimization process when a single simulation is considered. However, it can be a constraint when the parallel processing is used.

Voltage Stability Margin.
To determine the voltage stability margin, it is necessary to use the optimization formulation given in [19,21]: The maximum loadability will be given by the load increment parameter for each hour, and it is a parameter of interest for assessing voltage stability [19], and as described in (32) the loads are considered with a constant power model for the calculations. In (32), the limits in the generators capacity and the network power balance equations are included as problem constraints; these are presented in (15), (16), and (22) to (24). The calculation of was performed for hour 21, which is the peak hour and the hour with higher differences between the solutions with different objective functions.
Three scenarios were studied considering the use of single objective functions (operation cost and -index optimization) and the "best" solution for the multiobjective function. Table 4 shows the obtained solutions in period 21.
These results show that the proposed multiobjective energy resource scheduling methodology allows increasing the voltage stability margin effectively, and thus including the -index helps improving the solution when compared to the scheduling obtained only with the minimum operational cost (minimum 1 ).

Sensitivity Analysis.
In order to evaluate the proposed methodology behaviour in different operation scenarios, three different sensitivities analyses are performed. One of the key points of these analyses is the assumption that distribution network is connected to the transmission network. The transmission network can be represented by its Thevenin equivalent that is shown in Figure 11, in which 33-buss distribution network HV/MV connects to the transmission network from bus 0.
In the distribution network scheduling process, the slack bus is represented by the connection to the transmission network. However, this connection can be seen as a "virtual" slack bus due to the existence of a system slack bus. In practice, the "virtual" will impose the operation conditions of distribution network but depends on the Thevenin impedance of all system. Taking this aspect into account, the "virtual" slack bus can have different values of voltage magnitude and voltage angle according to the system operation scenario.
The analyses consist in the variation of the different parameters. In the first one (Figures 12 and 13), the "virtual" slack bus voltage angle is changed. In the second analysis, the "virtual" slack bus voltage magnitude is changed (Figures 14  and 15). In the third analysis, the consumers power demand is changed (Figures 16 and 17).
To perform these analyses some assumptions are considered: (i) The -index is computed considering the power system slack bus as reference.
(ii) The generation capacity has been increased two times.
(iii) The load consumption has been increased two times in analyses 1 and 2. (iv) The "virtual" slack bus voltage and angle in analysis 2 is of 50 ∘ .
(v) The lines thermal limits are increased to avoid violations.
By analysing Figures 12-17, it is possible to conclude that, in general, the minimization of the -index objective function provides better results concerning the -index value, and the minimization of the operation cost function provides better results in terms of operation cost. This is a logical and expected conclusion. Most interesting is the analysis of the "best" solution results. In fact, the results obtained using the multiobjective function give a solution for -index nearby the -index objective function and solutions nearby the operation cost objective function in terms of operation cost.
In Figures 12 and 13, it is possible to see that the -index increases with the voltage angle increasing, and the differences between the objective functions decrease. However, the use of the operation cost objective function ( 1 ) results in -index values higher than 1 for voltage angles higher than 55 ∘ , whereas the use of -index objective function ( 2 ), or the "best" solution objective function, keeps the values   of -index below 1. It is important to mention that in real operation, values of -index higher than 1 are not possible because the system would collapse. However, in the present analysis considering the day-ahead horizon this type of values can be obtained mainly because the operation conditions are significantly increased to analyse their limits. Figure 13 shows that the operation cost is more or less constant in the simulations. This means that the operation cost is independent from the voltage angle.
Regarding the voltage magnitude analysis, from Figures  14 and 15, it is possible to see that the -index is higher than 1 in three situations when the objective function 1 is used. In fact, the evolution is not constant. This happens because the objective function only considers the costs, and, in some cases, it can schedule resources at the same price, yet with different impact on the -index, as in this case. When the -index is included in the objective function, the values obtained for the -index have a constant evolution remaining below the limit. As in the voltage angle evaluation, the operation costs are more or less constant.
Analysing Figures 16 and 17, regarding the load sensitivity, it is possible to see that the use of objective function 2 limits the impact of the load increase in the -index. In this case,

L-index
Objective function 1 Objective function 2 "Best" solution  the electric vehicles charge and discharge are scheduled in different periods trying to minimize the impact of the load increase in -index. When the operation costs are included in the objective function, the -index increases, remaining below the limit. As expected, the cost increases with the load increase. The load increases four times, and the cost increases around two times. This fact can seem odd, but it can be justified by the total system demand which is composed of the load consumption and the electric vehicles charge. On the other hand, the PV panels have "take-or-pay" contracts with a price higher than the average. This means that the load increase is supplied by the generation units or external suppliers with a lower generation cost, reducing the impact in the operation cost.

Real Scenario Analysis.
Another important aspect to evaluate the proposed methodology is its application in a real scenario. The proposed methodology assumes a higher penetration of distributed generation and electric vehicles. This represents a future vision of power systems operation. In this way, the methodology was tested in a real network concerning the network characteristics, considering a scenario of DG and EVs penetration according to [70,71]. The considered network ( Figure 18) is a real 30 kV distribution network, supplied by one high voltage substation (60/30 kV) with 90 MVA of maximum power capacity distributed by 6 feeders, with a total number of 937 buses and 464 MV/LV power transformers [72]. This distribution network has already been used for many years and it has suffered many reformulations. It is partly composed of aluminium conductors and copper conductors, and the distribution is made by power lines and underground cables. The study results in 548 DG units most of them using photovoltaic panels and 464 aggregated loads in MV/LV power transformers. Figure 19 shows the obtained results using the weightedsum method considering 200 simulations with different weights.
The weighted-sum method found 70 nondominated solutions. The operation cost changes between 93,131 and 94,161 m.u. and the -index between 0.01293 and 0.9512 m.u. This means that the utopia point is defined by an operation cost of 93,131 m.u and an -index of 0.01293 m.u. The "best" solution ( Figure 20) is given by the weight factor ( ) of 0.69 resulting in an operation cost of 93,168 m.u. and an -index of 0.03506. Figure 21 shows the differences in the hourly operation costs.
The cost increases 1.87% with objective function 2 and increases only 0.04% when the "best" solution objective function is used. Regarding -index, the value obtained with objective function 2 is constant during the day and equal to 0.01242. When the objective function 1 is used, the -index is of 0.09512 imposed in the peak hour (period 12). In the "best" solution objective function, the -index is 0.03506. By analysing the differences in operation cost and in -index (see Figure 22), it is possible to see higher differences in theindex than in the operation costs. The differences in operation costs are mainly due to the electric vehicles scheduling. In fact, in some periods the EVs discharge happens to supply the charge of other EVs (see Figure 23). However, the charge and discharge allows reducing the power flow in some lines of the network, reducing the differences in the voltage magnitude and angle and consequently reducing the -index. Furthermore, the system uses a better reactive power scheduling to reduce the voltage differences between buses, reducing the -index values. The -index value is the same in most of periods (9 to 21), in the case of "best" solution. However, this value is not constant in the buses. In fact, theindex is imposed in each period by different buses according to the resources scheduling.
By analysing all the presented resources, it is possible to conclude that the proposed multiobjective energy resource scheduling and the two proposed weighted-sum and adapted goal programming methodologies can be successfully used in real networks with large penetration of distributed energy resources. The inclusion of -index objective function in the distribution networks can support significantly the management of all power system.

Conclusions
The future power systems will be operated considering a large range of different distributed energy resources connected in
A methodology to integrate the voltage stability in the scheduling of distributed energy resources in a distribution level was presented in this paper. The proposed formulation results in a multiobjective optimization problem considering the operation costs and the voltage stability. The voltage stability is "measured" by the load index value ( -index).
Two case studies are presented considering a 33-bus distribution network and a real Portuguese distribution network with 937 buses. In both cases, the obtained results show the advantage of the proposed methodologies mainly when a combination of the two objective functions is used. By using the "best" solution objective function, it is possible obtain a significantly better -index values with a short operation cost increase.
The voltage stability margin is evaluated in a peak period (in the 33-bus distribution network case study) considering the three objective functions, showing the effectiveness of the method. Additionally, a sensitivity analysis is presented considering extreme cases allowing evaluating the behaviour of the proposed method in complex operation scenarios.

1:
Operation cost function 2: Load index function Base: Base value BatMax: Battery energy capacity BatMin: Minimum stored energy to be guaranteed at the end of period Bus: Bus Ch: Charge process of the electric vehicle + : P o s i t i v ed e v i a t i o n − : N e g a t i v ed e v i a t i o n Dch: Discharge process of the electric vehicle DG: Distributed generation unit DGForecast: Forecast power of distributed generation unit in period EV: Electric vehicle : Objective function solution for the goal programming GCP: Generation curtailment power , : B u s and bus , Load: Load Max: Upper bound limit Min: Lower bound limit NSD: Nonsupplied demand Res Sche: Resource scheduled SP: External supplier Stored: Stored energy in the battery of the vehicle STP: Tap step TFR: Power transformer TFR HV MV: Transformer that connects from high voltage to medium voltage TFR MV LV: Transformer that connects from medium voltage to low voltage : L i n e .