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

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.

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,15,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,15,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 the N + 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,15,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:

• A solution of the system of N + 2 implicit equation as an optimization problem with restrictions, which can be used to explore other optimization methods to reduce the model’s solution time.
• Performance comparison of metaheuristic and deterministic optimization algorithms for solving the problem provides guidelines to continue exploring other optimization methods to reduce the computation time.
• A novel application of the algorithms of weighted differential evolution and search for symbiotic organisms to solve the implicit model of series-parallel photovoltaic arrays is presented.

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, Section 4 and Section 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 $N_s$ 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, $I_{ph}$ represents the current generated by the photovoltaic effect, the diode D models the nonlinearities of the PN junction, and the resistors $R_p$ and $R_s$ 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.

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 $I_{sat}$ and $\eta$ are the inverse saturation current, and ideality factor of D, respectively, $V_{t,d}$ is the diode thermal voltage and $I_{bd}$ is the current through BD. The thermal voltage is defined as $V_{t,d}=k·T/q$, where k is the Boltzmann constant, q is the electron charge and T is the cell’s temperature in K. Moreover, $I_{db}$ is defined in (2), where $I_{sat,bd}$ and ${\eta }_{bd}$ are the inverse saturation current and ideality factor of BD. It is worth noting that $V_{t,d}$ 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_{ph}-I_{sat}\left[exp\left(\frac{V+\left(I-I_{bd}\right)·R_s}{N_s·\eta ·V_{t,d}}\right)-1\right]-\frac{V+\left(I-I_{bd}\right)·R_s}{R_p}+I_{bd}=0.$$

(1)

$$I_{bd}=I_{sat,bd}·\left[exp\left(\frac{-V}{{{\eta }_{bd}·V}_{t,d}}\right)-1\right].$$

(2)

Although the SDM parameters (i.e., $I_{ph},I_{sat},V_{t,d},R_s$ and $R_h$) may vary with the irradiance and temperature [31,32,33], $I_{ph}$ 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 ($N_m$) and a blocking diode connected in series. In turn, each module is formed by one or more submodules connected in series ($N_{sm}$); then, each string is formed by N submodules, where $N=N_m·N_{sm}$. The generator is connected to a power converter, which sets the generator voltage ($V_{gen}$) to perform the MPP tracking; therefore, $V_{gen}$ is known for the model and the objective is to find: the generator current ($I_{gen}$), the current of each string, the voltage of each submodule, and the voltage of each blocking diode.

Considering the $V_{gen}$ 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 $\overrightarrow{V}$ is a vector formed by the voltages of the N submodules and the blocking diode ($\overrightarrow{V}=\left[V_1,V_2,\cdots ,V_N,V_{blk}\right]$) and $I_{str}$ is the string current.

$$F\left(\overrightarrow{V},I_{str}\right)=\left[\begin{array}{c}{f}_1\left(V_1,I_{str}\right)\\ f_2\left(V_2,I_{str}\right)\\ \begin{array}{c}⋮\\ f_N\left(V_N,I_{str}\right)\\ \begin{array}{c}{f}_{N+1}\left(V_{blk},I_{str}\right)\\ f_{N+2}\left(\overrightarrow{V}\right)=\sum _{i=1}^{N+1}{V}_i-V_{gen}\end{array}\end{array}\end{array}\right]=\overrightarrow{0}.$$

(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 $I_{sat,blk}$, and ${\eta }_{blk}$ are the inverse saturation current and the ideality factor of the blocking diode, respectively, and $V_{t,d}$ 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 $V_{gen}$.

$$f_{N+1}\left(V_{blk},I_{str}\right)=-I_{str}+I_{sat,blk}·\left(exp\left(\frac{-V_{blk}}{{\eta }_{blk}·V_{t,d}}\right)-1\right)=0.$$

(4)

All the electrical variables of the string are obtained by solving (3), i.e., $I_{str}$ and $\overrightarrow{V}$; hence, after solving all the strings, Igen can be calculated using (5).

$$I_{gen}=\sum _{j=1}^M{I}_{str,j}.$$

(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 $\mathbb{R}$ be the set of real numbers and $\mathbf{L}$ a non-empty subset of ${\mathbb{R}}^n$, i.e., $\mathbf{L}\subset {\mathbb{R}}^n$; then, the system of nonlinear equations given in (3) can be expressed as shown in (6), where $\overrightarrow{x}\in {\mathbb{R}}^n$ and $f_k:\mathbf{L}\to \mathbb{R},\forall k=1,2,\dots ,N+2$. The solution of (6) is given by all vectors $\overrightarrow{a}\in {\mathbb{R}}^n$, such that $f_k\left(\overrightarrow{a}\right)=0,\forall k=1,2,\dots ,N+2$.

$$\left\{\begin{array}{c}{f}_1\left(\overrightarrow{x}\right)=0\\ f_2\left(\overrightarrow{x}\right)=0\\ \begin{array}{c}⋮\\ f_{N+2}\left(\overrightarrow{x}\right)=0\end{array}\end{array}\right..$$

(6)

Now, let’s consider the function $F_{obj}:\mathbf{L}\to \mathbb{R}$ a function defined as shown in (7).

$$F_{obj}=\sum _{i=1}^{N+2}{f}_i{\left(\overrightarrow{x}\right)}^2,with\overrightarrow{x}\in \mathbf{L}.$$

(7)

Note that $F_{obj}$ is a well-defined function that guarantees the existence of a minimum, since $F_{obj}\left(\overrightarrow{x}\right)\ge 0$ for all $\overrightarrow{x}\in \mathbf{L}$ is a totally ordered set where exists the infimum of $F_{obj}\left(\mathbf{L}\right)\in \mathbb{R}$ and this infimum is non-negative, i.e., $inf\left\{F_{obj}\left(\mathbf{L}\right)\right\}\ge 0$. Therefore, the minimum of $F_{obj}$ over $\mathbf{L}$ exists and it is a non-negative minimum.

Theorem 1.

Suppose that the system of equations given in (3) has a solution $\mathbf{L}$ and $(\overrightarrow{a}\in \mathbf{L})$; hence, is a solution of (3) if and only if, $\overrightarrow{a}$ minimizes $F_{obj}$ given in (7).

Proof. 

If $\overrightarrow{a}$ is a solution of the system of equations given in (3), then $f_k\left(\overrightarrow{a}\right)=0,\forall k=1,2,\dots ,N+2$, therefore, $F_{obj}\left(\overrightarrow{a}\right)=0$. Moreover, considering that $F_{obj}\left(\overrightarrow{x}\right)\ge 0$ for all $\overrightarrow{x}\in \mathbf{L}$, then $\overrightarrow{a}$ is the minimum point of $F_{obj}$. Now, if $\overrightarrow{a}$ minimizes $F_{obj}$ but it is not a solution of (3), then $F_{obj}\left(\overrightarrow{a}\right)$ must be positive, since $F_{obj}\left(\overrightarrow{x}\right)\ge 0$ for all $\overrightarrow{x}\in \mathbf{L}$. Taking into account that (3) has a solution in $\mathbf{L}$, then there is a solution ${\overrightarrow{x}}^{\ast }\in \mathbf{L}$ such that $F_{obj}\left({\overrightarrow{x}}^{\ast }\right)=0$ and ${\overrightarrow{x}}^{\ast }\ne \overrightarrow{a}$. Hence, $F_{obj}\left({\overrightarrow{x}}^{\ast }\right)<F_{obj}\left(\overrightarrow{a}\right)$, which contradicts that $\overrightarrow{a}$ is a minimum point for $F_{obj}$. ▢

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 ($\overrightarrow{V}$) and the string current ($I_{str}$). However, to minimize (7) using numerical methods, it is necessary to establish the limits of $I_{str}$ and each unknown in $\overrightarrow{V}$ as shown in

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

Variable	Lower Limit	Upper Limit
$V_1\dots V_N$	$-V_{bd,max}$	$V_{OC,STC}$
$V_{blk}$	$-V_{blk,max}$	$0V$
$I_{str}$	$0A$	$I_{SC,STC}$

.

The lower limit of the string current is 0 A; while the upper limit is $I_{SC,STC}$, which corresponds to the submodule short-circuit current in Standard Test Conditions (STC), i.e., irradiance of 1 kW/m² and the submodule temperature is 25 °C (i.e., $T$ = 298.15 K). Moreover, the upper limit of the submodules voltages ($V_1\dots V_N$) is open-circuit voltage in STC ($V_{OC,STC}$), which corresponds to the submodule voltage when $I=0$ A. The lower limit of $V_1\dots V_N$ is the negative value of the maximum bypass diode voltage ($V_{bd,max}$), defined in (8), which corresponds to the bypass diode voltage when $I_{bd}=I_{SC,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 $I_{str}=0$ A; while the lower limit corresponds to the negative value of the maximum blocking diode voltage ($V_{blk,max}$) defined in (9).

$$V_{bd,max}={{\eta }_{bd}·V}_{t,d}·ln\left(\frac{I_{SC,STC}}{I_{sat,bd}}+1\right).$$

(8)

$$V_{blk,max}={{\eta }_{blk}·V}_{t,d}·ln\left(\frac{I_{SC,STC}}{I_{sat,blk}}+1\right).$$

(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\left(x\right)$, and searches the point $x^{\ast }$ that minimizes the function f through the iterative process indicated in (10), where $x_k$ is the actual position of the search agent and $p_k$ is the step of the k-th iteration.

$$x_{k+1}=x_k+p_k.$$

(10)

To determine the magnitude and direction of the vector $p_k$ that approaches $x_{k+1}$ to $x^{\ast }$ is posed as an optimization problem, as shown in (11), where $m_k\left(p\right)$ is an approximate function of $f\left(x\right)$ formed by the two first terms of the Taylor series in a region $Q_k$ that contains $x_k$.

$$\underset{p}{min}\left\{m_k\left(p\right),p\in Q_k\right\}.$$

(11)

If $m_k\left(p\right)$ is a good approximation of the objective function in a region $Q_k$, then this region is denominated thrust-interval and it is defined as $Q_k=\left\{\mathbf{p}\mid |\mathbf{p}|<{\Delta }_k\right\}$, where ${\Delta }_k$ is the radius of the thrust-region in the k-th iteration. Once the vector $p_k^{\ast }$ is found, i.e., the solution of (11), it is necessary to verify the fulfillment of the condition expressed in (12).

$$f\left(x_k+p_k^{\ast }\right)<f\left(x_k\right).$$

(12)

If condition (12) is fulfilled, the actual position is updated as $x_{k+1}=x_k+p_k^{\ast }$; then, the next step is to calculate a new model $m_{k+1}$ and a new trust-region $Q_{k+1}$ from the updated point $x_{k+1}$. Otherwise, if condition (12) is not fulfilled, the algorithm keeps the actual position, i.e., $x_{k+1}=x_k$, and the trust-region is reduced by reducing its radius (${\Delta }_{k+1}<{\Delta }_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 ${\mu }_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}{∥J\left(x_k\right)(x-x_k)+r{x}_k∥}_2^2.$$

(13)

$$s.t.∥x-x_k∥\le {\Delta }_k.$$

(14)

This is a least-squares problem, which is linear with restrictions because $r\left(x_k\right)$ is a linearized model valid in the trust-region defined by ${\Delta }_k$. The solution to the problem described in (13) and (14) is obtained by solving (15), where $s=x-x_k$, and the position of the next iteration is given by (16) [35].

$$(J{\left(x_k\right)}^{T}J\left(x_k\right)+{\mu }_{k}I)s=-J{\left(x_k\right)}^{T}r\left(x_k\right).$$

(15)

$$x_{k+1}=x_k-{(J{\left(x_k\right)}^{T}J\left(x_k\right)+\mu I)}^{-1}·J{\left(x_k\right)}^{T}r\left(x_k\right).$$

(16)

When $∥{\left(J{\left(x_k\right)}^{T}J\left(x_k\right)\right)}^{-1}J{\left(x_k\right)}^{T}r\left(x_k\right)∥\le {\Delta }_k$, then ${\mu }_k=0$ and (16) is solved by $s_k$; in this case, the best search direction is the same as the Gauss–Newton method. Otherwise, ${\mu }_k>0$ and the solution of (16) fulfills $s_k={\Delta }_k$; with this condition, it is possible to find the solution of (16) and the value of ${\mu }_k$. When ${\mu }_k\gg 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}Is=-J{\left(x_k\right)}^{T}r\left(x_k\right).$$

(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_{\left(i\right)}$ and $Sub{P}_{\left(i\right)}$). If $T_{\left(i\right)}$ is a better solution than $Sub{P}_{\left(i\right)}$, then $Sub{P}_{\left(i\right)}$ is updated ($Sub{P}_{\left(i\right)}$ = $T_{\left(i\right)}$), 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.

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, $X_i$ and $X_j$, where $X_i$ is the i-th organism into the ecosystem and $X_j$ 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 ($M_{vec}$) described in (18)

$$M_{vec}=\frac{X_i+X_j}{2}.$$

(18)

Then, the algorithm calculates two new candidates to replace $X_{inew}$ and $X_{jnew}$ 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 $X_{best}$ 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 $B{F}_1$ and $B{F}_2$ in (19) and (20), which are randomly assigned with 1 or 2. Once $X_{inew}$ and $X_{jnew}$ are calculated, the algorithm evaluates the objective function for both organisms. If $X_{inew}$ and $X_{jnew}$ are better solutions than $X_i$ and $X_j$, then $X_{inew}$ and $X_{jnew}$ replace $X_i$ and $X_j$ in the ecosystem.

$$X_{inew}=X_i+rand(0,1)·(X_{best}-M_{vec}·B{F}_1).$$

(19)

$$X_{jnew}=X_j+rand(0,1)·(X_{best}-M_{vec}·B{F}_2).$$

(20)

In the commensalism phase, $X_i$ is the same organism used in the previous phase, while $X_j$ is an organism randomly selected from the ecosystem, but it is different than $X_j$ used in the mutualism stage. The commensalism relationship between $X_i$ and $X_j$ is represented by (21), where $X_i$ benefits from $X_j$. If the objective function of the $X_{inew}$ is less than the one of the objective function of $X_i$, then $X_i$ is updated in the ecosystem.

$$X_{inew}=X_i+rand\left(-1,1\right)·\left(X_{best}-X_j\right).$$

(21)

The last symbiotic relationship is parasitism, where the organism $X_j$ is aleatory chosen from the ecosystem and $X_i$ is the same organism used in the commensalism phase, but with a random modification of one of its components. The modified $X_i$ is named parasite vector and if its objective function if better than $X_j$, then the parasite vector replaces $X_j$; otherwise, the parasite vector is discarded, and it is said that $X_j$ 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/m² ($G=1$ kW/m²), by following the procedure proposed in [31], obtaining the following values: $I_{ph}=9.311$ A, $I_{sat}=23.782\eta$A, $\eta =1.097,R_s=0.088\Omega$, and $R_p=246.670\Omega$. Afterward, the bypass diode parameters are calculated from the datasheet information as follows: $I_{sat,bd}=\mathrm{851,540}\mathsf{\mu }$A and ${\eta }_{bd}=1.634$. The irradiance condition of each submodule in the generator is represented by using the linear dependence on $I_{ph}$ with the submodule irradiance. Therefore, the irradiance reduction in a submodule produced by some kind of shading is represented by a coefficient $P_{irr}$ that multiplies $I_{ph}$ and varies between 0 and 1. Thus the reduction of $I_{ph}$ is proportional to the reduction of the submodule irradiance. The $P_{irr}$ coefficients of all the string submodules are grouped in a vector ($\overrightarrow{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., $\overrightarrow{P_{irrh}}=\left[111111\right]$) and partial shading (i.e., $\overrightarrow{P_{irrnh}}=\left[0.80.80.80.80.30.3\right]$) conditions are introduced in Figure 4 and Figure 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. Stop criteria for the optimization methods used for the small generator.

Stop Criteria	TRD	LMA	WDE	SOS
Number of particles	N/A	N/A	20	50
Maximun iteration number	$1\times {10}^3$	$1\times {10}^3$	$1\times {10}^4$	$5\times {10}^4$
Tolerance in the solution variation	$1\times {10}^{-3}$	$1\times {10}^{-9}$	N/A	N/A
Tolerance in the objective function	N/A	N/A	$1\times {10}^{-18}$	$1\times {10}^{-18}$

. Moreover, the evaluation criteria of the simulation results (see

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

Homogeneous Conditions
Algorithm	Computation Time (s)	RMSE (A)	${\mathit{F}}_{\mathit{eval}}$
EEC(REF)	4.9	0	0
TRD	1.1	0.0957	$2.03\times {10}^4$
LMA	1.2	0.0957	$2.32\times {10}^4$
WDE	308.0	0.0957	$1.36\times {10}^7$
SOS	164.0	0.0957	$7.48\times {10}^6$
Partial Shading Conditions
EEC(REF)	5.8	0	0
TRD	1.6	0.0387	$2.69\times {10}^4$
LMA	1.6	0.0387	$3.12\times {10}^4$
WDE	281.8	0.3488	$1.24\times {10}^7$
SOS	161.3	0.0388	$7.55\times {10}^6$

) are: the computation time, the root mean square error (RMSE), and the number of evaluations of the objective function ($F_{eval}$).

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 ($F_{eval}$) of the deterministic optimization methods (TRD and LMA) are significantly less than the computation time and $F_{eval}$ 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.

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. Stop criteria for the optimization methods used for the medium and large generators.

Stop Criteria	TRD	LMA	WDE	SOS
Number of particles	N/A	N/A	20	100
Maximun iteration number	$4\times {10}^3$	$4\times {10}^3$	$1\times {10}^4$	$1\times {10}^5$
Tolerance in the solution variation	$1\times {10}^{-5}$	$1\times {10}^{-9}$	N/A	N/A
Tolerance in the objective function	N/A	N/A	$1\times {10}^{-18}$	$1\times {10}^{-18}$

, to improve the performance of the algorithms.

The medium generator I–V curves for homogeneous and partial shading conditions are presented in Figure 6 and Figure 7, respectively; while the computation time, RMSE and $F_{eval}$ are introduced in

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

Homogeneous Conditions
Algorithm	Computation Time (s)	RMSE (A)	${\mathit{F}}_{\mathit{eval}}$
EEC(REF)	28.1	0	0
TRD	9.6	0.1034	$2.78\times {10}^5$
LMA	10.0	0.1034	$2.84\times {10}^5$
WDE	1898.8	0.1170	$5.24\times {10}^7$
SOS	1271.3	0.1660	$4.54\times {10}^7$
Partial Shading Conditions
EEC(REF)	28.3	0	0
TRD	15.9	0.0460	$4.56\times {10}^5$
LMA	20.3	0.0460	$5.07\times {10}^5$
WDE	1784.4	1.0273	$5.11\times {10}^7$
SOS	1296.9	0.7486	$4.54\times {10}^7$

. Moreover, the vectors that define the homogeneous ($\overrightarrow{P_{irrh}}$) and partial shading conditions ($\overrightarrow{P_{irrnh}}$) have 36 elements and are defined as follows: all elements in $\overrightarrow{P_{irrh}}$ are 1 and $\overrightarrow{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 Figure 6 and Figure 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 ($F_{eval}$). 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.

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 ($\overrightarrow{P_{irrh}}$) and partial shading ($\overrightarrow{P_{irrnh}}$) conditions have 72 elements defined as follows: all elements in ($\overrightarrow{P_{irrh}}$) are 1 and ($\overrightarrow{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 Figure 8 and Figure 9 for homogeneous and partial shading conditions, respectively; while the computation time, RMSE and $F_{eval}$ are presented in

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

Homogeneous Conditions
Algorithm	Computation Time (s)	RMSE (A)	${\mathit{F}}_{\mathit{eval}}$
EEC(REF)	61.8	0	0
TRD	42.0	0.1041	$1.13\times {10}^6$
LMA	36.6	0.1041	$8.25\times {10}^5$
WDE	5410.8	0.2415	$1.14\times {10}^8$
SOS	3047.3	0.8351	$9.10\times {10}^7$
Partial Shading Conditions
EEC(REF)	92.4	0	0
TRD	72.3	0.0404	$1.81\times {10}^6$
LMA	85.6	0.0408	$1.87\times {10}^6$
WDE	5586.4	1.3015	$1.13\times {10}^8$
SOS	3124.7	0.9886	$9.10\times {10}^7$

.

Once more, the results in Figure 8 and Figure 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 $F_{eval}$ 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 $F_{eval}$ (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 ($I_{ph}=9.316$ A, $I_{sat}=1.549$ nA, $\eta =0.972$, $R_s=0.175\Omega$ and $R_p=400.804\Omega$) and the bypass diodes ($I_{sat,bd}=\mathrm{851,540}\mathsf{\mu }$A and ${\eta }_{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$\Omega$. Therefore, such resistance was included in the $R_s$ parameter of each submodule; thus, the ohmic losses introduced by the cables are lumped in the $R_s$ of the submodules.

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 ($\overrightarrow{P_{irr,1}}$) = [0.4 0.4 0.4] and ($\overrightarrow{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. Evaluation criteria for a small generator operating under homogeneous and partial shading conditions from experimental results.

Scenario	Method	Computation Time (s)	RMSE (A)	${\mathit{F}}_{\mathbf{eval}}$
	1 EEC	3.04	0.0860	0
	2 TRD	0.67	0.0766	$1.18\times {10}^4$
C-1	3 LMA	0.78	0.0766	$1.33\times {10}^4$
	4 WDE	74.33	0.0766	$3.57\times {10}^6$
	5 SOS	98.17	0.0766	$4.75\times {10}^6$
	1 EEC	3.15	0.0898	0
	2 TRD	1.02	0.1012	$1.82\times {10}^4$
C-2	3 LMA	1.23	0.1012	$2.22\times {10}^4$
	4 WDE	83.40	0.3698	$3.87\times {10}^6$
	5 SOS	104.33	0.1012	$4.89\times {10}^6$

. 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 $F_{eval}$ 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 $F_{eval}$ 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.

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.