Optimal Sizing of Hybrid Microgrid in a Remote Island Considering Advanced Direct Load Control for Demand Response and Low Carbon Emission

Optimal sizing of the power system can drastically reduce the total cost, which is challenging due to the fluctuation in output power of RE (primarily wind and solar) and pollution from thermal generators. The main purpose of this study is to cope with this output power uncertainty of renewables by considering ADLC, residential PV, and BESS at the lowest cost and with the least amount of carbon emission, while putting less burden on consumers by minimizing the IL. This paper optimizes the cost and carbon emission function of a hybrid energy system comprising PV, WG, BESS, and DG at Aguni Island, Japan, using a multi-objective optimization model. To solve the proposed problem in the presence of ADLC, the e-constraint method and MILP are utilized. After obtaining all possible solutions, the FSM selects the best possible solution among all solutions. The result shows that while case 1 has a lower energy cost than the other cases, the quantity of IL is quite significant, putting customers in a burden. In case 2 and case 3, the total energy cost is 11.23% and 10% higher than case 1, respectively, but the sum of the IL is 99% and 95.96% lower than case 1 as the ADLC is applied only for the consumers who have residential PV and BESS, which can reflect the importance of residential PV and BESS. The total cost of case 3 is 1.72% lower than case 2, but IL is higher because sometimes home PV power will be used to charge the home BESS.


Introduction
Aguni Island is a remote island of Okinawa prefecture located 60 kilometres northwest of the Naha district on Okinawa Island in the East China Sea. It covers 7.64 km 2 and has a population of roughly 800 people. There is one bar, one cop, and a few restaurants in the area. Aside from the hotel, there are around ten guest homes that cater to scuba divers, and guests [1]. There is no power supply from the main island, and there is less possibility of it in the future. The power generation system of this island is the internal combustion of heavy oil, and the total rated output power is 1600 kW [2]. The cost of energy production in this way is very costly as Japan needs to import 99.7% oil from abroad [3]. On the other hand, greenhouse gas emissions are rising significantly due to massive electricity production from fossil fuels, such as oil, coal, and gas, which tend to change the climate, security, and price of natural resources. In FY2016, total greenhouse gas (GHG) emissions (preliminary figures) in Japan were 1322 Mt CO 2 eq. (0.2%, 6.2%, and 4.6% decreased as compared to FY2015, FY2013, and FY2005, respectively). The main reason for the lower emissions in • K-means clustering is used to select the representative days to represent the whole year data by representative days that can minimize the computational burden of simulating one-year data; • -constraint method with mixed-integer linear programming (MILP) is applied to solve the proposed multi-objective optimization problem and obtain multiple possible solutions. The fuzzy satisfying method (FSM) is then used to select the appropriate solution; • DR program is implemented to flatten the load curve for solving the important duck curve problem. Residential PV and BESS are introduced to minimize the amount of interruptible load (IL), making the system more reliable.
The rest of the paper is organized as follows: The literature review is described in Section 2. The proposed power system model and advanced direct load control are described in Section 3. The data selection method is discussed in Section 4. The objective function and constraints are explained in Section 5. The techniques of solving the multiobjective optimization are described in Section 6. In Section 7, the results of 3 different case studies (without home PV and BESS, with home PV and with home PV and battery) are compared to show the effectiveness of the proposed method. Finally, the conclusion and future modification of this paper are presented in Section 8.

Literature Review
Some other existing works on HRES are also reviewed for better understanding. In [21], the greenhouse gas emissions from current fossil fuel power facilities of Bangladesh is analyzed using hybrid optimization of multiple energy resources (HOMER). The results demonstrate that coal, diesel, and natural gas power plants release 0.90 kg, 0.76 kg, and 0.566 kg of CO2 per kWh, respectively, responsible for climate change.
In [22], the aim was to discover the optimal size of a hybrid energy storage system, consisting of a hydrogen fuel cell and a supercapacitor for a commercial load provided by solar panels. To examine the influence of the hydrogen cost on the cost of the system and the levelized cost of energy, a sensitivity analysis on the estimated costs of hydrogen storage is performed using HOMER Pro under Cape Town weather conditions. Although the cost of such a hybrid storage system is likely to fall in the coming years, it will still be too expensive to deploy for a commercial load. In [23], HOMER software generates three optimal configurations of hybrid renewable energy systems to meet the residential and agricultural power demand requirements of an energy-deprived village in India, with wind, PV, and battery-based HRES being the most cost-effective design for this specific area. The load is forecasted for a remote district in India, and HOMER software is used to optimize the design and conduct a techno-economic analysis of the proposed PV, WG, and bio-generator based HRES system in [24], which is a more cost-effective system than the conventional one. In [25], different combinations of the off-grid hybrid energy system (HES) are optimized for a rural hilly area of Bangladesh, considering minimizing the cost of energy, net present cost, and CO2 emissions using HOMER software where PV-Diesel-PHS based system is more cost-effective. In [26], HOMER software is used to optimize hybrid PV-DG-battery based systems for electrifying the rural area of Benin. According to the findings, the suggested method is more successful than the current approach in terms of lowering energy costs and carbon emissions.
In [27], intelligent flower pollination algorithm is used for the optimal design and energy management of the hybrid systems based on hydrogen storage, including PV, WG, and FC, to reduce the total net present cost of the northwest region of Iran that finds the optimal decision variables at the minimum cost, and with better reliability values in various reliability indices. In these papers, a demand response program was not considered, which may increase the cost of energy.
In [28], genetic algorithm (GA) and HOMER Pro Software are used to minimize the total system net present cost, cost of energy, unmet load, CO2 emissions of an off-grid HRES for supplying electricity to a group of three villages in India. A genetic algorithm is used in [29] to optimize the solar, wind, and storage based off-grid hybrid system for the remote island. In both papers, and the obtained results are compared with the HOMAR software, and the optimal result of GA is more cost-effective than HOMAR. For minimizing the total/net present cost of the HES encompassing solar PV, wind, diesel generator, and battery for electrifying rural areas in stand-alone applications, a suitable improved GA program has been developed in [30] utilizing the MATLAB toolbox. The size of gridconnected solar and battery systems for residential houses is optimized using a genetic algorithm in [31] to reduce the overall yearly cost of electricity. In this paper, 0% leakage and complete charge/discharge capabilities of the battery are considered, which is not realistic. In [32], GA is used to optimize the PV, wind, battery, and diesel-based hybrid system to satisfy the electricity demands of a remote village of northern Nigeria at the lowest possible cost and with the least amount of carbon emissions. The computational cost of GA is expensive as it takes a long time for convergence. In [33], the NSGA-II method is used for optimal sizing of a HES that includes PV, WG, BESS, combined cooling, heating, and power generation system, heat storage tank, gas boiler, and electric chiller to decrease economic and environmental consequences. A one year dataset is used for simulation, which makes slow convergence of the optimal result.
An improved multi-objective grey wolf optimizer is applied in [34] to minimize the annualized cost of the system and deficiency of power supply probability by determining the optimal size of a hybrid microgrid consisting of PV, WG, tidal current, battery, and diesel for an island. In this research, only the initial cost and running cost of installed equipment is considered that cannot be able to reflect the appropriate cost of energy.
A multi-objective mathematical model is established in [35] for optimizing the capacity of hybrid energy storage system for a grid-connected 99MW Caka wind farm in Qinghai Province, China, to maximize target satisfaction rate and minimize net present value. The only wind farm is considered with a storage system which is less efficient than multiple energy sources.
In [36], the Grasshopper Optimization Algorithm (GOA) is used to determine the optimal system configuration of an autonomous microgrid system that includes PV, WG, BESS, and DG to fulfil energy demand reliably and cost of energy. This system can fulfil the energy demand for only five residential houses.
A Firefly algorithm is applied in [37] to determine the optimal size of a solar, wind and battery storage based hybrid energy system considering the minimum cost of energy for electrifying remote villages in India. In this paper, only one day of summer and winter data are used for simulation, which is not enough for the whole year's representation.
In [38], a multi-objective crow search algorithm is proposed to reduce the total net present cost and loss of power supply probability of an off-grid hybrid system consisting of PV, FC, and DG to supply electric power in the south of Iran. According to the simulation result, the total cost of the HES can be reduced by the integration of hydrogen energy technology. A geographical information system module is used in [39], to choose the best site based on different factors. A hybrid optimization technique is then used to estimate the optimal capacity to satisfy the demand at the least cost, which is more accurate than other algorithms. One year of data have been used for simulation, which increased the computational burden.
RETScreen simulation software is deployed in [40] to estimate the cost and pollutant emission parameters of a PV, biomass, and additional storage based off-grid hybrid system in the remote areas of Ashuganj, Bangladesh, which is more reliable than the conventional kerosene-based system. No optimization method is used in this paper.
In [41], MILP optimization algorithm has been designed for a case study of a mountain hut in South Tyrol (Italy) with solar, wind, DG, and battery storage based HRES to design a tool capable of determining the optimal sizing of an HRES that can assist engineers in identifying the best trade-off between costs and energy deficiency during the planning stage. The reduction in carbon emission is not considered here.
For the Gobi Desert in China, the -constraint technique and elephant herd optimization algorithm are utilized in [42] to reduce the loss of load probability, CO2 emissions, and yearly cost of a solar, diesel, and battery-based hybrid system. According to the results, the suggested system emits less carbon than the PSO and HOMER-based systems. A significant portion of the energy generated by the PV throughout the day is lost for the limitation of battery capacity.
For optimization and sensitivity analysis of an autonomous HES consisting of PVdiesel-battery for a remote Saharan village in southern Algeria, particle swarm optimization and the -constraint method were proposed in [43] to reduce total system cost, unmet load, and CO2 emissions, which is more cost-effective compared to HOMER.
MILP and -constraint method were used to generate energy storage system dayahead scheduling in the context of wind farm uncertainty. According to the simulation result, the use of an energy storage unit lowers daily costs and emissions. Only one day of data have been used for simulation, which is insufficient for long term planning for optimal scheduling of energy storage system [44].
A multi-objective optimal scheduling model for CCHP microgrids integrated with RE, energy storage system, and incentive-based demand response is solved using MILP and augmented-constraint methods by minimizing pollutant gas emissions and lowering costs. By adjusting the peak of the exchange power curve, the IL and the battery may effectively adapt to peak load fluctuations. Residential PV and BESS can reduce the IL, which is not considered in this paper [45].
Among various optimization methods, the -constraint technique minimizes the computational cost of the system and becomes effective when the limits of objective functions are known. This approach is very good for finding convex and linear Pareto optimal front as it has the ability to find accurate Pareto front, rather than approximated solutions [46].
According to the limitation of above literature, the contributions of this study are as follows: • The advantages of minimizing cost-emissions function, as well as household PV and BESS, are examined, consider carbon emission as a constraint; • Advanced direct load control (ADLC) is used to flatten the load curve. This ADLC is exclusively applied to customers who possess residential PV and BESS in order to reduce the amount of IL on the system, making it more reliable by lowering the impact of power outages.

Power System Model
A decentralized power system model is designed and shown in Figure 1 in this paper for a remote island in Okinawa Prefecture named Aguni island with a peak load of 1 MW. It is considered that the power system is disconnected from the Okinawa main power system. As a hybrid power system, a small PV, WG, DG, and BESS are considered as the power generation source and energy storage. PV and WG power is used to meet the demand, and the battery is charged by the remaining power. If PV and WG power is not enough, BESS will be discharged. If the demand is not compensated by the PV, WG with BESS, DG will run, and DR will be applied.

Modeling of Installed Equipments
The load data used in this paper are collected from the Okinawa Electric Power website. Solar irradiation and wind speed data of one year are used to calculate the output power of PV and WG, respectively, which is taken from Japan Meteorological Agency. The output power generated by PV is calculated by Equation (1).
where E PV is the power produced by PV (kW), H is the quantity of solar irradiation (kW/m 2 ), A PV is the total solar panel area (m 2 ), and η PV is the solar panel yield or efficiency of the PV system which is considered as 12.3% in this paper [47]. The output power generated by WG is represented by Equations (2) and (3).
where E WG is the generated power by WG (kW), C p is the output power coefficient of performance of WG, ρ is the air density (kg/m 3 ), A WG is the swept area of blades of the WG (m 2 ). V hub is the speed of wind (m/s) at the target height of the hub of the WG, V re f is the speed of wind (m/s) at the reference height of the WG. Z hub and Z re f is the target hub height and reference hub height of the WG, respectively. The cost of installed equipment and other parameters for calculating the generated power of PV, WG, BESS and DG is listed in Table 1.

Modeling of Advanced Direct Load Control
DLC is a contract based DR signed by the consumers, which is provided by utilities. According to this contract between consumers and utilities, consumers permit utilities to control consumers' air conditioners and water heaters remotely. Utilities can shut down such appliances during peak-demand periods or power supply shortages or notify consumers about peak periods to shut down appliances. The consumers who participate in this contract obtain compensation by decreasing their electricity bill [48,49].
ADLC, an incentive-based DR, has been introduced in this paper to bridge the gap between supply and demand curves. We estimate that this contract will be made by 0%-50% of consumers. The consumers following the agreement accede to the power cut while the shortage of power. Each contracted household will face a maximum of 2 h of power cut per day, and for 1 kW of a power cut, they will receive JPY 10 as compensation.
where IL max is the maximum amount of IL per hour. • Compensation cost function in ADLC: In ADLC, the compensation cost is calculated by Equation (5) comp cost = n ∑ t=1 IL t (kW) * PR comp (Yen/kW) where comp cost is the compensation cost, IL t is the IL in hour t and PR comp is the compensation price per kW which is considered as JPY 10 here.
The advantages of this contract are listed below: • Electricity bill saving for consumers: For cutting of 1 kW of electricity from a household each time, it will receive a discount of JPY 10 from the electricity bill; • Lower electricity bill for other consumers: For using less energy by the contracted consumers, the wholesale market price of energy will not increase at the time of shortage; • Minimize the duck curve problem: Duck curve occurs during peak load with low generation. If the load is minimized by cutting off power, the duck curve will also be minimized.

Proposed Data Selection Method for Optimal Scheduling
In order to finalize the optimal power supply configuration with the best operation plan of the annual data, a large number of calculations are required. These data should be simplified to reduce the computational burden. The calculation cost can be minimized by choosing some qualitative days as representative days. K-means clustering is used to simplify the yearly data in this research, which is a simple unsupervised machine learning algorithm that classifies similar data into a number (k) of clusters. For k-means clustering, the optimal number of clusters is selected by the elbow method in this paper. K-means clustering is run by the elbow method based on the dataset for a range of values for k, and an average score is calculated for all clusters for each value of k. The distortion score is calculated by default from the sum of square distances between each point and its assigned center. After plotting these overall metrics, it is possible to visualize the best value for k. If the line graph looks like an arm, then the elbow point is the best value of k that can be either up or down [50]. In the power system model, all historical days which are classified into the same cluster will be illustrated by the same representative day.
The simplified working procedure of the elbow method for determining representative days for optimal operation in this paper is shown in Figure 2. The annual net demand dataset of clustering in this paper is shown in Figure 3, which depends on load demand, solar radiation and wind speed. Figure 4 shows the elbow diagram where number 18 is considered as the elbow point and the best value of k. Figure 5 shows the net load demand of the obtained dataset of the number of clusters [51]. The obtained dataset, actual load data, power generated from WG and PV is shown in Figure 6, ensuring that proper data selection has been accomplished here [52] which can be able to represent the whole year.

Time (hour)
Net load demand (kW)

Problem Formulation
In the assumed power system model of this paper, the number of installed WG, PV, and BESS can vary. The output power is proportional to the rising number of installed WG and PV. The minimization of the LCC of the system and carbon emission are the objective functions of this paper. The number WG, PV, and BESS, schedule of charging or discharging of BESS, and schedule of DG operation are determined by MILP so that the objective functions can be minimized.

Objective Function i
Cost function: In this paper, multi-objective optimization is considered with two objective function of minimization. The first objective function is to minimize the LCC which is defined by Equation (6).
where C F is the initial/fixed investment cost (JPY), Y is the total lifetime (Year) of the project which is considered 15 years in this paper. C O&M (y) operation and maintenance cost and C RP (y) is the replacement cost in the year y (JPY), respectively. D is the discount rate which is assumed to be 3% in this paper. RV is the residual value. The initial investment cost C F is the sum of the initial investment cost of installed equipment which is calculated by Equation (7) C where C F PV , C F WG and C F BESS are the fixed investment cost of PV, WG, and BESS, respectively. Similarly, the operation and maintenance cost C O&M (y) and replacement cost C RP (y) is calculated by Equations (8) and (9), respectively.
C O&M (y) = C O&M PV (y) + C O&M WG (y) + C O&M BESS (y) (8) where C O&M PV (y), C O&M WG (y), and C O&M BESS (y) are the the operation and maintenance cost of PV, WG, and BESS in the year y, respectively.
C RP (y) = C RP PV (y) + C RP WG (y) + C RP BESS (y) (9) where, C RP PV (y), C RP WG (y), and C RP BESS (y) are the the replacement cost of PV, WG, and BESS in the year y, respectively. ii Emission function: The second objective function is to minimize the carbon emission which can be calculated by Equation (10).
Min : where H c is the heating value of C fuel oil (GJ/kL), CF is the CO 2 conversion factor of CO 2 per calorific value (tCO 2 /GJ) and C coil is the consumption of C heavy oil that can be calculated by Equation (11).
where, H is the heating value (kJ/kW) and E DG is the total power generated by DG. The fuel cost of DG is shown by Equation (12).
E (DG i ) is the power generated by DG i at each time t (kW), P c is the rate of C fuel oil (JPY), and η i is the efficiency of DG.

Constraints
i. Constraints for charging and discharging of BESS are given below: where, B c and B d are the charging and discharging state of BESS, respectively. ii.
Maximum output power of BESS iii. Start/stop constraint of DG iv.
Power balance limit where, n B is the number of BESS; C B is the per unit capacity of BESS; η cB , η dB are charging and discharging efficiency which have taken by 80% in this paper; n PV and n WG are the number of installed PV and WG. E PV and E WG , are the power generated prom PV and WG, respectively. E d and E c are the discharging and charging power of BESS. E max d , E max c are the maximum discharging and charging power of BESS. E DG are generated power from DG; E sr is the surplus power; SOC i is the initial state of the SOC of BESS; E IL is the amount of IL; E L is the load demand. v.

-Constraint Method
The -constraint method is used to solve the proposed multi-objective optimization. One objective function is considered for optimization, and the rest of the objective functions are converted to inequality functions and considered as constraints. In this research, one objective function, φ 1 , is optimized, and another objective function, φ 2 , is considered as a constraint as follows.
Equations (21) and (22) describe the -constraint method where φ 2 is restricted by parameter. The minimum and maximum values of this parameter are set very carefully. Gradually the parameter varies from the minimum value to maximum value, and the solution of the modified single objective function is obtained for each parameter. The set of all obtained solutions from the variations of min to max are the Pareto optimal front of the multi-objective optimization problem.

Fuzzy Satisfying Method
To find the best possible solution from the Pareto optimal front of the multi-objective optimization problem, FSM is necessary. A fuzzy membership function with the interval (0, 1) is assigned to each solution in the Pareto front. The linear fuzzy membership functions can be obtained n th solution of i th objective function is defined as Equation (23) where φ n i shows the optimality degree of the n th solution of i th objective function.  (24) and (25).
The minimum value between φ n 1 and φ n 2 is determined. This value is called the membership function of the n th solution which can be calculated by Equation (26).
The solution with maximum lowest membership function can be selected as the best possible solution which is calculated by Equation (27).
The above solution method can be briefly explained by following steps

•
Step-1: Minimize φ 1 and φ 2 subject to all equal and unequal constraints using MILP; • Step-2: Calculate the maximum value of φ 1 and φ 2 ; • Step-3: A number of iterations are considered where in the first iteration, φ 2 is minimum and φ 1 is maximum. With increasing the number of iteration, φ 2 will increase, and in the last iteration, φ 2 is the maximum, and φ 1 is minimum; • Step-4: To select the best possible solution FSM is used. The solutions of the objective functions are converted to per unit values using Equations (24) and (25). Then the membership function and best possible solution can be calculated by Equations (26) and (27), respectively.

Result Analysis
Using one-year data for simulation is time-consuming. Selecting random days for simulation could not compute the optimal sizing of RE correctly because of the variation of output power and load demand. The sizing of the HRES could not meet the needs of the load demand if the proper representative days of the entire year were not selected. That is why K-means clustering is used to simplify the yearly data to minimize the computational burden.
In this work, a comparison of 3 case studies is considered to reflect the usefulness of the proposed method.
Case 1: Energy optimization with ADLC-In this case, if there is a shortage of power, cutting off the power supply will occur for all consumers who have a contract with ADLC.
Case 2: Energy optimization with ADLC and residential PV-In case 2, residential PV is considered for some houses. Since these PVs are not connected to the grid, reverse energy will not appear to the grid. Only the consumers who have the residential PV in their home would have the contract of ADLC and face the cut off power supply.
Case 3: Energy optimization with ADLC, residential PV and BESS-In this case, it is considered that some houses have residential PV and BESS. The consumer who owns the residential PV and BESS would have the contract of ADLC and face the cut off power supply.
The results of determining the optimal configuration of these cases are calculated by the selected data based on Section 3. 11 different solutions for each case with the combination of PV, WG, BESS, and DG are obtained in this paper with the variation of epsilon. The LCC, the fuel cost of DG and the carbon emission of each case is listed in, Table 2 shows that the LCC has decreased with the increase in carbon emission. The obtained Pareto set of each case is shown in Figure 7. Here, the best possible solution of each case is selected by the FSM discussed in Section 5 where the maximum value of constraint is set to 30% of the maximum value of .
From the Pareto set, the optimal configuration in case 1 can be achieved at the lowest cost, and in case 2 the LCC is the highest among these three cases.
The result of each case is listed in Table 3 where the maximum value of -constraint is set to 0% and 30% of the maximum value of .  Although the LCC, as well as total cost (TC), is low in case 1 in Table 3, the amount of IL is higher than in other cases. Additionally, in this case, contracted consumers cannot benefit from electricity during the cut off electricity because there is no source of electricity generation in the houses. At some point, it will become a burden on consumers.
In case 2, the LCC and TC are higher than in other cases, but IL is 99% low here. As only the houses that own residential PVs will face the cut off power during the shortage; they can use the power generated from PVs while the power interruption. This case can reduce the burden of consumers, but the duck curve phenomenon is occurred by PV generation in the daytime could not be solved.
In case 3, the LCC and TC are higher than case 1 and lower than case 2, and the IL is 95.96% lower than case 1 and higher than case 2. In this case, only the houses that own residential PVs and BESSs will face the cut off power during the shortage. They can use the power generated from PVs and battery storage while the time of power cut off. As there are residential BESS, sometimes the power generated from PVs is used for charging home BESSs. That is why the amount of IL is higher than case 2 here. However, it can slightly solve the duck curve problem. Figures 8-10 show the optimal operation of 2 days at 30% carbon emission, respectively.

Conclusions
In this paper, a multi-objective optimization model was proposed for optimal sizing of PV, WG, BESS, and DG based HES considering minimum cost and carbon emission utilizing ADLC, residential PV, and BESS program to electrify the remote Aguni Island in Japan. For inspiring consumers to install the residential PV and BESS, ADLC is modeled in a way that only the consumers with residential PV and BESS would have the contract of ADLC, which can significantly reduce the amount of IL. The proposed multi-objective model is solved by the -constraint method. MILP is applied to obtain all possible solutions of this model, and FSM is responsible for selecting the best solution. Three different case studies are considered here to validate the effectiveness of the DR program with residential PV and BESS. With the variation of epsilon, 11 different results are obtained for each case. FSM has selected the best solution where the value of epsilon is 30%. The best solution illustrates that:

•
In case 1, only ADLC is considered with no home PV and BESS. In this case, the energy demand is fulfilled at the lowest cost, but the amount of IL is very high, which will be a burden for consumers as they cannot be able to use electricity at that time; • In case 2, ADLC and home PV is considered where the total energy cost is 11.23% higher than case 1, but the amount of IL is almost 99% lower than case 1; • In case 3, ADLC, home PV and home BESS is considered where the total cost is 10% higher than case 1 and the amount of IL is 95.96% lower than case 1.
This paper does not take into account real-time uncertainty in load demand. Specifically, the proposed method considers 18-day sample data from 1 year of historical data. Calculating compensation for ADLC on historical data might have minor bugs due to real-time load uncertainty.