Estimation of Single-Diode Photovoltaic Model Using the Differential Evolution Algorithm with Adaptive Boundaries

Abstract

This study proposes a calculation methodology that determines the optimal boundary parameters of the single-diode photovoltaic model. It allows the calculation of the single-diode photovoltaic model when no reference parameter boundaries are available. The differential evolution algorithm, integrated with a step-by-step boundary definition module, is used to calculate the optimal parameters of the single-diode photovoltaic model, improving the performance of the classic algorithm compared with other studies. The solution is validated by comparing the results with well-established algorithms described in the state-of-the-art, and by estimating the five important points (cardinal points) of an IV curve, namely short-circuit, maximum power, and open circuit points, using a database composed of 100 solar photovoltaic modules. The results show that an optimal set of parameter boundaries enables the differential evolution algorithm to minimize the error of the estimated cardinal points. Moreover, the proposed calculus methodology is capable of producing high-performance response photovoltaic models for different technologies and rated powers.

1. Introduction

Different mono-facial solar photovoltaic module SPVM technologies have been developed e.g., mono-crystalline, poly-crystalline, and thin layer, among others. Their electrical behavior is usually described by the curve traced in the current and voltage plane (IV curve). Similarly to the mono-facial technology, SPVMs have recently incorporated bifacial technology, and unlike their mono-facial counterpart, bifacial SPVMs are capable of harnessing incident solar radiation from both faces of the photovoltaic PV module, increasing its production capacity [1]. Authors such as Singh et al. [2] and Liang et al. [3] claim that the behavior of bifacial solar technology can be described using the same method as mono-facial modules. In addition, Mujahed Al-Dhaifallah et al. [4] developed a fractional control system capable of tracking the maximum power point, which highly depends on the electrical model of the SPVM. These applications require modeling the electrical behavior of the SPVM with a high level of accuracy in order to improve their performance.

SPVM have a highly non-linear conduct, ranging from the ideal model behavior being described by three parameters to more complex models that use five or even seven parameters. Increasing the number of parameters included in the SPVMs electrical model increases the difficulty of obtaining these parameters; therefore, the electrical model must be carefully selected [5,6,7,8]. In many studies, the five parameter model displays sufficient accuracy to describe the operation of the SPVM, as it is capable of predicting the behavior of the module under different operating scenarios with high precision [9,10], in addition it includes the transition to different operating conditions starting from a particular state. The transition between two operating conditions was subsequently refined by Boyd et al. [11]. Different methodologies are used to determine the input parameters of the models, and the currently available tools are the explicit methods through Lambert’s W function [12], the iterative method, the heuristic method, and the metaheuristic approaches.

The calculation strategy is selected according to two main factors: (i) the input data, and (ii) the SPVM electrical model. On the one hand, the input data may be obtained from an IV curve, or from the SPVMs datasheet. The SPVMs datasheet is widely used as input parameter as all PV modules manufacturers provide this information. However, the datasheet data are limited to three PV operational points, i.e., short circuit, open circuit, maximum power point of both the standard test condition STC and normalized operating cell temperature NOCT. On the other hand, the electrical model is selected according to the required precision. Explicit and iterative methods usually use approximations to obtain the model parameters, while metaheuristic methods are based on stochastics and do not require any approximations to determine the output parameters of the selected model.

Heuristics refers to solving a problem subject to a series of restrictions. These restrictions allow determining an approximation of the solution with a higher computational speed when compared to conventional algorithms that calculate the exact solution. An improvement to this type of algorithm is the metaheuristic algorithm, based on a focused search of the solution. An example of this method is indicated by Elazab et al. [13], who present a Whale Optimization Algorithm WOA based on the whale hunting method, converging inwards as a spiral towards the approximate solution. Another metaheuristic algorithm is proposed by Ebrahimi et al. [14], which uses a Flexible Particle Swarm Optimization Algorithm FPSOA inspired by the social behavior of birds. The behavior of the FPSOA is based on the experience of each particle and of its neighbors. Another type of metaheuristic algorithm is based on the evolution of a given population. The metaheuristic algorithms are widely used in the industry to solve problems such as wire electric discharge machining [15], signal and image processing [16], economic dispatch [17], and spatial forest planning [18]; therefore, due to the versatility of this type of algorithm, the metaheuristic approach is selected to model the SPVM electrical behavior.

Ishaque et al. [19] developed a method to obtain a five parameter model based on the metaheuristic differential evolution algorithm DEA using only the information provided by a SPVMs datasheet. On the other hand, Biswas et al. [20] slightly modifies the standard DEA, obtaining equivalent five and seven parameter models by varying the objective function and introducing a dynamic population on a step-by-step basis. The objective function indicated by Biswas et al., provides a range of potential solutions. Such diverse solutions can be explained by the fact that the study only considers three operational points located on curve IV, leaving aside the behavior of the PV curve. As both studies are supported by the same base algorithm, they require that the bounds upon which the solution resides are defined. Ishaque et al., relies on the data provided by Villalva et al. [21], while Biswas et al., does not make any reference to the bounds used to determine these parameters. The information and methodology provided by these studies is limited, and the criteria for the selected limits used in the algorithms indicated by the state-of-the-art is studied by some authors like Toledo et al. [22] using an analytical approach. However, the boundaries calculus method addressed by Toledo et al., presents inconsistencies for certain SPVMs since it is not capable of providing a finite value for the maximum limit of the shunt resistance.

On the other hand, Abido et al. [23] uses the DEA as a base to determine a five parameter model that incorporates a correction on the thermal conversion of two of the five parameters, providing more precision to the model when changing from one operating condition to another. However, there are no specific criteria regarding the selection of the boundaries of the mentioned parameters. This is relevant as the quality of the solution is associated with the input boundaries provided to the algorithm. In addition, De Soto et al. [9] indicates that the value of these parameters is modified according to the operating condition at which they are calculated, and the PV specimen under study. This suggests a problem that must be solved. Several authors such as Villalva et al., and Yan et al. [24] experimentally determined bounds over which it is possible to find a solution for the PV equations; however, as mentioned above, these parameters are subject to specific operating conditions. Considering this background, as the DEA has shown adequate results and a high adaptability to solve different problems, this work proposes and autonomous DEA capable of determining electrical models of future solar photovoltaic technologies, such as bifacial or perovskite solar cells. Furthermore, this methodology can be used to estimate the electrical model of a solar power plant.

This study presents a methodology that enables the calculation of the SDM parameter boundaries. The SDM is calculated by means of the differential evolution algorithm. However, in order to improve the base DEA algorithm, an adaptive step-by-step limit module is incorporated. The calculation method improves the performance of the solution achieved by other authors, allowing the calculation of the SDM when no reference parameter boundaries are available. Two different operational conditions (for example, standard test conditions and normal operating cell temperature) are required as input data in order to calculate the SDM which minimizes the cardinal points error. These data are usually provided by the manufacturer on the SPVM datasheet.

The proposed parameter boundaries calculation method is studied in detail using the KC200GT SPVM. The performance of the proposed methodology is validated estimating the cardinal points for two and comparing with (i) results of with well-established algorithms for the Kyocera KC200GT SPVM, and (ii) 100 SPVM datasheets of different technologies and rated powers. The remaining publication is divided as follows: Section 2 indicates the SPVM electric model, Section 3 indicates the boundaries calculation method, Section 4 explains the differential evolution algorithm used, Section 5 presents the results of the simulations and the validation of the model and, finally, Section 6 presents the conclusions.

2. Photovoltaic Model

Given the constructive nature of SPVMs, the electrical model of the SPVM is characterized by having at least one ideal diode; therefore, its behavior is explained based on semiconductor theory [25,26,27].

2.1. Ideal Model

The ideal SPVM model relates the behavior of the PV current $I_{pv}$ and PV voltage $V_{pv}$ as a function of three parameters: the photovoltaic current source $I_{ph}$, the saturation current of the diode $I_o$, and the diode’s ideality factor A. Equations (1)–(5) indicate the relationships between the different variables present in the ideal PV model, indicated in Figure 1, where $V_D$ and $I_D$ represent the voltage and current of the diode. The thermal voltage of the diode $V_t$ is a function of Boltzmann’s constant $k_B$, the elemental charge of the electron q and the PV cell temperature T. However, SPVMs often have $N_s$ PV cells connected in series, so it is convenient to introduce the equivalent thermal voltage parameter, a, that concentrates the information of the thermal voltage, the number of cells connected in series, and the diode ideality factor. Usually, the temperature is measured in K and the Boltzmann constant is expressed in eV.

$$\begin{array}{ccc}\hfill V_D& =& V_{pv}\hfill \end{array}$$

(1)

$$\begin{array}{ccc}\hfill V_t& =& \frac{k_B·T}{q}\hfill \end{array}$$

(2)

$$\begin{array}{ccc}\hfill a& =& A·N_s·V_t\hfill \end{array}$$

(3)

$$\begin{array}{ccc}\hfill I_D& =& I_o·\{e^{\frac{V_D}{a}}-1\}\hfill \end{array}$$

(4)

$$\begin{array}{ccc}\hfill I_{ph}& =& I_D+I_{pv}\hfill \end{array}$$

(5)

The ideal model parameters can be determined according to the operation points indicated by the manufacturer of the SPVM, and these are called PV cardinal points: short-circuit current $I_{sc}$, open circuit voltage $V_{oc}$, maximum power current $I_{mpp}$, maximum power voltage $V_{mpp}$, and maximum power point $P_{mpp}$. These values are usually indicated for an operational condition called standard test condition STC, in which the irradiance level $S_{ref}$ and temperature $T_{ref}$ are 1000 W/m${}^2$ and 298.15 K. By replacing these parameters on Equations (1)–(5), the three parameters can be calculated by solving the system of non-linear equations indicated in Equations (6)–(8).

$$\begin{array}{ccc}\hfill I_{ph}& =& I_{sc}\hfill \end{array}$$

(6)

$$\begin{array}{ccc}\hfill I_{ph}& =& I_o·\{e^{\frac{V_{oc}}{a}}-1\}\hfill \end{array}$$

(7)

$$\begin{array}{ccc}\hfill I_{ph}& =& I_o·\{e^{\frac{V_{mpp}}{a}}-1\}+I_{mpp}\hfill \end{array}$$

(8)

2.2. Single-Diode Model

The single-diode model described by five parameters (SDM) is commonly used as it is able to express, with precision, the behavior of the SPVM under different operating conditions. This model is defined as a function of five parameters: the photovoltaic current source $I_{ph}$, the diode saturation current $I_o$, the equivalent thermal voltage a, the series resistance $R_s$, and the shunt resistance $R_{sh}$. The equations that govern the behavior of the SPVM are given in Equations (2)–(4) and (9)–(11), where the shunt current $I_{sh}$ is determined as a function of $V_D$ and $R_{sh}$, while the diode current $I_D$ remains as defined in Equation (4). The electrical model is presented in Figure 2.

$$\begin{array}{ccc}\hfill V_D& =& V_{pv}+R_s·I_{pv}\hfill \end{array}$$

(9)

$$\begin{array}{ccc}\hfill I_{sh}& =& \frac{V_D}{R_{sh}}\hfill \end{array}$$

(10)

$$\begin{array}{ccc}\hfill I_{ph}& =& I_D+I_{sh}+I_{pv}\hfill \end{array}$$

(11)

Equations (12)–(14) indicate the behavior of the five parameter model at the PV cardinal points. Unlike the ideal case, it is not possible to solve the equation system, as the number of variables exceeds the available equations. However, it is possible to add a condition which indicates the maximum PV power point, mathematically described in Equation (15), where the PV power function is dependent on the PV voltage. This additional condition requires knowledge of the derivative of $I_{pv}$ with respect to $V_{pv}$, given in Equation (16).

$$\begin{array}{ccc}\hfill I_{ph}& =& I_o·\{e^{\frac{R_s·I_{sc}}{a}}-1\}+\frac{R_s}{R_{sh}}·I_{sc}+I_{sc}\hfill \end{array}$$

(12)

$$\begin{array}{ccc}\hfill I_{ph}& =& I_o·\{e^{\frac{V_{oc}}{a}}-1\}+\frac{1}{R_{sh}}·V_{oc}\hfill \end{array}$$

(13)

$$I_{ph}=I_o·\{e^{\frac{V_{mpp}+R_s·I_{mpp}}{a}}-1\}+\frac{V_{mpp}+R_s·I_{mpp}}{R_{sh}}+I_{mpp}$$

(14)

$$\frac{\mathsf{\partial }{P}_{pv}}{\mathsf{\partial }{V}_{pv}}{|}_{mpp}=I_{mpp}+V_{mpp}·\frac{\mathsf{\partial }{I}_{pv}}{\mathsf{\partial }{V}_{pv}}{|}_{mpp}=0$$

(15)

$$\frac{\mathsf{\partial }{I}_{pv}}{\mathsf{\partial }{V}_{pv}}{|}_{mpp}=-\frac{R_{sh}·I_o·e^{\frac{V_{mpp}+R_s·I_{mpp}}{a}}+a}{R_s·\{R_{sh}·I_o·e^{\frac{V_{mpp}+R_s·I_{mpp}}{a}}+a\}+R_{sh}·a}$$

(16)

The fifth equation required to solve the system is introduced by several authors as an approximation for a given operational point, e.g., De Soto et al., indicates that it is possible to estimate the temperature coefficient for open circuit ${\beta }_{V_{oc}}$ conditions. On the other hand, Jadli et al. [28] assumes that the fifth condition is the slope of the PV current, determined according to the PV voltage in the short-circuit condition.

Although the system of equations indicated in Equations (12)–(15) is not sufficient to determine the five parameters of the model, it is enough to determine the cardinal points. The solution of this system of equations is obtained through multiple optimization techniques. However, in this study, the explicit expressions provided by Lambert’s W function $W\left(z\right)$ are used [29].

In order to determine the behavior of the IV curve under different operating conditions, it is necessary to transform the SDM from the reference condition to the desired condition. Such a method is exposed by De Soto et al., in Equations (17)–(21), where the sub index ref and op indicates the reference and the desired operational condition, respectively. ${\alpha }_{I_{sc,ref}}$ indicates the short-circuit current temperature coefficient at a known reference irradiance $S_{ref}$. The band gap energy $E_g$ is calculated using the Varshni equation applied to silicon materials, indicated in Equation (22) [30,31]. The parameters used by the Varshni equation are extracted from [32].

$$\begin{array}{ccc}\hfill I_{ph,op}& =& \frac{S_{op}}{S_{ref}}·\left[I_{ph,ref}+{\alpha }_{I_{sc,ref}}·(T_{c,op}-T_{c,ref})\right]\hfill \end{array}$$

(17)

$$\begin{array}{ccc}\hfill I_{o,op}& =& I_{o,ref}·{\left[\frac{T_{c,op}}{T_{c,ref}}\right]}^3·exp\left[\frac{1}{k_B}\left(\frac{E_{g,ref}}{T_{c,ref}}-\frac{E_{g,op}}{T_{c,op}}\right)\right]\hfill \end{array}$$

(18)

$$\begin{array}{ccc}\hfill a_{op}& =& \frac{T_{c,op}}{T_{c,ref}}·a_{ref}\hfill \end{array}$$

(19)

$$\begin{array}{ccc}\hfill R_{sh,op}& =& \frac{S_{ref}}{S_{op}}·R_{sh,ref}\hfill \end{array}$$

(20)

$$\begin{array}{ccc}\hfill R_{s,op}& =& R_{s,ref}\hfill \end{array}$$

(21)

$$E_g=1.166-\frac{4.73·{10}^{-4}·T_c^2}{636+T_c}$$

(22)

3. Parameter Boundaries Definition

This section proposes a methodology for the calculation of the boundaries on which the SDM parameters reside. However, it is only able to determine the boundaries for the following variables: (i) the equivalent thermal voltage, (ii) the series resistance, and (iii) the shunt resistance, since the remain parameters can be calculated according to Equations (23) and (24) [20].

$$\begin{array}{ccc}\hfill I_o& =& \frac{I_{sc}+\frac{R_s·I_{sc}}{R_{sh}-\frac{V_{oc}}{R_{sh}}}}{e^{\frac{V_{oc}}{a}}-e^{\frac{R_s·I_{sc}}{a}}}\hfill \end{array}$$

(23)

$$\begin{array}{ccc}\hfill I_{ph}& =& I_o·(e^{\frac{V_{oc}}{a}}-1)+\frac{V_{oc}}{R_{sh}}\hfill \end{array}$$

(24)

3.1. Equivalent Thermal Voltage Boundaries

The equivalent thermal voltage is directly correlated to the ideality factor of the diode in the electrical model. The impact of a on the PV power curve is displayed in Figure 3, and the shape of the curve is not significantly altered with respect to the ideal curve. However, as the value of a increases, the power and voltage of the PV curve increase in a near direct proportional manner with respect to a.

The range of the equivalent thermal voltage is determined by the values of the diodes ideality factor A, which usually ranges from 1 to 2. Therefore, the limits of a can be calculated by applying Equations (25) and (26). It should be noted that the calculation of a requires knowledge of the number of cells connected in series $N_s$ and the cell temperature $T_c$ for which the cardinal points were obtained.

$$\begin{array}{ccc}\hfill a_{min}& =& 1·N_s·V_t=1·N_s·\frac{k_B·T_c}{q}\hfill \end{array}$$

(25)

$$\begin{array}{ccc}\hfill a_{max}& =& 2·N_s·V_t=2·N_s·\frac{k_B·T_c}{q}\hfill \end{array}$$

(26)

3.2. Series Resistance Boundaries

The series resistance indicates the resistance of the cells material and resistive contacts that oppose the current flow [25]. The value of the series resistance has a significant impact on the output power between the open circuit voltage and maximum power point. As the magnitude of $R_s$ increases, the voltage produced by the SPVM tends to zero. Therefore, the available power that the SPVM can deliver to an external load is limited, reaching zero when $R_s$ tends to infinite. This translates into a degradation of the PV curve. The effect of the variation of $R_s$ on the PV curve is indicated in Figure 4, where $R_s$ varies between 0.5 and 50 $\mathsf{\Omega }$.

In order to calculate boundaries for $R_s$, geometric limitations are incorporated as photovoltaic power restrictions. The upper limit of $R_s$ is calculated by using the slope of the line that connects the maximum power point $(V_{mpp},I_{mpp})$ and the open circuit voltage $V_{oc}$, as indicated in Equation (27). On the other hand, the minimum value of $R_s$ corresponds to zero; however, it is useful to set an infimum value for series resistance $R_{s,inf}$ as different from zero. For this purpose, a detailed study on the variation of the maximum power point in function of the series resistance ${\mathsf{\Phi }}_s\left(R_s\right)$ is performed; nevertheless, the variation of the maximum power point is presented as the deviation of the ideal case as indicated in Equation (28). Here, the maximum power is calculated using the ideal photovoltaic model considering a series resistance values different from zero.

Figure 5 indicates the behavior of ${\mathsf{\Phi }}_s\left(R_s\right)$ for a test SPVM. In this figure, a minimum value of ${\mathsf{\Phi }}_s\left(R_s\right)$ (referred to as $\phi$) is defined in the $\left[0,{\mathsf{\Phi }}_s\left(R_{s,max}\right)\right]$ range. By setting $\phi$, it is possible to calculate $R_{s,inf}$ by means of Equation (29). This equation can be solved using an optimization tool, e.g., the midpoint method, as the the $R_{s,inf}$ bounds are known (defined between 0 and $R_{s,max}$). For the test module indicated in Figure 5, the maximum bound calculated is 867.20 m$\mathsf{\Omega }$, delivering a ${\mathsf{\Phi }}_s\left(R_{s,max}\right)$ value of 23.68%. Assuming that $\phi \ll {\mathsf{\Phi }}_s\left(R_{s,max}\right)$, according to Equation (30), a power margin ${\Delta }_P\approx {\mathsf{\Phi }}_s\left(R_{s,max}\right)=23.68$% is achieved.

$$R_{s,max}=\frac{V_{oc}-V_{mpp}}{I_{mpp}}$$

(27)

$${\mathsf{\Phi }}_s\left(R_s\right)=100\times \frac{P_{max,ideal}-P_{max}\left(R_s\right)}{P_{max,ideal}}\%$$

(28)

$$I_{mpp}\left(R_s\right)\times V_{mpp}\left(R_s\right)-P_{max,ideal}\times \left(1-\frac{\phi }{100}\right)=0$$

(29)

$${\Delta }_P=\frac{100}{P_{max,ideal}}\times (P_{max,R_{s,inf}}-P_{max,R_{s,max}})={\mathsf{\Phi }}_s\left(R_{s,max}\right)-\phi$$

(30)

3.3. Shunt Resistance Boundaries

The shunt resistance indicates the leakage of current through the cell around the edges of the device and between contacts of different polarities [25]. The shunt resistance $R_{sh}$ significantly impacts the PV curve as its value tends to 0 $\mathsf{\Omega }$, whereas for high values, the impact on the curve is minimal. This is because small values of $R_{sh}$ tend to short-circuit the current source of the five parameter model and, thus, the power is not able to leave the SPVM. The situation described is presented in Figure 6, where it can be appreciated that values in the order of 1 $\mathsf{\Omega }$ degrade the maximum power of the PV curves by more than 100 W (with respect to the ideal PV model).

A similar calculation method to the one presented for series resistance can be applied to determine the shunt resistance. Since a minimum value of the shunt resistance $R_{sh,min}$ is known, a finite maximum value needs to be defined. Then, a supremum value for the shunt resistance $R_{sh,sup}$ is calculated. In this case, Equation (31) presents the minimum value of $R_{sh}$ calculated using the short-circuit current $I_{sc}$ and the maximum power point $(I_{mpp},V_{mpp})$. To calculate $R_{sh,sup}$, the variation of the maximum power point ${\mathsf{\Phi }}_{sh}\left(R_{sh}\right)$ is studied as indicated in Equation (32). Here, the maximum power is calculated using the ideal photovoltaic model, considering a shunt resistance different from zero.

Figure 7 indicates the behavior of ${\mathsf{\Phi }}_{sh}\left(R_{sh}\right)$ for a test SPVM. In a similar manner to the $R_s$ scenario, $\phi$ indicates the minimum value expected for ${\mathsf{\Phi }}_{sh}\left(R_{sh}\right)$; however, a supremum value $R_{sh,sup}$ is achieved. $R_{sh,sup}$ can be calculated solving the Equation (33). For the test SPVM indicated in Figure 7, the value of the minimum shunt resistance is equal to 43.83 $\mathsf{\Omega }$, delivering a ${\mathsf{\Phi }}_{sh}\left(R_{sh,min}\right)$ value of 7.5%. Assuming that $\phi \ll {\mathsf{\Phi }}_{sh}\left(R_{sh,min}\right)$, according to Equation (34), a power margin ${\Delta }_{sh}\approx {\mathsf{\Phi }}_{sh}\left(R_{sh,min}\right)=7.5$% is achieved.

$$R_{sh,max}=\frac{V_{mpp}}{I_{sc}-I_{mpp}}$$

(31)

$${\mathsf{\Phi }}_{sh}\left(R_{sh}\right)=100\times \frac{P_{max,ideal}-P_{max}\left(R_{sh}\right)}{P_{max,ideal}}\%$$

(32)

$$I_{mpp}\left(R_{sh}\right)\times V_{mpp}\left(R_{sh}\right)-P_{max,ideal}\times \left(1-\frac{\phi }{100}\right)=0$$

(33)

$${\Delta }_P=\frac{100}{P_{max,ideal}}\times (P_{max,R_{sh,sup}}-P_{max,R_{sh,min}})={\mathsf{\Phi }}_s\left(R_{sh,sup}\right)-\phi$$

(34)

4. Differential Evolution Algorithm

The differential evolution algorithm DEA proposed by Price et al. [33], uses a solution-seeking mechanism with a mutation operation over a user-defined boundary space in a specific number of generations G. Each generation created by the algorithm has a population $P_G$ of fixed magnitude $NP$. Here, every individual item of the population ${\chi }_p$, is composed by D elements $x_d$, as indicated in Equations (35) and (36). The generations evolve until the $G_{max}$ iterations are reached.

In order to calculate the single-diode photovoltaic model, the following variables are used: equivalent thermal voltage, series resistance, and shunt resistance; therefore, D equals 3. The photovoltaic current and saturation current of the diode are calculated according to Equations (23) and (24).

$$\begin{array}{ccc}\hfill P_G=\left[{\chi }_1,\dots ,{\chi }_p,\dots ,{\chi }_{NP}\right]& & G=1,2,\dots ,G_{max}\hfill \end{array}$$

(35)

$$\begin{array}{ccc}\hfill {\chi }_p=\left[x_1,\dots ,x_d,\dots ,x_D\right]& & p=1,2,\dots ,NP\hfill \end{array}$$

(36)

4.1. Initialization

The algorithm begins with a population of magnitude $NP$, where $NP\in \mathbb{N}$. Each element of the population has a known dimension D and is constructed within lower and upper boundaries, ${\chi }_L$ and ${\chi }_H$. It is important to emphasize that the lower and upper boundaries have D dimension, as well as each element of the generated population. This is indicated in Equations (37) and (38). Once the boundaries are defined, a random vector ${\chi }_{rng}$ of dimension D is generated in the [0,1] domain, in order to determine the elements of the initial population, according to Equation (39). This scenario is illustrated in Figure 8a, where seven 2D random elements are generated within the boundaries, and represented by a dotted line.

$$\begin{array}{ccc}\hfill {\chi }_H& =& \left[x_{1,H},\dots ,x_{d,H},\dots ,x_{D,H}\right]\hfill \end{array}$$

(37)

$$\begin{array}{ccc}\hfill {\chi }_L& =& \left[x_{1,L},\dots ,x_{d,L},\dots ,x_{D,L}\right]\hfill \end{array}$$

(38)

$${\chi }_p={\chi }_L+{\chi }_{rng}·({\chi }_H-{\chi }_L)$$

(39)

4.2. Mutation

The mutation process creates a mutant vector $V_G$ that is used to modify $P_G$. This vector is defined as follows: three different elements belonging to $P_G$ are randomly selected: ${\chi }_{r1}$, ${\chi }_{r2}$, and ${\chi }_{r3}$. Then, $V_G$ is calculated subtracting the difference between vectors ${\chi }_{r2}$ and ${\chi }_{r3}$. This difference is then weighted by a mutant factor F and added to vector ${\chi }_{r1}$, as indicated in Equation (40). Note that F corresponds to a single number that lies in the range of $\left[0.1\right]$. Considering that $V_G$ is calculated from vectors of dimension D, it must also be composed of D elements $v_G$, as indicated by Equation (41). This calculus methodology is illustrated in Figure 8b, where ${\chi }_{r1}$, ${\chi }_{r2}$, and ${\chi }_{r3}$, are represented by vectors (9), (2), and (4). In this example, the mutant vector begins in the origin of the reference system and ends in the red dot.

$$\begin{array}{ccc}\hfill V_G& =& {\chi }_{r1}+F·({\chi }_{r2}-{\chi }_{r3})\hfill \end{array}$$

(40)

$$\begin{array}{ccc}\hfill V_G& =& \left[v_1,\dots ,v_d,\dots ,v_D\right]G=1,2,\dots ,G_{max}\hfill \end{array}$$

(41)

4.3. Crossover

Using the information provided by the mutant vector, it is possible to construct a new population ${\rho }_G$. The elements of ${\rho }_G$ have the possibility of replacing elements of the current population, in order to improve the overall performance. The population of ${\rho }_G$ is referred to as trial vectors ${\psi }_p$, and is composed by D elements $y_d$, as indicated in Equations (42) and (43). $y_d$ elements are chosen according to the binomial crossover strategy, which depends on three factors: (i) the crossover rate $CR$, (ii) the crossover random number $C{R}_{d,rng}$, and (iii) a random value $d_{rng}$, as indicated in Equation (44). Here, $CR$ remains the same for all generations, and usually is defined as a single number in the [0, 1] domain. On the other hand, each $y_d$ element has its own $C{R}_{d,rng}$ number in the [0, 1] domain. The third factor ensures that at least one chromosome of the mutant vector lives on the new population, since $d_{rng}$ is a natural number randomly selected from the [1, D] domain. However, it is possible that some of the $y_d$ elements breach the bounds set by $x_{d,L}$ and $x_{d,H}$; therefore, these values must be replaced. This is expressed in Equation (45), and referred to as the penalty function. Here, a new random number $y_{d,rng}$ is generated to calculate the value of the new element $y_d$.

$$\begin{array}{ccc}\hfill {\rho }_G=\left[{\psi }_1,\dots ,{\psi }_p,\dots ,{\psi }_{NP}\right]& & G=1,2,\dots ,G_{max}\hfill \end{array}$$

(42)

$$\begin{array}{ccc}\hfill {\psi }_p=\left[y_1,\dots ,y_d,\dots ,y_D\right]& & p=1,2,\dots ,NP\hfill \end{array}$$

(43)

$$y_d=\left\{\begin{array}{cc}{v}_d,\hfill & \mathit{if}\left(C{R}_{d,rng}\le CR\mathit{or}d=d_{\mathit{rng}}\right)\hfill \\ x_d,\hfill & \mathit{other}\mathit{case}\hfill \end{array}\right.$$

(44)

$$y_d=x_{d,L}+y_{d,rng}·\left(x_{d,H}-x_{d,L}\right)\mathit{if}(y_d<x_{d,L})\mathit{or}(y_d>x_{d,H})$$

(45)

4.4. Evaluation and Selection

The screening process for new generations $\left(G+1\right)$ is based on comparing the value of the objective function $J\left({\chi }_p\right)$ between G and $\left(G+1\right)$, maintaining the solution that minimizes the $J\left({\chi }_p\right)$ value. Ishaque et al., suggests replacing the ${\chi }_p$ element by the new ${\psi }_p$ element as long as this new element presents a better value for the objective function, i.e., the new elements will always provide a lower or at most equal value of $J\left({\chi }_p\right)$. Unlike the methodology proposed in [19], this paper adopts the method of replacing the element that produces the largest $J\left({\chi }_p\right)$ value of the entire population, as indicated in Equation (46). It can also be seen that the solution provided by the proposed algorithm continuously improves or, in the worst case, remains the same.

$J\left({\chi }_p\right)$ is calculated from the mean value of two errors: (i) the cardinal points from the first reference condition ${\xi }_{1,k}$, and (ii) the cardinal points of the second reference condition ${\xi }_{2,k}$, as indicated in Equation (47). Vectors ${\xi }_{p,k}\left({\chi }_p\right)$ and ${\xi }_{p,k}^{\prime }\left({\chi }_p\right)$ indicate the calculated cardinal points of the first and second reference conditions. It should be noted that the cardinal points vectors have 5 k-components. The same calculation method is valid to determine the ${\psi }_p$ element error.

$${\chi }_p=\left\{\begin{array}{cc}{\psi }_p\hfill & \mathit{if}J\left({\psi }_p\right)>min\left(J\left({\chi }_p\right)\right)\hfill \\ {\chi }_p\hfill & \mathit{other}\mathit{case}\hfill \end{array}\right.$$

(46)

$$J\left({\chi }_p\right)\%=100\times \frac{1}{2}\times \left(\sqrt{\frac{1}{5}\sum _{k=1}^5{\left|\frac{{\xi }_{1,k}-{\xi }_{p,k}\left({\chi }_p\right)}{{\xi }_{1,k}}\right|}^2}+\sqrt{\frac{1}{5}\sum _{k=1}^5{\left|\frac{{\xi }_{2,k}-{\xi }_{p,k}^{\prime }\left({\chi }_p\right)}{{\xi }_{2,k}}\right|}^2}\right)$$

(47)

4.5. Adaptive Boundaries

When compared to the conventional algorithm, the proposed differential evolution algorithm considers the inclusion of two additional sections. These sections are (A) automatic initial definition of boundaries based on the information of the cardinal points and (B) adaptive optimization of the parameters based on statistical information.

In order to find the best solution, the DEA is executed $K_{max}$ times, generating a vector of dimension $K_{max}$. It is important to note that the values of the population K depend on $G_{max}$, $NP$, D, F, $CR$, ${\chi }_L$, and ${\chi }_H$. Once the population K has been obtained, a, $R_s$, and $R_{sh}$ upper and lower bounds must be rewritten according to the statistical maximum and minimum of the elements that comprise their populations $\overline{x_p}$. Equations (48) and (49) indicate the mathematical calculation of the above, where $q_{n,\overline{x_p}}$ represents the $n-$th quartile of the $\overline{x_p}$ elements. It is important to consider that the statistical maximum and minimum values cannot be above or below the smallest element in the population. To solve this problem, the results should be verified according to Equations (51) and (50).

This procedure must then be repeated B times until the solution reaches the required tolerance $J_{tol}$, set by a user-defined threshold, or if the algorithm meets the maximum number of iterations $B_{max}$. The performance of the solution is determined as indicated in Equation (47). Figure 9 indicates the proposed DEA flowchart.

$$\begin{array}{ccc}\hfill x_{d,H}& =& q_{2,\overline{x_p}}+1.5·(q_{3,\overline{x_p}}-q_{1,\overline{x_p}})\hfill \end{array}$$

(48)

$$\begin{array}{ccc}\hfill x_{d,L}& =& q_{2,\overline{x_p}}-1.5·(q_{3,\overline{x_p}}-q_{1,\overline{x_p}})\hfill \end{array}$$

(49)

$$x_{d,H}=\left\{\begin{array}{cc}max\overline{x_p}\hfill & \mathit{if}{x}_{d,H}>max\overline{x_p}\hfill \\ x_{d,H}\hfill & \mathit{other}\mathit{case}\hfill \end{array}\right.$$

(50)

$$x_{d,L}=\left\{\begin{array}{cc}min\overline{x_p}\hfill & \mathit{if}{x}_{d,L}<min\overline{x_p}\hfill \\ x_{d,L}\hfill & \mathit{other}\mathit{case}\hfill \end{array}\right.$$

(51)

5. Results and Analysis

This section is divided as follows (i) definition of the algorithm parameters, and (ii) model validation and discussion under different scenarios.

5.1. Algorithm Parameters Definition

The behavior of the proposed algorithm is simulated from two different perspectives: (i) impact on the solution when the adaptive evolution limit module is included, and (ii) examination of the algorithm performance according to different values of $\phi$. The differential algorithm parameters used are $NP=30$, $F=0.4$, and $CR=0.4$, adopted from [19]. The maximum values of the iterations B, K, and G are assumed equal to $B_{max}=15$, $K_{max}=9$, and $G_{max}=6$. $J_{tol}$ is set to zero and the KC200GT SPVM is used to execute these experiments. Note that since the optimized variables are a, $R_s$, and $R_{sh}$, D equals three.

5.1.1. Adaptive Evolution of Limits

The performance of the adaptive limit module is studied by comparing the DEA $J\left({\chi }_p\right)$ % value on three trials. For this purpose, the SDM is used to execute the optimization, assuming $\phi$ equal to 1%. The initial limits for the proposed $\phi$ value are indicated in

Table 1. Initial parameter limits for the KC200GT SPVM using $\phi =1$%.

Bounds	a V	${\mathit{R}}_{\mathit{s}}$$\mathbf{m}\mathbf{\Omega }$	${\mathit{R}}_{\mathit{sh}}$$\mathbf{\Omega }$
$X_{L,B=1}$	1.3874	0.8673	43.833
$X_{H,B=1}$	2.7748	867.20	341.42

. It can be seen from Figure 10 that the error of the classical DEA algorithm decreases stochastically as improved solutions are found within the initial boundaries specified. On the other hand, the solution found by the proposed DEA algorithm improves in each B iteration, until the fifth iteration. This is explained by the fact that in the fifth iteration, the values of the upper and lower limits are similar; therefore, the algorithm has converged to a final value.

Table 2. Results for the different tests conducted in the experiment.

Parameter	Classic	Proposed
$I_{ph}$ A	8.2235	8.2236
$I_o$ nA	1.5173	1.6605
a V	1.4693	1.4752
$R_s$$\mathrm{m}\mathsf{\Omega }$	315.58	313.11
$R_{sh}$$\mathsf{\Omega }$	192.57	188.63
$J\left({\chi }_p\right)$ %	0.3652	$0.3558$

indicates the best SDM solution achieved for both the classic and proposed DEA. It can be seen that the error of the solution found by the proposed algorithm is approximately 2.58% better than the result obtained by the classic method. However, this does not mean that the classic algorithm cannot reach the same solution, as this value is contained within the search range of the classic DEA.

The evolution of the boundaries of a, $R_s$, and $R_{sh}$ can be seen in Figure 11. This figure indicates that the final value of the optimized variables is very similar for each trial. These results suggest that the proposed DEA is capable of delivering a single high performance solution in few iterations every time it is executed. Therefore, the adaptive boundary module will be considered in every B iteration. Additionally, the algorithm is run three times selecting the best solution.

5.1.2. Selection of $\phi$

This test allows to study the performance of the proposed DEA for $\phi$ values contained in the $\left[{10}^{-3}\%,20\%\right]$ range. It is important to note that $\phi$ only affects $R_{s,inf}$ and $R_{sh,sup}$; therefore, the remaining variables involved in the SDM are not considered in this sub-section.

Figure 12 shows the evolution of the objective function, where each point in the figure represents the mean value of 100 executions of the proposed DEA. It can be seen from this figure that as $\phi$ increases from ${10}^{-3}\%$, the value of $J\left({\chi }_p\right)$ decreases until it reaches a minimum value which resides in the range of [1%, 2%]. From this minimum value, as $\phi$ continues to increase, the objective function starts to increase. This performance degradation is caused by the fact that the value of the upper and lower limits of the series and shunt resistances have similar values; therefore, the DEA can only modify the equivalent thermal voltage in order to find the optimal solution.

Figure 13a,b indicate the evolution of the optimal series and shunt resistances as well as their boundaries. In the same manner to the results presented in Figure 12, the algorithm is executed 100 times, obtaining the mean value of the DEA. It can be seen that for a given $\phi$ critical value, the optimal resistances equals its upper (or lower) boundary. This fact suggest that $\phi$ must be lower than the critical value. For the KC200GT SPVM, a critical value near $2\%$ is achieved. This work assumes an appropriate $\phi$ value of 1%, as a wide range of solutions can be obtained from this value.

5.2. Model Validation

This sub-section is divided as follows: (i) comparison of the electrical parameters obtained for the Kyocera KC200GT SPVM, and (ii) comparison of the cardinal points estimated using the database added in Appendix A. The cardinal points considered as input data correspond to (i) STC (reference condition), and (ii) NOCT (operating condition), both available in the datasheet of each module.

5.2.1. Kyocera KC200GT Comparison

The performance of the SDM is compared to the information presented in the state-of-the-art. In order to calculate the boundaries, a value of $\phi$ equal to 1% is used according to the results indicated in Section 5.1.2. In addition, the DEA is executed three times, selecting the best solution. The datasheet information is displayed in

Table 3. Error of the calculated parameters for the Kyocera KC200GT SPVM. The table displays the cardinal points at the STC and NOCT conditions.

	Datasheet	Acarino  et al. [34]	N. Eddine et al. [35]	Jadli et al. [28]	Elazab et al. [13]	Ebrahimi et al. [14]	Biswas et al. [20]	Proposed Method
Solution type	-	Explicit	Explicit	Iterative	Heuristic	Heuristic	Heuristic	Heuristic
$I_{ph,STC}$ A	-	8.2100	8.2233	8.2119	8.2800	8.2186	8.2197	8.2236
$I_{o,STC}$ nA	-	2.1546	2.1524	196.06	85.580	1.4360	68.000	1.6784
$a_{STC}$ V	-	1.4921	1.4926	1.87656	1.7897	1.4641	1.7702	1.4759
$R_{s,STC}$ m$\mathsf{\Omega }$	-	284.40	308.00	210.89	281.5	240.94	191.10	313.06
$R_{sh,STC}$$\mathsf{\Omega }$	-	157.54	193.05	895.80	424.22	130.28	161.74	189.38
$I_{sc,STC}$ A	8.21	8.1952	8.2102	8.2100	8.2745	8.2034	8.2100	8.2100
$V_{oc,STC}$ V	32.9	32.879	32.901	32.926	32.892	32.849	32.900	32.900
$I_{mpp,STC}$ A	7.61	7.5728	7.6087	7.6058	7.6436	7.5662	7.5279	7.6103
$V_{mpp,STC}$ V	26.3	26.449	26.305	26.340	25.968	26.762	26.613	26.299
$P_{mpp,STC}$ W	200.14	200.29	200.15	200.34	198.49	202.49	200.34	200.14
$I_{sc,NOCT}$ A	6.62	6.6144	6.6262	6.6242	6.6764	6.6211	6.6255	6.6261
$V_{oc,NOCT}$ V	29.9	29.436	29.457	27.970	28.276	29.516	28.366	29.557
$I_{mpp,NOCT}$ A	6.13	6.0569	6.0819	6.0169	6.0614	6.0584	5.9899	6.0857
$V_{mpp,NOCT}$ V	23.2	23.404	23.295	21.828	21.869	23.753	22.442	23.390
$P_{mpp,NOCT}$ W	142.22	141.75	141.68	131.34	132.56	143.90	134.42	142.34
$J\left({\chi }_p\right)$ %	-	0.6501	0.3815	2.6112	2.6796	1.2281	2.2014	0.3563

.

Figure 14 indicates the evolution of the differential evolution algorithms error. It can be observed that the final solutions displays an error under 0.4%. To complement this result, Table 3 indicates different SDMs that describe the behavior of the KC200GT SPVM as well as the datasheet information. These parameters are obtained from the state-of-the-art and represent several methods for determining the SDM model of the SPVM under study, i.e., E.S. explicit solution, I.S. iterative solution, and H.S. heuristic solution. Within the heuristic approaches, three methods have been presented (i) whale optimization algorithm WOA, (ii) flexible particle swarm algorithm FPSOA, and (iii) differential evolution algorithm DEA. The last row summarizes the performance values of the different groups of parameters assessed, as indicated in Equation (47). Here, the three best-performing SDMs are highlighted in bold-type.

According to Table 3, the method proposed by N. Eddine et al., has the best performance, with a $J\left({\chi }_p\right)$ value of $0.3815$%. In contrast, the worst solution is indicated by [13], with a $J\left({\chi }_p\right)$ value of $2.6796$%. The proposed method reaches a $J\left({\chi }_p\right)$ value equal to 0.3563%, representing an improved solution when compared with the results indicated in the state-of-the-art. This set of parameters allow predicting the electrical behavior of the SPVM for both STC and NOCT conditions with high precision. In addition, since the proposed algorithm defines the boundaries of the variables automatically, it can be executed by any user without prior knowledge being required. Figure 15a,b indicate IV curves and PV curves for the three best SDMs solutions highlighted in Table 3. It can be seen that there are no appreciable differences between the models presented; however, the proposed DEA has the advantage of being adaptable and is capable of redetermining other electrical models, such as the double-diode model.

5.2.2. Database Comparison

In order to study the applicability of the proposed method on different solar module technologies, a database composed of 100 SPVM datasheets is used, and included in Appendix A). The SPVMs database includes mono-Si, Si poly-Si, and thin-film module technologies. The rated power ranges from 36 to 455 W.

Figure 16 indicates the behavior of the $R_{s,inf}$, and $R_{sh,sup}$ parameters using a $\phi$ value equal to 1%. It can be seen that $R_{sh,sup}$ remains in the range of ${10}^2$ $\mathsf{\Omega }$ to ${10}^3$ $\mathsf{\Omega }$ for the monocrystalline and polycrystalline SPVMs. In a similar way, the values of $R_{s,inf}$ remain in the range of ${10}^{-2}$ $\mathsf{\Omega }$ and ${10}^{-1}$ $\mathsf{\Omega }$. However, when calculating the optimum SDM for the thin-film SPVMs, the value of the boundaries increases approximately by an order of magnitude. This is explained by the fact that the shape of the IV curve of thin-film technologies is significantly different when compared to monocrystalline and polycristalline technologies.

On the other hand, Figure 17 resumes the performance of the different collected SPVMs. The configuration of the optimization algorithm remains the same as indicated in the previous section. It can be seen that the error remains under 3% for all SPVMs tested. This result indicates that the use of a $\phi$ value equal to 1% is capable of achieving a high quality response, even for the thin-film SPVMs. The mean value of $J\left({\chi }_p\right)$ obtained from the database is 0.77%.

6. Conclusions

This work presents a methodology able to determine the optimal boundaries of the parameters, which describe the single-diode photovoltaic model SDM. The boundaries are calculated by evaluating the deviation of the power at the maximum power point for different values of series and shunt resistors, defined by the $\phi$ parameter. This study uses the differential evolution algorithm DEA with step-by-step boundary definition module, improving the performance of the classic DEA algorithm. The proposed calculation method uses cardinal points provided from two different operating conditions, e.g., standard test conditions, and normal operating cell temperature, data usually included in the technical datasheet. The methodology was tested using two separate tests: (i) by comparing the electrical parameters obtained with the model for the Kyocera KC200GT SPVM with the results of well-established algorithms, and (ii) by comparing the algorithms predictions with 100 SPVM datasheets of different technologies and rated powers.

A value of $\phi$ equal to 1% is established from the results of the KC200GT SPVM comparison. This study also demonstrates that the changes introduced in the optimization algorithm reach a superior solution compared to the values achieved by traditional algorithms presented in the state-of-the-art. The second section of the results indicates that the methodology presented is capable of determining an optimal SDM for different solar module technologies and rated powers, reaching a maximum error near 3%.

Future work will focus on the following points: (i) extension of the methodology to support complex models, such as the double-diode model, and (ii) scalability of the proposed algorithm applied to a group of PV panels. The first point considers a different number of parameters to be optimized; therefore, the dimension D must be redefined in the same manner as the initial definition of the boundaries. On the other hand, the second point considers the application of the proposed algorithm for high accuracy estimation of electrical power in PV plants.