Efficient Unbalanced Three-Phase Network Modelling for Optimal PV Inverter Control

High penetration levels of renewable energy generation in the distribution network require voltage regulation to avoid excessive voltage at generating nodes. To effectively control the network and optimize network hosting capacity, the distribution system operator must have an efficient model for power flow analysis. This paper presents the formulas and steps to express the power flow analysis equations of an unbalanced 3-phase network in matrix form suited to programmed solutions. A benchmark MATLAB/Simulink network with unbalanced distribution lines, photovoltaic inverters, and loads is built to verify the matrix model. To demonstrate the application of the model, the control of reverse energy flow from the photovoltaic inverters to keep the voltage in the network below the regulated level is simulated. Two decentralized control algorithms are applied in the network, including an on/off and a multi-objective constrained optimization controller. The detailed construction of the optimization problem for the 3-phase network in matrix form, which is consistent with the power flow calculation, is described. Simulation with the control methods over a day shows that the total active power of the on/off and optimized controllers deliver 41.92% and 99.39% of the available solar power, respectively, while maintaining the network node voltages within limits.


Introduction
The increasing penetration of photovoltaic (PV) distributed generation (DG) in the low voltage (LV) distribution system raises difficulties for the distribution system operator (DSO) in maintaining voltage stability. Reverse power flow from PV sources causes a rise in node voltage and this is a key factor limiting the PV hosting capacity of a network [1]. The introduction of emerging smart-grid technology with advanced system control should enable the widespread integration of renewable energies. The smart grid with smart meters and high-speed broadband communication allows for the use of decentralized controls to address these distribution system problems [2,3]. There are different methods to determine the network voltage by Gauss-Seidel [4], Newton-Raphson [4][5][6][7][8][9], Fast Decoupled Load Flow (FDLF), Backward/Forward sweep power flow (BFS) [10][11][12] and a linear approach [13][14][15]. In this paper, the Newton-Raphson method is selected for its robustness and quadratic convergence [16,17]. In the Newton-Raphson method, the components can be expressed in three coordinates of polar, Cartesian and complex forms for both the power and current mismatch formulations [6].
The power flow analysis problem is developed mainly for constant power P-Q buses with one constant voltage swing bus. The power flow is also applied for constant power-voltage P-V buses [18] and different load models, such as constant impedance (Z), constant current (I) [8] and polynomial or combined constant impedance (Z), constant current (I) and constant power (P) ZIP loads [19]. There are power flow investigations of the transformer tap change [20], the FACTS (Flexible AC Transmission System) controllers, including IPFC (Interline Power Flow Controller) [21], UPFC (Unified Power Flow Controller) [22][23][24], STATCOM (Static Synchronous Compensator) [18,20]. These power flow calculations cover devices and components in the single-phase and balanced 3-phase distribution network but an extension is needed for the unbalanced 3-phase network.
In practice, the loads and DG in the distribution network are not all equally connected between three phases and the components of resistance and inductance of the 3-phase distribution conductors are not the same. Therefore, in the case of unbalanced 3-phase [25][26][27], there is the need to formulate an unbalanced 3-phase power flow analysis for the network with different voltage magnitude in each phase. Similar to the balanced network, there are studies to address the unbalanced 3-phase, including P-V buses [9,14,15] and ZIP load [14,15]. The effect of transformers and loads with different wye and delta connections is considered in the 3-phase network [14,28]. The Newton-Raphson power flow calculation methods are also presented in an unbalanced network as in [7,9]. These studies show the power flow calculation of the network and the steady-state results of the network analysis. The next step is the application of the power flow calculation into a network with a time-scale and control strategies.
Some previous research applied the power flow for the optimization of the PV inverters in the network [29], including the single-objective [30] and multi-objective optimization [31][32][33]. However, these applications are limited to the power calculation of the unbalanced 3-phase network. In this paper, the unbalanced power flow analysis is applied to the distribution network integrating residential loads and PV inverters for the optimization in matrix form. The control algorithms applied in this work are based on a decentralized control where all PV inverters in the network are controllable with the reactive power support to regulate the bus voltage and power factor within the prescribed limits.
In this paper, the formulas are all presented in matrix form which shows advantages in coding and presentation. This is similar to using the complex form whose amplitude and angle or the real and imaginary parts are presented in a single complex number. There are previous efforts to reduce the complexity of the formulas and calculations of power flow in a different presentation [34] and in complex form [6,35], which cut down the dimensions of Jacobian, voltage and power matrices by half. Below are the typical formulas of Newton-Raphson current mismatch in a specific phase p at bus I, where p is presented for phase a, b and c [7].
The presentation of phase in an index listed format as in the above equations is also in [9,15]. The matrix form of complete power flow formulas and the optimization construction is presented in Sections 2 and 3. The 3-phase matrix presentation shows the same goal to the research of complex power flow analysis. Some papers represent the power flow equation in matrix form [13,14]; however, they are applied for the linear approach and the optimization control is not presented. This paper put an effort to present all elements in the 3-phase matrices and treat them in formulas as the balanced system and reduce the amount of formula by a third by combining the three specific phases into a matrix. This presentation also avoids the confusion of switching between the scalar and matrix forms in the same formula. Besides, this reduces the mistakes or errors in the programming and reduces the number of coding loops. After the calculation of the power flow of the unbalanced 3-phase network, the optimization control is applied to the network and the construction of the control is also presented in the matrix form which is consistent with the power flow calculation. This paper includes two main parts. The first part develops a process for network description and power flow analysis in a matrix form as seen in Section 2. The process is applied to describe a benchmark 20-bus distribution network model in Section 3. Power flow analysis is conducted with the model under a measured irradiation profile and this results in overvoltage conditions with node voltages as high 1.3 pu. The second part, Sections 4 and 5, employs the same network model and applies two different decentralized control methods to eliminate overvoltage conditions. The decentralized algorithms use smart grid communication to provide knowledge of the load and potential PV power at each node to a central computer where the algorithm employs power flow analysis to identify overvoltage conditions and modifies the PV power output to maintain the voltage within limits. Section 4 describes the two control algorithms: a simple overvoltage control and a multi-objective optimization algorithm. The detailed construction of the constrained optimization problem in matrix form, consistent with the power flow calculation, is described. Section 5 presents a comparison of the performance results of the control methods. The discussion of the contribution and comparison to other research as well as the limit scope of this paper is presented in Section 6.

Power Flow Matrix Formation of an Unbalanced 3-Phase Network
In a network of N buses with known distribution line parameters and scheduled power of P sh and Q sh , the voltage magnitude |V| and voltage angle δ can be calculated using the Newton-Raphson method [4,5] with the current-mismatch function ∆I in the polar-coordinate format [6,7].
The Newton-Raphson method uses an iteration of calculations to find the results. The initial values of angle and magnitude of all buses are set at δ (0) and |V| (0) as in (3). In Equation (4), the next estimates of δ (h+1) and |V| (h+1) are computed by adding the previous values to the differences of ∆δ and ∆|V|, which are determined in (5).
Bus 1 is the swing bus, the voltage magnitude |V 1 | and angle δ 1 of which are known. The values of the remaining (N − 1) buses need to be determined using the Jacobian matrices, J 1 to J 4 , that calculate the partial derivatives of current with respect to voltage angles and magnitudes.
Each element inside matrices ∆δ, ∆|V| and ∆I is a 3 × 1 matrix. The matrices J 1 to J 4 have the same dimension as J, and each element of J is a matrix whose dimension is 3 × 3. The formulas for all elements of the matrices in (5) are presented in Table 1. Matrix Y is the system admittance whose conductance G and susceptance B is calculated from the netlist of the distribution network, as shown in (6) in Section 3. The values P sh and Q sh are the scheduled powers set for all buses in the network. Where there is no load or generator, the scheduled active and reactive powers are set to zero.
The flowchart in Figure 1 summaries the steps to follow for the calculation of the power flow voltage magnitude and angle. Steps 1 to 5 are presented in Section 3 with a sample benchmark network. As seen in Figure 1, the iteration is stopped when all elements of ∆I are smaller than a specified accuracy of 10 −6 or the iteration count h exceeds 50.
The operations , , and •2 are the matrix element-wise multiplication, division and power. The diag(A 3×1 ) creates a 3 × 3 matrix with the elements of A on the main diagonal.
The operations ⊙, ⊘, and °2 are the matrix element-wise multiplication, division and power. The diag(A3×1) creates a 3 × 3 matrix with the elements of A on the main diagonal.
The flowchart in Figure 1 summaries the steps to follow for the calculation of the power flow voltage magnitude and angle. Steps 1 to 5 are presented in Section 3 with a sample benchmark network. As seen in Figure 1, the iteration is stopped when all elements of ∆ are smaller than a specified accuracy of 10 -6 or the iteration count h exceeds 50.

Network Simulation Using Matrix Formation
A sample network, as shown in Figure 2, is used for model testing in the MATLAB/Simulink environment. The 20-bus network, C1 to C20, is a radial, 3-phase LV distribution system with 400 V line-to-line voltage [27]. This LV distribution network benchmark configuration presents a test environment to study the integration of renewable energy and smart grids.

Network Simulation Using Matrix Formation
A sample network, as shown in Figure 2, is used for model testing in the MATLAB/Simulink environment. The 20-bus network, C1 to C20, is a radial, 3-phase LV distribution system with 400 V line-to-line voltage [27]. This LV distribution network benchmark configuration presents a test environment to study the integration of renewable energy and smart grids. Bus C1 is the swing or slack bus which has a fixed voltage of 1.03 pu for each phase. The MATLAB/Simulink is simulated in phasor mode at 50 Hz. The three types of distribution lines, OH1 to OH3, have the same length of 30 m and their impedances, to , are provided as seen in Table 2 or they can be computed as in [5]. The MATLAB generalized mutual inductance block [36] is used in the simulation to represent the distribution conductors, including the self and mutual resistance and inductance. The netlist of branches or distribution lines connecting all buses is presented in Table 3. In this 20-bus network, there are 19 branches. There is an impedance value between two buses and the admittance is derived by taking the inverse of the impedance matrix.
Bus C1 is the swing or slack bus which has a fixed voltage of 1.03 pu for each phase. The MATLAB/Simulink is simulated in phasor mode at 50 Hz.
The three types of distribution lines, OH1 to OH3, have the same length of 30 m and their impedances, z OH1 to z OH3 , are provided as seen in Table 2 or they can be computed as in [5]. The MATLAB generalized mutual inductance block [36] is used in the simulation to represent the distribution conductors, including the self and mutual resistance and inductance.
The netlist of branches or distribution lines connecting all buses is presented in Table 3. In this 20-bus network, there are 19 branches. There is an impedance value between two buses and the admittance is derived by taking the inverse of the impedance matrix.
Energies 2020, 13, 3011 6 of 14 From the netlist in Table 3, the admittance y ik between bus i and bus k is known. Then, the admittance matrix Y [4] can be calculated. G and B are the conductance and susceptance matrix, respectively, derived from Y whose elements are 3 × 3 matrices.
To demonstrate the matrix-based power flow analysis with the benchmark network, the maximum power point profile P MPP in Figure 3a, which is taken from measured irradiation sensor data, is applied to each of the three PV inverters, PV1 to PV3, which are all connected in a 3-phase connection. The four loads, Load1 to Load4, with single-phase, 2-phase and 3-phase connections, have the profile of each phase with the active power P Load [27] and reactive power Q Load , as seen in Figure 3b. The reactive power Q Load is calculated from P Load to meet the power factor (pf ) of 0.95 lagging. The ratings of the solar power and load profiles are chosen for the system to see the impact of voltage increase and the effectiveness of the control methods which are presented in the following section.
Energies 2020, 13, x FOR PEER REVIEW 6 of 15 From the netlist in Table 3, the admittance between bus i and bus k is known. Then, the admittance matrix [4] can be calculated. and are the conductance and susceptance matrix, respectively, derived from whose elements are 3 × 3 matrices.
To demonstrate the matrix-based power flow analysis with the benchmark network, the maximum power point profile P in Figure 3a, which is taken from measured irradiation sensor data, is applied to each of the three PV inverters, PV1 to PV3, which are all connected in a 3-phase connection. The four loads, Load1 to Load4, with single-phase, 2-phase and 3-phase connections, have the profile of each phase with the active power P [27] and reactive power Q , as seen in Figure 3b. The reactive power Q is calculated from P to meet the power factor (pf) of 0.95 lagging. The ratings of the solar power and load profiles are chosen for the system to see the impact of voltage increase and the effectiveness of the control methods which are presented in the following section. The scheduled active and reactive powers of the network, as well as the initial voltage magnitudes and angles of all buses, are presented in Table 4. The powers at PV buses are adjustable and are determined by the decentralized control method.   The scheduled active and reactive powers of the network, as well as the initial voltage magnitudes and angles of all buses, are presented in Table 4. The powers at PV buses are adjustable and are determined by the decentralized control method.
To investigate the impact of voltage increase, the PV inverters are applied to run at full active power profile, as in Figure 3a without any control. Figure 4 shows the results of per-unit voltages at To investigate the impact of voltage increase, the PV inverters are applied to run at full active power profile, as in Figure 3a without any control. Figure 4 shows the results of per-unit voltages at the three PV connected buses. The voltage levels without applying any control method reach as high as 1.3 pu. The result illustrates the need for control algorithms to keep the voltage within limits during periods of excess generation capacity.

Network Control Application
Two decentralized control strategies are investigated to solve the overvoltage problem and demonstrate the application of the matrix power flow analysis. The first is a simple on/off control for PV inverters when there is a predicted overvoltage. The second method is to apply a multi-objective optimization to get maximum energy from the PV generation, while keeping the voltage under the limit. For this analysis, the voltage in the network is kept to an upper limit of 10% higher than the rated value [37].

Decentralized Overvoltage Control (DOVC)
The decentralized overvoltage control receives P data from each inverter and load data from each node and after power flow analysis simply turns off any inverter where the voltage at the bus is 10% over the rated level. The flowchart of the control is shown in Figure 5a. In this control, only active power is controlled. Initially, the available active power from the PV system is divided equally across the three phases. Then, the power flow analysis is conducted to determine the bus voltage as in Section 2. When there is an overvoltage on a PV bus, the PV power at this bus would be turned off until the voltage is under the limit, when it will be turned on.

Network Control Application
Two decentralized control strategies are investigated to solve the overvoltage problem and demonstrate the application of the matrix power flow analysis. The first is a simple on/off control for PV inverters when there is a predicted overvoltage. The second method is to apply a multi-objective optimization to get maximum energy from the PV generation, while keeping the voltage under the limit. For this analysis, the voltage in the network is kept to an upper limit of 10% higher than the rated value [37].

Decentralized Overvoltage Control (DOVC)
The decentralized overvoltage control receives P MPP data from each inverter and load data from each node and after power flow analysis simply turns off any inverter where the voltage at the bus is 10% over the rated level. The flowchart of the control is shown in Figure 5a. In this control, only active power is controlled. Initially, the available active power from the PV system is divided equally across the three phases. Then, the power flow analysis is conducted to determine the bus voltage as in Section 2. When there is an overvoltage on a PV bus, the PV power at this bus would be turned off until the voltage is under the limit, when it will be turned on.

Decentralized Power Dispatch Control (DPDC)
The flowchart of this control is shown in Figure 5b. The initial steps of the DPDC are the same as the DOVC, the active power is divided the 3-phase equally without reactive power. The controller receives P data from each inverter and load data from each node and a power flow calculation as in Section 2 is taken to get the voltage magnitude values. As seen in the flowchart, when no

Decentralized Power Dispatch Control (DPDC)
The flowchart of this control is shown in Figure 5b. The initial steps of the DPDC are the same as the DOVC, the active power is divided the 3-phase equally without reactive power. The controller receives P MPP data from each inverter and load data from each node and a power flow calculation as in Section 2 is taken to get the voltage magnitude values. As seen in the flowchart, when no overvoltage issues exist, the PV inverters output equally on the three phases, which is similar to the DOVC. In the case of overvoltage, multi-objective optimization is applied. The purpose of the optimization is to find the active and reactive powers of PV inverters to maximize the dispatched active power and minimize the reactive power absorbed while maintaining the limit of the voltage.
The power factor of the PV inverter is kept higher than p f min or 0.9 [38] so that the relation between active and reactive power is shown as Equation (7).
The optimization uses the fmincon function on MATLAB [39]. The construction of the optimization includes the specification of the objective function, boundaries, initial points, and constraints. The objective function is shown in (8). In the system of N PV inverters, the total active power delivered from the PV inverters to their buses is maximized. The reactive power is negative to reduce the voltage level. Thus, there are two objectives of this control, including maximizing active power and minimizing reactive power.
The lower and upper bounds of active and reactive power are shown in (9) and (10), and the linear inequality conditions at PV bus i are shown in Equations (11) and (12).
The non-linear inequality condition is shown in (13). The voltage magnitude |V PV | in (13) is a function of P PV and Q PV , and it is computed using the Newton-Raphson power flow as mentioned in Section 2. This step of power flow calculation for |V PV | is also described in the second block of Figure 5b.
Energies 2020, 13, 3011 9 of 14 The initial or starting points of the active and reactive powers are shown in (14). The active power of each phase is equal to an amount of 50% the initial setting. In this network configuration, the initial points are chosen to satisfy all constraints, especially the voltage limit.
The control sequence is summarized as follows. The power profiles are applied to the loads and PV inverters. There is no reactive power applied for the PV inverters and they operate at unity power factor. The power flow is carried out to calculate the network bus voltages. These bus voltages are constantly checked and when there is no overvoltage, the system keeps running with the pre-set power profiles. When any calculated voltage reaches a limit, the optimization control is applied to the PV inverters. Reactive power is consumed by the PV inverters to reduce the voltage level. When the reactive power reaches the maximum capacity set by the power factor limit, the active power from the inverters is curtailed. The active power curtailment is the last option of the control when the reactive power capacity is at maximum level or power factor is at minimum acceptable value.

Simulation Results
The simulation results for the two decentralized control methods are shown in Figures 6-8. The results are placed together for comparison. As seen in Figure 6, at around 11 a.m., the voltage increases to 1.1 pu and the DOVC turns the PV off and turns it on again at around 4 p.m. The DPDC performs an optimization at each time step and sets reactive and active power levels to keep the voltage at or below 1.1 pu.
Energies 2020, 13, x FOR PEER REVIEW 9 of 15 The non-linear inequality condition is shown in (13). The voltage magnitude | | in (13) is a function of and , and it is computed using the Newton-Raphson power flow as mentioned in Section 2. This step of power flow calculation for | | is also described in the second block of Figure 5b.
The initial or starting points of the active and reactive powers are shown in (14). The active power of each phase is equal to an amount of 50% the initial setting. In this network configuration, the initial points are chosen to satisfy all constraints, especially the voltage limit.
The control sequence is summarized as follows. The power profiles are applied to the loads and PV inverters. There is no reactive power applied for the PV inverters and they operate at unity power factor. The power flow is carried out to calculate the network bus voltages. These bus voltages are constantly checked and when there is no overvoltage, the system keeps running with the pre-set power profiles. When any calculated voltage reaches a limit, the optimization control is applied to the PV inverters. Reactive power is consumed by the PV inverters to reduce the voltage level. When the reactive power reaches the maximum capacity set by the power factor limit, the active power from the inverters is curtailed. The active power curtailment is the last option of the control when the reactive power capacity is at maximum level or power factor is at minimum acceptable value.

Simulation Results
The simulation results for the two decentralized control methods are shown in Figures 6-8. The results are placed together for comparison. As seen in Figure 6, at around 11 a.m., the voltage increases to 1.1 pu and the DOVC turns the PV off and turns it on again at around 4 p.m. The DPDC performs an optimization at each time step and sets reactive and active power levels to keep the voltage at or below 1.1 pu. In Figure 7, the DOVC strategy has PVs turned off and there is no active power delivered from 11 a.m. to 4 p.m. due to voltage constraints. For the DPDC, the algorithm initially responds to overvoltage conditions by absorbing reactive power while the active power keeps close to the P profile. Around noon, the limit on reactive power absorption results in the curtailed active power at PV2 to maintain the voltage level under 1.1 pu.  In Figure 7, the DOVC strategy has PVs turned off and there is no active power delivered from 11 a.m. to 4 p.m. due to voltage constraints. For the DPDC, the algorithm initially responds to overvoltage conditions by absorbing reactive power while the active power keeps close to the P profile. Around noon, the limit on reactive power absorption results in the curtailed active power at PV2 to maintain the voltage level under 1.1 pu.  The reactive power and power factor of DPDC are shown in Figure 8. The reactive power is negative or being absorbed by the PVs and the power factor of the PV inverters is kept above its minimum threshold of 0.9. The maximum energy of the solar power profile over the simulated 24-h period is 1.494 MWh. The total energy delivered using the DOVC is 0.626 MWh or 41.92% of the available total. The energy generated from the DPDC is 1.484 MWh or 99.39% of the maximum available total. The use of control strategies, as outlined, could greatly increase the capacity to host PV systems within a distribution network. Combined with demand-side management and energy storage, the generation peaks In Figure 7, the DOVC strategy has PVs turned off and there is no active power delivered from 11 a.m. to 4 p.m. due to voltage constraints. For the DPDC, the algorithm initially responds to overvoltage conditions by absorbing reactive power while the active power keeps close to the P MPP profile. Around noon, the limit on reactive power absorption results in the curtailed active power at PV2 to maintain the voltage level under 1.1 pu.
The reactive power and power factor of DPDC are shown in Figure 8. The reactive power is negative or being absorbed by the PVs and the power factor of the PV inverters is kept above its minimum threshold of 0.9.
The maximum energy of the solar power profile over the simulated 24-h period is 1.494 MWh. The total energy delivered using the DOVC is 0.626 MWh or 41.92% of the available total. The energy generated from the DPDC is 1.484 MWh or 99.39% of the maximum available total. The use of control strategies, as outlined, could greatly increase the capacity to host PV systems within a distribution network. Combined with demand-side management and energy storage, the generation peaks resulting from high PV penetration could be prevented from disrupting DSO voltage regulation requirements while delivering maximum renewable energy benefits to end-users.

Discussion
The Newton-Raphson method for power flow analysis of unbalanced 3-phase systems in the format of current mismatch in polar coordinates has been presented in detail. After the construction of the power flow formulas in matrix form, they are applied to an unbalanced 3-phase benchmark network. The network operation is simulated over a full day with residential loads and PV inverters in the network. The effect of voltage rise due to reverse power flow from the PV inverters is investigated and control strategies based on ideal decentralized communication are applied for the overvoltage scenarios. Initially, reactive power support levels from each of the PV inverters to reduce the node voltages are controlled by the optimization method. The power factors of the PV inverters are considered when applying the reactive powers of the PV inverters in relation to the active power supplied. The construction of the optimization problem is also presented in the matrix form, which is consistent with the power flow formulas.
The formulas of power flow are reconstructed from the expression of a single-phase into a matrix form where each matrix consists of three phase components. This gives advantages of reducing the complexity of the formulas and this compact matrix form is constructed from the perspective of a programmer. For comparison to the conventional formulas, Equations (1) and (2) are rewritten in Equations (15) and (16). These two are put "side-by-side" with the two equations from Table 1 as seen in Equations (17) and (18) In Equations (17) and (18), the phase index (p = a, b, c) is omitted. This means that Equation (17) presents the current mismatch ∆I r for 3-phase components while in Equation (15), the current ∆I p r presents only a single phase. Furthermore, in Equation (15), there are two sum loops for the k (k = 1 to k = N) and q (q = a, b, c), while there is only one sum loop in Equation (17). The matrix expression, therefore, reduces the program loops and reduces the code verification burden. The Hadamard or element-wise operations , , •2 and diag() are practical in coding.
The matrix presentation is also used in other papers of power flow calculation but there are no complete formulas in the specific Newton-Raphson of the current mismatch in polar coordinates. Moreover, the consistency of using the matrix form for both the power flow and optimization construction is mentioned. There are several papers on power flow analysis methods and these methods are verified by static power flow results [6][7][8][9][10][11][12][13][14][15][16][17][18][19][20][21][22][23][24]. Other research shows optimization controls, over a time period, applied to the network but the power flow calculations are briefly introduced or neglected [29][30][31][32][33]. This paper presents a complete power flow calculation in matrix form, from the construction of the admittance matrix illustrated in an unbalanced benchmark network over a full day and the optimization control strategy in matrix form and effects of the solar PV inverters in the network.
In this paper, the decentralized approach is applied, which ignores practical issues of communication and device delay. It is the first step to further investigation of other control methods and scenarios including battery storage, demand-side management, the effect of communication errors, lack of bus information and fault analysis. The power flow formulas are applied for the constant power P-Q buses with a slack bus. From this base formulation, it is possible to expand the matrix expression for a constant power-voltage P-V bus, transformer with tap change, FACTS and other load types, such as constant impedance Z, constant current I and combined impedance-current-power ZIP with both wye and delta connections. The formula may also be presented in an integrated matrix and complex form. This research is an environment for further investigation of high penetration PV energy into the distribution network and future smart microgrids and optimization of renewable energy management in a stable grid with dynamic changes of residential loads and weather conditions.

Conclusions
The research has presented a detailed formulation and process to perform power flow analysis in an unbalanced 3-phase distribution network using the Newton-Raphson method in the format of current mismatch in polar coordinates. The network is described in a matrix format suited to the development of a programmed solution. A benchmark 20-node network with unbalanced distribution lines, PV inverters, and loads was built using MATLAB/Simulink to verify the matrix model. Power flow analysis on this network for a sunny day showed that node voltages without control could rise to 1.3 pu.
To address this overvoltage issue, two decentralized control methods, using a power flow analysis with the matrix model, were applied to the network. A decentralized overvoltage controller (DOVC) simply turns off any inverter where the voltage at the bus is 10% over the rated level. A decentralized power dispatch controller (DPDC) performs a multi-objective optimization to find the active and reactive powers of PV inverters to maximize the dispatched active power and minimize the reactive power absorbed while maintaining the voltage and power factor within limits. The detailed construction of the optimization problem for the 3-phase network in matrix form which is consistent with the power flow calculation is also described. Simulation with the control methods over a day shows that the total active power of the on/off and optimized controllers deliver 41.92% and 99.39% of the available solar power, respectively, while maintaining the network node voltages within limits. Funding: This research received partial funding from an EU Interreg project whose main purpose is the creation of Energy Community Co-operatives (ECCO).