Optimization Method for Operation Schedule of Microgrids Considering Uncertainty in Available Data

: Operation scheduling in electric power grids is one of the most practical optimization problems as it sets a target for the efﬁcient management of the electric power supply and demand. Advancement of a method to solve this issue is crucially required, especially in microgrids. This is because the operational capability of microgrids is generally lower than that of conventional bulk power grids, and therefore, it is extremely important to develop an appropriate, coordinated operation schedule of the microgrid components. Although various techniques have been developed to solve the problem, there is no established solution. The authors propose a problem framework and a solution method that ﬁnds the optimal operation schedule of the microgrid components considering the uncertainty in the available data. In the authors’ proposal, the objective function of the target problem is formulated as the expected cost of the microgrid’s operations. Since the risk of imbalance in the power supply and demand is evaluated as a part of the objective function, the necessary operational reserve power is automatically calculated. The usefulness of the proposed problem framework and its solution method was veriﬁed through numerical simulations and the results are discussed.


Introduction
Microgrids are a framework of smart power grids that manage a localized group of electrical power sources and loads. They can be operated in both connected and disconnected to the conventional bulk power grid [1][2][3]. Based on the growth in the market penetration of renewable energy sources (RESs), microgrids are now regarded as one of the most practical, sustainable power grids. In fact, many studies and developments in the field of microgrid-related technologies are currently being carried out [4,5], as well as demonstrative field tests [6,7].
There are two types of microgrid components: one is the controllable component, and the other is the uncontrollable one. Controllable power generation systems (CGs), energy storage systems (ESSs), and controllable loads (CLs) are examples of the former. The attributes of the controllable components usually depend on the judgement of whether the power grids accept the reverse power flow from them or not. Electric vehicles (EVs) are a typical example (accepted: ESS, not accepted: CL). In contrast, the latter type consists of electrical appliances on the consumer side and variable renewable energy-based generation systems (VREGs). Since it is difficult for microgrid operators to control electrical appliances, this paper does not include them as a type of CL. In microgrids, the VREGs take a significant portion of the electrical power source, and thus, they create difficulties in the management of the power supply and demand. This is because VREGs, whose outputs are strongly dependent on the weather conditions, introduce additional uncertainty into the readily

Problem Framework
As shown in Figure 1, microgrids generally consist of CGs, ESSs, CLs, electrical loads, and VREGs. Components 1-3 are controllable, and components 4 and 5 are uncontrollable. We can aggregate the latter into one uncontrollable component. This is the net load, and its values are set by the sum of the values of electrical loads and VREG outputs in each time interval. Based on the predicted values of the net load, operation scheduling of the microgrids is represented as the problem of determining a set of start-up/shut-down times and output shares of the CGs, charging/discharging states of the ESSs, and charging states of the CLs.
Energies 2021, 14, 2487 3 of 13 electrical resources at the operation scheduling stage, e.g., one day ahead, can reduce the difficulties in rescheduling the management of the power supply and demand. In other words, appropriate scheduling of the electricity trade contributes to improving the economic efficiency, reliability, and quality of the power management. Under these circumstances, trading electricity is regarded as an additional controllable factor in this paper, although microgrid operators cannot control components of the bulk power grids.

Problem Formulation
Microgrid operators design an operation schedule for all controllable components in advance. This problem is formulated as an MIP problem, and its optimization variables are defined as , ∈ , , for ∀ , ∀ , , ∈ , , for ∀ , ∈ , where is time ( = 1, ⋯ , ); is the number assigned to the CGs ( = 1, ⋯ , ); , is the ON/OFF state variable of CG (ON: 1, OFF: 0), which is an element of vectors and ; , is the output of CG , which is an element of vectors and ; and are the maximum and the minimum outputs of CG ; is the number assigned to the ESSs ( = 1, ⋯ , ); , is the output of ESS , which is an element of vectors and ; and are the maximum and the minimum capable outputs of ESS ( ≤ 0 ≤ ); is the set of periods that the ESSs are available; is the trading electricity, which is an element of vector .
In this paper, the CLs are regarded as a type of ESS, and their states are treated by Equation (3) by substituting with zero. The traditional problem frameworks often require accurate data for the uncontrollable components, which is impractical in the actual operation scheduling of the microgrids. The only available data are the predicted values of these components with the uncertainty. Here, we define the net loads as the stochastic variables and set their likely ranges as  If the microgrid operators cannot appropriately manage the balance of the power supply and demand, an extra payment is required for additional electricity trading to compensate for the resulting power surplus or shortage; this is the imbalance penalty. The imbalance penalty is normally very expensive, and therefore, microgrid operators secure reserve power to prevent any unexpected payment. On the other hand, the reservation of electrical resources at the operation scheduling stage, e.g., one day ahead, can reduce the difficulties in rescheduling the management of the power supply and demand. In other words, appropriate scheduling of the electricity trade contributes to improving the economic efficiency, reliability, and quality of the power management. Under these circumstances, trading electricity is regarded as an additional controllable factor in this paper, although microgrid operators cannot control components of the bulk power grids.

Problem Formulation
Microgrid operators design an operation schedule for all controllable components in advance. This problem is formulated as an MIP problem, and its optimization variables are defined as u i,t ∈ {0, 1}, for ∀i, ∀t, where t is time (t = 1, · · · , T); i is the number assigned to the CGs (i = 1, · · · , NG); u i,t is the ON/OFF state variable of CG i (ON: 1, OFF: 0), which is an element of vectors u t and u; g i,t is the output of CG i, which is an element of vectors g t and g; G max i and G min i are the maximum and the minimum outputs of CG i; j is the number assigned to the ESSs (j = 1, · · · , NS); s j,t is the output of ESS j, which is an element of vectors s t and s; S max j and S min j are the maximum and the minimum capable outputs of ESS j (S min ; T j is the set of periods that the ESSs are available; e t is the trading electricity, which is an element of vector e. In this paper, the CLs are regarded as a type of ESS, and their states are treated by Equation (3) by substituting S max j with zero. The traditional problem frameworks often require accurate data for the uncontrollable components, which is impractical in the actual operation scheduling of the microgrids. The only available data are the predicted values of these components with the uncertainty. Here, we define the net loads as the stochastic variables and set their likely ranges as where d t is the predicted net loads, which is an element of vector d; d max t and d min t are the maximum or the minimum values of the assumed net loads.
The operation scheduling problem is to determine the set of (u, g, s, e) in terms of minimizing the sum of the operational costs of the microgrid components. Fuel costs of the CGs are usually approximated as quadratic functions by means of their outputs, and these can be reduced by steady control of the CG outputs. As opposed to the CGs, the ESSs do not directly consume fuel, and this brings new challenges to the evaluation of the ESS operations. Based on Equation (5), the objective function is represented as min u,g,s,e E f u, g, s, e, d , where d † t is the realization values of net loads; g † i,t is the CG outputs calculated according to d † t ; p( ) is the probability density function; A i , B i , and C i are the coefficients of fuel cost of CG i; D i is the start-up cost of CG i; M t is the price of electricity trade in operation scheduling stage; V IO t is the price of the imbalance penalty.
As shown in Equation (10), ESSs contribute to reducing the operational cost of CGs and the payment for trading electricity.
The microgrid components must satisfy the following constraints: • State of duration of the CGs where u on i,t and u off i,t are the consecutive operating and suspended duration of CG i, and UT i and DT i are the minimum operating and suspended duration of CG i, respectively. • Ramp rate of the CGs where ∆G • State of charge (SOC) of the ESSs , for ∀j, t ∈ T j , If s j,t > 0 then q j,t = q j,t−1 − 1 η j s j,t ; if s j,t < 0 then q j,t = q j,t−1 − η j s j,t .
Energies 2021, 14, 2487 where q j,t is the SOC level of the ESS j; Q max j and Q min j are the maximum and the minimum SOC levels of the ESS j; η j is the overall efficiency of the ESS j.

•
Output range of the ESSs If the time interval ∆t, the ramp-up rate ∆G up i , and the ramp-down rate ∆G down i satisfy conditions (16) and (17), the constraints (11) and (12) become inactive [20,25]. Similarly, evaluation of the variables g max i,t and g min i,t is unnecessary in Equation (13). In this case, g max i,t and g min i,t become equal to G max i and G min i , respectively, and this constraint can be integrated into Equation (2).
Since these conditions are satisfied in the operation scheduling stage of the microgrids, the authors ignore the constraints (11) and (12) and treat the constraint (13) as the same as Equation (2).

Characteristics of Proposed Problem Framework
The formulated problem is to determine the optimal operation schedule, that is the set of (u * , g * , s * , e * ). There are two types of optimization variables: u is the set of discrete variables, and the others are sets of continuous variables. Moreover, we cannot omit the influence of the ESSs and the CLs from the microgrid operations because their capacities are large in regard to their share of electrical power sources. These are the reasons why the target problem is complicated as compared to that of the bulk power grids.
As shown in Equations (6)-(10), the objective function is formulated as the expected operational cost of the microgrids. Equations (7) and (10), evaluate the impact of the unexpected imbalance as a potential risk of operational cost increases. In the process, the reserve power is evaluated corresponding to the risk of a cost increase. Hence, we can remove the traditional constraints for the balance of power supply and demand and the operational margin from the problem formulation. These are integrated into the objective function.

Solution Method
In conventional approaches [29,30], the calculation of Equation (7) is often replaced with an iterative calculation using discrete probability distributions of the past VREG outputs. The results of this calculation strongly depends on the number and the variety of samples in the past record. Too small a sample number produces inaccuracy in the calculation results, while a larger number increases the costs of the iterative calculation. The authors propose a strategy enabling us to treat p d † t as a continuous function, and then apply PSO to the target problem framework.

Approximation of Outputs of Controllable Generation Systems
As shown in Equation (8), shape of the fuel cost function depends on the operating CGs, and economic output share of the operating CGs is determined by the fuel cost coefficients of each CG. In the proposed solution method, first, possible UC candidates are enumerated exactly. As opposed to bulk power grids, the number of CGs in microgrids is sufficiently small, and thus we can identify all possible UC candidates. Next, the optimal output share of the operating CGs is calculated in association with the enumerated UC candidates. Since the ELD problems become a quadratic programming problem by fixing the values of u i,t , quadratic programming solvers are applicable in this step. Finally, by nonlinear regression, each set of optimal output shares are approximated as the aggregated Energies 2021, 14, 2487 6 of 13 fuel cost functions of the UC candidates. Through this process, we can replace Equation (13) with (18), and then Equation (8) with (19).
where h t is the sum of outputs of the operating CGs (= ∑ NG i=1 g i,t u i,t ); h † t is the sum of operating CG outputs calculated according to d † t ; A u t , B u t , and C u t are the coefficients of approximated fuel cost associated with the vector u t . Table 1 summarizes the specifications of sample CGs, and Figure 2 displays their fuel cost functions. These were formulated by referring to [18,20,25,30]. When we focus on microgrids with three CGs, the possible UC candidates are (0, 0, 0), (0, 0, 1), (0, 1, 0), (1, 0, 0), (0, 1, 1), (1, 0, 1), (1, 1, 0), and (1, 1, 1). However, we can remove (0, 0, 0)-(1, 0, 0) from the approximation target because their optimal output shares are obvious. Quadratic programming calculates the optimal output shares for each of the UC candidates (0, 1, 1)-(1, 1, 1), and they are independently approximated by the nonlinear regression. Figure 3 illustrates the fuel cost functions of the aggregated operating CGs. As a result, we can obtain the coefficients A u t , B u t , and C u t as summarized in Table 2. Table 1. Specifications of CGs ("¥" means "any currency unit is applicable").
where ℎ is the sum of outputs of the operating CGs (= ∑ , , ); ℎ is the sum of operating CG outputs calculated according to ; , , and are the coefficients of approximated fuel cost associated with the vector . Table 1 summarizes the specifications of sample CGs, and Figure 2 displays their fuel cost functions. These were formulated by referring to [18,20,25,30]. When we focus on microgrids with three CGs, the possible UC candidates are (0,0,0) , (0,0,1) , (0,1,0) , (1,0,0), (0,1,1), (1,0,1), (1,1,0), and (1,1,1). However, we can remove (0,0,0)- (1,0,0) from the approximation target because their optimal output shares are obvious. Quadratic programming calculates the optimal output shares for each of the UC candidates (0,1,1)-(1,1,1), and they are independently approximated by the nonlinear regression. Figure 3 illustrates the fuel cost functions of the aggregated operating CGs. As a result, we can obtain the coefficients , , and as summarized in Table 2.  Figure 3 ("¥" means "any currency unit is applicable").     Figure 3 ("¥" means "any currency unit is applicable"). The proposed approximation can be completed as a preprocessing of the operation scheduling problem using readily available information, but it is applicable only when the microgrid does not have so many CGs. For example, with 10 CGs, the number of possible UC candidates is 2 10 (= 1024). In addition, its results are strongly dependent on the specifications of the CGs. However, there was no significant difference in the precalculation results of the widely used power grid models, e.g., IEEE 118 bus model [31,32]. Through simple verification, the authors concluded that this approximation is normally applicable to microgrids.

Reformulation of Objective Function
Before reformulating the objective function, the authors introduced a new variable ε t and replaced the stochastic variables as Then, based on Equations (18)-(21), function (7) is reformulated as If If In the case where we set s t and e t by any means, the probability density function p ε † t in Equation (22) is a horizontally shifted function of p d † t in (7). In the original formulation, the function ( f 1 + f 3 ) is determined by the set of u i,t , g i,t , s j,t , e t , and it is impossible to specify its shape in the range of ε min t , ε max t . Therefore, the integral calculus was replaced with an iterative calculation using discrete probability distributions. On the other hand, in Equation (22), the function f 1 + f 3 becomes a set of continuous function depending on u i,t , s j,t , e t once the process of 3.1 is completed. As a result, we can directly treat the calculation of Equation (22) as an integral calculus.

Application of Particle Swarm Optimization
By applying the process described in Sections 3.1 and 3.2, the CG outputs g i,t (or h t ) are removed from the optimization variables. However, the target problem is still extremely difficult to solve exactly. In this paper, a PSO was selected as the basis of the solution method. The standard PSO is the typical population-based stochastic algorithm that iteratively searches the solution space while improving a given measure of solution quality [33], which is the fitness function. An initial set of randomly generated solutions, which is the initial swarm, propagates in the search space towards the globally optimal solution until a termination condition is met. During the process, all members of the swarm fit and share the information of the search space. Each particle k (k = 1, 2, . . . , K) has a position x n k and a velocity v n k in the iteration n (n = 1, 2, . . . , N), and these are updated as where ω is the inertia weight factor that controls the iteration size; θ 1 and θ 2 are the cognitive factors that represent the trust for each particle and the swarm; r 1 and r 2 are the random numbers in the range of [0, 1]; x * k is the personal best for the particle k (pbest); min k x * k is the best in the swarm (gbest). That is, the particles fly through the search space based on Equations (26) and (27) to find the globally best location.
According to the problem formulation, the fitness function is represented by Equation (22), and the position of the particles is defined as x n k = (u, s, e) n k , for ∀k, ∀n.
The standard PSO is successful in solving many continuous problems; however, there are still difficulties in the treatment of discrete optimization problems [34]. Therefore, a strategy of binary particle swarm optimization (BPSO) was applied to the proposed solution method as shown below: If 0.5 < 1 1 + e −u i,t then u i,t = 1, else u i,t = 0, for ∀i, ∀t. (29) In (29), the CG i becomes ON at time t if the value of the sigmoid limiting transformation function is greater than or equal to 0.5, otherwise it becomes OFF.

Numerical Simulations
With a view to verifying the validity of the authors' proposal, numerical simulations were carried out on the microgrid model illustrated in Figure 1. The microgrid components were set as three CGs, two aggregated ESSs (one aggregated stationary ESS and one aggregated CL), one aggregated electrical load, and one aggregated VREG. In the numerical simulations, photovoltaic generation systems (PVs) were selected as the VREG. As for the probability density function, the authors used a Laplace distribution whose location parameter was set to the predicted net load, and the scale parameter was set to 1 [35]. The specifications of the CGs are already summarized in Table 1 and Figure 2. Table 3 shows the specifications of the aggregated ESSs. Profiles of the predicted net load and the price for traded electricity are displayed in Figure 4. These were made with reference to [20,25,30].   The time interval ∆t was set to 1 h, and the daily operation schedule (t = 1, 2, . . . , 24) was set as the optimization target. Here, the target period was set from 8:00 (t = 1) to 7.00 of the next day (t = 24) considering operation of the CLs. The available duration of the aggregated CL was set from 21:00 (t = 14) to 7:00 (t = 24). The initial SOC level of the aggregated ESS was set to 50% of its capacity, and it had to be returned to the original level until the end of the scheduling period. On the other hand, the initial SOC level of the aggregated CL was set to 50%, and it had to be 100% at the end of the scheduling period. In other words, the optimal operation schedule must satisfy the following additional constraints: According to the preliminary trial and errors, the parameters of the PSO were set as the following: the total number of particles K was 40; the maximum iteration N was 300; the inertia weight factor was 0.9, and the cognitive factors θ 1 and θ 2 were 2.0 in both.
Under these conditions, the authors determined the optimal operation schedule of the microgrid. For comparison, the traditional operation schedule was also determined under the assumption that the microgrid operators trusted the predicted data as the basis of the operation scheduling. In the traditional problem framework, the operational cost based on the predicted net load was used as the objective function. Furthermore, the authors assumed that the reserve power must compensate for a deviation within 5% of the predicted net load because the traditional framework requires a constraint for the operational margin.  Table 4.  The time interval ∆ was set to 1 h, and the daily operation schedule ( = 1,2, … ,24) was set as the optimization target. Here, the target period was set from 8:00 ( = 1) to 7.00 of the next day ( = 24) considering operation of the CLs. The available duration of the aggregated CL was set from 21:00 ( = 14) to 7:00 ( = 24). The initial SOC level of the aggregated ESS was set to 50% of its capacity, and it had to be returned to the original level until the end of the scheduling period. On the other hand, the initial SOC level of the aggregated CL was set to 50%, and it had to be 100% at the end of the scheduling period. In other words, the optimal operation schedule must satisfy the following additional constraints:

Numerical Simulation Results
According to the preliminary trial and errors, the parameters of the PSO were set as the following: the total number of particles was 40; the maximum iteration was 300; the inertia weight factor was 0.9, and the cognitive factors and were 2.0 in both. Under these conditions, the authors determined the optimal operation schedule of the microgrid. For comparison, the traditional operation schedule was also determined under the assumption that the microgrid operators trusted the predicted data as the basis of the operation scheduling. In the traditional problem framework, the operational cost based on the predicted net load was used as the objective function. Furthermore, the authors assumed that the reserve power must compensate for a deviation within 5% of the predicted net load because the traditional framework requires a constraint for the operational margin.  Table 4.  As shown in both Figures 5 and 6, the balance of power supply and demand on the predicted net load was maintained by the coordinated operation of the CGs, the aggregated ESSs, and the electricity trade. As shown in Table 4, if the accuracy of the net load prediction is sufficiently high, the optimal schedule according to the traditional problem framework is better than that of the proposed one. However, this condition is impractical in actual microgrid operations. In contrast, the proposed framework obtained a better solution from the viewpoint of the expected operational cost. Moreover, the difference between the cost using the predicted net load and the expected cost was remarkably small compared to that in the traditional one. This indicates that the optimal operation schedule of the proposed framework generated resilience in the microgrid operations.

Numerical Simulation Results
These results were obtained within a few minutes on Intel Core i7-8700K @ 3.70 GHz. Therefore, the authors concluded that the proposed solution method is applicable to actual operation scheduling.

Discussion on Numerical Simulation Results
As described in Section 2.2, the necessary reserve power was automatically calculated in the proposed problem framework. Figure 7 shows the differences between the operational margins in Figures 5 and 6. Figure 8 illustrates the transitions of the objective functions.  As shown in both Figures 5 and 6, the balance of power supply and demand on the predicted net load was maintained by the coordinated operation of the CGs, the aggregated ESSs, and the electricity trade. As shown in Table 4, if the accuracy of the net load prediction is sufficiently high, the optimal schedule according to the traditional problem framework is better than that of the proposed one. However, this condition is impractical in actual microgrid operations. In contrast, the proposed framework obtained a better solution from the viewpoint of the expected operational cost. Moreover, the difference between the cost using the predicted net load and the expected cost was remarkably small compared to that in the traditional one. This indicates that the optimal operation schedule of the proposed framework generated resilience in the microgrid operations.
These results were obtained within a few minutes on Intel Core i7-8700K @ 3.70 GHz. Therefore, the authors concluded that the proposed solution method is applicable to actual operation scheduling.

Discussion on Numerical Simulation Results
As described in Section 2.2, the necessary reserve power was automatically calculated in the proposed problem framework. Figure 7 shows the differences between the operational margins in Figures 5 and 6. Figure 8 illustrates the transitions of the objective functions.  (a) (b) Figure 8. Comparison of costs in optimal solutions: (a) Proposed framework; (b) Traditional framework ("¥" means "any currency unit is applicable").
In Figure 7, the proposed framework enhanced the upper-side margin remarkably from 13:00 to 15:00 and after 22:00. With reference to Figure 7, it can be confirmed that neighboring periods include the peak of the PV output or rapid change in the net load. In In Figure 7, the proposed framework enhanced the upper-side margin remarkably from 13:00 to 15:00 and after 22:00. With reference to Figure 7, it can be confirmed that neighboring periods include the peak of the PV output or rapid change in the net load. In general, these conditions lead to difficulty in balancing the operation of power supply and demand. As a result, as shown in Figure 8b, the expected cost was large in comparison with the cost using the predicted net load during the periods. Besides, the cost transitions in Figure 8a were similar, while the transitions in Figure 8b were different. Figure 8b shows that the transition of the expected cost was usually larger than that of the cost using the predicted net load. These results show that the necessary reserve power was automatically secured (in Figure 8), corresponding to the risk of increases in operational costs. Therefore, we can conclude that our proposal functioned well.

Conclusions
The authors proposed a problem framework and a solution method that finds the optimal coordinated operation schedule of the microgrid components. In the problem formulation, the objective function was represented as the expected cost of the microgrid operations. Since the operational margin has a large impact on the expected cost, the necessary reserve power is automatically secured in the obtained operation schedule. As a result, we can remove the traditional constraints for the balance of power supply and demand and the operational margin from the problem formulation. In the numerical simulation results, it was confirmed that the proposed objective function made it possible to evaluate the uncertainty in the prediction of the VREG output as the potential risk of cost increases. Meanwhile, the proposed solution method enabled us to treat the trend in VREG output prediction as a continuous probability density function. This indicates that there is no need to be concerned about the number and the variety of samples of the past VREG output in the proposed solution method. Through the numerical simulations and discussion of their results, we were able to conclude that the proposed strategy did not have any significant influence on the process of the solution evaluation and it provided a useful operation schedule.
In future work, the authors will improve the solution method, including a discussion on the appropriate selection of its basis. In addition, the influence ofapproximating the CG outputs will be analyzed in more detail.