Comparative Study on Optimization Solvers for Implementation of a Two-Stage Economic Dispatch Strategy in a Microgrid Energy Management System

: A microgrid energy management system (MEMS) optimally schedules the operation of dispatchable distributed energy resources to minimize the operation costs of microgrids (MGs) via an economic dispatch (ED). Actual ED implementation in the MEMS relies on an optimization software package called an optimization solver. This paper presents a comparative study of optimization solvers to investigate their suitability for ED implementation in the MEMS. Four optimization solvers, including commercial as well as open-source-based ones, were compared in terms of their computational capability and optimization results for ED. Two-stage scheduling was applied for the ED strategy, whereby a mixed-integer programming problem was solved to yield the optimal operation schedule of battery-based energy storage systems. In the first stage, the optimal schedule is identified one day before the operating day; in the second stage, the optimal schedule is updated every 5 min during actual operation to compensate for operational uncertainties. A modularized programming strategy was also introduced to allow for a comparison between the optimization solvers and efficient writing of codes. Comparative simulation case studies were conducted on three test-bed MGs to evaluate the optimization results and computation times of the compared optimization solvers .


Introduction
Microgrids (MGs) have gained much attention as an important building block for future power systems. MGs provide an effective means to accommodate the high penetration of renewable generators and minimize power transmission losses by supplying local loads using distributed energy resources (DERs). A microgrid energy management system (MEMS) is a hierarchical supervisory control system responsible for the reliable, secure, and economical operation of MGs [1,2]. For the economical operation, the MEMS optimally coordinates the operation schedules of dispatchable DERs (e.g., micro turbines, batteries, and controllable loads), such that the MG's daily operating cost is minimized while considering the forecasted load demand and renewable generation. This optimization process is referred to as economic dispatch (ED), which is performed as the secondary control in the three-level hierarchical control actions conducted by the MEMS [3]. Actual ED implementation in the MEMS relies on an optimization solver. The optimization solver is a software package developed for solving optimization problems according to pre-programmed solution algorithms. Various off-the-shelf optimization solvers, including the commercial and open-Energies 2020, 13, 1096 2 of 22 source-based ones, are available for ED applications. Commercial optimization solvers generally have superior computational capability compared with their open-source counterparts. The ED implementation environment is then programmed either in general-purpose programming languages (GPLs), e.g., Python and Julia, or in algebraic modeling languages (AMLs). AML is a highlevel programming language dedicated to optimization applications. Popular AMLs include GAMS and AMPL. Table 1. Summary of previous studies on the implementation of economic dispatch (ED) strategies.

Ref
Optimization Solver Types (2) Language Types (3) Demonstration Methods (  Numerous studies on ED have been conducted from the perspectives of its optimization algorithms [4][5][6][7][8][9][10] and actual implementation methods [11][12][13][14][15]. Table 1 summarizes the optimization algorithms, optimal solvers, programming languages, and demonstration methods used in these studies. The optimization algorithms to the ED problem have been well studied. Particularly, in references [4][5][6][7][8][9][10], various analytic-based optimization methods, such as linear programming (LP), mixed-integer programming (MIP), and convex optimization, have been investigated for the application to ED. In reference [4], day-ahead ED was formulated as an LP problem to determine the hourly planned operation set points (i.e., charging and discharging power references) of batterybased energy storage systems (ESSs) and electric vehicles (EVs). In reference [5], convex optimization method was used to determine the optimal charging and discharging schedules of an ESS. In reference [6], an online battery power control method based on an MIP formulation was used over a rolling horizon window. The ED strategies proposed in [4][5][6] planned the operation schedule only once and the schedule is followed without any modifications during operation. However, one-time scheduling strategies are vulnerable to uncertainty in the operational conditions, such as forecasting errors in load demand and renewable generation. To address this problem, a two-stage ED strategy proposed in [7][8][9][10] was employed. In the first stage, the optimal operation schedule is developed one day before the operating day (i.e., day-ahead scheduling). In the second stage, the operation schedule obtained from the first-stage ED is updated repeatedly during actual operation (i.e., real-time scheduling). In references [7,8], the two-stage ED strategy was applied for combined cooling, heating and power (CCHP) MG, and networked MGs. Two-stage ED strategy determines the daily charging and discharging power profiles of batteries in reference [9], and the schedules of ESSs and controllable loads in reference [10].
With regard to research on the actual implementation of ED, in references [11][12][13][14][15], field demonstrations were conducted to validate ED strategies under real MG operation. Specifically, all these studies considered battery-based ESSs as their dispatchable DERs for ED and most of them [11][12][13][14] established their ED environment using commercial optimization solvers, as shown in Table 1. In references [11][12][13], the ED implementation environment was implemented using GAMS and an optimal schedule was determined using a commercial optimization solver. In reference [13], an ESS scheduling model called a power generation-side strategy, which was defined as a MIP problem, was implemented to schedule the 2.2 kVA ESS under forecasting errors in load demand and photovoltaic (PV) generation. In reference [14], based on the power generation-side strategy used in [13], an MIPbased ED strategy was implemented considering Spanish self-consumption regulatory constraints, which were introduced to facilitate the high penetration of renewable energy resources into MGs, via a commercial AML of AIMMS and its add-in optimization solver AIMMS Outer Approximation (AOA). Due to the guaranteed high computational capabilities of commercial optimization solvers, they have been commonly adopted for an optimal scheduling application of power systems. They are particularly effective for the optimization of complex problems involving a substantial number of decision variables. However, such extreme computational capabilities provided by commercial solvers may not be necessary for the application to ED of MGs, which generally involves determining the operation schedules of a small number and simple combination of energy devices. The licensing fees to use the commercial solvers can lead to significant increase in the setup costs for MGs, which can aggravate the economic feasibility of the MG operation [16].
In this context, open-source-based optimization for power systems has recently received much attention [17,18]. For example, in reference [15], the ED was implemented in a completely opensource-based environment established using Python, with the incorporation of a coin-or branch and cut (CBC) open-source optimization solver. In this case, a campus MG scheduling was on time, and the daily operational costs were reduced by 58% using open-source optimization. In reference [19], a completely open-source-based ED environment was established using the open-source programming language Julia and its corresponding open-source optimization solver ECOS; computational capability was verified over a 1-year period for a massive operation dataset measured from a German transmission network. The results demonstrated that the open-source environment is suitable for hourly ED of a large-scale grid.
However, for an open-source-based optimization solver to be further widely adopted for ED of MGs, the suitability of their computational capabilities for ED application should be examined in detail. In particular, to lay a firm basis to use open-source-based solvers, the computation time and optimal results obtained using open-source-based solvers must be compared with those accomplished using the commercial ones. A comparison of the commercial and open-source optimization solvers under various operating conditions will help MG operators select an appropriate optimization solver that satisfies their technical as well as budgetary requirements. The differences in the computational capabilities of open-source-based solvers and commercial ones were examined for several test optimization problems in references [16,20], and the open-source-based solvers showed comparable computational capability as those of the commercial solvers in some cases.
Being motivated by the previous comparative studies on optimization solvers and the recent increase in the needs for open-source-based ED environment, this paper presents a comparative study on various optimization solvers to investigate their suitability for ED in MEMS environments. The main contributions of this manuscript are as follows:


The computational capabilities of the four widely employed optimization solvers, including CPLEX and Gurobi as commercial solvers and GNU linear programing kit (GLPK) and CBC as open-source solvers, are compared in terms of optimal cost, scheduling results, and computational time. The results will be helpful for MG operators, who seek to find a costeffective optimization solver that best fits their technical requirements and budgetary constraints.  To examine the applicability to both day-ahead and real-time scheduling, the two-stage scheduling algorithm was adopted as an ED strategy, whereby operation set points of ESSs and the power injected from the main-grid are coordinated to minimize the operating costs of gridconnected MGs.  To compare the optimization solvers considering actual operating conditions, a 6-bus campus MG with actual renewable generation and load data [15] was used as a simulation test-bed. For a more comprehensive investigation, additional simulation case studies were conducted on IEEE 33-bus and 123-bus test systems [21], which were modified slightly to emulate mid-and largescale MGs, respectively.  A modularized programming strategy is presented for a fair comparison between the optimization solvers. The modularized programming strategy provides an effective way to avoid redundant usage of functions and variable definitions. It also allows the overall architecture of the code to be developed easily with enhanced readability, which is critical for managing and debugging the programming code.
The rest of the manuscript is organized as follows. Section 2 presents an overview of the ED implementation environment. Section 3 reviews the two-stage ED algorithm and explains its implementation using modular programming architecture. Section 4 discusses simulation case study results, and Section 5 concludes the paper. Figure 1 presents an overview of the implementation environment for the ED in the MEMS. The implementation environment can be programmed either in GPLs or in AMLs. Compared with GPLs, AMLs have the advantage of having a syntax similar to that of the mathematical notation of optimization problems, which allows for intuitive writing and understanding of the code. In the GPLbased implementation, this difference can be complemented by adopting a modeling package. The modeling package provides a convenient programming interface similar to that of AMLs, so it is widely employed with optimization solvers when an optimization environment is established using GPLs. However, the overall computational speed slows significantly as additional processes are added due to the modeling package. This will be discussed in detail, along with the simulation results, in Section 4. Optimization solvers consist of functional libraries, a set of functions dedicated to optimization. Functional libraries are often classified into two types according to their tasks: modeling libraries and optimizer libraries. The modeling library formulates the optimization problem for the ED using input data, such as the predicted load demands and renewable generation, MG parameters and energy resources, and electricity price. The optimizer library is responsible for finding the optimal solution according to a pre-programmed solution algorithm, which is often based on analytic methods involving LP, MIP, and convex optimization.

Configuration of the ED Implementation Enviroment
When an optimization solver is used within the environment established in AMLs or employed with a modeling package, the modeling library is not used. Instead, the optimization problem is first formulated by the AML or modeling package and is then passed on to the optimizer library. An optimization solver is compatible with certain kinds of modeling packages. Table 2 lists the optimization solvers based on their type and compatible modeling packages. Each GPL supports different modeling packages. For example, PuLP and Pyomo are the modeling packages supported in Python environments, and JuMP is the modeling package available with Julia.

Modularized Programming Architecture
For a fair comparison of computational capability, the optimization solvers should be executed within the ED environment programmed with the same architecture. In this study, this was achieved by modularizing the entire code according to the tasks. Figure 2 presents the modularized architecture for implementing the two-stage ED strategy. Module 1 (M1) is responsible for loading the data, e.g., load demands, renewable generation forecasts, parameters of the energy resources, and electricity price, into the ED environment. Such information can be obtained from the forecasting system or the supervisory control and data acquisition (SCADA) system established in the MEMS [22]. Based on the acquired information, M2 and M3 formulate the A and b matrices, representing the equality and inequality constraints of the ED, respectively. The equality constraints include equality operating conditions, as well as physical models of the MGs and devices (e.g., power balance equations and battery energy conservation). The inequality constraints comprise inequality operating conditions for reliable MG operation, such as the state-of-charge (SOC) limits and charging and discharging power limits of ESSs; Section 3 presents the details of the constraints considered in this investigation. In M4, A and b matrices are provided as inputs to the modeling package or optimization solver to formulate and solve the optimization problem. Figure 2 presents the first and second scheduling stages of the two-stage ED strategy, executed using an identical program architecture. Therefore, redundant usage of functions and variable definitions that occur when programming the two scheduling stages separately is effectively avoided. Additionally, the overall structure of the code can be developed easily with enhanced readability to make the writing and debugging of the programming code more efficient.  Figure 3 presents a schematic diagram of a grid-connected MG, where PV generators and Li-ion battery-based ESSs are considered as DERs. Battery-based ESSs have high energy densities and rapid dynamic responses, which make them suitable for high-power and high-frequency cyclic operation [23]. The operation set points (i.e., charging and discharging power references) of the ESSs were dispatched from the MEMS via communication links. The operating cost of an MG is given by the total daily electricity fee charged to the imported power from the main grid (i.e., utility power). To minimize operating costs, the MEMS implements the two-stage ED strategy, while considering the forecasted load demand and renewable generation. Note that the performance of the optimization module in the MEMS depends on the accuracy of the forecasting module. Even though the forecasting error is inevitable, it is assumed that the forecasting module provides accurate forecasts of load demand and renewable generation output power. This is not a radical assumption, and already there have been many studies that show high accuracy with many forecasting algorithms such as linear regression, clustering, and support vector machine (SVM) [24,25]. When the amounts of accumulated data in the database increase over time, the accuracy of forecasts will improve. In another approach to deal with the variation and uncertainties, robust ED strategies based on stochastic optimization were proposed in [26,27]. This approach can provide robust and efficient solutions, but modeling with scenarios can increase the computational difficulty of the optimization problem. To address this problem, a problem modeling approach proposed in [28] and a scenario selection algorithm in [29] were proposed. However, we did not consider the representative scenarios since we applied not the stochastic optimization method but the deterministic optimization approach to the proposed ED strategy. For brevity, this stochastic optimization approach is not further discussed in this paper, where instead we focus on the deterministic one. Furthermore, it is possible to get more accurate forecasts as real-time operation becomes more imminent [30]. Thus, we proposed the two-stage scheduling to mitigate forecast errors and improve the accuracy of scheduling.

Two-Stage ED Strategy
In the two-stage ED strategy, the operation schedule of the first stage is such that it minimizes the total daily operating cost of the MG. First-stage scheduling is conducted only once, one day prior to the actual day of operation (i.e., day-ahead scheduling). The single scheduling time interval of the first stage is 1 h. The scheduled profiles do not change during an hour. The second stage is implemented to mitigate the difference between the scheduled utility power and the second-stage scheduled utility power, which arises due to variabilities in load demand and renewable generation. Figure 4 presents the timeframe of the ED strategy. In second-stage scheduling, the operating set points of ESSs and the utility power scheduled in the first stage are updated every 5 min during actual operation, with a unit scheduling time step of 5 min. The second stage is carried out in an hour-ahead scheduling mode, with a unit scheduling time step of 5 min. The time horizon of the second stage diminishes by 5 min at each execution, as shown in Figure 4, and recovers over the subsequent hour. For the first 5-min interval of every hour, the operating set points of the ESSs are given as scheduled in the first stage. After the first 5-min interval, the operation schedule is updated every 5 min by the second-stage scheduling during the remaining 55 min, with a diminishing time horizon. Thus, the initial scheduling of every hour has more time steps to schedule (i.e., 11 time steps) compared with subsequent scheduling within the hour. In this way, second-stage scheduling effectively compensates for the uncertainties inherent in dayahead scheduling. 2. Second-stage scheduling: In the second stage (i.e., hour-ahead scheduling), the optimal schedule obtained from the first stage is updated every 5 min during the day of operation to compensate for the uncertainties in load demands and PV output power. It should be noted that hour-ahead scheduling is not performed at the beginning of every 5-min interval, but rather follows the dayahead scheduled profile. In other words, the second stage runs 11 times per an hour. Only the first interval of each run set serves as a final decision and the rest of the intervals are for reference only.
The following subsections present the detailed objective functions, constraints, and model for implementing the two-stage ED strategy; most of these have been employed in the author's previous work [15].

Objective Function
The objective function of the first stage is given as where ct is the hourly electricity price of the main grid and Pu,t d is the imported power from the utility at time t; α1 is the weighting coefficient for the penalty function and U1,t is the continuous decision variable determining the penalty function at time t; Pcont is the contracted power and Td is set of hourly periods. Equation (1) represents the objective to minimize daily operation costs, which consist of the charged electricity fee on the imported power from the utility side and the penalty function to maintain the utility power within the contracted power. For the hourly varying electricity fee, timeof-use (TOU) rates were considered in this study. Equation (2) describes the penalty term. The objective of second-stage scheduling is to mitigate the difference between scheduled utility power and second-stage scheduled power due to uncertainties in load demands and renewable generation. The objective function of the second stage is given as where Pu,t h is the hour-ahead scheduled utility power at time t; α2 and α3 are the weighting coefficients for penalty functions; U3,t and U4,t are the continuous decision variables determining the penalty function at time t, which consider the exceeded SOC value; SOCn,t h is the state-of-charge of nth ESS in the second stage at time t; SOCn min and SOCn max are the minimum and maximum SOC limit of the nth ESS; Th is the set of second stage scheduling time periods. In Equation (3), the second and third terms are penalty functions to prevent the overcharge and discharge of ESSs due to parameter uncertainties and nonlinear operation characteristics, such as natural battery discharge.

Equality Constraints
The equality constraints consist of the physical model of the MG and ESSs, as well as equality operating conditions to assure the uniform and continuous participation of ESSs in the ED. The physical model of MG and ESSs are represented by the power balance equation and energy convergence in the battery, respectively. The power balance between net power production and consumption in an MG is represented by the linearized DistFlow formulation [31]. The power balance equation used in day-ahead and hour-ahead scheduling is given as where Pj,i,t and Pi,k,t are the power flow from bus j to bus i and from bus i to bus k at time t, respectively; PLi,t f and PPV,t f are the forecasted load demand and PV generation at bus i at time t, respectively; Pch,n,t and Pdch,n,t are the scheduled charging and discharging power of nth ESS at time t, respectively; B is the set of buses; and L is the set of power lines in an MG. The battery-based ESSs are modeled by the energy convergence stored in a battery, which is represented by the variation in SOC level with respect to the input power, given by , , where SOCn,t is the SOC level of the nth ESS at time t; ƞch,n and ƞdch,n are the charging and discharging efficiency of the nth ESS; CAPn max is the maximum capacity of the nth ESS. Generally, charging and discharging of batteries is strictly determined by cost-effective operation constraints. However, without using constraint (8), the final SOC level might be extremely low, especially when the PV output power is deeply low and load demand is high. To ensure the continuous participation of ESSs with uniform performance in next-day ED, the SOC levels at the initial and final time steps of first-stage scheduling should be equal [32]. The energy stored at the last time of scheduling should be set to its initial value, as where SOCn,init and SOCn,final are the SOC levels at the initial and final time steps of a day, respectively. The equality constraint (8) may be relaxed to inequality constraints that represent the acceptable variation in the SOC with respect to the initial SOC level, to allow for more flexible operation of ESSs [33]. Meanwhile, the objective functions given in Equations (1) and (3) should be linearized to be applied to the two-stage MIP strategy. Here, the piecewise linearization method is employed. The linearization introduces additional equality constraints and variables, as , h n t t t t h SOC =U +U +U +r t T   (10) where r0 and r3 are the parameters for objective function linearization and U2,t and U5,t are the continuous variables for objective function linearization at time t.

Inequality Constraints
Inequality constraints include charging/discharging power capacity and SOC limits of ESSs, as well as contracted power penalty functions. The SOC limit constraint is given by Charging/discharging power limits of an ESS is given as  (13) where un,t is a binary decision variable to determine the operation status of the nth ESS at time t (i.e., un,t = 1: charging and un,t = 0: discharging).
To prevent the utility power from exceeding the contracted power in real-time scheduling, a surcharging system is considered in the second scheduling stage. The electric meter installed by Korea Electric Power Corporation (KEPCO), South Korea, estimates the peak power over a 15-min period [34]. When the peak power that updated from the electric meter exceeds the contracted power, KEPCO imposes extra charges for overuse. The surcharging system requires that the utility power averaged over a 15-min period be maintained within the contracted power to stabilize the supply. This requirement is imposed by Additionally, linearizing the first-stage objective function (1) using the piecewise linearization method yields the inequality constraints where r1, r2 are parameters for objective function linearization and S1,t is a binary variable for linearization of objective function at time t.
Similarly, linearizing the second-stage objective function (9) gives additional inequality constraints expressed by where r4 and r5 are parameters for objective function linearization and S2,t and S3,t are binary variables for objective function linearization at time t. Finally, nonlinear objective functions (1) and (3) are linearized to a set of equality constraints (9), (10), and inequality constraints (15)- (19), which are in the form of an MIP problem tractable for analytic method-based optimization solvers. Figure 5 presents the coupling of the two-stage ED strategy using a modular program architecture, as discussed in Section 2. Four modules were deployed in both the first and second scheduling stages. The detailed tasks of each module are described below.

Module (1) Loading the input data for each scheduling stage as follows:
-First-stage scheduling: MG topology; PV generation and load demand forecasts; electricity price; initial SOC level of ESSs; ESS parameters.
-Second-stage scheduling: First-stage scheduled utility power; MG topology; 5-min sampled measurements of PV generation, load demand, SOC levels, and utility power; electricity price; initial SOC level of ESSs; ESS parameters.

Module (2)
Creating the A and b matrices that represent the equality constraints of an optimization problem using the following Equations: -First-stage scheduling: Equations (6)-(9).

Module (4)
Formulating the optimization problem from the A and b matrices obtained from M2 and M3 using the modeling package or modeling library of the optimization solver in a GPL-based platform or using the inherent interface provided in the AML-based environment; solving the optimization problem by executing the optimizer library of the optimization solver.

Simulation Conditions
Four widely adopted optimization solvers, CPLEX, Gurobi, GLPK, and CBC, were applied for the two-stage ED strategy in simulation case studies, to compare computational capabilities. Figure  6 presents three grid-connected MGs considered for the test-beds containing loads, PV generators, and ESSs. The 6-bus MG shown in Figure 6a represents a small-scale campus MG in operation [15]. The IEEE 33-bus and 123-bus radial test systems [21] shown in Figure 6b and Figure 6c, respectively, were also examined, with slight modifications to emulate mid-and large-scale MG operations. All ESSs in the test-bed MGs were assumed to have the same parameters listed in Table 3, which is an actual parameter set provided by an ESS manufacturer. Table 4 lists the three-level (i.e., offpeak, mid-, and peak-loads) TOU rates used in the simulation case studies; these are the rates used by KEPCO, South Korea [35]. Figure 7 presents the hourly sampled actual load demand and PV generation forecasting data provided in [36]. For simplicity, the PV generation and load demand information of each bus are assumed to have the same forecasting profiles. The real-time variations in load demands and PV generations were emulated by applying the 5-min sampled stochastic variations to the forecasting data shown in Figure 7, which are uniformly distributed within the range of ±5% of them. The contracted utility power was set to 2 MW. Table 5 lists the specifications of the implementation environment. The compared optimization solvers were applied with the opensource modeling package PuLP in Python, using an Intel Xenon Silver 4114 CPU with 192.0 GB RAM. The computational capability of the optimization solvers was evaluated by investigating the optimal costs, scheduling results, and computation time. Simulations were performed 100 times for each optimization solver, and the average computation times were compared.

Comparison of Optimal Costs and Scheduled Profiles
All the compared optimization solvers identically determined the minimal operating costs of the three MGs as USD 86.20, USD 1698, and USD 8610; thus, the same optimal results were achieved. Figure 8 presents the scheduled power profiles of the ESS and utility power for the three MGs. The ESS power profiles shown in Figure 8b and Figure 8c are those of arbitrarily selected ESSs in the 33bus and 123-bus systems (Figure 6b and Figure 6c, respectively). Note that positive power represents discharging of the ESS. In contrast, the ESS power profiles differed from each other, as shown in Figure 8a-c. However, they have a common tendency; the ESSs are charged when the TOU is low (e.g., from 18:00 to 24:00) and discharged when the TOU is high (e.g., from 10:00 to 18:00). Note that all optimization solvers scheduled the ESS power to be maintained well within its limit (i.e., 50 kW). Similarly, the scheduled utility power profiles shown in Figure 8d-f have a common tendency to import the high power from the main grid when the TOU is inexpensive (e.g., from 18:00 to 24:00) to enhance MG operation efficiency by lowering operating costs. This tendency was due to not only the low TOU price but also the increased load demand at that time, as shown in Figure 7. Figure 9 presents the variations in the scheduled SOC levels corresponding to the ESS power profiles of Figure 8a-c. All of the scheduled SOC levels were well maintained within the limits for secure operation (i.e., 20% to 90%). Additionally, all of the compared optimization solvers succeeded in making the SOC levels at the final time step equal to the initial SOC levels, assuring continuous participation in the next-day ED with uniform contribution.  Table 6 lists the computation times of the optimization solvers for the three simulation cases. To investigate the computation time without a modeling package, additional tests were conducted with CPLEX and Gurobi, which provide application programming interfaces (APIs) for the Python environment. Using APIs, the optimization solver can be interfaced directly with the ED environment without a modeling package. In first-stage scheduling, all of the compared optimization solvers yielded a suitable computation time for the day-ahead scheduling application in all cases.

Comparison of Computation Time
Specifically, for all three cases, the open-source optimization solver CBC yielded the maximum computation time (i.e., 1. 32, 12.42, and 113.64 s, respectively). For the second-stage scheduling, the maximum and average computation times of the initial scheduling of each hour were investigated. This is because the initial scheduling has the most time steps to schedule (i.e., 11 time steps) as shown in Figure 4. Thus, the longest computation time will occur during the initial scheduling, and subsequent scheduling will have shorter computation times. For all optimization solvers, the computation time for the second-stage scheduling also yielded sound results for the hour-ahead or real-time scheduling application, with the maximum and the average computation time for initial scheduling of all the optimization solvers completed within 300 s (i.e., 5 min). Figure 10 presents the computation times obtained with the modeling package. In both firststage and second-stage scheduling, the open-source optimization solvers (i.e., CBC and GLPK) and the commercial optimization solvers (i.e., Gurobi and CPLEX) yielded similar computation times. The main reason for the similar computation times from the different optimization solvers is due to the additional processes by the modeling package: converting the input data syntax given in highlevel mathematical notation into low-level syntax that can be processed by optimization solvers took up most of the total computation time. Therefore, when the optimization solvers were executed using their APIs (i.e., Gurobi (API) and CPLEX(API)), they yield much shorter computation times; however, as noted above, the prolonged computation time due to the modeling package is a trade-off to convenient interfaces. Another factor that determines the computation times is a programming language, and the study that compared programming languages with respect to the computational time, readability, accessibility, and strengths/weaknesses was examined for general-purpose computational problems in reference [43]. However, this is not further discussed in this paper, where we instead focus on the optimization solvers and modeling package.  Table 7 lists the detailed computation times spent by the modeling package with respect to its performing of tasks (i.e., formulating equality and inequality constraints). The time ratio represents the portion they occupy in the total computation times listed in Table 6. In all cases, the modeling package accounts for more than half of the total computation time. For further investigation, the computation time data are visualized in Figure 11. In both scheduling stages, the computation time for formulating constraints becomes more dominant than the other processes for optimization in the total computation time as the scale of an MG increases. This result means that the computation speed of an ED environment is largely affected by a modeling package as well as an optimization solver. Furthermore, adopting an efficient modeling package becomes the more important issue, as opposed to which type of optimization solver to use, when the MG involves a large number of buses and DERs.

Conclusions
This paper presents a comparative study on optimization solvers to investigate their suitability for ED applications. Four widely employed optimization solvers, both commercial and open-sourcebased ones, were compared in terms of optimal cost, scheduling results, and computation time. A two-stage scheduling strategy was applied for the ED algorithm. Equations representing the operating conditions and models of MGs and DERs were formulated and combined into a set of constraints to constitute the MIP problem, which is tractable for off-the-shelf optimization solvers adopting analytic solution algorithms. A modularized programming strategy was also introduced to allow a fair comparison between the optimization solvers and the efficient establishment of the twostage ED environment. Simulation case studies were conducted on three MGs of varying size. The simulation results revealed that all of the compared optimization solvers were achieved the same optimal costs and operated within the time required for each step of the two-stage scheduling strategy. With regard to computational capability, a much shorter computation time was achieved when the optimization solvers were used with their APIs compared to the cases where the modeling package was employed together with the optimization solvers. However, the delay in computation time due to the modeling package can be compensated for by the convenient interfaces that the modeling package provides. Additionally, for large-scale MG applications, the modeling package is a critical factor that determines the computation speed of the ED environment.
Further work is required to investigate how the computational efficiency of an ED environment varies under different combinations of modeling packages, optimization solvers, and programming languages. Specifically, PuLP was adopted as a modeling package in this study, though other modeling packages, such as Pyomo or JuMP, need to be employed and analyzed in subsequent studies. The relative computational efficiency of AML-based ED environments (e.g., GAMS and AMPL) compared with those that are GPL-based needs to be examined, because the simulation results obtained in this study revealed that formulating the constraints accounts for a large portion of the computation time.