Co-Planning of Demand Response and Distributed Generators in an Active Distribution Network

: The integration of renewables is fast-growing, in light of smart grid technology development. As a result, the uncertain nature of renewables and load demand poses signiﬁcant technical challenges to distribution network (DN) daily operation. To alleviate such issues, price-sensitive demand response and distributed generators can be coordinated to accommodate the renewable energy. However, the investment cost for demand response facilities, i.e., load control switch and advanced metering infrastructure, cannot be ignored, especially when the responsive demand is large. In this paper, an optimal coordinated investment for distributed generator and demand response facilities is proposed, based on a linearized, price-elastic demand response model. To hedge against the uncertainties of renewables and load demand, a two-stage robust investment scheme is proposed, where the investment decisions are optimized in the ﬁrst stage, and the demand response participation with the coordination of distributed generators is adjusted in the second stage. Simulations on the modiﬁed IEEE 33-node and 123-node DN demonstrate the effectiveness of the proposed model.


Introduction
Recently, the penetration of renewables increases constantly in the modern distribution network (DN), due to their worldwide availability and sustainability [1][2][3]. In addition to the benefits of conventional distributed generations, such as energy loss reduction and infrastructure upgrade deferment [4], renewable energy sources also play significant roles in carbon emission reduction, energy conservation, and the promotion of fossil fuel alternatives. Nevertheless, due to their intermittent and variable nature, renewables, such as wind power and solar energy, are known as non-dispatchable and uncontrollable sources. The expanding integration of renewables will introduce tremendous difficulties in order to maintain power system reliability [5][6][7]. The bidirectional power flow that emerges with renewables installation leads to either a rise or a drop in voltage, depending on the output fluctuation of renewables. In addition, a large injection of renewables may lead to substations being adversely overloaded and line congestion.
In order to ease these negative impacts, dispatchable distributed generator (DG) is usually considered as an energy supplement and reserve, by reason of its controllability and small ramping limit [4]. Demand response (DR) has been commonly known as an efficient approach in peak load reduction [8]. Generally, the DR program is designed to adjust terminal customers' electricity consumption by offering financial incentives, or providing time-varying prices, in order to decrease

Demand Response Model
In the proposed demand response model, we assume the load at each DR bus is divided into inelastic and elastic parts. Considering the basic inelastic electricity consumption, only the part of demand that is flexible in electricity usage is selected as the DR program by investing DRFs. Other electricity demand, such as in schools and hospitals, is not affected by electricity price. A linearized price-elastic demand response (PEDR) model is implemented in the second stage to adjust the firststage investment decisions. The power demand is usually modeled using a price-elastic curve [13,14,22], the expression of which is given as below: where P d j,t and Prj,t denote power demand and electricity price, respectively, A is the proportion and index constant parameters given in [22], and γ is the price elasticity which describes the sensitivity of load demand as to the change of price. For a given A and γ, the DN revenue from electricity consumption is equal to the integral of the constant inelastic demand and the price-elastic curve (present as r d j,t ), as shown in Figure 2. Unfortunately, (1a) is non-convex. Thus, to remove nonlinear terms, the stepwise linearization approximation (SLA) is adopted to estimate the grid revenue contributed by DR program, which is modelled with the following auxiliary variables and constraints: where h d j,t in (3a) donates the actual demand at DR node which is composed of constant inelastic demand and elastic demand. Constraint (3b) gives the upper and lower bound of the elastic demand.

Demand Response Model
In the proposed demand response model, we assume the load at each DR bus is divided into inelastic and elastic parts. Considering the basic inelastic electricity consumption, only the part of demand that is flexible in electricity usage is selected as the DR program by investing DRFs. Other electricity demand, such as in schools and hospitals, is not affected by electricity price. A linearized price-elastic demand response (PEDR) model is implemented in the second stage to adjust the first-stage investment decisions. The power demand is usually modeled using a price-elastic curve [13,14,22], the expression of which is given as below: where P d j,t and Pr j , t denote power demand and electricity price, respectively, A is the proportion and index constant parameters given in [22], and γ is the price elasticity which describes the sensitivity of load demand as to the change of price. For a given A and γ, the DN revenue from electricity consumption is equal to the integral of the constant inelastic demand and the price-elastic curve (present as r d j,t ), as shown in Figure 2. Unfortunately, (1a) is non-convex. Thus, to remove nonlinear terms, the stepwise linearization approximation (SLA) is adopted to estimate the grid revenue contributed by DR program, which is modelled with the following auxiliary variables and constraints: where d k j,t is the auxiliary variable recommended for demand at step k, l k j,t is the kth step length of the stepwise function, and Pr k j,t is the corresponding electricity price. It is noted that Pr k j,t stringently decreases with k, so that the DR program constraints are given as where h d j,t in (3a) donates the actual demand at DR node which is composed of constant inelastic demand and elastic demand. Constraint (3b) gives the upper and lower bound of the elastic demand.  Figure 2.
Step-wise function approximation of PEDR curve.

Objective Function
The investment objective is to minimize the total cost, so as to find out the optimal allocation location of the DRF and DG. The objective function is to optimize the total costs and revenues during the whole planning horizon, which is formulated as: Overall, the function includes five parts, which are investment costs (Cinv), operation and maintenance costs for DGs (Cope), energy transaction cost (Cadn), DN revenue (Cearn), and penalty costs for operation violation (Cpena).
Herein, the total investment cost is defined in (4b), which consists of the costs of DRFs and DGs, respectively. The investment cost of the DRF defined in (4c) contains two terms. The first term is the procurement and installation cost of LCS and AMI. The second term contains the education cost and financial incentive cost paid to electricity users for DR promotion. The investment cost of DG is defined in (4d). Notice that only the procurement and installation cost of LCS and DG has a linear relationship with their capacity. The capital recovery factor ε in (4c) and (4d) is utilized to transform the total investment cost into daily cost [23]. Equations (4e)-(4l) give the detailed expression of the operation costs and revenues, which contain the operation and maintenance cost of the DG, wind turbine (WT), and photovoltaic generator(PV), respectively; the fuel cost of DG; the energy transaction cost with the main-grid; and Step-wise function approximation of PEDR curve.

Objective Function
The investment objective is to minimize the total cost, so as to find out the optimal allocation location of the DRF and DG. The objective function is to optimize the total costs and revenues during the whole planning horizon, which is formulated as: Overall, the function includes five parts, which are investment costs (C inv ), operation and maintenance costs for DGs (C ope ), energy transaction cost (C adn ), DN revenue (C earn ), and penalty costs for operation violation (C pena ).
Herein, the total investment cost is defined in (4b), which consists of the costs of DRFs and DGs, respectively. The investment cost of the DRF defined in (4c) contains two terms. The first term is the procurement and installation cost of LCS and AMI. The second term contains the education cost and financial incentive cost paid to electricity users for DR promotion. The investment cost of DG is defined in (4d). Notice that only the procurement and installation cost of LCS and DG has a linear relationship with their capacity. The capital recovery factor ε in (4c) and (4d) is utilized to transform the total investment cost into daily cost [23].
Equations (4e)-(4l) give the detailed expression of the operation costs and revenues, which contain the operation and maintenance cost of the DG, wind turbine (WT), and photovoltaic generator(PV), respectively; the fuel cost of DG; the energy transaction cost with the main-grid; and the revenue by selling electricity to customers. The operation and maintenance cost for WT and PV in (4f) are normally considered as the regular refurbishment and maintenance cost arisen from the depreciation with daily operation [24]. Additionally, the fuel cost of DG is dependent on the total energy generation, as illustrated with equation (4g).
The electricity exchange between DN and the main grid is presented as (4h). Considering the uncertain load demand, the DN revenue in (4i) is divided into three parts: the revenue from DR-enabled demand in (4j), the expected revenue from DR-not-enabled demand in (4k), and the uncertain revenue that is violated from its expected value in (4l).
In Equations (4m)-(4p), the penalty costs are expressed as constraint violations of voltage fluctuation, line congestion, and substation overload, which are computed by multiplying slackness variables with their penalty fees.

Constraint
The economic and technical constraints describing the investment and daily operation are listed in this subsection: (1) The physical and technical constraints that describe the limitation of investment decision. The limitation of maximum number of DRF and DG investment locations is given in (5a) and (5b). The limitations of DG and DR capacities under a DN node are given in (5c) and (5d). In order to reduce the non-linearity and computing burden of the problem, the bilinear terms α dr j β dr j and α dg j β dg j,b are replaced by the auxiliary variables η dr j,t and γ j,b , which are defined and limited by constraints (5e)-(5h). In addition, the capital recovery factor is defined in (5j).
(2) The physical and technical constraints describing the DR-enabled DN operation. The branch flow model (BFM) proposed in [25,26] is further simplified for distribution network modeling [27], which is described with constraints (6a)-(6c). Equation (6d) fixes the voltage at the substation as V 0 . Constraints (6e)-(6j) are the linearized PEDR model, in which the uncertainty of the load demand is considered by introducing the deviation rate R d j,t from expected value. Constraints (6k) and (6l) guarantee the energy consumption with DR investment is no smaller than original consumption, and the electricity bill for consumers with DR investment is no larger than the original bill. The constraints (6m)-(6n) express the energy exchange between DN and main grid, where P in represent power purchase from the main grid, and P out means DN sells energy to the main grid. Constraint (6o) gives the substation capacity limitation. The output limitation of DG and local reactive compensation devices are presented in constraints (6p)-(6q). Constraint (6r) is to guarantee that the power supply can satisfy the total demand with reserve factor θ. Constraint (6s) gives the voltage fluctuation limitations, and constraint (6t) guarantees the active power flow within the line capacity. It is worth noting that the operation constraints in (6o), (6s)-(6t) are relaxed, with positive variables S uv j,t , S lv j,t , S usub t , S lsub t , S ucap ii,t , and S lcap ii,t , so that the penalty cost for operation violation can be computed.

Uncertainty Set
In this paper, three polyhedral uncertainty sets are defined for the wind power, solar energy, and the load demand, respectively, given as follows, µ wt up,n,t P wt mean,n,t ≤ P wt n,t ≤ µ wt low,n,t P wt mean,n,t , n ∈ N wt , ∀t} µ pv up,n,t P pv mean,n,t ≤ P pv n,t ≤ µ pv low,n,t P pv mean,n,t , n ∈ N pv , ∀t} In (7a) and (7b), we assume that the wind power and solar energy outputs are within an interval where P wt,pv mean,n,t represents the mean value of their outputs. This interval can generated through is allowed to be any value in the given range [28].
The uncertainty set for load demand in (7c) is expressed as the deviation rate R d j,t from expected value. Likewise, the upper and lower bounds µ d low,j,t /µ d up,j,t indicate the lower and upper interval of the uncertain load demand for each single hour. Furthermore, the uncertainty budgets are enforced to all the uncertainty set in (7a)-(7c), with Γ low /Γ up to quantify temporally and spatially aggregated values. These budgets are to control the robustness of the proposed model. They can be adjusted according to the system planner's attitude towards practice application.
With the consideration of uncertain variables, the two-stage optimization objective (4a) can be rewritten as follows: min In (8), the first term is the minimum investment cost of DRF and DG, while the second term is the minimum the other costs and revenues from DN operation. The worst-case scenario is first formed in search of the uncertain variables in the uncertainty set, and then the operation variables are optimized within the operation constraints to obtain minimal corresponding costs in this scenario, which makes the "max-min" bi-level problem in (8). From this feature of robust mathematic model, the first-stage variables serve as "here and now" decisions; meanwhile, the operation variables work as adjustable "wait and see" decisions. Hence, the optimal solution is able to hedge against any possible scenario, as well as ensure every operation constraint is satisfied in all possible uncertainty case.

Compact Formulation and Duality
The compact form of the proposed model is written as below: According to their specialty, all the binary variables are grouped in x (i.e., x = {α j , β j }). Moreover, the second-stage variables are divided into operation variable y and slack variable s. Operation variables include DGs and Static Var Compensators (SVCs )outputs P dg and Q svc , adjusted DR program participation d k , and uncontrollable BFM variables P, Q, V. The uncertain variables, including wind power, PV outputs P wt,pv n,t , and deviation rate R d j,t of load demand are expressed as u. Otherwise, the constraints (5a)-(5h) and (6t) for x can be grouped into (9b), and the operation constraints (9c)-(9e) can account for (6a)-(6s) and (6u)-(6w). It is noted the I(u) and K(u) are the coefficient matrix of x, and include the uncertainty variables according to (6a)-(6q). The uncertainty sets (7a)-(7c) are sorted in (9f).
The optimization problem (9a)-(9f) can be divided into two stages. Mathematically, with a given u * generated from sub-problem, the master problem can be driven from (9) and become a mix-integer linear program (MILP), which can be directly solved by a commercial solver. Hence, an optimal solution (x * , λ * ) can be obtained, in which λ * denotes the value of second-stage objective function value. With these decision variables x * , the sub-problem becomes a "max-min" bi-level program, which is Energies 2018, 11, 354 9 of 18 noted as O(u, x * ) in this paper. Due to its strong duality, O(u, x*) is equivalent to the dual problem as follows: where {φ 1 , φ 2 , φ 3 } are the dual variable vectors for constraints (9c)-(9e).

Algorithm
The CCG algorithm is an efficient decomposition algorithm for two-stage robust optimization [29]. Thus, a two-level algorithm is proposed in this subsection, in which the CCG algorithm is adopted in the outer level and the inner level uses an OA algorithm [30] to solve the bilinear problem in (10). The procedure of the proposed algorithm is given in following steps, with a defined convergence index ε: The outer CCG algorithm is given as: (1) Initialize LB ccg = −∞, UB ccg = +∞, g = 0 and set the sub-problem solution set soa = ∅.

Algorithm
The CCG algorithm is an efficient decomposition algorithm for two-stage robust optimization [29]. Thus, a two-level algorithm is proposed in this subsection, in which the CCG algorithm is adopted in the outer level and the inner level uses an OA algorithm [30] to solve the bilinear problem in (10). The procedure of the proposed algorithm is given in following steps, with a defined convergence index ε: The outer CCG algorithm is given as: (1). Initialize LBccg = −∞, UBccg = +∞, g = 0 and set the sub-problem solution set soa = ∅.
(2). Solve the master problem, obtain the optimal solution (x (4). Solve the OA master problem, which is the linearized version of the second stage problem, defined as below:  The general framework of the algorithm is described in Figure 3.
(4) Solve the OA master problem, which is the linearized version of the second stage problem, defined as below: Energies 2017, 10, x noted as О(u, x * ) in this paper. Due to its strong duality, О(u, x*) is equivalent to the dual pr follows: where {φ1, φ2, φ3} are the dual variable vectors for constraints (9c)-(9e).

Algorithm
The CCG algorithm is an efficient decomposition algorithm for two-stage robust optimiza Thus, a two-level algorithm is proposed in this subsection, in which the CCG algorithm is in the outer level and the inner level uses an OA algorithm [30] to solve the bilinear problem The procedure of the proposed algorithm is given in following steps, with a defined convergenc The outer CCG algorithm is given as: (1). Initialize LBccg = −∞, UBccg = +∞, g = 0 and set the sub-problem solution set soa = ∅.
The general framework of the algorithm is described in Figure 3.

Case Studies
In this section, the proposed framework and methodology were tested on the modified IEEE 33node distribution system. All the computational tasks were implemented with a Yalmip toolbox in MATLAB (R2016a, MathWorks, Natick, MA, USA) on a 2.4 GHz personal computer with 8 GB of RAM, using Gurobi as the MILP solver (7.0.1, Gurobi Optimization, Inc., Houston, TX, USA). The optimality gap was set as 10 −3 .

Modified IEEE 33-Node Distribution Network
The dataset of the modified test system can be obtained from [31].The total peak load demand and capacity of lines were five times of their original given values. The required voltage violation range ε was set as 0.05 (p.u.). The 24-h load demand profile (percentages of daily peak load, shown in Figure 4) were obtained from [32], and they were multiplied by the peak load to determine the expected daily load demand profile. In addition, the electricity price (Prj,t) was assumed to be timevarying, following the change of a 24-h load demand profile, with its peak value set as 50$/MWh. The price to purchase energy from the main grid (Prin) was assumed to be equal to Prj,t, and the price to sell energy to the main grid (Prout) was set as 20% of Prin. Other parameters for describing the renewables and SVCs are listed in Table 1.

Case Studies
In this section, the proposed framework and methodology were tested on the modified IEEE 33-node distribution system. All the computational tasks were implemented with a Yalmip toolbox in MATLAB (R2016a, MathWorks, Natick, MA, USA) on a 2.4 GHz personal computer with 8 GB of RAM, using Gurobi as the MILP solver (7.0.1, Gurobi Optimization, Inc., Houston, TX, USA). The optimality gap was set as 10 −3 .

Modified IEEE 33-Node Distribution Network
The dataset of the modified test system can be obtained from [31].The total peak load demand and capacity of lines were five times of their original given values. The required voltage violation range ε was set as 0.05 (p.u.). The 24-h load demand profile (percentages of daily peak load, shown in Figure 4) were obtained from [32], and they were multiplied by the peak load to determine the expected daily load demand profile. In addition, the electricity price (Pr j , t ) was assumed to be time-varying, following the change of a 24-h load demand profile, with its peak value set as 50$/MWh. The price to purchase energy from the main grid (Pr in ) was assumed to be equal to Pr j , t , and the price to sell energy to the main grid (Pr out ) was set as 20% of Pr in . Other parameters for describing the renewables and SVCs are listed in Table 1.

Case Studies
In this section, the proposed framework and methodology were tested on the modified IEEE 33node distribution system. All the computational tasks were implemented with a Yalmip toolbox in MATLAB (R2016a, MathWorks, Natick, MA, USA) on a 2.4 GHz personal computer with 8 GB of RAM, using Gurobi as the MILP solver (7.0.1, Gurobi Optimization, Inc., Houston, TX, USA). The optimality gap was set as 10 −3 .

Modified IEEE 33-Node Distribution Network
The dataset of the modified test system can be obtained from [31].The total peak load demand and capacity of lines were five times of their original given values. The required voltage violation range ε was set as 0.05 (p.u.). The 24-h load demand profile (percentages of daily peak load, shown in Figure 4) were obtained from [32], and they were multiplied by the peak load to determine the expected daily load demand profile. In addition, the electricity price (Prj,t) was assumed to be timevarying, following the change of a 24-h load demand profile, with its peak value set as 50$/MWh. The price to purchase energy from the main grid (Prin) was assumed to be equal to Prj,t, and the price to sell energy to the main grid (Prout) was set as 20% of Prin. Other parameters for describing the renewables and SVCs are listed in Table 1.

Location
Capacity WT (MW) PV (MW) SVC (MVar) 17 3.675 1.125 3.250  In this paper, electricity price was linearized stepwise, with ten-piece segments used in [13] to approximate the price elastic demand curve in (2a). Furthermore, the inelastic demand at each node was set as 80% of the actual demand, and the upper bound of total demand was equal to 120%. Accordingly, Table 2 shows the parameters for various costs in this case study, which were obtained from [16,33]. The discrete increment of DG was set to be 10 kVA in this case study. The interest rate z was set as 3%, and the planning horizon was set as 20 years, respectively.

First-Stage Co-Investment Scheme in a Different Uncertainty Set
We set Case 1 as the base case, in which the uncertainty budget is given in Table 3. To comprehensively illustrate the co-investment scheme in different uncertain budgets of renewable energy, three cases were defined with different Γ wt,pv low /Γ wt,pv up . Table 4 depicts the investment of DG in each case. It can be observed that the locations of DGs are similar, but the sizes installed in the particular locations are relatively different. That is mainly because the DG location depends on the topology, which remains unchanged for all cases, but the DG size depends on the severity of DN operation violation, which is further affected by the uncertainty budget. Figure 5 shows the results of DR investment in three robust optimization-based cases, and one deterministic optimization-based case. It can be observed the DR program investment increases when the uncertainty budget is broadened. Furthermore, Table 5 presents the corresponding investment cost of each case. It is worth noticing that the investment costs for DGs are almost invariant, while the DR investment cost is enlarged with different renewables' uncertainty budgets. However, considering that DR investment holds a small proportion of the co-planning cost, a slight increase of total investment cost can help the system to hedge against variation of renewable energy sources, which reveals the economic benefit of the DR program investment.   In this paper, electricity price was linearized stepwise, with ten-piece segments used in [13] to approximate the price elastic demand curve in (2a). Furthermore, the inelastic demand at each node was set as 80% of the actual demand, and the upper bound of total demand was equal to 120%. Accordingly, Table 2 shows the parameters for various costs in this case study, which were obtained from [16,33]. The discrete increment of DG was set to be 10 kVA in this case study. The interest rate z was set as 3%, and the planning horizon was set as 20 years, respectively.

First-Stage Co-Investment Scheme in a Different Uncertainty Set
We set Case 1 as the base case, in which the uncertainty budget is given in Table 3. To comprehensively illustrate the co-investment scheme in different uncertain budgets of renewable energy, three cases were defined with different Γ wt,pv low /Γ wt,pv up . Table 4 depicts the investment of DG in each case. It can be observed that the locations of DGs are similar, but the sizes installed in the particular locations are relatively different. That is mainly because the DG location depends on the topology, which remains unchanged for all cases, but the DG size depends on the severity of DN operation violation, which is further affected by the uncertainty budget. Figure 5 shows the results of DR investment in three robust optimization-based cases, and one deterministic optimization-based case. It can be observed the DR program investment increases when the uncertainty budget is broadened. Furthermore, Table 5 presents the corresponding investment cost of each case. It is worth noticing that the investment costs for DGs are almost invariant, while the DR investment cost is enlarged with different renewables' uncertainty budgets. However, considering that DR investment holds a small proportion of the co-planning cost, a slight increase of total investment cost can help the system to hedge against variation of renewable energy sources, which reveals the economic benefit of the DR program investment.

Second-Stage Operation Result
After the first-stage optimal DR and DG allocation location is given, the worst-case scenario is found based on the second-stage operation strategy. The total hourly DG output and the hourly energy transaction with the main grid in Case 1 are presented in Figure 6. It can be clearly seen that the DG outputs are high during period from 8:00 to 9:00 and from 19:00 to 23:00. According to the daily load profile shown in Figure 4, this is mainly because the electricity consumption is high and renewable energy is low at these time periods. Specifically, the generation of DGs reach their maximum output between 19:00 and 20:00, when the load demand hits its peak value and renewables reach their lower bound. In this situation, the energy provided by all the DGs and the energy reduced by DR program combined is not sufficient to fulfill the peak load, and thus the DN needs to buy energy from external grid. Besides, during the period from 24:00 to 6:00 and from 13:00 to 16:00, when the wind and solar energy is supposed to be abundant and the load demand is low and stable, the DN sells energy to the main grid and the DGs are not committed on. Figure 7 displays the load shifting under DR-enabled nodes in Case 1, where the load shifting is defined as the difference between DR load and expected demand. It can be seen that the load demand is enlarged between 14:00 and 15:00, when electricity price is low but the wind power and solar energy are more than enough to supply the low power demand. Similarly, the DR is also fitted up to augment the electricity consumption between 24:00 and 6:00, due to the abundant wind power and low energy price. When the wind energy generation is insufficient, and the energy purchased from external grid is expensive, the electricity consumption is reduced from 8:00 to 9:00 and 18:00 to 22:00, respectively. Thus, it can be found that the coordinated operation of DR program and DG can effectively accommodate the uncertainties of renewables and load. respectively. Thus, it can be found that the coordinated operation of DR program and DG can effectively accommodate the uncertainties of renewables and load.

Voltage Profile with Different Demand Response Ability
To verify the robustness in maintaining system security, the worst-case scenario obtained in Case 1 is implemented, using both the proposed model and the deterministic model. Figure 8 depicts the 24-h voltage profit for the critical node in the test system. It can be observed that severe voltage issues, including voltage over-drop and over-rise, emerge in deterministic results at 15:00 and from 19:00 to 21:00, which leads a large penalty cost. Nevertheless, with the co-investment decision in Case 1 and the proposed operation strategy, the voltage fluctuations in the worst-case scenario are maintained within the acceptable range and the voltage violations are significantly diminished. Again, the proposed co-planning scheme and operation strategy is proved to be able to accommodate the renewable energy both economically and robustly. respectively. Thus, it can be found that the coordinated operation of DR program and DG can effectively accommodate the uncertainties of renewables and load.

Voltage Profile with Different Demand Response Ability
To verify the robustness in maintaining system security, the worst-case scenario obtained in Case 1 is implemented, using both the proposed model and the deterministic model. Figure 8 depicts the 24-h voltage profit for the critical node in the test system. It can be observed that severe voltage issues, including voltage over-drop and over-rise, emerge in deterministic results at 15:00 and from 19:00 to 21:00, which leads a large penalty cost. Nevertheless, with the co-investment decision in Case 1 and the proposed operation strategy, the voltage fluctuations in the worst-case scenario are maintained within the acceptable range and the voltage violations are significantly diminished. Again, the proposed co-planning scheme and operation strategy is proved to be able to accommodate the renewable energy both economically and robustly.

Voltage Profile with Different Demand Response Ability
To verify the robustness in maintaining system security, the worst-case scenario obtained in Case 1 is implemented, using both the proposed model and the deterministic model. Figure 8 depicts the 24-h voltage profit for the critical node in the test system. It can be observed that severe voltage issues, including voltage over-drop and over-rise, emerge in deterministic results at 15:00 and from 19:00 to 21:00, which leads a large penalty cost. Nevertheless, with the co-investment decision in Case 1 and the proposed operation strategy, the voltage fluctuations in the worst-case scenario are maintained within the acceptable range and the voltage violations are significantly diminished. Again, the proposed co-planning scheme and operation strategy is proved to be able to accommodate the renewable energy both economically and robustly. issues, including voltage over-drop and over-rise, emerge in deterministic results at 15:00 and from 19:00 to 21:00, which leads a large penalty cost. Nevertheless, with the co-investment decision in Case 1 and the proposed operation strategy, the voltage fluctuations in the worst-case scenario are maintained within the acceptable range and the voltage violations are significantly diminished. Again, the proposed co-planning scheme and operation strategy is proved to be able to accommodate the renewable energy both economically and robustly.

Statistical Feasibility Check
To verify the robustness in addressing various scenarios, the proposed model was tested with 1079 scenarios, randomly generated from a three-year dataset of wind power, PV, and load demand

Statistical Feasibility Check
To verify the robustness in addressing various scenarios, the proposed model was tested with 1079 scenarios, randomly generated from a three-year dataset of wind power, PV, and load demand in [32,34,35]. The deterministic model was also implemented as benchmark. Table 6 shows the statistic results of the simulations. As seen from the results, the merits and the defects of the proposed model were validated from different aspects. The average operation cost of the proposed model at $28,420.20 was larger than that of deterministic model, at $26,478.56. Using the proposed model, more DGs were invested, and thus the fuel cost increased accordingly. However, the average energy transaction cost was considerably reduced. The DN revenue in the proposed model slightly decreased as more demand load joined the DR program. Besides, the average maximal voltage violation rate was significantly reduced, from 7.546% to 0.006%. As a result, the average penalty cost was also considerably reduced. It should be noticed that the probability that the proposed model outperforms the deterministic model in total second-stage cost and voltage security is 95.85% and 98.99% respectively, which indicates the robustness of proposed model in maintaining the DN operations.

Modified IEEE 123-Node Distribution Network
To demonstrate the scalability and effectivity of the proposed framework, a modified IEEE 123-node DN in [36] was also used as a test system, which additionally included eight 1.1 MW WT, 0.5 MW PV, and 1.5 MVAR SVC in the same locations, which were connected at nodes 6, 37, 42, 56, 62, 69, 86 and 114. The 24-h load demand profile was the same as the 33-node system, and the total peak load demand and capacity of lines were five times of their original given values. The other parameters were the same as Section 5.1. Table 7 shows the co-planning scheme in three robust optimization-based cases and one deterministic optimization-based case, in which the uncertainty sets were the same as in Section 5.2. Furthermore, Table 8 presents the corresponding investment cost of each case. In general, it could be observed that results of DGs and DR investment showed the same trend with the results in the 33-node test system, which verifies the efficiency and scalability of the proposed framework.

Conclusions
This paper proposes a DRF and DG co-investment model to determine the optimal capacity and location of demand response and DG. A robust, optimization-based two-stage model was adopted to accommodate the uncertainties of renewables and load demand. The price-sensitive demand was modeled with modified PEDR, where both inelastic and elastic demand was taken into account. In the proposed model, the overall cost was optimized, taking into consideration the first-stage investment decision, as well as the hourly DG and DR coordinated operation in second-stage operation. The customer bills, total energy supply, and operation regulations were considered to motivate the electricity end users to participate in the DR program. Finally, a modified CCG/OA algorithm was employed to solve the proposed two-level programming problem. Case studies demonstrate the effectiveness of the proposed model from aspects of economic and security operations. With the aid of the proposed model, not only is the operation violation of the DN alleviated, but the system cost is also minimized. Integer variable indicating bth increments in size of DG at j P j,t /Q j,t Active/reactive power flow at j,t V j,t Voltage at j,t