**Estimation of Single-Diode and Two-Diode Solar Cell Parameters by Using a Chaotic Optimization Approach**

#### **Martin Calasan ´ 1,\*, Dražen Jovanovi´c 2, Vesna Rubeži´c 1, Saša Mujovi´c <sup>1</sup> and Slobodan Ðukanovi´c <sup>1</sup>**


Received: 3 October 2019; Accepted: 29 October 2019; Published: 4 November 2019

**Abstract:** Estimation of single-diode and two-diode solar cell parameters by using chaotic optimization approach (COA) is addressed. The proposed approach is based on the use of experimentally determined current-voltage (*I-V*) characteristics. It outperforms a large number of other techniques in terms of average error between the measured and the estimated *I-V* values, as well as of time complexity. Implementation of the proposed approach on the *I-V* curves measured in laboratory environment for different values of solar irradiation and temperature prove its applicability in terms of accuracy, effectiveness and the ease of implementation for a wide range of practical environment conditions. The COA-based parameter estimation is, therefore, useful for PV power converter designers who require fast and accurate model for PV cell/module.

**Keywords:** Solar cell parameters; single-diode model; two-diode model; COA

#### **1. Introduction**

The contribution of solar energy in total electric energy production is growing constantly. As the price of solar inverters and solar panels constantly decreases, most countries are basing their energy policy on higher use of solar energy. Studies on energy networks, and especially testing of the integration of solar energy sources into power networks, requires accurate calculation of the solar output power, as well as accurate modeling of solar cells. For that reason, modeling of solar cells (corresponding equivalent circuit and accurate parameters value) represents a very popular research field.

In the literature, two basic models of the equivalent circuits of solar cell can be found, namely the single-diode model (SDM) [1] and the double-diode model (DDM) [2]. DDM considers the composite effect of the neutral region of the junction, and, therefore, models the solar cells more accurately [3]. However, it is characterized by seven unknown parameters. Because of the complexity of DDM, some authors reduce the number of unknown parameters [3,4], which can greatly affect the model accuracy [5]. In this work, we focus on both SDM and DDM, without neglecting any of the model parameters of the solar cell.

For the estimation of solar cell parameters, two main sets of "input data" and corresponding estimations can be found, namely


The former uses the datasheet information (open circuit voltage, short circuit current, voltage and current value at maximum power point characteristics) provided by photovoltaic (PV) manufacturers under standard test conditions. However, recent research [7] on the usage of datasheet values for solar cell parameter estimation shows that current-voltage characteristic is not unique when designers focus on three datasheet points (open circuit, short circuit and maximum power). It is shown that by observing only three points, we can have multiple *I-V* characteristics, although in reality a solar cell has a single defined *I-V* characteristic corresponding to a specific set of cell parameters. To obtain unique and accurate *I-V* characteristic of PV cell, experimental data on more than three major points are necessary [7]. Research [7] has also implied that only approaches based on experimental data generate accurate models.

Evaluating the performance of solar cells (or PV panels) requires as accurate an estimation of the equivalent circuit solar cell parameters as possible. The approaches used for this purpose can be categorized as follows:


Analytical techniques provide mathematical expressions for solving equivalent circuit parameters based on some input data (manufacturer data or/and data obtained from measurements). A review and comparative assessment of non-iterative methods for the extraction of the single-diode model parameters of photovoltaic modules is given in [11]. In general, analytical techniques provide rapid solution. On the other hand, these techniques give erroneous results when the estimated and measure solar cell output characteristics are compared [12].

Numerical techniques are based on curve fitting, usually via iterative methods. However, the application of curve fitting to nonlinear diode equations is quite complex, making numerical determination of solar cell parameters unpopular [14].

Recently, meta-heuristic algorithms for solar cell parameter estimation have been proposed [14–58]. They impose no restrictions on the problem formulation, they are excellent in dealing with nonlinear equations, and they can be applied for different numbers of unknown parameters.

Of all the proposed techniques, none excels in terms of accuracy and efficiency with respect to others. This was our main incentive for doing research in this field. We propose both accurate and efficient parameter optimization of solar cell SDM and DDM through chaotic optimization approach (COA).

Recently, COA has been used in solving various optimization problems: parameter identification of Jiles-Atherton hysteresis model [59], single-phase transformer parameter estimation [60], design of PID parameters for automatic voltage regulation of synchronous machine [61], antenna array radiation pattern synthesis [62,63]. The main advantages of COA over other optimization techniques are easy implementation and short execution time [64]. It should be noted that different versions of chaotic algorithm have also been used in solar cell parameters estimation, namely chaotic heterogeneous comprehensive learning particle swarm optimizer variants [16], chaotic asexual reproduction optimization [33], mutative-scale parallel chaos optimization algorithm [41], chaos-embedded gravitational search algorithm [65], chaotic improved artificial bee colony algorithm [66], improved chaotic whale optimization algorithm [67], etc. Unlike methods proposed in [16,33,41,65–67], this paper will use COA based on Logistic map for solar cell parameter estimation. Logistic map has a simple form with one variable and one control parameter, and it can produce chaotic behavior similar to more complex chaotic systems. [68]. The power of this optimization approach is demonstrated in [60], where it was used for the estimations of the transformer's parameters.

The effectiveness of the proposed approach will be evaluated on different solar cells (different with respect to solar cell voltage and current level) found in the literature and the laboratory environment. Furthermore, COA-based parameter estimation will be compared with 50 various literature techniques

for SDM, and with 12 various techniques for DDM. Also, we will compare parameters obtained by using the proposed method on the measured data with analytically and numerically obtained parameters. Finally, we will apply COA on the measured *I-V* characteristics using the laboratory environment.

The paper is organized as follows. SDM and DDM for solar cells are described in Section 2. In Section 3, COA and its implementation for solar cell parameters estimation are described. The comparison of solar cell parameters estimation accuracy obtained by using the COA-based and other methods, for one solar cell and one solar module, is presented in Section 4. The experimental setup for measuring *I-V* curves is presented in Section 5, along with the COA-based parameter estimation results. The concluding remarks are given in Section 6.

#### **2. Mathematical Modeling of Single and Double Solar Cells**

SDM is commonly used model for solar cell representation [1], and the equivalent circuit is shown in Figure 1a. The *I-V* relationship for this model can be described by the following equation:

$$I = \left. I\_{\overline{P}^{\overline{v}}} - I\_o \left( e^{\frac{V + IR\_o}{n \cdot V\_{\overline{h}}}} - 1 \right) - \frac{V + IR\_s}{R\_p} \right. \tag{1}$$

where *Ipv* is the photo-generated current, *Rs* the series parasitic resistance, *Rp* the parallel parasitic resistance, *I*<sup>0</sup> the saturation current, *n* is ideality factor and *Vth* = *k*B*T*/*q* is the thermal voltage (*k*<sup>B</sup> is Boltzmann constant equal to 1.38 <sup>×</sup> 10−<sup>23</sup> J/K, *T* the temperature and *q* the electron charge equal to 1.602 <sup>×</sup> <sup>10</sup>−<sup>19</sup> C).

**Figure 1.** (**a**) SDM, and (**b**) DDM of a solar cell.

The equivalent circuit with DDM for the solar cell is shown in Figure 1b. Therefore, unlike the SDM model, the DDM model of the solar cell, in addition to the rectifying diode, includes one more diode to consider the space charge recombination current [48]. The *I-V* characteristic of DDM is given as

$$I = \ I\_{pv} - I\_{o1} \left( e^{\frac{V + IR\_{\rm{th}}}{n\_1 - V\_{\rm{th}}}} - 1 \right) - I\_{o2} \left( e^{\frac{V + IR\_{\rm{th}}}{n\_2 - V\_{\rm{th}}}} - 1 \right) - \frac{V + IR\_{\rm{th}}}{R\_{\rm{P}}} \tag{2}$$

where *Io*<sup>1</sup> and *Io*<sup>2</sup> are the diffusion and saturation currents, whereas *n*<sup>1</sup> and *n*<sup>2</sup> are the diffusion and recombination diode ideality factors [48]. The ideality factor is discussed in [69,70], whereas [71] presents a method for ideality factor calculation.

#### **3. COA and Objective Function**

COA is a very powerful optimization technique that has found numerous scientific applications [59–63]. This approach is based on the theory of chaos, which is, in a mathematical sense, described by ordinary differential equations or by an iterative map [64].

Different chaotic systems, including the logistic map, lozi map, tent map and Lorenz system, can be found in the literature. In this paper, we will base COA on the logistic map [59–64].

The task of COA is to estimate a set of unknown parameters X which minimizes the objective function (*OF*). In our case, for SDM, X = [*Rs*, *Rp*, *Ipv*, *Io*, *n*], and for DDM, X = [*Rs*, *Rp*, *Ipv*, *Io*1, *Io*2, *n*1, *n*2]. Therefore, in general, vector X = [x1, x2, ... x*n*] contains variables limited to the lower (*LV*) and upper (*UV*) permitted value, i.e, *xi* ∈ [*Li*, *Ui*]. On the other side, the *OF* for SDM is

$$OF \;= \sum\_{t=-1}^{P} \left( I\_{pv} - I\_o \left( e^{\frac{V\_t + I\_t R\_o}{n \cdot V\_{th}}} - 1 \right) - \frac{V\_t + I\_t R\_s}{R\_P} - I\_t \right) \tag{3}$$

whereas for DDM it reads

$$OF = \sum\_{t=-1}^{P} \left( I\_{\mathbb{P}^{\mathcal{U}}} - I\_{\mathbb{G}1} \left( e^{\frac{V\_{t} + l\_{\mathbb{P}^{R\_{\mathcal{S}}}}{n\_{\mathcal{I}} + V\_{\mathbb{R}^{\mathcal{S}}}}}{n\_{\mathcal{I}} + V\_{\mathbb{G}1}} - 1 \right) - I\_{\mathbb{G}2} \left( e^{\frac{V\_{t} + l\_{\mathbb{P}^{R\_{\mathcal{S}}}}{n\_{\mathcal{I}} + V\_{\mathbb{R}^{\mathcal{S}}}}}} - 1 \right) - \frac{V\_{\mathbb{I}} + I\_{\mathbb{I}}R\_{\mathbb{S}}}{R\_{\mathcal{P}}} - I\_{\mathbb{I}} \right) \tag{4}$$

where *P* is the number of measured *I-V* pairs from the *I-V* characteristics, and *Vt* and *It* represent the voltage and current value of pair *t*.

Figure 2 presents the search procedure, i.e., the COA flowchart. The detailed description of COA flowchart can be found in [59].

**Figure 2.** COA flowchart.

In this paper, the following COA parameters were used: *M* = 1000, *N* = 50,000. The COA-based estimation is compared with other approaches through the root mean square error (*RMSE*), defined as follows:

$$RMSE = \sqrt{\frac{\sum\_{k=1}^{P} \left(I\_{\text{cst},k} - I\_{\text{meas},k}\right)}{P}} \tag{5}$$

where *Iest,k* and *Imeas,k* represent the estimated and the measured values of solar output current in point *k*, respectively.

#### **4. Simulation Results**

To evaluate COA for solar cell parameters estimation, we first applied the proposed method to an experimental current-voltage characteristic extracted from the manufacturer's datasheets of a well-known R.T.C. France solar cell operating under standard test conditions.

The values of parameters obtained by using COA for the R.T.C. France solar cell are summarized, by year of publication, in Table 1, for SDM and Table 2 for DDM. These values are compared with the values of parameters published in recent papers (column Reference) for the same experimental data. During the estimation process, the parameter ranges for SDM estimation were *Rs*(Ω) ∈ [0.02, 0.05], *Ipv*(*A*) ∈ [0.74, 0.78], *Io*(μ*A*) ∈ [0.2, 0.4], *Rp*(Ω) ∈ [50, 55] and *n* ∈ [1.35, 1.6], whereas for DDM, they were *Rs*(Ω) ∈ [0.02, 0.04], *Rp*(Ω) ∈ [54, 58], *n*<sup>1</sup> ∈ [1.4, 1.5], *n*<sup>2</sup> ∈ [1.95, 2], *Ipv*(*A*) ∈ [0.75, 0.77] *Io*1(μ*A*) ∈ [0.2, 0.25] and *Io*2(μ*A*) ∈ [0.7, 0.8].

**Table 1.** Calculated SDM parameters for the R.T.C France solar cell.


\* for this method, a real RMSE are given [56].


**Table 2.** Calculated DDM parameters for the R.T.C France solar cell.

Tables 1 and 2 report parameters as they appear in the cited papers with no modification. However, in some papers in the Energy Conversion and Management journal (in Table 1 marked by \*), inaccuracies occurred in parameter estimation of the PV cell using metaheuristic techniques. Namely, the results proposed in [15–18,21] do not correspond to the objective function [56].

The presented results, especially the value of RMSE, show that COA offers solar characteristics closer to the measured characteristics than the other existing methods, i.e., it outperforms other methods in terms of accuracy. In addition, by observing Tables 1 and 2, it is also evident that DDM characterizes solar cells more accurately than SDM, which supports the conclusion regarding DDM accuracy noted in [3].

It can be seen that COA outperforms several other techniques, such as evaporation rate-based water cycle algorithm (ER-WPA) [19] and cat swarm optimization (CSO) [24] for SDM, and with the generalized opposition-flower pollination algorithm-nelder-mead simplex method (GOFPANM) [53], by a small margin. However, the implementation of COA is simpler than implementation of ER-WPA, CSO and GOFPANM. Furthermore, COA is computationally less demanding than CSO since CSO requires changing the operation mode during the estimation process [24]. On the other side, GOFPANM is a hybrid algorithm which combines local and global search as well as different algorithms during estimation [53]. In general, most evolutionary algorithms have the complexity of O((*np* + *C*of *p*)*N*i), where O is the big O notation, *n* is the dimension of the parameter space, *p* is the population size, *N*<sup>i</sup> is the number of iterations and *C*of is the complexity of the *OF*. The complexity of COA is O(*QC*of), where *Q* is the number of points in the parameter space in which the OF is calculated. Therefore, the proposed COA-based estimation has significantly lower computational complexity than evolutionary algorithms.

To show the additional advantage of COA over other techniques, we conducted a comparison in terms of required time for one iteration. In that sense, in MATLAB 2015 (MathWorks, Natick, MA, USA) we have implemented the following algorithms for solar cell parameter estimation: evaporation rate-based water cycle algorithm (ER-WCA) [19], cuckoo search (CS) [46] and harmony search (HS) [48]. ER-WCA algorithm has a very good accuracy, very close to that obtained by the proposed method (see Table 1). On the other hand, HS and CS also have a good accuracy (~10−4). All computer simulations were carried out on a PC with Intel(R) Core (TM) i3-7020U CPU @ 2.30 GHz and 4 GB RAM. The obtained results, i.e., the mean, maximal and minimal required time per one iteration, obtained over 20 runs, are presented in Table 3. Clearly, the COA-based algorithm is the most efficient method, as it is characterized by the lowest value of required time per iteration. Note, in order to draw a fair comparison between the considered algorithms, MATLAB implementation follows the same rules for each algorithm (e.g., avoiding loops and using array operations such as dot product and matrix product whenever possible).


**Table 3.** Time per iteration comparison.

The measured *I-V* and *P-V* characteristics and the corresponding simulated characteristics, for parameters obtained by using COA, are shown in Figure 3. Very good agreement can be seen between the measured and estimated curves. Also, the difference between the DDM and SDM simulated curves is small but consistent and always in favor of DDM. In addition, in Table 4, we presented the estimated value of the unknown DDM parameters of the BPSolar MSX-60 module. These parameters are obtained by using COA as well as by using analytical, numerical, iteration and Newton methods presented in the literature. In the COA-based estimation, the ranges of parameters were *Rs*(Ω) <sup>∈</sup> [0.2 0.4], *Rp*(Ω) <sup>∈</sup> [150, 300], *<sup>n</sup>*<sup>1</sup> <sup>∈</sup> [0.5, 1.5], *Ipv*(*A*) <sup>∈</sup> [3.5, 4], *Io*1(*A*) <sup>∈</sup> - 10<sup>−</sup>10, 10−<sup>6</sup> . , *Io*2(*A*) ∈ - 10<sup>−</sup>10, 10−<sup>6</sup> . , and *n*<sup>2</sup> ∈ [1.5, 2]. From the presented results, it is clear that COA outperforms the considered non-metaheuristic methods for solar cell parameters determination in terms of accuracy.

**Figure 3.** *I-V* and *P-V* characteristics of R.T.C. France solar cell.

**Table 4.** Calculated DDM parameters for the BPSolar MSX-60 module.


The measured and corresponding simulated *I-V* and *P-V* characteristics, for parameters obtained by using COA and other methods, are shown in Figures 4 and 5, respectively. It is evident that COA outperforms the other methods in terms of approaching the measured characteristics.

**Figure 4.** *I-V* characteristics of BPSolar MSX-60 module.

**Figure 5.** *P-V* characteristics of BPSolar MSX-60 module.

Based on all of the presented results, it can be concluded that COA can precisely estimate the solar cell/module circuit parameters, outperforming the other metaheuristics as well as analytical or numerical methods in terms of estimation accuracy.

#### **5. Experimental Results and Analysis**

To check the applicability and efficiency of COA for solar cell parameter estimation, we also observed solar cells from the Clean Energy Trainer setup. The main motivation to use these solar modules is that this setup enables adjustable solar insolation, USB data monitoring for PC-supported data acquisition and analysis, as well as highly advanced didactic software for system control and real-time data plotting.

The observed system contains of:


The experimental setup, installed in Laboratory for Automatics, at the Faculty of Electrical engineering, University of Montenegro, is presented in Figure 6.

**Figure 6.** Experimental setup.

Firstly, we measured the *I-V* characteristics for insolation of 1285 W/m2 and temperature of 42 ◦C. For the measured *I-V* pairs, we determined single and double diode solar cell parameters (see Table 5). The parameter ranges for SDM estimation were *Rs*(Ω) <sup>∈</sup> [0.1, 0.4], *Ipv*(*A*) <sup>∈</sup> [0.2, 0.4], *Io*(*A*) <sup>∈</sup> - <sup>5</sup> <sup>×</sup> <sup>10</sup><sup>−</sup>8, 15×10−<sup>8</sup> . , *Rp*(Ω) ∈ [200, 600] and *n* ∈ [0.2, 1], whereas for DDM were *Rs*(Ω) ∈ [0.1, 0.4], *Rp*(Ω) <sup>∈</sup> [600, 900], *<sup>n</sup>*<sup>1</sup> <sup>∈</sup> [0.2, 1], *<sup>n</sup>*<sup>2</sup> <sup>∈</sup> [1.95, 2], *Ipv*(*A*) <sup>∈</sup> [0.2, 0.4], *Io*1(*A*) <sup>∈</sup> - <sup>5</sup> <sup>×</sup> <sup>10</sup><sup>−</sup>8, 15×10−<sup>8</sup> . , and *Io*2(*A*) ∈ - <sup>5</sup> <sup>×</sup> <sup>10</sup><sup>−</sup>8, 15×10−<sup>8</sup> . . Then we measured the *I-V* and *P-V* characteristics for different values of insolation and temperature. The corresponding simulated characteristics were determined by taking into account the change of parameters with insolation and temperature (see [13]). The measured and estimated *I-V* and *P-V* characteristics for different values of insolation and temperature are presented in Figures 7–10. The agreement between the measured and estimated characteristics is evident (see zoomed parts in these figures). Finally, we repeated the estimation procedure on all measured *I-V* characteristics. The estimated values of parameters were in range of ±4% of the initially observed, which confirms that we can use any of the measured characteristics for parameter estimation. On the other hand, by observing the data provided in Table 5, even for this module, it is evident that DDM is more accurate than SDM.

**Table 5.** Estimated value of experimentally tested solar module parameters.


**Figure 7.** Current-voltage characteristics for two different insolation values and for temperature T = 42 ◦C.

**Figure 8.** Power-voltage characteristics for two different insolation values and for temperature T = 42 ◦C.

**Figure 9.** Current-voltage characteristics for two different temperatures and insolation values.

**Figure 10.** Power-voltage characteristics for two different temperatures and insolation values.

#### **6. Conclusions**

Modeling of solar cells is a very popular research direction, which is supported by numerous recent contributions in the literature. This paper proposes COA as a very successful approach for this purpose.

The proposed method is verified using practical data from various manufacturers. Its accuracy is confirmed by comparing its RMSE with numerous metaheuristics and non-metaheuristics methods for different solar cells. Experimental testing of COA applicability for parameter estimation is also implemented in laboratory environment. In all considered scenarios, a high level of accuracy is demonstrated. Apart from this, excellent matching of the simulated *I-V* and *P-V* curves with the measured characteristics additionally confirms the COA accuracy and its applicability for parameter estimation.

In future work, our attention will be focused on the usage of COA for estimation of solar cell parameters when solar cell output current is represented through the Lambert W function.

**Author Contributions:** Conceptualization, M.C. and S. ´ Đ.; Methodology, M.C.; software, D.J. and V.R.; validation, ´ V.R and M.C., formal analysis, S.M.; investigation, M. ´ C., V.R and S. ´ Đ.; resources, S.M.; writing—original draft preparation, M.C. and S. ´ Đ.; writing—review and editing, S.Đ.; visualization, D.J.; supervision, M.C. ´

**Funding:** This work has been supported through European Union's Horizon 2020 research and innovation program under project CROSSBOW-CROSS BOrder management of variable renewable energies and storage units enabling a transnational Wholesale market (Grant No. 773430).

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

### *Article* **A Solution of Implicit Model of Series-Parallel Photovoltaic Arrays by Using Deterministic and Metaheuristic Global Optimization Algorithms**

#### **Luis Miguel Pérez Archila 1, Juan David Bastidas-Rodríguez 2, Rodrigo Correa 1, Luz Adriana Trejos Grisales 3,\* and Daniel Gonzalez-Montoya <sup>4</sup>**


Received: 19 December 2019; Accepted: 30 January 2020; Published: 12 February 2020

**Abstract:** The implicit model of photovoltaic (PV) arrays in series-parallel (SP) configuration does not require the LambertW function, since it uses the single-diode model, to represent each submodule, and the implicit current-voltage relationship to construct systems of nonlinear equations that describe the electrical behavior of a PV generator. However, the implicit model does not analyze different solution methods to reduce computation time. This paper formulates the solution of the implicit model of SP arrays as an optimization problem with restrictions for all the variables, i.e., submodules voltages, blocking diode voltage, and strings currents. Such an optimization problem is solved by using two deterministic (Trust-Region Dogleg and Levenberg Marquard) and two metaheuristics (Weighted Differential Evolution and Symbiotic Organism Search) optimization algorithms to reproduce the current–voltage (I–V) curves of small, medium, and large generators operating under homogeneous and non-homogeneous conditions. The performance of all optimization algorithms is evaluated with simulations and experiments. Simulation results indicate that both deterministic optimization algorithms correctly reproduce I–V curves in all the cases; nevertheless, the two metaheuristic optimization methods only reproduce the I–V curves for small generators, but not for medium and large generators. Finally, experimental results confirm the simulation results for small arrays and validate the reference model used in the simulations.

**Keywords:** implicit model solution; photovoltaic array; series–parallel; global optimization; partial shading; deterministic optimization algorithm; metaheuristic optimization algorithm

#### **1. Introduction**

The continuous effects of climate change and environmental contamination have motivated a growing interest in the use of renewable energies to replace conventional energy sources based on fossil fuels. Within renewable energies, photovoltaic (PV) energy offers interesting advantages: long life cycle, lack of mobile parts, low maintenance costs, modularity, among others [1]; that is why the installed PV capacity has been growing in the last years [2]. In a PV system, one of the most important elements is the generator; therefore, it is important to have models of this element to reproduce its electrical behavior when they operate under homogeneous and non-homogeneous conditions [3]. These models are used in different applications like the sizing of PV generators, the power, and energy production estimation, the analysis, and evaluation of maximum power point tracking techniques, model-based diagnosis, and other applications [4]. A PV generator can be connected in different configurations, nonetheless, series-parallel (SP) is one of the most used configurations. In an SP generator, two or more modules are connected in series to form a string and two or more strings are connected in parallel to form the generator. As a consequence, all the strings have the same voltage and can be analyzed independently [5]. It is worth noting that each module is formed by one or more submodules connected in series and each submodule has a bypass diode connected in antiparallel. When the irradiance and temperature of all the submodules in the generator are the same, it can be said that the generator operates under homogeneous conditions [6]. In these conditions, the whole generator can be represented by the single-diode model (SDM) scaled in voltage, according to the number of modules in the string, and scaled in current, according to the number of strings in parallel [7]. In this case, the current vs. voltage (I–V) curve has a single knee and the power vs. voltage (P-V) curve has a single maximum power point (MPP), which can be easily tracked. In real applications, the PV generator is commonly shaded, or partially shaded, by surrounding objects, passing clouds, or other objects. Hence, the operating conditions (i.e., irradiance and temperature) of the shaded submodules are different than the operating conditions of the rest of the submodules [8]. In this case, it is said that the generator operates in nonhomogeneous conditions and the generator's I–V curve may present multiple knees, which is translated into multiple MPPs in the P-V curve [9]. To model a PV generator operating in nonhomogeneous conditions, it is necessary to consider that each submodule in a module is represented by the SDM. In this way, it is possible to include the effects of the nonhomogeneous conditions in the electrical behavior of the generator [10]. Therefore, the generator can be analyzed as an equivalent circuit that is obtained by connecting strings in parallel, where each string is formed by a set of submodules connected in series [11,12].

In the literature, most of the models use the SDM to represent each submodule and they use the LambertW function [13] to obtain an explicit equation of the submodule's current as a function of its voltage [7,14–16]. Then, each string formed by N submodules and a blocking diode is modeled independently by a system of N + 1 nonlinear equations, where the unknowns are the voltages of the blocking diode and the N submodules [7,14–16]. Solving the system of nonlinear equations of each string, it is simple to calculate the string current and finally the generator's current. However, the evaluation of the LambertW function implies a high computational cost for the solution to the system of nonlinear equations associated with each string. For example, in [7] the modeling of a PV array based on the use of the Lambert W-function to obtain an explicit relationship between the voltage and current of any PV module is presented. In addition, the non-linear equations that describe the PV array are solved by means of explicit symbolic calculation of the inverse of the Jacobian matrix. A similar solution was published in [14], where a new formulation of the one-diode model of the PV module is solved using the LambertW function with the terminal voltage expressed as an explicit function of the current, this approximation resulting in a reduction of simulation calculation time. Another solution was presented in [16], where a model of the PV field is presented by means of a set non-linear equations characterized by a sparse Jacobian matrix, which requires a moderate computational burden, both in terms of memory use and processor speed. Optimization methods for estimating the PV module parameters are presented as a suitable option to overcome the drawbacks of deterministic and iterative methods. In addition, optimization methods are not only used in PV area for estimating parameters, since the reconfiguration of PV modules to mitigate the effect of partial shading conditions in PV arrays is currently one of the most studied topics [17]. The model proposed in [18] uses the implicit equation that describes the current-voltage relationship in the SDM. From such an equation, is it possible to construct a system of N + 2 implicit equations to model each string, where the unknowns are the string current and theN+1 voltages of the N submodules and the blocking diode. Then, solving each string it is possible to calculate the generator's current as in the other models. In [18] the system of N + 2 implicit equations is solved by using the Trust-Region Dogleg (TRD) algorithm, which is a deterministic optimization method. However, the authors in [18] do not

formulate the solution of the system of N + 2 implicit equations as an optimization problem with restrictions; moreover, the authors do not provide any justification for the selection of TRD for solving the optimization problem nor evaluate other deterministic and metaheuristic optimization algorithms. Finally, those works [7,14–16,18] are focused on the model, i.e., the system of nonlinear equations, and not in the optimization algorithms used to solve such models. Therefore, from those papers it is not clear which type of optimization methods should be used to reduce the calculation time.

This paper formulates the solution of the system of N + 2 implicit equations, associated with each string, as an optimization problem with restriction for the N + 2 unknowns, i.e., string current and the voltages of the N modules and the blocking diode. Moreover, the paper evaluates four optimization algorithms, two deterministic and two metaheuristics, to solve the problem and generate the I–V and P-V curves of generators with different sizes. The two deterministic algorithms are TRD and Levenberg-Marquardt (LMA), which are widely used to solve optimization problems in different areas; while the two metaheuristic optimization algorithms are Weighted Differential Evolution (WDE) and Symbiotic Organism Search (SOS), which have been recently used for PV applications. The evaluation of the optimization algorithms is performed with simulations of small, medium, and large generators in homogeneous and non-homogeneous conditions. Finally, experimental results are used to evaluate the different optimization algorithms for a small PV generator. In light of the previous analysis, the main advantages of the proposed procedure are listed as follows:


The rest of the paper is organized as follows: Section 2 describes the implicit model used in this paper, Section 3 introduces the optimization problem and the operating principle of the optimization algorithms, Sections 4 and 5 show the simulation and experimental results, respectively, and Section 6 closes the paper with the conclusions.

#### **2. Implicit Model of a Series-Parallel Generator**

This section shows a brief description of the implicit model of SP generators proposed in [18], which includes the SDM used to represent each submodule, the system of N + 2 implicit equation of each string, and the calculation of the array current.

#### *2.1. Submodule Model*

The SDM (see Figure 1) has been widely used in the literature to represent the electrical behavior of PV submodules [10,19,20], which can be defined as a set of *Ns* cells connected in series. Such a model is able to represent mono-crystalline and poly-crystalline PV cells technologies [20]; in this paper PV modules formed by poly-crystalline cells have been used in simulation and experimental validations. In the SDM, *Iph* represents the current generated by the photovoltaic effect, the diode D models the nonlinearities of the PN junction, and the resistors *Rp* and *Rs* represent the leakage currents and ohmic losses, respectively. The bypass diode (BD) is an additional connected in antiparallel to the submodule (see Figure 1) to limit the power dissipated in the cells in mismatching conditions. The calculation of the PV module parameters has been widely reported in literature, not only for the SDM model [21] but also for the two-diode [22] and three-diode model [23]. Nevertheless, the estimating parameter techniques have been more widely studied for the SDM. Such an issue has been addressed by means of different approaches which can be grouped into: analytical, deterministic and metaheuristic [24,25]. The first two groups present drawbacks concerning accuracy

and convergence. Therefore, metaheuristic algorithms have been introduced as an attractive alternative for estimating the PV cells and modules parameters [25]. Such approaches may include particle swarm optimization [26], genetic algorithms [27], differential evolution [28], cuckoo search [29], among others.

**Figure 1.** Submodule model including the bypass diode.

From the circuit shown in Figure 1, it is possible to write an implicit equation that describes the relationship between the submodule's output voltage (V) and current (I), as shown in (1), where *Isat* and *η* are the inverse saturation current, and ideality factor of D, respectively, *Vt*,*<sup>d</sup>* is the diode thermal voltage and *Ibd* is the current through BD. The thermal voltage is defined as *Vt*,*<sup>d</sup>* = *k* · *T*/*q*, where *k* is the Boltzmann constant, *q* is the electron charge and *T* is the cell's temperature in K. Moreover, *Idb* is defined in (2), where *Isat*,*bd* and *ηbd* are the inverse saturation current and ideality factor of BD. It is worth noting that *Vt*,*<sup>d</sup>* is the same for diodes D and BD since it can be assumed that the cell's temperature and the bypass diode temperature are the same if there isn't a fast change in the irradiance [30].

$$f\left(V, I\right) = -I + I\_{pl} - I\_{sat}\left[\exp\left(\frac{V + (I - I\_{bd}) \cdot R\_s}{N\_s \cdot \eta \cdot V\_{t,d}}\right) - 1\right] - \frac{V + (I - I\_{bd}) \cdot R\_s}{R\_p} + I\_{bd} = 0. \tag{1}$$

$$I\_{bd} = I\_{sat,hd} \cdot \left[ \exp\left(\frac{-V}{\eta\_{bd} \cdot V\_{t,d}}\right) - 1 \right]. \tag{2}$$

Although the SDM parameters (i.e., *Iph*, *Isat*, *Vt*,*d*, *Rs* and *Rh*) may vary with the irradiance and temperature [31–33], *Iph* is directly proportional to the irradiance; therefore, it can be used as an indication of the incident irradiance on the submodule.

#### *2.2. Generator's Model*

In general, an SP generator is formed by M parallel-connected strings (see Figure 2), where each string is formed by a set of PV modules (*Nm*) and a blocking diode connected in series. In turn, each module is formed by one or more submodules connected in series (*Nsm*); then, each string is formed by N submodules, where *N* = *Nm* · *Nsm*. The generator is connected to a power converter, which sets the generator voltage (*Vgen*) to perform the MPP tracking; therefore, *Vgen* is known for the model and the objective is to find: the generator current (*Igen*), the current of each string, the voltage of each submodule, and the voltage of each blocking diode.

Considering the *Vgen* is common for all the strings, each string can be analyzed and solved independently. Then, a string formed by N submodules and a blocking diode connected in series can be represented by a system of N + 2 nonlinear and implicit equations, as shown in (3) [18]; where *V* is a vector formed by the voltages of the N submodules and the blocking diode (*V*- = [*V*1, *V*2, ··· , *VN*, *Vblk*]) and *Istr* is the string current.

**Figure 2.** Series–parallel photovoltaic generator with M strings, of N submodules each, connected in parallel, where each string includes a blocking diode and each module has two submodules.

$$F\left(\vec{\mathcal{V}}, I\_{\rm str}\right) = \begin{bmatrix} f\_1\left(V\_{1r}, I\_{\rm str}\right) \\ f\_2\left(V\_{2r}, I\_{\rm str}\right) \\ \vdots \\ f\_N\left(V\_{Nr}, I\_{\rm str}\right) \\ f\_{N+1}\left(V\_{\rm blk}, I\_{\rm str}\right) \\ f\_{N+2}\left(\vec{\mathcal{V}}\right) = \sum\_{i=1}^{N+1} V\_i - V\_{\rm gen} \end{bmatrix} = \vec{0}. \tag{3}$$

The first N equations of (3), correspond the voltage-current relationship in each submodule given in (1), while equation N + 1 describes the model of the blocking diode and it is described in (4), where *Isat*,*blk*, and *ηblk* are the inverse saturation current and the ideality factor of the blocking diode, respectively, and *Vt*,*<sup>d</sup>* is the diode thermal voltage defined in Section 2.1, since the blocking diode temperature is assumed the same of the submodules. Finally, equation N + 2 in (3) is obtained from the Kirchhoff voltage law in the loop formed with the string and *Vgen*.

$$f\_{N+1}\left(V\_{blk}, I\_{str}\right) = -I\_{str} + I\_{sat, blk} \cdot \left(\exp\left(\frac{-V\_{blk}}{\eta\_{blk} \cdot V\_{t,d}}\right) - 1\right) = 0. \tag{4}$$

All the electrical variables of the string are obtained by solving (3), i.e., *Istr* and *V*- ; hence, after solving all the strings, Igen can be calculated using (5).

$$I\_{\rm gcn} = \sum\_{j=1}^{M} I\_{\rm str,j}.\tag{5}$$

#### **3. Optimization Problem Formulation and Solution Methods**

This Section defines the solution of (3) as an optimization problem with its restrictions and presents a brief explanation of the optimization algorithms used in the paper to solve the optimization problem.

#### *3.1. Optimization Problem and Restrictions*

Let <sup>R</sup> be the set of real numbers and **<sup>L</sup>** a non-empty subset of <sup>R</sup>*n*, i.e., **<sup>L</sup>** <sup>⊂</sup> <sup>R</sup>*n*; then, the system of nonlinear equations given in (3) can be expressed as shown in (6), where *<sup>x</sup>* <sup>∈</sup> <sup>R</sup>*<sup>n</sup>*

and *fk* : **<sup>L</sup>** <sup>→</sup> <sup>R</sup>, <sup>∀</sup> *<sup>k</sup>* <sup>=</sup> 1, 2, ... , *<sup>N</sup>* <sup>+</sup> 2. The solution of (6) is given by all vectors *<sup>a</sup>* <sup>∈</sup> <sup>R</sup>*n*, such that *fk* (*a*) = 0, ∀ *k* = 1, 2, . . . , *N* + 2.

$$\begin{cases} \quad f\_1\left(\vec{x}\right) = 0\\ \quad f\_2\left(\vec{x}\right) = 0\\ \quad \vdots\\ \quad f\_{N+2}\left(\vec{x}\right) = 0 \end{cases} \tag{6}$$

Now, let's consider the function *Fobj* : **<sup>L</sup>** <sup>→</sup> <sup>R</sup> a function defined as shown in (7).

$$F\_{obj} = \sum\_{i=1}^{N+2} f\_i \left(\vec{\mathbf{x}}\right)^2 \; \; with \; \vec{\mathbf{x}} \in \mathbf{L}. \tag{7}$$

Note that *Fobj* is a well-defined function that guarantees the existence of a minimum, since *Fobj*(*x*) ≥ 0 for all *<sup>x</sup>* <sup>∈</sup> **<sup>L</sup>** is a totally ordered set where exists the infimum of *Fobj*(**L**) <sup>∈</sup> <sup>R</sup> and this infimum is non-negative, i.e., *inf Fobj*(**L**) ≥ 0. Therefore, the minimum of *Fobj* over **L** exists and it is a non-negative minimum.

**Theorem 1.** *Suppose that the system of equations given in* (3) *has a solution* **L** *and* (*a* ∈ **L**)*; hence, is a solution of* (3) *if and only if,a minimizes Fobj given in* (7)*.*

**Proof.** If*a* is a solution of the system of equations given in (3), then *fk* (*a*) = 0, ∀ *k* = 1, 2, ... , *N* + 2, therefore, *Fobj*(*a*) = 0. Moreover, considering that *Fobj*(*x*) ≥ 0 for all *x* ∈ **L**, then*a* is the minimum point of *Fobj*. Now, if*a* minimizes *Fobj* but it is not a solution of (3), then *Fobj*(*a*) must be positive, since *Fobj*(*x*) ≥ 0 for all *x* ∈ **L**. Taking into account that (3) has a solution in **L**, then there is a solution *x*<sup>∗</sup> ∈ **L** such that *Fobj*(*x*∗) = 0 and *x*<sup>∗</sup> = *a*. Hence, *Fobj*(*x*∗) < *Fobj*(*a*), which contradicts that*a* is a minimum point for *Fobj*.

The last paragraph proves that the minimum point of the objective function given in (7) is the set of values that solve the system of equations shown in (3). Such a solution contains the unknown voltages in the string (*V*- ) and the string current (*Istr*). However, to minimize (7) using numerical methods, it is necessary to establish the limits of *Istr* and each unknown in *V*as shown in Table 1.


**Table 1.** Upper and lower limits of the electrical variables in a string.

The lower limit of the string current is 0 A; while the upper limit is *ISC*,*STC*, which corresponds to the submodule short-circuit current in Standard Test Conditions (STC), i.e., irradiance of 1 kW/m2 and the submodule temperature is 25 ◦C (i.e., *T* = 298.15 K). Moreover, the upper limit of the submodules voltages (*V*<sup>1</sup> ... *VN*) is open-circuit voltage in STC (*VOC*,*STC*), which corresponds to the submodule voltage when *I* = 0 A. The lower limit of *V*<sup>1</sup> ... *VN* is the negative value of the maximum bypass diode voltage (*Vbd*,*max*), defined in (8), which corresponds to the bypass diode voltage when *Ibd* = *ISC*,*STC*. This current corresponds to the worst-case scenario where the submodule is completely shaded, and the string current is maximum. Finally, the upper limit of the blocking diode voltage is 0 V, which corresponds to *Istr* = 0 A; while the lower limit corresponds to the negative value of the maximum blocking diode voltage (*Vblk*,*max*) defined in (9).

$$\left| V\_{\rm bd,max} = \eta\_{\rm bd} \cdot V\_{\rm t,d} \cdot \ln \left( \frac{I\_{\rm SC,STC}}{I\_{\rm sat,hd}} + 1 \right) \right. \tag{8}$$

$$V\_{blk,max} = \eta\_{blk} \cdot V\_{t,d} \cdot \ln\left(\frac{I\_{SC,STC}}{I\_{sat,blk}} + 1\right). \tag{9}$$

#### *3.2. Optimization Methods*

This section presents a brief description of the optimization algorithms used to minimize the objective function introduced in (7) and to find the solution to the system of nonlinear equations associated with the string (3).

#### 3.2.1. Trust-Region Dogleg (TRD)

This is a deterministic algorithm with a single search agent. It starts from an initial point *x*0, that belongs to the domain of *f*(*x*), and searches the point *x*∗ that minimizes the function *f* through the iterative process indicated in (10), where *xk* is the actual position of the search agent and *pk* is the step of the *k*-th iteration.

$$
\mathbf{x}\_{k+1} = \mathbf{x}\_k + p\_k.\tag{10}
$$

To determine the magnitude and direction of the vector *pk* that approaches *xk*<sup>+</sup><sup>1</sup> to *x*<sup>∗</sup> is posed as an optimization problem, as shown in (11), where *mk*(*p*) is an approximate function of *f*(*x*) formed by the two first terms of the Taylor series in a region *Qk* that contains *xk*.

$$\min\_{p} \left\{ m\_k \left( p \right), \, p \in Q\_k \right\}. \tag{11}$$

If *mk*(*p*) is a good approximation of the objective function in a region *Qk*, then this region is denominated thrust-interval and it is defined as *Qk* = {**p** | |**p** | < Δ*k*}, where Δ*<sup>k</sup>* is the radius of the thrust-region in the *k*-th iteration. Once the vector *p*∗ *<sup>k</sup>* is found, i.e., the solution of (11), it is necessary to verify the fulfillment of the condition expressed in (12).

$$f\left(\mathbf{x}\_{k} + p\_{k}^{\*}\right) < f\left(\mathbf{x}\_{k}\right). \tag{12}$$

If condition (12) is fulfilled, the actual position is updated as *xk*<sup>+</sup><sup>1</sup> = *xk* + *p*<sup>∗</sup> *<sup>k</sup>* ; then, the next step is to calculate a new model *mk*<sup>+</sup><sup>1</sup> and a new trust-region *Qk*<sup>+</sup><sup>1</sup> from the updated point *xk*+1. Otherwise, if condition (12) is not fulfilled, the algorithm keeps the actual position, i.e., *xk*<sup>+</sup><sup>1</sup> = *xk*, and the trust-region is reduced by reducing its radius (Δ*k*+<sup>1</sup> < Δ*k*). The details of this algorithm can be found in [34] along with different variants of this method.

#### 3.2.2. Levenberg-Marquardt (LMA)

This is a deterministic optimization algorithm that requires the calculation of the jacobian matrix (**J**) and it is especially effective for solving nonlinear problems, where methods that used linear models, like Gauss–Newton, are not effective for the entire solution region [35]. The essence of this method is to select the search direction of the next step between the one given by the Gauss–Newton method and a direction close to the one provided by the gradient descent method. That direction search is selected by a criterion denominated *μk*.

Let's consider an optimization problem that has the structure of a trust-interval, as shown in (13) and (14) [35].

$$\min \ \frac{1}{2} \left\| f(\mathbf{x}\_k)(\mathbf{x} - \mathbf{x}\_k) + r\mathbf{x}\_k \right\|\_2^2. \tag{13}$$

$$\text{s.t. } \|\mathbf{x} - \mathbf{x}\_k\| \le \Delta\_k. \tag{14}$$

This is a least-squares problem, which is linear with restrictions because *r*(*xk*) is a linearized model valid in the trust-region defined by Δ*k*. The solution to the problem described in (13) and (14) is obtained by solving (15), where *s* = *x* − *xk*, and the position of the next iteration is given by (16) [35].

$$(f(\mathbf{x}\_k)^T f(\mathbf{x}\_k) + \mu\_k I)\mathbf{s} = -f(\mathbf{x}\_k)^T r(\mathbf{x}\_k). \tag{15}$$

$$\mathbf{x}\_{k+1} = \mathbf{x}\_k - \left(f(\mathbf{x}\_k)^T f(\mathbf{x}\_k) + \mu I\right)^{-1} \cdot f(\mathbf{x}\_k)^T r(\mathbf{x}\_k). \tag{16}$$

When (*J*(*xk*)*<sup>T</sup> <sup>J</sup>*(*xk*))−<sup>1</sup> *<sup>J</sup>*(*xk*)*Tr*(*xk*) <sup>≤</sup> <sup>Δ</sup>*k*, then *<sup>μ</sup><sup>k</sup>* <sup>=</sup> 0 and (16) is solved by *sk*; in this case, the best search direction is the same as the Gauss–Newton method. Otherwise, *μ<sup>k</sup>* > 0 and the solution of (16) fulfills *sk* = Δ*k*; with this condition, it is possible to find the solution of (16) and the value of *μk*. When *μ<sup>k</sup>* 0, the search direction is close to the one of the gradient descent method and (15) approaches to (17). Further details of the Levenberg–Marquardt algorithm are provided in [35].

$$
\mu\_k I \mathbf{s} = -I(\mathbf{x}\_k)^T r(\mathbf{x}\_k). \tag{17}
$$

#### 3.2.3. Weighted Differential Evolution (WDE)

This metaheuristic optimization method uses the evolution concept to solve numerical optimization problems and it is characterized by calculating the search direction and the search steps in an evolutionary way [36]. The general process of the WDE algorithm is introduced in Figure 3, where the first step is to generate two populations denominated initial population (P) and subpopulation (SubP), both randomly created with a normal distribution. Then, SubP is mutated to create TempP and vector F and binary map M are generated, by using random variables. Population SubP, F and M are used in a mutation and crossover proves to generate the new population (T), which requires verification that all the individuals fulfill the problem restrictions (boundaries control). The individuals that do not fulfill the restrictions are adjusted to be within the problem borders. Now, the objective function is evaluated for all the individuals of T and SubP and the best individuals in both populations are identified (*T*(*i*) and *SubP*(*i*)). If *T*(*i*) is a better solution than *SubP*(*i*), then *SubP*(*i*) is updated (*SubP*(*i*) = *T*(*i*)), and population P is updated from SubP; otherwise, TempP, F and M are calculated to generate a new T. Finally, the best individual of P is evaluated to verify if it fulfills the stop criteria; else, a new algorithm iteration starts.

**Figure 3.** Flow chart of the weighted differential evolution algorithm.

#### 3.2.4. Symbiotic Organism Search (SOS)

This is a metaheuristic optimization algorithm inspired by the different interactions between the living beings of different species that share an ecosystem, named symbiotic relationships. The first phase of the algorithm is to create an initial population of particles that represents an ecosystem, and each particle represents an organism. Then, it is necessary to define some parameters, like the solution boundaries and the stop criteria, before starting the iterative process, which is based on three types of relationships: mutualism, commensalism, and parasitism. In the mutualism phase, two organisms participate, *Xi* and *Xj*, where *Xi* is the *i*-th organism into the ecosystem and *Xj* is a randomly chosen organism within the ecosystem. In this phase, both organisms could be benefited from each other through the creation of a new vector (*Mvec*) described in (18)

$$M\_{\text{rec}} = \frac{X\_i + X\_j}{2}.\tag{18}$$

Then, the algorithm calculates two new candidates to replace *Xinew* and *Xjnew* by using (19) and (20), which model the mutualist symbiotic relationship between two organisms. In (19) and (20) *rand*(0, 1) is a vector with random numbers between 0 and 1, and *Xbest* is the organism with the lowest value of the cost function, also known as the organism with the best adaptation in the ecosystem. Moreover, in mutualist relationships some species benefit more than the others; this phenomenon is represented by the factors *BF*<sup>1</sup> and *BF*<sup>2</sup> in (19) and (20), which are randomly assigned with 1 or 2. Once *Xinew* and *Xjnew* are calculated, the algorithm evaluates the objective function for both organisms. If *Xinew* and *Xjnew* are better solutions than *Xi* and *Xj*, then *Xinew* and *Xjnew* replace *Xi* and *Xj* in the ecosystem.

$$X\_{invw} = X\_i + rand(0, 1) \cdot (X\_{best} - M\_{\text{rev}} \cdot BF\_1). \tag{19}$$

$$X\_{jucw} = X\_j + rand(0, 1) \cdot (X\_{\text{best}} - M\_{\text{rec}} \cdot BF\_2). \tag{20}$$

In the commensalism phase, *Xi* is the same organism used in the previous phase, while *Xj* is an organism randomly selected from the ecosystem, but it is different than *Xj* used in the mutualism stage. The commensalism relationship between *Xi* and *Xj* is represented by (21), where *Xi* benefits from *Xj*. If the objective function of the *Xinew* is less than the one of the objective function of *Xi*, then *Xi* is updated in the ecosystem.

$$X\_{inrw} = X\_i + rand\left(-1, 1\right) \cdot \left(X\_{best} - X\_j\right). \tag{21}$$

The last symbiotic relationship is parasitism, where the organism *Xj* is aleatory chosen from the ecosystem and *Xi* is the same organism used in the commensalism phase, but with a random modification of one of its components. The modified *Xi* is named parasite vector and if its objective function if better than *Xj*, then the parasite vector replaces *Xj*; otherwise, the parasite vector is discarded, and it is said that *Xj* developed immunity to the parasite. Once the three phases end for the *i*-th organism, the algorithm moves to the next organism to repeat the process until all the organisms in the ecosystem have been processed. At this point, one iteration if finished and the stop criteria are evaluated, if no stop condition is met, then the algorithm starts a new iteration. The detailed explanation of this algorithm is shown in [37]. Keep in mind that the evaluated algorithms (many of them) are sensitive to the assumed values for their multiple parameters. It is so sensitive that if they are not chosen with caution, they may not even converge to the expected solution even for PV arrays operating under uniform conditions. This aspect is a notable disadvantage of this metaheuristic solution strategy. Until now, there is no solid mathematical basis to predict the convergence and influence of these parameters. It is certainly possible that the results presented could have been obtained in shorter computation times, having selected other values for the algorithm parameters. However, it is worth noting that we tune the optimization algorithms in order to obtain the best results

for all the algorithms. Particularly, the parameters of the metaheuristic methods selected for this work were defined randomly.

#### **4. Simulations Results**

This section introduces the simulation results for small, medium and large generators, i.e., PV generators formed by a single string with 6, 36 and 72 series-connected submodules, respectively. All the simulated generators are formed by Trina Solar TMS-PD05 de 270 W [38] modules, which are composed of three series-connected submodules. Moreover, the bypass diodes of the three submodules in this section are GF3045T [39]. It is worth noting that the generators considered in this section only have one string. This is because for generators with M strings each string is independently solved, which is equivalent to solve one string M times and then to sum the currents of all the strings to calculate the array current. The SDM parameters in STC are calculated, with the procedure proposed in [40], from the module datasheet information. Then, the parameters are adjusted for a module temperature of 44 ◦C (*T* = 44 ◦C) and an irradiance of 1 kW/m2 (*G* = 1 kW/m2), by following the procedure proposed in [31], obtaining the following values: *Iph* = 9.311 A, *Isat* = 23.782 *η*A, *η* = 1.097, *Rs* = 0.088 Ω, and *Rp* = 246.670 Ω. Afterward, the bypass diode parameters are calculated from the datasheet information as follows: *Isat*,*bd* = 851, 540 μA and *ηbd* = 1.634. The irradiance condition of each submodule in the generator is represented by using the linear dependence on *Iph* with the submodule irradiance. Therefore, the irradiance reduction in a submodule produced by some kind of shading is represented by a coefficient *Pirr* that multiplies *Iph* and varies between 0 and 1. Thus the reduction of *Iph* is proportional to the reduction of the submodule irradiance. The *Pirr* coefficients of all the string submodules are grouped in a vector (*P*- *irr*) that represents the shading conditions of the string. The solutions obtained with the four optimization methods used in this paper are compared with the solution of the equivalent electrical circuit (EEC) of each generator. Those EECs for the small, medium, and large generators are implemented in Simulink of Matlab; hence, the errors calculated for the optimization algorithms are calculated concerning the EEC solution.

#### *4.1. Small Generator*

This generator is formed by two modules connected in series that correspond to six submodules. This low power generator may be used in grid-connected applications, by using a microinverter, or in stand-alone applications to charge a battery or to feed a set of lights. The simulation results of this generator operating under homogeneous (i.e., *P*- *irrh* = [111111]) and partial shading (i.e., *P*- *irrnh* = [0.8 0.8 0.8 0.8 0.3 0.3]) conditions are introduced in Figures 4 and 5, respectively. Those figures show the I–V curves and the error in the current calculation for all the implemented optimization methods, which use the stop criteria shown in Table 2. Moreover, the evaluation criteria of the simulation results (see Table 3) are: the computation time, the root mean square error (RMSE), and the number of evaluations of the objective function (*Feval*).

**Table 2.** Stop criteria for the optimization methods used for the small generator.


(**b**) **Figure 4.** Simulation results of the small generator (six submodules) under homogeneous conditions. (**a**) Current vs. voltage (I–V) curve. (**b**) Current calculation error.

**Figure 5.** Simulation results of the small generator (six submodules) under partial shading conditions. (**a**) Current vs. voltage (I–V) curve. (**b**) Current calculation error.

In Figure 4a it can be observed that all the implemented optimization methods can reproduce the I–V curve of a small string (six submodules) under homogeneous conditions. Additionally, the errors in the string current calculation are similar for the different methods, as shown in Figure 4b and Table 3. When the generator operates in partial shading conditions, three of the four algorithms (TRD, LMA, and SOS) can solve the system of implicit equations for all the voltages to reproduce the I–V with similar errors (see Figure 5 and Table 3). However, the WDE algorithm presents significant errors in the current calculation for voltages between 20 V and 55 V, as it is evidenced in Figure 5 and Table 3. Therefore, WDE is not a suitable algorithm to solve the model of a small generator in partial shading conditions. Note that under homogeneous and partial shading conditions the computation time and number of cost function evaluations (*Feval*) of the deterministic optimization methods (TRD and LMA) are significantly less than the computation time and *Feval* of the metaheuristic algorithms (WDE and SOS), as shown in Table 3. Moreover, the RMSE values of the methods able to solve the optimization problem provide similar errors in the I–V curve reproduction.


**Table 3.** Evaluation criteria for the small generator operating under homogeneous and partial shading conditions.

#### *4.2. Medium Generator*

In this case, the generator is formed by one string with 12 modules (36 submodules) connected in series, which is a typical configuration for a residential application with a string inverter (e.g., ABB UNO-7.6-TL-OUTD). For this generator, the stop criteria of the optimization algorithms are modified concerning the previous subsection, as shown in Table 4, to improve the performance of the algorithms.

**Table 4.** Stop criteria for the optimization methods used for the medium and large generators.


The medium generator I–V curves for homogeneous and partial shading conditions are presented in Figures 6 and 7, respectively; while the computation time, RMSE and *Feval* are introduced in Table 5. Moreover, the vectors that define the homogeneous (*P*- *irrh*) and partial shading conditions (*P*- *irrnh*) have 36 elements and are defined as follows: all elements in *P*- *irrh* are 1 and *P*- *irrnh* has 24 elements equal to 0.8, 6 elements equal to 0.6, and 6 elements equal to 0.2.

The results in Figures 6 and 7 and Table 5 show that the deterministic methods (TRD and LMA) solve the system of equations of the string for homogeneous and partial shading conditions. Additionally, the RMSE values for the I–V curves are the same for both deterministic methods (less than 0.1 A); nevertheless, computation times of TRD are less than the computation times of LMA for both homogeneous and partial shading conditions. In general, metaheuristic methods present greater errors and computation times than the deterministic methods in the I–V curve's calculation. For homogeneous conditions, the current errors are concentrated in medium and high generator voltages (from 250 V to 400 V), as shown in Figure 6; those errors generate RMSE values between 13% and 60% greater than the ones of the deterministic methods. Whilst for partial shading conditions, the RMSE values of the metaheuristic methods between 15.3 and 21.3 times greater than the RMSE values of the deterministic methods, and the errors occur in the entire range of the generator voltage (see Figure 7). As observed, the computation times of the metaheuristic algorithms are two orders of magnitude greater than the computation times of the deterministic methods, which is also reflected in the number of cost function evaluations (*Feval*). Finally, the results in Table 5 indicate that the implemented metaheuristic optimization methods are not suitable for solving the implicit model of a medium PV generator.

**Figure 6.** Simulation results of the medium generator (36 submodules) under homogeneous conditions. (**a**) Current vs. voltage (I–V) curve. (**b**) Current calculation error.

(**b**) **Figure 7.** Simulation results of the medium generator (36 submodules) under partial shading conditions. (**a**) Current vs. voltage (I–V) curve. (**b**) Current calculation error.

**Table 5.** Evaluation criteria for the medium generator operating under homogeneous and partial shading conditions.


#### *4.3. Large Generator*

This generator is formed by a single string with 24 modules (72 submodules) connected in series, which can be found in industrial applications or high-power generators that use inverters whit maximum input voltages of 1 kV (e.g., Sunny Tripower 20000TL). The stop criteria are the same shown in Table 4 and the vectors that define homogeneous (*P*- *irrh*) and partial shading (*P*- *irrnh*) conditions have 72 elements defined as follows: all elements in (*P*- *irrh*) are 1 and (*P*- *irrnh*) has 30 elements equal to 0.8, 30 elements equal to 0.6, and 12 elements equal to 0.2. The generator I–V curves and the errors in the current calculation are shown in Figures 8 and 9 for homogeneous and partial shading conditions, respectively; while the computation time, RMSE and *Feval* are presented in Table 6.

(**b**)

**Figure 8.** Simulation results of the large generator (72 submodules) under homogeneous conditions. (**a**) Current vs. voltage (I–V) curve. (**b**) Current calculation error.

(**b**) **Figure 9.** Simulation results of the large generator (72 submodules) under partial shading conditions. (**a**) Current vs. voltage (I–V) curve. (**b**) Current calculation error.

**Table 6.** Evaluation criteria for the large generator operating under homogeneous and partial shading conditions.


Once more, the results in Figures 8 and 9 and Table 6 indicate that the deterministic optimization methods (TRD and LMA) are capable to reproduce the I–V curves for the evaluated conditions with RMSE less than 0.10 A. In homogeneous conditions, the computation time of LMA is 12.9% less than the one of TRD; while in partial shading conditions the computation time of TRD is 15.5% less than the computation time of LMA. Moreover, the values of *Feval* for TRD and LMA have similar behavior to the computation time. The metaheuristic algorithms show errors and computation time greater than the deterministic algorithms in the reproduction of the I–V curves in the evaluated operating conditions. For homogeneous conditions (see Figure 8) the biggest errors are present in high voltages; as consequence, the RMSE values of the metaheuristic methods are significantly greater (between 132% and 700%) than the deterministic algorithms. In partial shading conditions (see Figure 9) the errors are in the whole voltage range, making the RMSE values of the metaheuristic methods between 24.2 and 31.9 times the RMSE values of the deterministic methods. Additionally, the computation time of the metaheuristic algorithms is two orders of magnitude greater than the one of the deterministic algorithms, which is similar to the differences in *Feval* (see Table 6). After the evaluation of the four optimization algorithms for small, medium and large PV arrays, it seems that metaheuristic optimization methods do not worth their evaluations because the evaluated deterministic optimization methods outperform the evaluated metaheuristic methods. However, those results could not be obtained if the comparison is not performed; therefore, we consider that the evaluation is proposed in the paper is valid, since it is necessary to explore different optimization methods to reduce the computation time of the model. Such a reduction is important for applications like reconfiguration, energy production estimation, evaluation of MPPT techniques, among others.

#### **5. Experimental Results**

The performance of the optimization methods was validated with experimental data, which were obtained with the test bench shown in Figure 10 This test bench if formed by a PV with Matlab, a programmable electronic load BK Precision 8500, and a Trina Solar TMS-PD05 de 270 W [38] PV module. The PV module is the same one used in the simulations (formed by 3 submodules); however, the SDM parameters of the submodules (*Iph* = 9.316 A, *Isat* = 1.549 nA, *η* = 0.972, *Rs* = 0.175 Ω and *Rp* = 400.804 Ω) and the bypass diodes (*Isat*,*bd* = 851, 540 μA and *ηbd* = 0.2261) are calculated from experimental I-V curves. The experimental validation considers two cables AWG 12 of 28.96 m each, which introduce losses that can be represented by resistance in series with the module of 313 mΩ. Therefore, such resistance was included in the *Rs* parameter of each submodule; thus, the ohmic losses introduced by the cables are lumped in the *Rs* of the submodules.

**Figure 10.** Test bench used to obtain the experimental data. (**a**) PV module with three submodules. (**b**) Electronic load. (**c**) PC with Matlab.

The small generator used for the experimental validation was formed by one string of three submodules, which was evaluated for homogeneous (C-1) and partial shading (C-2) defined by the vectors (*P*- *irr*,1) = [0.4 0.4 0.4] and (*P*- *irr*,2) = [0.620 0.598 0.210], respectively. Moreover, the stop criteria used for the optimization algorithms were the same ones defined in Table 2. The experimental results for conditions C-1 and C-2 are summarized in Figure 11 and Table 7. Figure 11 introduces the I–V curves obtained with the experimental test bench, the four optimization algorithms (TRD, LMA, WDE, and SOS), and the equivalent electrical circuit (EEC) for conditions C-1 and C-2; while Table 7 shows the computation time, RMSE values, and *Feval* for the different methods.

In general, the experimental results agree with the simulation results obtained for small generators (Section 4.1) regarding computation time and current errors. On the one hand, in homogeneous conditions (C-1), the four optimization algorithms generate I–V curves that fit the experimental data with practically the same RMSE values. Moreover, the error in the current calculation for condition C-1 increases for high voltages (circles in Figure 11b), which agrees with the behavior shown in the simulation results for small generators (see Figure 4b). On the other hand, in partial shading conditions (C-2) three algorithms (TRD, LMA and SOS) reproduce the experimental I–V curves with the same RMSE values. However, the WDE presents significant errors in the current calculation for medium voltages (blue x in Figure 11b), which agree with the simulation results in Figure 5 of Section 4.1 (small generators). The behavior of the computation time in the experimental validation is similar to the one presented in simulation results. The computation time of the metaheuristic algorithms are, approximately, two orders of magnitude greater than the deterministic algorithms. Additionally, the computation time of the deterministic methods was less than the EEC in a proportion similar to the one presented above. The computation time of the four optimization algorithms is reflected in the *Feval* values, which are similar to the values in the simulation results. Finally, Table 7 shows that the smallest errors were obtained with the EEC model, which verifies the usefulness of this model to be used as the reference in the simulation results.


**Table 7.** Evaluation criteria for a small generator operating under homogeneous and partial shading conditions from experimental results.

(**b**)

**Figure 11.** Experimental results in a small generator (3 submodules) under homogeneous and partial shading conditions. (**a**) Current vs. voltage (I–V) curve. (**b**) Current calculation error.

#### **6. Conclusions**

The solution of the implicit model of a PV generator in SP configuration as an optimization problem has been proposed and solved with deterministic and metaheuristic algorithms. The implicit model of an NxM SP generator (M strings with N submodules each) is formed by M systems of N+2 implicit equations, where each system of equations corresponds to one string and it can be solved independently. Moreover, in a system of N + 2 implicit equations, the unknowns correspond to the N submodules voltages, the blocking diode voltage, and the string current. To solve a system of N+2 implicit equations, the paper proposes an objective function and defines the upper and lower restrictions for the unknowns. Finally, the paper evaluates two deterministic optimization algorithms (TRD and LMA) and two metaheuristic algorithms (WDE and SOS) to solve the optimization problem and generate the I–V curves of small, medium, and large generators.

The proposed objective function and the performance of the four evaluated methods were validated with simulation and experimental results. The simulations show that the deterministic optimization methods can solve the optimization problem for small, medium and large generators operating in homogeneous and partial shading conditions; therefore, with these algorithms, it is possible to reproduce the generator I–V curves. Moreover, simulation results show that the

metaheuristic algorithms are not able to solve the optimization problem in all cases. On the one hand, SOS correctly reproduces the I–V curves for small generators under homogeneous and non-homogeneous conditions and medium generators in homogeneous conditions. Nevertheless, for the other simulation scenarios, the errors were significant. On the other hand, WDE only solved the optimization problem for small generators under homogeneous conditions; while the errors were large for the other cases concerning the deterministic methods. We believe that metaheuristic algorithms performed worse because due to the random selection of parameter values instead of the nature of the selected algorithms (weighted differential evolution and search for symbiotic organisms). This could even be explained based on the No Free Lunch Theorems (NFL) for Optimization. Additionally, experimental results confirm simulation results for small generators in homogeneous and non-homogeneous conditions, where all the methods correctly reproduce the I–V curves except for WDE under partial shading conditions. Moreover, experimental results illustrate the usefulness of the EEC model as a reference for simulation results.

Finally, the results suggest that the solution of the implicit model of SP generators should be performed with deterministic algorithms, which suggests the necessity of evaluating other metaheuristic optimization methods to solve the proposed problem and reduce the computation time. Those presented results should be used to evaluate different solutions for energy yield calculations or economic analysis of a PV system; hence as future work, an addition to the presented algorithms with energy predictions or economic calculations will be considered.

**Author Contributions:** Conceptualization and methodology was performed by J.D.B.-R. and R.C.; software was developed by L.M.P.A. and J.D.B.-R.; simulation were performed by L.M.P.A.; experiments were developed by L.A.T.G. and D.G.-M.; analyses were performed by L.M.P.A., J.D.B.-R., R.C., L.A.T.G. and D.G.-M.; writing of the original paper draft was performed by L.M.P.A., J.D.B.-R., R.C.; writing of the reviewed draft was made by R.C., L.A.T.G. and D.G.-M. All authors have read and agree to the published version of the manuscript.

**Funding:** This work was funded by Universidad Industrial de Santander under the scholarship "Crédito Condonable para Estudiantes de Posgrado" and Universidad Nacional de Colombia Sede Manizales under de project "Desarrollo de un micro-inversor solar multietapa de inyección a red con capacidad de seguimiento del punto de máxima potencia" code HERMES 43049.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


c 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

#### *Article*
