Warm-start piecewise linear approximation-based solution for load pick-up problem in electrical distribution system

: As the core sub-problem of both network reconfiguration and service restoration of the electrical distribution system (EDS), the load pick-up (LPP) problem in EDS searches the optimal configuration of the EDS, aiming to minimise the power loss or provide as much power as possible to load ends. The piecewise linearisation (PWL) approximation method can be used to tackle the network power flow constraints’ non-linearity in the LPP problem model, and transform the LPP model into a mixed- integer linear programming model (LPP-MILP model). The errors in the PWL approximation of the network power flow constraints may affect the feasibility of the solving results of the LPP-MILP model. The single method to reduce the PWL approximation errors by increasing the number of discretisations in PWL function is not stable. Moreover, the solution efficiency of the LPP-MILP model is sacrificed severely. In this study, a warm-start PWL approximation-based solution for the LPP problem is proposed. The variable upper bounds in the PWL approximation functions are renewed dynamically in the warm-start solution procedure to reduce the PWL approximation errors with higher computational efficiency. Modified IEEE 33-bus radial distribution test system and a real and large distribution system, 1066-bus system, are used to test and verify the effectiveness and robustness of the proposed methodology.


Introduction
The load pick-up (LPP) problem in the electrical distribution system (EDS) searches the optimal configuration of the EDS, aiming to minimise the power loss or provide as much power as possible to load ends. The LPP problem is the core sub-problem of both network optimisation reconfiguration and service restoration (reconfiguration after fault or blackout) of the EDS. One fundamental model of the LPP problem for network optimisation reconfiguration is presented in [1]. And based on [1], the model of the LPP problem for service restoration is proposed in [2]. For limited generating capacity in the EDS, the optimisation models of the LPP problems for network optimisation reconfiguration (minimise power loss) and service restoration (maximise outage load restored) are mathematically equivalent [3]. The LPP problem belongs to the category of mixed-integer problem since the operational status of feeder or bus is defined strictly as 0 or 1. Along with other continuous variables, objective function and constraints, the LPP problem can be formulated as a mixed-integer non-linear programming problem. Referring to models in [4], network power flow constraints are an essential part of the constraints in the LPP problem model (LPP model in short). And the network power flow constraints make up the majority of non-linear constraints in the LPP problem.
To tackle the network power flow constraints' non-linearity, based on the DistFlow equations for radial EDS network [5], some major power flow constraint formulations (e.g. DC power flow constraints [6], second-order cone programming (SOCP) based power flow constraints [7,8], relaxed semi-definite programmingbased power flow constraints [9] and piecewise linearisation (PWL) approximation-based power flow constraints [10]) are constructed and used in the LPP model. In [7], the single quadratic equality constraint in the DistFlow equations is relaxed to secondorder conic constraint, and the original non-convex three-phase optimisation model was converted to a mixed-integer SOCP model, which can be solved by commercial mathematical programming solvers (e.g. CPLEX and Gurobi). In [10], to further reduce the computational complexity, quadratic variables in the DistFlow equations are approximated by PWL functions, and then the LPP model is transformed into a mixed-integer linear programming (MILP) model, and can be solved by commercial solvers effectively.
However, there are some disadvantages existing in the candidate power flow constraints in the LPP model. The DC power flow constraints cannot account for the magnitudes of bus voltage and the distribution of reactive power in the EDS, which may generate an impractical solution. For the SOCP-based power flow constraints, the rigorous SOCP relaxation optimality preconditions to guarantee the SOCP exactness may hinder the applicability of the SOCP-based power flow constraints. What is more, as indicated in [11], the recovery of AC feasibility from the solution results of SOCP model may be affected due to too large relaxation error incurred by SOCP relaxation. For the PWL approximationbased power flow constraints, quadratic variables in the DistFlow equations are approximated by PWL functions, the approximation errors and the accuracy analysis must be researched especially. Previous works of literature hardly concentrate on these topics. This paper utilises the PWL function approximation method (PWL approximation method in short) to linearise the network power flow constraints of the LPP model and analyses the accuracy and efficiency of the application of this method systematically.
The concept of 'warm-start' has been mentioned in previous literature [12]. The main idea of the warm-start solution is based on the multi-step solution procedure, which includes optimal solution obtained from the first step and then using the first step results for an iterative solution in the next step. Take the two-step solution procedure in [13] for example. In the first step, the voltage magnitude variables in the MILP model are set to the system nominal voltage magnitude (constant), and then the MILP model is solved. In the second step, voltage magnitude variables are renewed by the quadratic root of voltage magnitude square variables results obtained from the first step, and the MILP model is then resolved for a more accurate optimisation result. Some other current published research work related to warm-start can be found in [14][15][16]. Yang et al. [14] propose a warm-start network model adopting the initial points obtained from the DC OPF model, which can further enhance the performance of the OPF model. In [15], the authors investigate using learning approaches to find warm-start points with high quality to solve OPF iteratively. In [16], to achieve efficient management of distribution generations (DGs) and smart control of capacitors in EDSs, a quadratic programming (QP) model built on linearised quadratically constrained programming (QCP) model is proposed. To solve the QP model, estimated values for bus voltage magnitudes and line power flows are solved from the QCP model and then used for initialising the solution of the QP model.
As mentioned before, the accuracy and efficiency of the application of the PWL approximation method to linearise the power flow constraints of the LPP model in the EDS must be considered. In this paper, to improve the performance of the PWL approximation method, a warm-start PWL approximation-based solution for the LPP problem in the EDS is proposed. Apart from the network power flow constraints, other essential constraints of the LPP problem in the EDS are presented, and make up the entire LPP model. The PWL approximation method is introduced to linearise the network power flow constraints, and transform the LPP model into an MILP model (LPP-MILP model). And a multistep warm-start solution procedure of the LPP-MILP model is established. Error indices are constructed to evaluate the accuracy of the proposed warm-start PWL approximation-based solution for the LPP problem in the EDS.
For the rest of the paper, the mathematical formulas of the LPP model are presented in Section 2. The detailed warm-start PWL approximation-based solution for the LPP model are analysed in Section 3. The LPP-MILP model complexity assessment and error indices to evaluate the accuracy of the proposed methodology are presented in Section 4. Case studies to demonstrate the performance of the proposed methodology are included in Section 5. Eventually, Section 6 concludes the paper.

Objective function
For the two main different application scenarios for the LPP problem in the EDS (i.e. LPP problem for network optimisation reconfiguration and LPP problem for service restoration), different objective functions are as follows: For the LPP problem for network optimisation reconfiguration, minimising the power loss of the EDS is considered as the optimisation goal in the LPP model, and the corresponding objective function is shown in (1) [17]. For the LPP problem for service restoration, shown in (2), maximising the restored load is set as the objective function in the LPP model [2]. As indicated in [3], for deterministic generation output in the EDS (i.e. ∑ n ∈ S bus P G, n = Constant), objective functions in (1) and (2) are mathematically equivalent.

Constraints
To formulate the mathematical formulas of the LPP model, assumptions in [4] are considered in this paper: (a) a balanced three-phase AC power flow is considered, and the EDS is represented by a single-phase equivalent; (b) the radial topology should be maintained for the EDS network; (c) the operation limits, such as bus voltage magnitude limits and feeder current capacity limits must be satisfied. It is worthy to point out that, the proposed optimisation model is formulated statically. In other words, it focuses on the optimisation of a single period. Consequently, the temporal-relating constraints are not considered in the paper. Once the optimal solution is obtained, the optimal topology and the decision of switching of the EDS is decided. In this part, constraints considered in the LPP model are presented in detail, including the network power flow constraints, the radial topology constraints, and other essential constraints. Fig. 1 is the schematic diagram of a partial radial EDS network. For radial EDS network, the widely used DistFlow network power flow constraints can be described as follows [4]:

Network power flow constraints:
The traditional equations of active and reactive power balance for each bus are represented as (3) and (4). The voltage difference across the feeder ij can be obtained from (5). The relationship among the bus voltage magnitude, feeder current magnitude, and feeder active/reactive power is illustrated in (6), which is the only non-linear constraint in the DistFlow network power flow constraints. Limits for bus voltage magnitudes and feeder current magnitudes are shown in (7) and (8), respectively. Note that only quadratic terms of the magnitudes of bus voltage and feeder current (i.e. V i sqr and I i j sqr ) are used in (3)- (8), hence they shall not be treated as quadratic variables By introducing the state variables of buses and feeders, constraints (5)-(8) can be reformulated in (9)- (12). The state variables assumes a value of 1 if the bus or feeder is in the use state, and 0 if the bus or feeder is in the non-use state. Referring to [4], by introducing a sufficiently big constant number M, constraints in (9)- (12) can be relaxed if the corresponding bus or feeder is in the non-use state. In the LPP model, (3)-(4) and (9)-(12) make up the integrated network power flow constraints.

Radial topology constraints:
Incorporating DGs into EDS planning and operation will be gainful for strengthening the efficiency and tackling some technical issues [18]. In recent years, with the large-scale DGs integrated in the EDS, the island is an important operation mode of the EDS with DGs [19]. On basis of the idea of the fictitious network and fictitious power in [20], a directed fictitious network-based radial topology constraint is constructed in [6], and can be suitable for the EDS with an indeterminate number of islands Fig. 2 shows a schematic diagram of a directed fictitious network of the EDS. In Fig. 2, a fictitious bus (bus 0 ) with fictitious feeders connecting bus 0 and all DG buses are added to the original network. The bus and feeder sets of the fictitious network are defined in (13) and (14), respectively. All feeders in the fictitious network are transformed into directed feeders. The directed feeder set is defined in (15). Directed feeder ij means feeder with start point i and end point j The formulas of the radial constraints based on the directed fictitious network are presented in (16)-(22). The name of the fictitious bus (bus 0 ) is simplified as 0 for clearer expression in the formulas. In (16), every bus (except bus 0 ) in use state in the fictitious network has a unit of fictitious load. All fictitious power transmitted in the fictitious network is generated from bus 0 . The fictitious power balance equations in (16) ensures the connectivity of the fictitious network (all buses and feeders in use state). It is worthy to note that H ki and H i j takes a discrete integral value and v i f is integer variable which takes value 0 or 1. Therefore, (16) can be satisfied. The constraints of the state variables of buses and feeders in the directed fictitious network are shown in (18)- (20). The mapping relationships for state variables in the fictitious network and the original network are in (18) and (19). The default values of bus 0 and fictitious feeders pointed to bus 0 are set in (20) and (21), respectively. For a connected network including all buses and feeders in use state, constraint (22) makes every bus (except bus 0 ) in use state capable of having one and only parent bus which ensures that the connected network is a tree-like branching network. For the constraints in (16)-(22), the radial topology of the fictitious network can be guaranteed, and accordingly the radial topology of all islands in the EDS with DGs (original network) can be guaranteed.

Other constraints:
Other essential constraints in the LPP model are listed as follows: The relationship between bus and feeder state variables are indicated in (23). When the feeder ij is in use state, the corresponding endpoints of the feeder v i and v j are accordingly in use state. Restricted by the state variables, the limits of P G, i ,Q G, i , P L, i ,Q L, i , P i j and Q i j (i.e. security operational constraints) are shown in (24)-(29), respectively. It is worthy to point out that the LPP problem is the core subproblem of network reconfiguration and service restoration problems. We do not intend to solve the network reconfiguration and service restoration problems simultaneously. Herein we formulate their respective objective functions and common constraints in the previous sections and summarise them in Table 1. Constraints (28) and (29) are the network load restoration constraints in service restoration problem, which indicate that nodal load curtailment is considered due to the total generating capacity of DGs in the EDS may be less than the total active loads. In contrast, in the network reconfiguration problem, the load must be satisfied 3 Warm-start PWL approximation-based solution for the LPP model

PWL approximation of network power flow constraints
In the LPP model presented in Section 2, only constraint (6) or (10) (considering the state variables) contains quadratic variables is non-linear in the network power flow constraints. As a standard procedure in electrical system optimisation problems, the PWL approximation method is used to process the quadratic variables in (6) or (10), and to linearise the non-linear network power flow constraint [10]. To make the presentation clearer, constraint (6) is used to illustrate the PWL approximation procedure: Firstly, the term V j sqr in the left side of (6) is approximated by (V bus norm ) 2 (shown in (30)), and constraint (6) is reformulated as (31). Due to the narrow voltage magnitude interval [V bus min , V bus max ] in the EDS, the simplification of (30) has a relatively low approximation error [10]. Secondly, the quadratic variables P i j 2 and Q i j 2 in (30) are approximated by PWL functions f (y, y¯, Λ) (shown in (32) and (33)), and then constraint (6) can be further expressed as (34). Similarly, considering the state variables of the EDS, constraint (10) can be formulated as (35) The PWL function f (y, y¯, Λ) is defined in (36)-(41). Considering the paper length, the schematic diagram of the PWL function approximates the quadratic curve and the detailed presentation of the PWL function can be referred to [10] f (y, y¯, Λ) = ∑ 0 ≤ Δ y, λ ≤ y¯/Λ∀λ = 1, 2, . . . , Λ ϕ y, λ = (2λ − 1)y¯/Λ∀λ = 1, 2, . . . , Λ In the subsequent sections, by selecting (2) in Section 2 as the objective function of the LPP model, by linearisation of the network power flow constraints, we choose the service restoration problem as the example to demonstrate the effectiveness of the proposed warm-start PWL approximation-based solution.

Warm-start solution procedure of the LPP-MILP model
The errors in the PWL approximation of the network power flow constraints may affect the feasibility of the solving results of the LPP-MILP model. In general, when implementing the PWL approximation method, the bigger value of Λ contributes to better fitting of the quadratic curves of P i j 2 and Q i j 2 , and therefore the approximation and the solution will be more accurate. However, the solution efficiency may be meanwhile sacrificed severely due to the large-scale variables generated by the PWL functions, which sharply increases the dimensions of the corresponding optimisation problem Except for the number of discretisation Λ, the set value of upper limit (e.g. P i j max or Q i j max in (35)) in the PWL function is another important factor influencing the approximation accuracy of the quadratic variable. For given V bus max and I i j max , P i j max and Q i j max can be calculated by (42) and (43), respectively. For a certain I i j max for an EDS (e.g. I i j max = 200 A for a 44-bus distribution system in [8]), the actual active/reactive power through feeder P i j and Q i j may be far less than P i j max and Q i j max in some cases, especially in the service restoration problem due to the limited generating capacity in the EDS. The big gaps between the actual values and the upper limit set values will affect the accuracy of the PWL approximation method heavily.
In terms of the solution of the LPP-MILP model, for a small set value of Λ, if the upper limit set values can be rated closer to the actual values, the accuracy of the PWL approximation method can be improved, meanwhile the computational complexity is not increased. On basis of this idea and the concept of 'warm-start' presented in the introduction, a multi-step warm-start solution procedure of the LPP-MILP model is proposed. The flow chart of the warm-start solution procedure of the LPP-MILP model is shown in Fig. 3.
In the warm-start solution flow chart, the initial generated LPP-MILP model is solved firstly. And then the active/reactive power upper limit set values of feeders in use state (w i j = 1) are renewed by the results solved in the first step. The renewal principles of P i j max and Q i j max designed in Fig. 3 can ensure the feasibility of the regenerated LPP-MILP model. After error analysis, if errors between P i j 2 and f (P i j , P i j max , Λ) or Q i j 2 and f (Q i j , Q i j max , Λ) exceed the error threshold and within the maximal iteration number, the warmstart solution step can be reduplicated. Table 2 shows the complexity assessment of the LPP-MILP model in this paper. To assess the complexity of the LPP-MILP model, the counts of total variables and constraints in the LPP-MILP model are analysed and listed in the table.

LPP-MILP model complexity assessment and error indices construction
The PWL approximation method of network power flow constraints has been detailed presented in Section 3.1. The approximations in (30), (32) and (33) make up the core steps of the PWL approximation method. As the simplification of (30) has a relatively low approximation error [10], the errors caused by approximations in (32) and (33) mainly determine the accuracy of the PWL approximation method in the LPP-MILP model

Case study
In this section, the proposed multi-step solution procedure for the LPP problem is tested using two EDSs, i.e. the modified IEEE 33 bus system and a real 1066-bus distribution system with 44 DGs in Shandong province of China [21,22]. The LPP-MILP model and the warm-start solution procedure are implemented in MATLAB   with YALMIP toolbox, and CPLEX 12.5 is the MILP solver. The case studies are performed on a 2.2-GHz dual-core PC with 4 GB RAM.

Case I: modified IEEE 33-bus system
Considering different DGs allocation plans, two cases based on the modified IEEE 33-bus radial distribution test systems are studied to demonstrate the performance of the proposed methodology. The IEEE 33-bus radial distribution test system are with 33 buses and 37 feeders (including 5 tie lines). The detailed network structure, bus load parameters, and feeder parameters of the test system can be referred to [23]. The system base voltage is 12.66 kV and the base capacity is 10 MVA. Referring to [10], the bus voltage and feeder current limits are set as follows: V bus min = 0.95V bus norm , V bus max = 1.05V bus norm , and I i j max = 200 A i j ∈ S feeder . : Fig. 4 shows the topology of a modified IEEE 33-bus radial distribution test system with 3 DGs. The DG parameters of the modified test system are listed in Table 3.

Modified IEEE 33-bus system with 3 DGs
For the modified test system with 3 DGs, Table 4 shows the solving results of the LPP-MILP model without the warm-start solution procedure (i.e. warm-start iteration procedure not selected in Fig. 3) under different set values of Λ. Take the No. 1 scenario (Λ = 10) in Table 4 for example, the big errors in the approximation of the network power flow constraints (E p m > 100% and E q m > 1000%) will affect the feasibility of the solving results of the LPP-MILP model. Referring to Table 4, the bigger set value of Λ can lead to lower error index values, and promote better accuracy of the PWL approximation method basically. The complexity assessment of the LPP-MILP model has been analysed and listed in Table 2, with the growth of the set value of Λ, the counts of variables and constraints of the LPP-MILP model meet dramatically increase. Accordingly, as observed and demonstrated in Table 4, the modelling and solving time of the LPP-MILP model increases significantly with the growth of the set value of Λ.
By further comparing the solving results under different set values of Λ in Table 4, we find the single method to reduce the error index values by increasing the set value of Λ is not stable. The detailed LPP results of the modified test system under Λ = 50 (No. 5 scenario in Table 4) are shown in Fig. 5. According to Table 4, for Λ = 50, the solving results miss the optimal objective function value, and the error index values meet relative increment.
By referring to [10] and selecting a more stringent set value of Λ, relaxing the error threshold and the maximum iteration number limit in Fig. 3, Table 5 shows the solving results of the LPP-MILP model by the warm-start solution procedure under different iteration number with fixed set value of Λ (Λ = 10). In Table 5, the   No. 1 scenario (iteration number = 0) is just the solving results without the warm-start solution procedure in Table 4 under Λ = 10.

Fig. 4 Modified IEEE 33-bus radial distribution test system with 3 DGs
The detailed LPP results of the modified test system by the warm-start solution procedure (Λ = 10, iteration number = 1, No. 2 scenario in Table 5) are shown in Fig. 6. After one warm-start iteration process, the error index values have a significant decline, and the accuracy performance is better than the No. 9 scenario (Λ = 90) in Table 4 and with a much shorter modelling and solving time. After two warm-start iteration processes (No. 3 scenario in  Table 5 shows that with the increasing of iteration numbers, the error index values reveal a trend of gradual decrease, which can demonstrate the effectiveness and robustness of the proposed warm-start solution procedure of the LPP-MILP model. By comparing the modelling and solving time in Table 4, the accumulated modelling and solving time in Table 5 show the advantage in solution efficiency of the proposed methodology. Fig. 7 shows the topology of a modified IEEE 33-bus radial distribution test system with 2 DGs. The DG parameters of the modified test system are listed in Table 6.

Modified IEEE 33-bus system with 2 DGs:
To further test the effectiveness and robustness of the proposed methodology under smaller set value of Λ, the number of discretisation Λ is set as 5 in this case. For the modified test system with 2 DGs, solving results of the LPP-MILP model without the warm-start solution procedure are listed in Table 7. For warm-start solution procedure in Fig. 3, set the error thresholds as E p w − m ≤ 0.1% and E q w − m ≤ 0.1%, and set the maximal iteration number as 5. Solving results of the LPP-MILP model by the warmstart solution procedure are listed in Table 8. The iteration number of the solving results by the warm-start solution procedure is 3.  For the solving results without the warm-start solution procedure in Table 7, under small set value of Λ (Λ = 5), the error index values are very large. And compared to solving results in Table 8, the optimality of the objective function value has also been affected. For the solving results by the warm-start solution procedure in Table 8, after three warm-start iteration processes, error indices are effectively decreased within the error thresholds. Solving results in Table 8 demonstrates that under smaller set value of Λ, the proposed warm-start solution procedure of the LPP-MILP model can also keep its effectiveness and robustness.

Case II: 1066-bus system with 44 DGs
To demonstrate the effectiveness of the proposed method when applying to real and large EDS, in this case, the 1066-bus system in Shandong province of China is selected as the test system.
For 1066-bus system, the limit of bus voltage magnitude, feeder current and the number of discretisation in the PWL function are set as below: V bus min = 0.95V bus norm , V bus max = 1.05V bus norm , I i j max = 500 A i j ∈ S feeder , and Λ = 10. For the solution procedure in Fig. 3, set the error thresholds as E p m ≤ 0.1% and E q m ≤ 0.1%, the maximal iteration number as 5, and the optimality gap for MILP solver as 0.1%.
For the solving results by the warm-start solution procedure (Λ = 10) in Table 9, after 3 iteration processes, the error indices of the solving results are effectively reduced within the error thresholds, and the multi-step solution procedures are completed.
The number of restored feeders, buses and self-healing islands of the 1066-bus EDS by the multi-step solution procedure are listed in Table 10. As the iteration number increases, the number of restored feeders and buses gradually converges to the optimal results with progressively higher accuracy. From Table 10, the solution results with relatively high approximation errors may lead to impractical or sub-optimal solving results of the service restoration optimisation scheme.

Conclusion
This paper presents a warm-start PWL approximation-based solution for LPP problem in EDSs. For the LPP-MILP model constructed through linearisation of the network power flow constraints by the PWL approximation method, the errors in the PWL approximation method may affect the feasibility of the solving results. The accuracy of the PWL approximation method and the solution efficiency of the LPP-MILP model are considered and analysed systematically in this paper. A warm-start solution procedure of the LPP-MILP model is established in this paper to reduce the PWL approximation errors with higher computational efficiency.
Comparisons of solving results in the case studies of IEEE-33 bus system show that for the direct solution procedure (without the warm-start solution procedure) of the LPP-MILP model, the single method to improve the accuracy of the PWL approximation method by increasing the number of discretisation in PWL function is not stable, meanwhile, the solution efficiency is sacrificed severely. The proposed warm-start solution procedure of the LPP-MILP model can effectively reduce the PWL approximation errors with an advantage in solution efficiency. Considering different set values of the number of discretisation in the PWL function, the robustness of the proposed methodology is tested and demonstrated in the case studies. To stress the issue of application in relative large EDS, the proposed multi-step solution procedure for the LPP problem is further tested on a 1066-bus real distribution system in Shandong province of China.
Future work in this area will include the improvement of the linearised network PF constraints aiming at application into threephase unbalanced distribution system.