Fast calculation of the maximum power point of photovoltaic generators under partial shading

This paper presents a method to calculate the energy production of photovoltaic generators considering partial shading or mismatched conditions. The proposed method is based on the complete one-diode model including the bypass diode in its exponential form, where the current and voltage values of the modules composing the photovoltaic panel array are calculated without using the Lambert-W function. In addition, the method introduces a procedure to calculate the vicinity of the maximum power points, which enables the reduction of the operations required to obtain the global maximum. The proposed method provides short simulation times and high accuracy. On the other hand, since the method does not require complex mathematical functions, it can be implemented straightforwardly on known software packages and development languages such as C and C++. Those characteristics make this method a useful tool to evaluate the economic viability and return-of-investment time of photovoltaic installations. Simulation results and comparisons with a classical procedure confirm the good performance of the proposed method in terms of execution time and accuracy.


Introduction
PV systems are commonly analyzed by means of software packages, but due to the complexity of the PV circuits, simulation times are excessively long, which is a significant drawback if large PV systems must be analyzed.In addition, some computational tools do not enable to simulate actual operational conditions such as the mismatching caused by shadows.Some reported methods (Patel & Agarwal, 2008) that aimed to correct partially these drawbacks are based on commercial simulation environments that limit its use due to license issues.Others are based on obtaining systems of non-linear equations, but their solution require using complex functions such as the Lambert-W function, which still produces long processing times (Petrone, Spagnuolo, & Vitelli, 2007).Finally, other methods are based on simplifying the model of the PV modules and bypass diodes, which reduces the simulation time but also reduces the accuracy.On the other hand, most of the reported energy production analysis methods consider the reconstruction of the Current vs. Voltage (I-V) and Power vs. Voltage (P-V) curves on all the voltage range, which introduces long simulation times due to the significant amount of points simulated.Figure 1 presents the typical structure of a commercial PV system rated under 5 kW, in which several PV modules are connected in series (forming a string) in order to provide the voltage level required by commercial PV inverters.Such inverters include an implementation of a Maximum Power Point Tracking (MPPT), grid-connection functions and electrical protections.However, shades covering even a small part of the modules significantly change the Maximum Power Point (MPP) of the PV generator (Herrmann, Wiesner, & Vaanen, 1997).Hence, to accurately predict (i.e.simulate) the power production of such commercial systems it is required to account for partial shading conditions.
(partial shading of the string) or due to faults in some of the PV cells (Garcia, Hernandez & Jurado, 2012), which produce a different electrical characteristic in comparison with the non-damaged cells.
Concerning shading conditions, the main problem is the operation of the PV module as a load.Figures 2a and 2b, present the Voltage vs. Current (V-I) and Power vs. Current (P-I) curves of a ERDM-85 PV (ERDM Solar) module operating under two different solar irradiance conditions, 943 W/m 2 (I ph = 5 A) and 566 W/m 2 (I ph = 3 A), where a polarity change in PV voltages will force the modules to consume power.Such a polarity change condition appears due to the ohmic components of the PV module when the module current is higher than its maximum current production.From the model in Figure 1 it can be observed that when i s > i ph , without the presence of the bypass diode D by , the current in excess of i s over i ph will flow through R h producing a polarity change in v d voltage.The operation of the PV module as a load reduces the life-time of the module (Silvestre & Chouder, 2007), or could even produce a hotspot condition destroying the module (Herrmann, Wiesner, & Vaanen, 1997).In order to avoid such a detrimental operation, a bypass diode is connected in anti-parallel configuration with the module as illustrated in Figure 1 (Hernandez, Garcia & Jurado, 2012).The current in excess of i s over i ph flows through the bypass diode D by and not by R h .The diode activation imposes a small negative voltage to the module, forcing it to consume a small amount of power which is considered acceptable by commercial standards (Orozco-Gutierrez M. , Ramirez-Scarpetta, Spagnuolo, & Ramos-Paja, 2014).Such a condition is observed in Figures 2c and 2d, which present the V-I and P-I curves of the string formed by two ERDM-85 PV modules.The simulation shows that, for string currents lower than 3 A, both modules produce power; while for string currents between 3 A and 5 A only the highly irradiated module produces power.Such a condition puts in evidence the activation of the bypass diode associated to the less irradiated module for i s > 3 A.This paper presents a method to estimate the energy production of commercial PV systems under shading conditions based on a non-simplified system model avoiding the use of the Lambert-W function.The method introduces a procedure to approximate the location of the maximum power points, which reduces the number of evaluation points in order to obtain the global maximum power.By means of this approach, the method provides a significant reduction of the simulation time with high accuracy.The effectiveness of the proposed method is presented by means of simulations and comparisons with a classical procedure.

Mismatching phenomenon and classical model
The mismatching phenomenon occurs when several units of series-connected PV modules, i.e.PV strings, are subjected to different operation conditions.Such mismatched operation occurs due to different irradiation  From the results shown in Figure 2, it is noted that the PV string operates with both modules, or with a single one, depending on the operation conditions.This behavior holds for any number of modules, producing a system with N possible conditions, where N represents the number of PV modules.Such a simulation also shows that each PV module has an optimal operation condition, i.e. the MPP, in which the module produces the maximum power.Hence, if both modules have the same parameters, and operation conditions the MPPs are the same, then the MPP voltage (vMPP) of the string is calculated by multiplying the MPP voltage of any module by N. Similarly, the MPP power (pMPP) of the string is the MPP power of any module multiplied by N (Femia, Petrone, Spagnuolo, & Vitelli, 2012).However, under mismatched conditions the MPPs of the modules are different, as illustrated in Figure 2c and 2d.Moreover, the differences between the modules MPPs produce multiple local MPP, known as LMPP, from which the one with highest power is the global MPP, known as GMPP.The maximum number of LMPPs is equal to N.
Figure 2 illustrates such a multiple-maximum condition.In addition, since the modules produce different voltages for the same current, the MPP current (iMPP) of the string does not match the MPP current of the modules.Therefore, a non-trivial mathematical model is required to calculate the string maximum power.

Mathematical model of a PV module
The single-diode model, presented in Figure 1, is widely adopted in literature to represent the photovoltaic module behavior (Petrone, Spagnuolo, & Vitelli, 2007).In such a model, the current source I ph represents the photo-induced current, the diode D j represents the junction behavior, while R h and R s represent the parallel and series ohmic losses.The junction current i d is given by (1), the parallel current i h and module current i pv are given by (2), the module voltage v pv is given by (3), the bypass diode current i by and the string current i s are given by (4).
In these expressions i sat , v td , R h and R s depend on the PV module construction.Similarly, the parameters of the bypass diode i sat,by and v td,by depend on the device construction and i ph is proportional to the irradiance powering the module.
The main problem of this model concerns the implicit relation of the PV current, which is evident from ( 1)-( 4): i pv depends on v pv , which in turn depends on i pv .Hence, the PV current imposed by a PV voltage, or the PV voltage imposed by a PV current, can be calculated by solving the equations system (1)-( 4) using a numerical method.
Traditionally, the Lambert-W function has been used in order to obtain an explicit solution for such a system as given in ( 5) and ( 6) (Petrone, Spagnuolo, & Vitelli, 2007), where W is the Lambert-W function.
The main drawback of such a solution is the high computational load required to compute the Lambert-W with high accuracy (Veberič, 2012).In fact, the most widely used implementations of the Lambert-W function are available in Matlab and in the GNU Scientific Library GSL (Galassi, et al., 2013).To provide high accuracy, Matlab processes the Lambert-W function using the Symbolic Math Toolbox, which requires long calculation times in comparison with classical numerical solutions based on series or explicit expressions.In contrast, the work presented by Veberič (2012) reports that GSL uses a recursive solver (similar to the Newton-Raphson method) to compute the Lambert-W function, which also requires long processing times in comparison with explicit solutions.

Classical procedure to calculate the string power
Taking into account that PV strings are a series-connection of PV modules, an additional diode is added in series to avoid the injection of negative currents into the string.Such a diode is named blocking diode, as illustrated in Figure 1.
It is noted that the current of the blocking diode, given in (7), is the string current, but this diode exhibits a negative voltage (in relation with the modules voltage), hence it consumes power.
Applying the Kirchhoff laws to the electrical scheme of the PV string in Figure 1 drives to the following relations: the modules and blocking diode currents (i pv and i bk ) are equal to the string current (i s ), while the string voltage (V s ) is equal to the sum of the modules and blocking diode voltages (V pv,i and V bk ).Such electrical relations form the non-linear system F(x) given in ( 8), where the PV current of each module is calculated using the non-linear expressions ( 5) and ( 6).In such a system There are several papers dealing with the solution of (8) (Petrone, Spagnuolo, & Vitelli, 2007; Orozco-Gutierrez, Ramirez-Scarpetta, Spagnuolo, & Ramos-Paja, 2013; Bastidas, Franco, Petrone, Ramos-Paja, & Spagnuolo, 2013), with two main approaches: solving (8) using the exact system Jacobian (Petrone, Spagnuolo, & Vitelli, 2007;Orozco-Gutierrez, Ramirez-Scarpetta, Spagnuolo, & Ramos-Paja, 2013), or simplifying the Jacobian (e.g. the bypass diode equation) to speed-up the calculation (Bastidas, Franco, Petrone, Ramos-Paja, & Spagnuolo, 2013).In both cases a recursive solver is required, e.g.Newton-Raphson or Trust-region.Hence, the solution of (8) requires a large amount of Lambert-W calculations, which in turn require a large amount of time.The solution of ( 8) is based on the Jacobian J of F(x) given in ( 9), which requires the derivative of ( 5) and ( 6) with respect to the PV voltage, involving also the Lambert-W function.Finally, the recursive solver used to find the string current must invert J to approximate, successively, the system solution.

Calculation of the polarization curve without solving non-linear systems or Lambert-W functions
In order to avoid the use of the Lambert-W function to calculate each module voltage and current, i.e. not using ( 5) and ( 6), this paper proposes to perform a sweep over v d to calculate i pv from ( 1) and ( 2), which are explicit equations, i.e. no Lambert-W needed.Then, the module voltage v pv is calculated from (3) and the string current is calculated from (4), again without involving the Lambert-W function.With this procedure all the modules are characterized in terms of the string current without any simplification to the model, using two ordered vectors: PV voltage and string current vectors.However, since the modules are in series, they operate at the same string current i s .Therefore, the modules' electrical characteristics must be interpolated in order to calculate the modules' voltage at the same current i s .
Moreover, the voltage drop introduced by the blocking diode must be also taken into account: since the blocking diode current is explicit (7), the blocking diode voltage is calculated as v bk = v td,bk × ln (1 + i s / i sat,bk ).Finally, all the modules' and blocking diode voltages are added to obtain the string voltage, which enables to calculate the string power.The previous procedure has two main parts: first, the characterization of the modules, which must be done one time for a given irradiance condition S independent of the string current to be tested; and second, the calculation of the string voltage, which depends on the string current.The flowchart in Figure 3 summarizes the string voltage calculation procedure.
Calculating the whole power curve to detect the GMPP, as in Petrone, Spagnuolo, & Vitelli (2007), Orozco-Gutierrez, Ramirez-Scarpetta, Spagnuolo, & Ramos-Paja (2013) and Bastidas, Franco, Petrone, Ramos-Paja, & Spagnuolo (2013), requires executing a single time the first part of the procedure, while the second part must be executed for each string current value.In both cases no Lambert-W functions and recursive solvers are used.Figure 4 presents the performance comparison between the proposed fast method and the classical method, in calculating the complete P-I curve for strings with different number of modules (2 to 20).The results presented in Figure 4a show that the fast method requires processing times between 4 and 5 orders of magnitude shorter.For example, with 2 PV modules the classical method requires 2,47 × 10 6 % more time than the proposed fast method, while with 20 modules the classical method requires 2,03 × 10 7 % more time.In fact, for PV systems up to 20 modules, the proposed fast method never requires more than 1 ms to obtain the solution; instead, the classical method requires from 2,63 minutes (2 modules) up to 4,58 hours (20 modules).In order to validate the method's accuracy, three characteristics were evaluated: string voltages v s and power p s for each current value and the global maximum power p GMPP .Since the PV voltage and power are zero at the short-circuit current, a suitable error formula must be used to avoid divisions by zero.
The Range Average Absolute Error (RAAE) is adopted (Saavedra-Montes, Ramirez-Scarpetta, Ramos-Paja, & Malik, 2007), which normalizes the difference between the signals over the maximum variation range of the reference data to provide a fair comparison between the different operation conditions.The RAAE expression is given in ( 10), where h is the number of data, while y* k and y are the under-tests and reference signal vectors.Figure 4b presents the errors between the fast and classical methods for v s , p s and p GMPP , where the differences between the results of both methods are under 0,0001 %, which is a negligible difference.

RAAE
Then, those results put in evidence the high accuracy and fast processing time provided by the fast solution.
However, the calculation of the p MPP is based on performing a complete sweep to the string power curve, which requires calculating a large number of unnecessary points far from the GMPP.Avoiding such unnecessary calculations will enable to improve, even more, the performance of the proposed algorithm in terms of speed and required memory.Therefore, the following section proposes a method to estimate the neighborhood of the GMPP to reduce the number of P-I points to be calculated.

Estimation of the maximum power point vicinity
Since PV strings are formed by modules connected in series, the current in which the MPP of each module occurs is near to the string currents in which the LMPPs (and GMPP).This condition is due to the large change in the module power curve derivative present at the module MPP as reported in Orozco-Gutierrez, Ramirez-Scarpetta, Spagnuolo, & Ramos-Paja (2013): previous to the MPP the slope is positive and after the MPP the slope is negative having a higher amplitude.Therefore, the power slopes around the module MPP cause a change in the sign of the power slope of the string that produces a LMPP (or GMPP) as depicted in Figure 2: the MPP of each module in Figure 2b produces a MPP, at almost the same current, in the string power curve depicted in Figure 2d.Using that information, this section is devoted to estimate the vicinity of the string LMPPs (and GMPP) currents by estimating the modules' MPP currents.Then, starting from those currents, a hill-climbing algorithm is used to reach the LMPP without evaluating the complete power curve, thus avoiding the evaluation of an unnecessary large number of points.The first approximation used to find the vicinity of the MPP current is to neglect the effect of the series resistance R s and bypass diode, i.e. p pv ≈ p d = i pv • v d as given in (11), which leads to an explicit relation.Then, at the maximum p d value the derivative of p d with respect to v d is zero as given in (12).Hence, Equation ( 12) must be solved for each panel to find the v d  values near the MPP.However, Equation ( 12) is an implicit expression that cannot be solved explicitly even using the Lambert-W function.Therefore, some simplifications must be introduced into (12) to avoid the use of recursive solvers: first, since the MPP voltage of commercial PV modules is commonly larger than 16 V, and the thermal voltages v td of commercial modules is around 1,1 V (Eicker, 2003) Figure 5 presents the simulation of a PV string formed by 3 PV modules with the following parameters: i sat = 1,5415 × 10 -8 A, v td = 1,1088 V, R s = 0,0045 W, R h = 109,495 Ω, v oc,STC = 21,78 V, i sat,by = 1 × 10 -6 A, v td,by = 0,015 V, i sat,bk = 1 × 10 -6 A and v td,bk = 0,15 V. approximated MPP currents obtained with (13), which are very close to the exact MPP currents.With the previous results, the search of each LMPP can be started from the v d,MPP values, which strongly reduces the number of points of the power curve to be calculated in comparison with a full sweep of the curve as performed in traditional approaches.

Calculation of the maximum power with a reduced number of P-I points
Taking into account that the v d,MPP values are near to the LMPP conditions, a simple hill-climbing algorithm is used to detect the LMPP power: increase (or decrease) the string current while the power increases.Figure 6 presents the flowchart of the proposed algorithm, which is divided in modular blocks to evaluate all the LMPP possible conditions.The first block, named initialization block, characterizes the PV modules in the same way as the algorithm presented in Figure 3, where a v d sweep to each module is performed.
The initialization block also calculates the LMPP vicinity in terms of v d,MAX values as described in the previous section.Finally, the control variables N (number of modules), k (LMPP in evaluation) and p GMPP (maximum power) are initialized.
After detecting the zones of the LMPPs, the algorithm evaluates each zone to calculate the corresponding LMPP with the LMPP block, where the string current is calculated in the same way previously proposed in the flowchart of Figure 3: the modules' voltages are interpolated at the same string current value.This process starts at the approximated LMPP currents I MPP,approx, increasing the string current while the power increases.However, in the case that the first iteration produces a power reduction, it means the LMPP is at a lower current, hence the current is decreased while the power increases.
When the LMPP is detected it is contrasted with the other LMPPs to update the GMPP.This process is repeated until all the LMPPs are evaluated.To illustrate this procedure, a PV string formed with four modules is simulated, where the parameters are the same ones described in the previous section, but the irradiances of the modules are 940 W / m 2 , 600 W / m 2 , 400 W / m 2 , and 200 W / m 2 .Figure 7 illustrates the number of P-I points calculated using the sweep and proposed (reduced) methods, where the reduced number of calculations can be appreciated.
In fact, for the string made of four modules the sweep method calculates 950 P-I points, while the reduced method only calculates 58 points.Since the number of P-I points calculated by the reduced method depends on the approximation of I MPP,approx, multiple simulations were performed to provide an average relation between the number of points calculated by both methods to reach the same result: the reduced method calculates an average of 6,39 % of the P-I points calculated by the sweep method.Therefore, the proposed GMPP calculation algorithm, i.e. Figure 6, improves significantly the performance in comparison with classical solutions based on implicit For the sake of simplicity, the three modules have the same parameters, which does not introduce any distortion to the calculation results.To force the mismatching condition, the three modules are considered under different irradiance levels (i.e.partial shading): 500 W/m 2 , 400 W/m 2 and 200 W/ m 2 , hence three LMPP exist.The simulation also shows the expressions for i pv recursive solvers and sweeps for calculating the string power.
to process a 20 modules PV system, the proposed method requires less than 1 ms.On the other hand, the comparison between a classical procedure and the proposed method results in a RAAE lower than 0,0001 %, which is a negligible difference.Since the proposed method enables to obtain the GMPP in a fast and accurate way, large PV systems can be analyzed in terms of economic viability accounting for partial shading conditions, which was not possible with classical methods.Such analysis could include Monte Carlo simulations to quantify the impact of shading on the power production of the PV installation to select the best location for the PV array.In addition, this method can be used in the dynamic reconfiguration of the PV system to reduce the impact of partial shading in the power production.Finally, the simplicity of the method makes it suitable for implementation in software packages like Matlab and development languages such as C and C++.

Conclusions
The method proposed in this paper is a suitable tool to estimate the energy production of commercial PV systems.This method significantly reduces the use of the Lambert-W function and only calculates the string power within the neighborhood of the LMPPs.Therefore, this proposed method provides a significant reduction in the simulation time: while a classical method requires hours

Figure 1 .
Figure 1.Structure of commercial single-string photovoltaic systems.

Figure 2 .
Figure 2. Electrical behavior of PV strings under mismatched conditions.

Figure 3 .
Figure 3. Fast procedure to calculate the string voltage.

Figure 4 .
Figure 4. Electrical behavior of PV strings under mismatched conditions.

Figure 5 .
Figure 5. Approximation of v d values at the LMPPs.

Figure 6 .
Figure 6.Complete algorithm for calculating the maximum power.

Figure 7 .
Figure 7. Calculation of the GMPP with a reduced number of P-I points.