# Rethinking solar photovoltaic parameter estimation: global optimality
analysis and a simple efficient differential evolution method

Shuhua Gao111Department of Electrical and Computer Engineering, National
University of Singapore, 119077 Singapore<EMAIL_ADDRESS>Cheng
Xiang222Department of Electrical and Computer Engineering, National University
of Singapore, 119077 Singapore<EMAIL_ADDRESS>Yu Ming333Power Automation Pte
Ltd, 438B Alexandra Road, Alexandra TechnoPark, 119968 Singapore
<EMAIL_ADDRESS>Tan Kuan Tak444Engineering Cluster, Singapore Institute of
Technology, 10 Dover Drive, 138683 Singapore<EMAIL_ADDRESS>Tong Heng Lee555Department of Electrical and Computer Engineering, National
University of Singapore, 119077 Singapore<EMAIL_ADDRESS>
###### Abstract

A large variety of metaheuristics have been proposed for photovoltaic (PV)
parameter extraction. Our aim here is not to develop another new metaheuristic
method but to investigate two important yet rarely studied issues. We focus on
the two most widely used benchmark datasets and try to answer (i) whether the
globally optimal parameters have already been found; and (ii) whether a
significantly simpler metaheuristic, in contrast to currently sophisticated
ones, can achieve equally good performance. We address the former using a
branch and bound algorithm, which certifies the global minimum rigorously for
the single diode model (SDM) and locates a fairly tight upper bound for the
double diode model (DDM) on both datasets. These obtained values can serve as
useful references for evaluation and design of metaheuristic methods. However,
this algorithm is excessively slow and unsuitable for time-sensitive
applications (despite the great insights that it yields). Next, extensive
examination and comparison reveals that, perhaps surprisingly, a classic and
remarkably simple differential evolution (DE) algorithm can consistently
achieve the certified global minimum for the SDM and obtain the best-known
result for the DDM on both datasets. Moreover, the simple DE algorithm takes
only a fraction of the runtime required by other contemporary metaheuristics
and is thus preferable in real-time scenarios. This novel, unusual, and
notable finding also indicates that the employment of increasingly complicated
metaheuristics might be somewhat overkill for regular PV parameter estimation.
Finally, we discuss the implications of these results and suggest promising
directions for future development.

###### keywords:

Photovoltaic modeling, Parameter identification , Metaheuristic algorithms,
Global optimization , Differential evolution

††journal: Journal

##  1 Introduction

Solar energy is one of the most important renewable energy resources
characterized by its emission-free cleanness, ubiquitous abundance, high
sustainability, favorable power density, and almost cost-free exploitation
except the initial investment (Navabi et al., 2015; Jordehi, 2016a; Guo et
al., 2016). The major application of solar energy is photovoltaic (PV) power
generation through a solar PV system composed of many PV cells connected in
series or parallel (Li et al., 2019; Li and Xu, 2018). A PV cell is a
fundamental unit that converts sunlight to electricity, which is basically a
semiconductor diode with a P-N junction exposed to the light (Villalva et al.,
2009). Despite the potential advantages of solar energy, current PV systems
still suffer from low efficiency and large generation fluctuations, whose
power primarily depends on the varying environmental factors like irradiation
and temperature (Navabi et al., 2015; Jordehi, 2016a). Accurate modeling of PV
cells is thus necessary for the design, simulation, power forecasting, and
optimal control of solar PV systems (Jordehi, 2016a; Navabi et al., 2015; Gong
and Cai, 2013; Li et al., 2019; Guo et al., 2016; Yang et al., 2020).

The dominating method to describe solar PV systems depends on an analogous
electrical circuit model (Villalva et al., 2009), which, based on the number
of diodes, has been further specialized to the single-, double-, or even
three-diode models (Chenouard and El-Sehiemy, 2020; Li et al., 2019; Yang et
al., 2020; Jordehi, 2016a; Chin et al., 2015; Qais et al., 2019; Zaimi et al.,
2019). Despite the intuitiveness of these circuit models, the main difficulty
lies in the accurate determination of unknown parameters in a model. To
address this challenge, an increasing number of studies have been conducted in
recent years, yielding various algorithms to identify PV model parameters to
high precision (see, e.g., (Chin et al., 2015; Yang et al., 2020; Jordehi,
2016a) for detailed reviews). There are generally two kinds of data sources
for PV parameter estimation: datasheets provided by manufacturers and
experimentally measured I-V (current-voltage) curves (Chin et al., 2015;
Yousri et al., 2020; Ishaque and Salam, 2011). In contrast to the few
essential data points restricted to the standard test condition in the
datasheet (Almonacid et al., 2010), an experimental I-V curve generally
contains more data points measured in any practical operation condition.

Parameter estimation methods are broadly divided into two classes, analytical
methods and nonlinear optimization algorithms, which correspond roughly to the
two data sources mentioned above (Yang et al., 2020; Li et al., 2019; Yousri
et al., 2020). Analytical methods are distinguished by their simplicity and
rapid estimates through extracting information at the essential points to
formulate a set of equations to identify unknown parameters (Yousri et al.,
2020; Li et al., 2019; Navabi et al., 2015; Tong and Pora, 2016; Zaimi et al.,
2019). Nonetheless, to allow complete parameter estimation with limited
information extracted from merely essential points, some possibly restrictive
assumptions and simplifications have to be made, which may compromise the
accuracy and physical reliability of the results (Li et al., 2019; Yousri et
al., 2020). We refer readers to (Chin et al., 2015) and (Cotfas et al., 2013)
for a more detailed review of analytical methods.

The second class of methods usually formulates PV parameter estimation as a
nonlinear optimization problem from the perspective of I-V curve fitting
(Jordehi, 2016b), which has drawn a great deal of research interest in the
past decade. Due to the nonlinearity and nonconvexity present in the problem,
a large number of nature-inspired metaheuristics have thus been employed to
tackle this challenge in view of their global search ability (Yang, 2020). The
employment of metaheuristics for PV parameter estimation is also the focus of
our study in this paper.

Most metaheuristics are population-based by exploiting a swarm of interacting
agents to search the solution space efficiently, e.g., the classic genetic
algorithm (Dizqah et al., 2014) and particle swarm optimization (Nunes et al.,
2018). Since such nature-inspired metaheuristic algorithms are generally not
problem-specific, any metaheuristic optimizer may be applied to parameter
estimation of PV cells in principle. As reported in a recent survey (Yang,
2020), there are more than 100 different nature-inspired algorithms and
variants. It is thus unsurprising that numerous metaheuristic methods have
been employed for PV parameter estimation. Some recent examples include cat
swarm optimization (Guo et al., 2016), artificial been swarm optimization
(Askarzadeh and Rezazadeh, 2013), moth-flame optimizer (Sheng et al., 2019),
improved JAYA optimization (Yu et al., 2017b, 2019), sunflower optimization
(Qais et al., 2019), marine predators algorithm (Yousri et al., 2020), self-
adaptive ensemble-based differential evolution (Liang et al., 2020a), multiple
learning backtracking search (Yu et al., 2018), teaching-learning-based
optimization (Li et al., 2019; Chen et al., 2018), chaotic whale optimization
(Oliva et al., 2017), repaired adaptive differential evolution (Gong and Cai,
2013), and variants of these algorithms, among many others. Interested readers
may refer to the latest survey (Yang et al., 2020) for a comprehensive review
of metaheuristics in PV parameter identification. Finally, we note that there
also exist studies that try to hybridize traditional numerical optimization
methods with metaheuristic algorithms in the hope of bringing benefits from
both worlds, for example, the combination of the Nelder-Mead simplex algorithm
with the artificial bee colony algorithm (Chen et al., 2016) and with the moth
flame optimization (Zhang et al., 2020).

Despite the great expansion in interest of such metaheuristics, it should be
noted that, to our best knowledge, none of them can guarantee the discovery of
global optima; even though almost all of them have designed specific
mechanisms to help get away from local minima (Yang, 2020) (e.g., the well-
known probabilistic acceptance of worse solutions during iterations in
simulated annealing (El-Naggar et al., 2012)). In the literature of PV
parameter estimation, the performance of a metaheuristic algorithm is
typically justified by examining it empirically on some benchmark I-V curve
datasets, like the widely used RTC France solar cell data and the Photowatt-
PWP201 PV module data originally introduced in (Easwarakhanthan et al., 1986).
However, it can be observed in recent studies that, on the two representative
benchmark datasets, the difference in performance measured by root mean square
error (RMSE) among many metaheuristic methods appears negligible. For example,
it is only at the level of 1E-6 that the majority of methods exhibit
perceptible difference in terms of RMSE using either the single diode or the
double diode model on the RTC France dataset (see (Li et al., 2019, Table 3),
(Gong and Cai, 2013, Table 9), (Guo et al., 2016, Table 3), and Table 8 in
this paper). Thus, one may wonder naturally whether a significant improvement
is still possible by introducing more advanced and typically more complex
metaheuristics. Also, since a variety of metaheuristics characterized by
different complexity and distinct working principles are all able to get the
same RMSE value, a natural query is how sophisticated a metaheuristic has to
be, e.g., with extra parameter adaptation (Liang et al., 2020a) and strategic
hybridization (Chen et al., 2016), in order to estimate PV parameters
effectively. To put it more formally, we have at least the following two open
problems at present:

  1. 1.

Have the existing methods arrived at the minimum possible RMSE in respect of
the two most widely used benchmark datasets, since a bottleneck (i.e., no
further reduction) seems to emerge for methods proposed in the last five
years?

  2. 2.

If the above global minimum has already been attained by a few sophisticated
metaheuristic methods in the literature, can a simple and computationally
cheap metaheuristic get equally good results?

We believe that answers to the above questions are valuable to both the
research community and industrial practitioners. In one aspect, if the best-
known minimum of the objective function that the current literature has
achieved can be shown to be the global minimum (in the sense of
$\epsilon$-optimality (Floudas, 2013)), future researchers should not expect
to get even better solutions on this benchmark by striving for novel and
convoluted metaheuristics. Instead, researchers may put more efforts on other
factors of PV parameter identification aside from accuracy, e.g., performance
robustness, usability, and time efficiency. In another aspect, a practitioner
in the solar industry may get unfortunately overwhelmed by the abundance of
metaheuristics in the current literature and cannot figure out quickly which
method is ideal, if any, for a practical task. Instead of a blind and tedious
examination of existing algorithms, it is more meaningful and of practical
importance to start with a simple yet promising algorithm without time-
consuming parameter tuning. Such a simple, efficient, yet effective algorithm
will be recommended in this study, which answers the second question above and
provides valuable insights for forthcoming studies.

In this study, we investigate the above two problems using the two most
broadly adopted benchmark datasets (Easwarakhanthan et al., 1986). The main
contributions of this paper are summarized as follows.

  * 1.

We consider two popular PV cell models, including the five-parameter single
diode model (SDM) and the seven-parameter double diode model (DDM). An
interval arithmetic based branch and bound method is used to get the certified
global minimum for the SDM and a useful upper bound of the global minimum for
the DDM in terms of RMSE. These values can serve as useful references for the
evaluation and further development of metaheuristics. The obtained results
reveal that a few current methods have already obtained the $\epsilon$-global
minimum for the SDM.

  * 2.

We show, for the first time, that a simple differential evolution (DE)
algorithm can adequately attain the global minimum for the SDM and achieve
equally good accuracy for the DDM compared with a variety of sophisticated
state-of-the-art metaheuristics. Besides, the DE algorithm is characterized by
high stability and incomparable time efficiency. Overall, the DE algorithm
turns out to be a favorable alternative for PV parameter estimation despite
its simplicity.

  * 3.

All obtained results are analyzed, compared, and validated with other
contemporary algorithms. Based on our findings, we recommend the simple DE to
solar industry practitioners as the first choice in practical applications,
especially time-sensitive ones. Besides, we provide useful suggestions for
future studies to refresh viewpoints on PV parameter estimation and to refrain
from possible over-engineering caused by overcomplicated metaheuristics.

The remainder of this paper is organized as follows. Common PV models are
introduced in Section 2. The parameter estimation is formulated as an
optimization problem in Section 3. After that, the two optimization methods we
use in this study are described in Section 4. We then present results obtained
on two benchmark datasets for both models and conduct a detailed comparison in
Section 5. Finally, in Section 6, we give concluding remarks on our findings
and suggest possible future work.

##  2 Mathematical modeling of PV systems

In this section, we introduce the two most widely used circuit models for PV
systems, namely the single diode model (SDM) and the double diode model (DDM).
The physical rationale behind both models is that an ideal PV cell acts like a
semiconductor diode whose P-N junction is exposed to the light to produce an
electrical current through the photovoltaic effect (Villalva et al., 2009).

###  2.1 Single diode model (SDM)

The SDM is arguably the most popular model in the literature, especially for
parameter estimation, due to their least number of parameters and yet provable
effectiveness. The electrical circuit corresponding to the SDM is shown in
Fig. 1. Specifically, the circuit contains a current source $I_{ph}$, which
refers to the photocurrent generated by the PV cell, a diode flowing current
$I_{d}$, and two resistors with resistance $R_{p}$ and $R_{s}$, respectively.
The output current $I$ is obviously

$I=I_{ph}-I_{d}-I_{p}.$ (1)

where $I_{p}$ is the current through the parallel resistance $R_{p}$ (also
known as shunt resistance (Jordehi, 2016a; Zhang et al., 2020)). Supposing the
output voltage is $V$, the current $I_{p}$ is obviously

$I_{p}=\frac{V+IR_{s}}{R_{p}}.$ (2)

Figure 1: Equivalent circuit of the single diode model.

Next, we can calculate the diode current $I_{d}$ using the Shockley equation
as follows,

$I_{d}=I_{0}\left[\exp\left(\frac{q(V+IR_{s})}{nkT}\right)-1\right],$ (3)

where $I_{0}$ is the reverse saturation current of the diode, $n$ is the diode
ideal factor, $T$ is the temperature in Kelvin, and $V$ is the output voltage
of the cell. The other terms are just physical constants: the electron charge
$q=1.60217646\times 10^{-19}$ C and the Boltzmann constant $k=1.3806503\times
10^{-23}$ J/K.

By inserting Eq. (2) and (3) into Eq. (1), we get the complete equation of the
SDM:

$I=I_{ph}-I_{0}\left[\exp\left(\frac{q(V+IR_{s})}{nkT}\right)-1\right]-\frac{V+IR_{s}}{R_{p}}.$
(4)

There are five unknown parameters in (4), which are collected into a vector
$\bm{\theta}_{S}=[I_{ph},I_{0},n,R_{s},R_{p}]$.

###  2.2 Double diode model (DDM)

Despite the simplicity and usefulness of the above SDM, it does not consider
the effect of recombination current loss in the depletion region (Jordehi,
2016a; Chin et al., 2015; Zhang et al., 2020). An additional diode can be
introduced into the circuit to compensate for this specific loss to attain
higher accuracy. The equivalent circuit of the DDM is illustrated in Fig. 2,
in which there are two parallel diodes with current $I_{d1}$ and $I_{d2}$
respectively. One diode plays the role of a rectifier, while the other
emulates the recombination effect (Zhang et al., 2020; Yousri et al., 2020).

Figure 2: Equivalent circuit of the double diode model.

The only difference between the SDM and the DDM from the perspective of
mathematical equations is the inclusion of an additional current term. In
analogy to the SDM (4), the DDM is derived straightforwardly as follows:

$\displaystyle I$ $\displaystyle=$ $\displaystyle I_{ph}-I_{d1}-I_{d2}-I_{p}$
$\displaystyle=$ $\displaystyle
I_{ph}-I_{01}\left[\exp\left(\frac{q(V+IR_{s})}{n_{1}kT}\right)-1\right]-$
$\displaystyle
I_{02}\left[\exp\left(\frac{q(V+IR_{s})}{n_{2}kT}\right)-1\right]-\frac{V+IR_{s}}{R_{p}},$
(5)

where $I_{01}$ and $I_{02}$ are the reverse saturation current of the two
diodes, and $n_{1}$ and $n_{2}$ denote the ideality factor of the two diodes,
respectively. In contrast to the SDM discussed above, the DDM has seven
parameters in total, denoted by
$\bm{\theta}_{D}=[I_{ph},I_{01},I_{02},n_{1},n_{2},R_{s},R_{p}]$.

We note that, though in theory the DDM may better capture fine-grained
electrical characteristics of PV cells especially at low irradiance
conditions, the inclusion of two additional parameters poses more difficulty
for parameter estimation. In particular, the DDM is seldom used in datasheet-
based analytical methods because the accurate determination of seven
parameters requires at least seven equations, but at most five equations can
be formed with the manufacturer’s datasheet in general (Chin et al., 2015;
Tong and Pora, 2016). To further refine the circuit-based model, another diode
may be inserted to form the three diode model (Yousri et al., 2020; Qais et
al., 2019), which has in total of nine parameters. However, there is no solid
evidence justifying the usage of such complex models: they not only challenge
the effective estimation of PV parameters but also increase the computation
time, whose application in PV system simulations is consequently quite limited
(Jordehi, 2016a; Villalva et al., 2009).

###  2.3 PV module model

A practical PV module usually contains multiple PV cells that are connected in
series or parallel. Supposing $N_{s}$ is the number of series PV cells and
$N_{p}$ is the number for parallel connection, the overall model of a PV
module based on the SDM of an individual cell is (Li et al., 2019; Zhang et
al., 2020; Yu et al., 2019)

$\displaystyle I=$ $\displaystyle
I_{ph}N_{p}-I_{0}N_{p}\left[\exp\left(\frac{q(V/N_{s}+IR_{s}/N_{p})}{nkT}\right)-1\right]-$
$\displaystyle\frac{V/N_{s}+IR_{s}/N_{p}}{R_{p}/N_{p}}.$ (6)

Note, however, that Eq. (2.3) holds only if all PV cells in the module are
identical, i.e., they share the same parameter values. This assumption is
unlikely to be true in reality due to the intrinsic heterogeneity among cells.
For this reason, it is standard in the literature to fit the SDM (4) to PV
module data directly by lumping all cells into a single, functionally
equivalent cell (Easwarakhanthan et al., 1986; Yu et al., 2019). Of course,
the corresponding parameter values of a PV module (2.3) are supposed to be
much larger than those of an actual cell. The extension of the DDM (2.2) to a
PV module is straightforwardly analogous to (2.3) and omitted here.

##  3 Problem formulation

We focus on the I-V curve fitting methods in this study, which typically
formulate parameter estimation into a nonlinear optimization problem. The core
principle is to find appropriate parameter values such that the current values
calculated with either the SDM (4) or the DDM (2.2) match the measurement
values closely for a set of current-voltage data points. To simplify
notations, we shorten the SDM (4) and the DDM (2.2) as the following two
equations:

$\displaystyle I$ $\displaystyle=f_{S}(V,I;\bm{\theta}_{S})$ (SDM) (7a)
$\displaystyle I$ $\displaystyle=f_{D}(V,I;\bm{\theta}_{D})$ (DDM) (7b)

Without loss of generality, we continue our discussion towards problem
formulation using the SDM $f_{S}$ in (7a) (with $\bm{\theta}_{S}$ omitted for
conciseness) in the following part, whose principle can be transplanted to the
DDM seamlessly.

From a traditional perspective of curve fitting, e.g., our previous work (Gao
et al., 2018), we may view the voltage as the input and the current as the
output in (7a). However, the peculiarity of PV parameter estimation is that
the current is related to the voltage only implicitly (Cotfas et al., 2013),
that is, we cannot write down a simple closed-form solution $I=f_{S}^{-1}(V)$
for the model $f_{S}$ to compute $I$ given $V$ (Nunes et al., 2018).
Consequently, given a voltage value $V$, the corresponding current value $I$
generally has to be computed iteratively using a numerical root-finding
algorithm like the Newton-Raphson method to solve the nonlinear equation
$I-f_{S}(V,I)=0$ (Yousri et al., 2020; Jordehi, 2016b; Nunes et al., 2018;
Kler et al., 2019). It deserves to be mentioned that in some recent studies an
analytical solution to $f_{S}^{-1}$ has been obtained using the Lambert W
function (Yousri et al., 2020; Ćalasan et al., 2020). However, this method
only applies to the 5-parameter SDM (Ćalasan et al., 2020), and the
computation of the Lambert W function itself is iterative (Abdulrazzaq et al.,
2020) (see equations (6)-(11) in (Ćalasan et al., 2020)). Another potential
drawback of such iterative methods is the increased computational burden,
especially in the case of curve fitting based parameter estimation, since
thousands of function evaluations usually need to be performed. Hence, given a
measurement $(V^{m},I^{m})$ and a parameter vector $\bm{\theta}_{S}$ (not
necessarily the ground truth), the majority of metaheuristic-based studies in
the literature compute the predicted current with the model approximately but
computationally economically as

$I\approx f_{S}(V^{m},I^{m};\bm{\theta}_{S}),$ (8)

and try to reduce the deviation between $I$ and $I^{m}$ by adjusting
$\bm{\theta}_{S}$ (see (Chen et al., 2016), (Li et al., 2019), (Nunes et al.,
2018), (Yu et al., 2019), and (Liang et al., 2020b) among others). It is
evident that, if we can find a parameter vector $\hat{\bm{\theta}}_{S}$ such
that the computed $I$ in (8) is sufficiently close to $I^{m}$ for each
measurement, then $\hat{\bm{\theta}}_{S}$ is a good estimation of the ground
truth ${\bm{\theta}}_{S}^{*}$, which in turn justifies the rationality of the
approximation in (8).

Two metrics are frequently used in the literature to quantify the difference
between a computed current value $I$ and the true measurement value $I^{m}$
(Yang et al., 2020; Jordehi, 2016a), namely, the mean absolute error (MAE)
(Chenouard and El-Sehiemy, 2020; Guo et al., 2016) and the root mean square
error (RMSE) (Liang et al., 2020a; Zhang et al., 2020; Yousri et al., 2020;
Gnetchejo et al., 2019). In line with the majority of the literature and for
purposes of comparison, we adopt the RMSE metric in this study as follows:

$\textnormal{RMSE}=\sqrt{\frac{1}{N}\sum_{i=1}^{N}(I_{i}-I^{m}_{i})^{2}},$ (9)

where $N$ is number of measured data points.

A well-known trick to slightly promote computational efficiency of (9) is to
avoid the square root calculation and use an equivalent metric called the sum
of squared error (SSE), whose optimal solution is exactly the same as the one
of RMSE. In summary, we get the following constrained optimization problem,
which is known widely as nonlinear least squares regression in the literature.

minimize $\displaystyle
J(\bm{\theta})=\sum_{i=1}^{N}\left(f(V^{m}_{i},I^{m}_{i};\bm{\theta})-I^{m}_{i}\right)^{2},$
(10a) subject to $\displaystyle\bm{\theta}\in\Theta.$ (10b)

In the above equation (10), $f$ refers to either the SDM $f_{S}$ (4) or the
DDM $f_{D}$ (2.2), and $\bm{\theta}$ is the corresponding parameter vector
$\bm{\theta}_{S}$ or $\bm{\theta}_{D}$. $(V^{m}_{i},I^{m}_{i})$ is the $i$-th
data point in measurement. $\Theta$ denotes the bound constraints of
$\bm{\theta}$ that take physical reality into consideration (see the parameter
ranges listed in Table 1 for examples).

##  4 Optimization methods

To address the open problems we listed in Section 1, we attempt to (i)
determine rigorously the global minimum or at least a reasonably tight bound
of the global minimum using numerical optimization algorithms and (ii)
demonstrate that an intentionally simple metaheuristic can achieve the global
minimum or at least attain the best minimum known in the literature. To our
best knowledge, no existing metaheuristic algorithms, belonging to stochastic
global optimization family, can certify the global minimum or quantify the
optimality gap for the PV parameter estimation problem (10) (Chenouard and El-
Sehiemy, 2020; Yang, 2020). To evaluate the optimality rigorously in task (i),
we apply a branch and bound (B&B) based deterministic global optimization
technique. For task (ii), we choose a simple and widely used evolutionary
algorithm, differential evolution (DE). We detail the mathematical principles
of both methods in the remainder of this section.

###  4.1 Deterministic global optimization with an interval arithmetic based
branch-and-bound algorithm

The fundamental task of deterministic global optimization (DGO) is to
determine rigorously (i.e., with theoretical guarantees) the global minimum of
an objective function $f$ subject to a set of constraints $C$ (Floudas and
Gounaris, 2008; Floudas, 2013). DGO is notoriously tricky when handling
nonconvex optimization problems: even to check the feasibility of a given
problem may be out of reach (Floudas and Gounaris, 2008; Burer and Letchford,
2012). In fact, a general nonconvex optimization problem like the nonlinear
least squares in (10) has been proved to be NP-hard (Floudas, 2013). From a
practical point of view, we are usually more interested in identifying a
solution sufficiently close to the true global minimum, called the
$\epsilon$-global minimum, which is defined formally below.

###### Definition 1

(Floudas, 2013) Consider the problem of minimizing $f(\bm{x}),\bm{x}\in S$. A
feasible solution $\tilde{\bm{x}}\in S$ is an $\epsilon$-global minimum if
$f(\bm{x})\geq f(\tilde{\bm{x}})-\epsilon,\forall\bm{x}\in S$, where
$\epsilon\geq 0$ is a small prescribed tolerance.

The most popular algorithmic framework of DGO to handle such general NLP
problems is arguably the branch-and-bound (B&B) method and its variants like
branch-and-reduce (Burer and Letchford, 2012; Floudas and Gounaris, 2008;
Floudas, 2013). The B&B framework for global optimization is usually
attributed to McCormick (McCormick, 1976), whose general principle is
intuitive. The search space is divided recursively into smaller subspaces and
forms accordingly a tree structure of subproblems. One crucial factor to
success lies in the determination of proper bounds of each subproblem, usually
via convex or concave relaxation. The consequential pruning of search space is
performed by eliminating subproblems whose lower bound is no better than the
best upper bound found so far (Burer and Letchford, 2012; McCormick, 1976). In
the worst case, if only a fraction of subproblems can be pruned, the
performance of B&B is close to exhaustive search. Fortunately, interval
analysis is often a useful tool for efficient estimation of the lower and
upper bounds of regions/branches of the search space.

The definitive characteristic of interval analysis (also known as interval
arithmetic) is to represent a single number $x$ as an interval $[a,b]$ that
contains $x$ in order to track numerical errors such as rounding errors and
measurement errors. Apart from their numerical reliability, interval methods
have important applications in rigorous global optimization and location of
all solutions to nonlinear equations (Hansen and Walster, 2003). Common
arithmetic operators and mathematical functions have been extended to
intervals. For example, given any binary operator $\circ$ among
$(+,-,\times,\div)$, we have

$[x_{1},x_{2}]\circ[x_{2},y_{2}]=\\{x\circ
y|x\in[x_{1},x_{2}],y\in[y_{1},y_{2}]\\}.$

As mentioned above, one critical operation in B&B for global optimization is
to reduce efficiently the hierarchical subproblems by pruning the search tree;
otherwise, B&B degenerates to brute-force search that is intractable for most
nontrivial problems. Recall that pruning efficiency relies significantly on
the tightness of estimated bounds (especially the lower bound) of the
objective function for a specific subproblem. Since applying a function to an
interval yields another interval, interval analysis turns out to be valuable
in this particular aspect due to its capacity to compute rigorous yet sharp
bounds for a function $f$ over an interval vector $[\bm{x}]$ (also known as an
interval box, whose components are intervals) using techniques like hull
consistency, Taylor expansion, and interval Newton method (Hansen and Walster,
2003, Chapter 7). Note that our problem (10) has only bound constraints, which
is viewed as an unconstrained optimization problem with an initial interval
box defined by the bound constraints in interval analysis. Once the function
bounds over a subregion (i.e., a subbox) are available, we can detect and
discard an unpromising subregion, whose lower bound is greater than a known
minimum value at hand, from the subsequent search for a global minimum.

The general B&B framework with interval arithmetic is depicted in Algorithm 1,
where $[\bm{x}]$ denotes an interval box and
$[\underline{f}([\bm{x}]),\bar{f}([\bm{x}])]$ denotes the interval value of
applying function $f$ to $[\bm{x}]$. Each iteration is composed of three main
components: box selection (Line 7), box contracting (Line 8), and box
splitting (Line 14). Starting from the initial box $[\bm{x}]^{0}$
corresponding to the bound constraints $X$, the list $L$ tends to accumulate
many boxes over the time. There is no standard way in Line 7 to choose the box
$[\bm{x}]$ from $L$ to be processed next, but a simple heuristic is to pick
the box that has the smallest lower bound $\underline{f}([\bm{x}])$, which is
more likely to contain the global minimizer $\bm{x}^{*}$. The next key step is
to contract (i.e., prune) the chosen box $[\bm{x}]$. The purpose of
contracting is to delete subboxes inside $[\bm{x}]$ that cannot contain a
globally optimal solution. There are multiple techniques to perform
contraction, e.g., to delete subboxes wherein the gradient of $f$ is nonzero
or wherein $f$ is not convex. The rationale of the latter is that $f$ must be
convex in some neighborhood of a global minimum point $\bm{x}^{*}$ (Hansen and
Walster, 2003, Chapter 12). Another important technique to prune $[\bm{x}]$ is
to remove a subbox $[\bm{x}^{\prime}]\subseteq[\bm{x}]$ with
$\underline{f}([\bm{x}^{\prime}])\geq\bar{f}$, where $\bar{f}$ is the minimum
value of $f$ we have evaluated so far, which serves as an upper bound of the
global minimum value $f^{*}$. In general, such a box pruning operation is
conceptualized as a contractor, which is basically an operator that transforms
a box to a smaller (included) box by safely removing points without affecting
the global minimum. Similarly to contracting, there is no single best strategy
for box splitting in Line 14. Typical choices include splitting along the
widest edge of $[\bm{x}]$ or over a coordinate in which $f$ varies mostly.
After multiple iterations of box contracting and splitting, a small box may be
added to the candidate solution list $L_{S}$ as a tentative solution, i.e.,
one that may contain the global minimum $\bm{x}^{*}$. In order to be included
in $L_{S}$, a box $[\bm{x}]$ must satisfy two conditions that are checked in
Line 11:

$w([\bm{x}])\leq\epsilon_{x},\ w(f([\bm{x}]))\leq\epsilon_{f},$ (11)

where $w(\cdot)$ denotes the width of a box defined by its largest diameter
among all components. $\epsilon_{x}$ and $\epsilon_{f}$ are two tolerance
parameters provided by the user, often known as the precision. Obviously,
smaller values of $\epsilon_{x}$ and $\epsilon_{f}$ demand more iterations and
prolongs the execution time. After the main loop finishes, we post-process the
solution list $L_{S}$ in Line 17 to discard boxes which cannot contain the
global minimum $\bm{x}^{*}$ according to the latest knowledge of $\bar{f}$.

At the end of Algorithm 1, we get the list of boxes $L_{S}$ containing the
global minimum $\bm{x}^{*}$ and the (usually very tight) bounds of the global
minimum value $f^{*}\in[\underline{f},\bar{f}]$. Note that
$f([\bm{x}])-f^{*}\leq 2\epsilon_{f},\forall[\bm{x}]\in L_{S}$, is guaranteed
(Hansen and Walster, 2003). Recall Definition 1, and we see that any $\bm{x}$
inside the remaining boxes in $L_{S}$ becomes an $\epsilon$-global minimum
with $\epsilon=2\epsilon_{f}$ in this case. Though the exact global minimum
$\bm{x}^{*}$ and $f^{*}$ are still unknown in general, a reasonably tight
bound acquired by setting small $\epsilon_{f}$ and $\epsilon_{x}$ in Algorithm
1 is usually enough in practice.

Algorithm 1 Interval arithmetic based branch-and-bound algorithm for
deterministic global optimization

1:objective function $f:\mathbb{R}^{n}\rightarrow\mathbb{R}$ and bound
constraints $X\subset\mathbb{R}^{n}$, precision parameters $\epsilon_{f}$ and
$\epsilon_{x}$

2:bounds of the global minimum $[\underline{f},\bar{f}]$ and a list of boxes
$L_{S}$ that contain all global minima

3:initialize a list $L\leftarrow\\{[\bm{x}]\\}$ with $[\bm{x}]$ corresponding
to $X$

4:initialize an empty candidate solution list $L_{S}$

5:$\bar{f}\leftarrow\infty$ $\triangleright$ Upper bound of $f^{*}$

6:while $L\neq\emptyset$ do

7: choose $[\bm{x}]\in L$ and remove $[\bm{x}]$ from $L$

8: contract $[\bm{x}]$

9: evaluate $f$ at the center of $[\bm{x}]$ and get value $f_{c}$

10: Update $\bar{f}$ by $\bar{f}\leftarrow\min\\{\bar{f},f_{c}\\}$

11: if $[\bm{x}]$ satisfies criteria (11) then

12: append $[\bm{x}]$ to $L_{S}$

13: else

14: split $[\bm{x}]$ into subboxes and add them to $L$

15: end if

16:end while

17:remove any box $[\bm{x}]\in L_{S}$ from $L_{S}$ with
$\underline{f}([\bm{x}])>\bar{f}$

18:$\underline{f}\leftarrow\min_{[\bm{x}]\in L_{S}}{\underline{f}([\bm{x}])}$

Note that Algorithm 1 only sketches out the basic skeleton of interval B&B
algorithms. There are many practical issues and improvements to be considered
for high-quality implementation. We thus resort to a well-developed open-
source numerical library called ibexopt666http://www.ibex-lib.org (v2.8),
which is based on interval analysis and contractor programming, to perform
interval B&B optimization in PV parameter estimation. We refer interested
readers to the monograph (Hansen and Walster, 2003) and the documentation of
ibxopt for more details.

###### Remark 1

The global optimization solver ibexopt has been employed in a recent study
(Chenouard and El-Sehiemy, 2020). Nonetheless, the results obtained therein
have insufficient precision due to the timeout setting of 20000 s and,
consequently, cannot yield sufficiently tight optimality bounds. In this
study, we use the interval B&B algorithm (implemented by ibexopt) only to
validate global optimality but do not recommend it as a routine optimizer for
PV parameter identification due to its excessively long running time. Thus, we
are free to allow it more time and set smaller tolerances for promoted
precision. We will give more details in Section 5.2.

###  4.2 Stochastic global optimization with a simple differential evolution
algorithm

As we have discussed above, though an interval B&B algorithm has the capacity
to ascertain the global optimum rigorously in theory, it is generally much
more computationally expensive and runs considerably slower than stochastic
metaheuristic algorithms, which can avoid local minima to a certain extent but
with no guarantee (Yang, 2020; Yang et al., 2020). The ability to find
satisfactorily high-quality solutions in reasonable computation time motivates
the broad application of stochastic search based metaheuristics in PV
parameter estimation (Chen et al., 2016; Yang et al., 2020; Nunes et al.,
2018). In this section, we take a remarkably simple evolutionary algorithm on
purpose, namely the differential evolution (DE), to approach the parameter
estimation problem (10) in order to address the second concern listed in
Section 1.

DE is a simple, efficient, yet powerful evolutionary algorithm (EA) originally
proposed by Price and Storn (Price et al., 2006), which has got proved success
in a variety of optimization problems arising from diverse domains of science
and engineering (Das and Suganthan, 2011). The exceptional popularity of DE is
partly attributed to its simplicity, e.g., the main body of the algorithm
takes only four to five lines of code (see Algorithm 2), as well as the small
number of control parameters, only four in classical DE. Similar to other EAs,
DE follows the Darwinian principle to evolve a population of solutions (called
vectors in the DE community), and each iteration of DE comprises three key
steps: selection, crossover, and mutation. The distinguishing feature of DE is
its mutation with difference vectors (Price et al., 2006; Das and Suganthan,
2011). In addition, by making a parent vector compete with a child vector for
survival into the next generation, DE’s selection operator enforces elitism
and ensures that the highest quality vector ever found will never get lost
during long evolution.

Suppose that we want to minimize a scalar function
$f:\mathbb{R}^{n}\rightarrow\mathbb{R}$ with bound constraints
$X\subset\mathbb{R}^{n}$. A candidate solution is represented by a vector
$\bm{x}\in\mathbb{R}^{n}$ in DE. The procedure of classical (and also the most
simple) DE is outlined in Algorithm 2. The initial population $P$ comprises
$N_{p}$ vectors, and each vector $\bm{x}_{i}^{0},i\in[1,N_{p}]$ is generated
randomly as follows,

$x_{i,j}^{0}=\underline{b}_{j}+\text{rand}(0,1)\cdot(\bar{b}_{j}-\underline{b}_{j}),$
(12)

where $x_{i,j}^{0}$ denotes the $j$-th component of $\bm{x}_{i}^{0}$,
$j\in[1,n]$, $\underline{b}_{j}$ and $\bar{b}_{j}$ represent the lower and
upper bound of the $j$-th variable respectively, and $\text{rand}(0,1)$
generates a random number between 0 and 1, which is called anew independently
for each pair of $i$ and $j$.

Algorithm 2 Differential evolution with bound constraints handled by bounce-
back

1:objective function $f:\mathbb{R}^{n}\rightarrow\mathbb{R}$ and bound
constraints $X\subset\mathbb{R}^{n}$, control parameters $N_{p},C_{r},F,G$

2:the best vector discovered $\hat{x}$ and the corresponding function value
$\hat{f}$

3:generate randomly an initial population
$P^{0}\leftarrow\\{\bm{x}_{i}^{0}\\},i=1,2,\cdots,N_{p}$ with (12)

4:for $g$ from $0$ to $G-1$ do

5: for each vector $\bm{x}_{i}^{g}\in P^{g}$ do

6: generate a donor vector $\bm{v}_{i}^{g}$ by (13) $\triangleright$ mutation

7: $\bm{v}_{i}^{g}\leftarrow\text{bounce-back}(\bm{v}_{i}^{g})$ by (14)

8: $\bm{u}_{i}^{g}\leftarrow\text{crossover}(\bm{v}_{i}^{g},\bm{x}_{i}^{g})$
by (15)

9: $\bm{x}_{i}^{g+1}\leftarrow\text{select}(\bm{u}_{i}^{g},\bm{x}_{i}^{g})$ by
(16)

10: insert $\bm{x}_{i}^{g+1}$ into the new population $P^{g+1}$

11: end for

12:end for

13:$\hat{x}\leftarrow$ the best vector in $P^{G}$

14:$\hat{f}\leftarrow f(\hat{x})$

The mutation operation in DE is characterized by its usage of difference
vectors. Several mutation strategies have been developed in the literature
(Das and Suganthan, 2011; Price et al., 2006). Here we adopt the most commonly
used one called the “DE/rand/1” scheme. For each vector in the $g$-th
iteration, an auxiliary vector $\bm{v}_{i}^{g},i\in[1,N_{p}],$ is generated by

$\bm{v}_{i}^{g}=\bm{x}_{a}^{g}+F(\bm{x}_{b}^{g}-\bm{x}_{c}^{g}),\ a\neq b\neq
c\neq i,$ (13)

where $a,b,c\in[1,N_{p}]$ are randomly chosen and $F$ is a control parameter,
called the scaling factor, which typically lies in the range $[0.4,1]$ (Das
and Suganthan, 2011). The resultant $\bm{v}_{i}^{g}$ is known as the donor
vector in contrast to the target vector $\bm{x}_{i}^{g}$.

Note, however, that the donor vector $\bm{v}_{i}^{g}$ produced in (13) may lie
outside the feasible region if the problem has additional constraints, like
bound constraints $X$ (see Table 1 for PV parameter estimation). We adapt a
simple strategy called bounce-back from (Das and Suganthan, 2011; Price et
al., 2006) to handle bound constraints in Line 7, which relocates each
infeasible component between the bound it exceeds and the corresponding value
of the target vector, as follows.

$v_{i,j}^{g}=\begin{cases}\underline{b}_{j}+\text{rand}(0,1)\cdot(x_{i,j}^{g}-\underline{b}_{j})\
\text{if }v_{i,j}^{g}<\underline{b}_{j}\\\
\bar{b}_{j}-\text{rand}(0,1)\cdot(\bar{b}_{j}-x_{i,j}^{g})\ \text{if
}v_{i,j}^{g}>\bar{b}_{j}\end{cases}$ (14)

As is common in EAs, a crossover operation is applied in DE to enhance the
diversity of the population. The fundamental idea is to exchange components
between the donor vector $\bm{v}_{i}^{g}$ and the target vector
$\bm{x}_{i}^{g}$ so as to produce a new vector $\bm{u}_{i}^{g}$ named the
trial vector. An intuitive and widely used crossover scheme is the binomial
crossover, which is performed as follows:

$u_{i,j}^{g}=\begin{cases}v_{i,j}^{g}\ \text{if }\text{rand}(0,1)\leq
C_{r}\text{ or }j=\beta\\\ x_{i,j}^{g}\ \text{otherwise}\end{cases}$ (15)

where $\beta\in[1,n]$ is a random integer that is generated anew for each $i$.
Besides, $C_{r}$ is the user provided crossover rate. The role of $\beta$ is
to ensure that at least one component of the trial vector $\bm{u}_{i}^{g}$
comes from $\bm{v}_{i}^{g}$ (Das and Suganthan, 2011).

Finally, differing from other EAs like genetic algorithm, DE imposes elitism
by selecting the better one between the target vector $\bm{x}_{i}^{g}$ and the
trial vector $\bm{u}_{i}^{g}$ as the $i$-th vector in the next generation.
That is,

$\bm{x}_{i}^{g+1}=\begin{cases}\bm{u}_{i}^{g}\ \text{if }f(\bm{u}_{i}^{g})\leq
f(\bm{x}_{i}^{g})\\\ \bm{x}_{i}^{g}\ \text{otherwise}\end{cases}.$ (16)

One interesting observation is that $\bm{u}_{i}^{g}$ replaces $\bm{x}_{i}^{g}$
even if they share the same objective value, a useful feature that allows DE
to move over flat fitness landscapes for efficient exploration of the search
space (Das and Suganthan, 2011). After the selection with (16) in Line 9, the
population size remains $N_{p}$. Note that the termination criterion of DE can
be defined in various forms, e.g., when the minimum function value obtained so
far has not changed for a few consecutive generations (early stopping), apart
from specifying a fixed number of generations, $G$, in Algorithm 2.

###### Remark 2

Many variants of DE have been proposed (see the latest survey (Bilal et al.,
2020) for an extensive review), which include advanced mutation and/or
crossover operators and control parameter adaptation. For example, the
repaired adaptive DE (Gong and Cai, 2013) and the self-adaptive ensemble-based
differential evolution (SEDE) (Liang et al., 2020a) has been used to extract
solar cell parameters. We refer readers to Section 4.1.2 of (Yang et al.,
2020) for a review of advanced DE in PV parameter identification. Algorithm 2
represents the simplest form of classic DE, which is chosen intentionally here
to validate our hypothesis: the increasingly complicated metaheuristics
pervading the literature may not be necessary for pragmatic PV parameter
estimation.

##  5 Results and discussion

In this section, we apply the optimization methods introduced in Section 4 to
two widely used benchmark datasets and try to answer the two questions raised
in Section 1 by analyzing and comparing the results.

###  5.1 Datasets and experimental settings

The experimental data are acquired from (Easwarakhanthan et al., 1986), which
serve as the de facto benchmark in the PV parameter estimation literature
(Yousri et al., 2020; Yang et al., 2020). The first dataset contains 26 data
points measured for an RTC France solar PV cell of 57 mm diameter that
operates at 33 ∘C and under irradiance 1000 W/m2. The second dataset is
measured from a solar module called Photowatt-PWP201 that contains 36
polycrystalline silicon cells in series and operates at 45 ∘C temperature with
an irradiance of 1000 W/m2. Though the original dataset for Photowatt-PWP201
in (Easwarakhanthan et al., 1986, Table II) includes 26 data points, most
studies discard the first data point, e.g., (Gong and Cai, 2013; Li et al.,
2019; Chen et al., 2016; Liang et al., 2020b; Jordehi, 2016b). We follow this
convention in our study for fair comparison and compatible analysis. Note
that, the methodology developed in this study can be applied seamlessly to any
other dataset acquired under other conditions in practical applications.

Regarding the bound constraints in (10b), we adopt the parameter search range
used in the majority of literature for fair comparison purposes (see, e.g.,
(Chenouard and El-Sehiemy, 2020; Chen et al., 2016; Nunes et al., 2018; Li et
al., 2019; Liang et al., 2020b), among others). The parameter search ranges
for both datasets are listed in Table 1. Since the reverse saturation current
of a diode is extremely small, the unit $\mu$A is adopted to permit similar
scales among all parameters for numerical stability.

Algorithm 1 is implemented in the open-source software ibexopt (v2.8). We code
algorithm 2 with the Julia programming language for its flexibility and
efficiency. The results presented below were obtained on a desktop PC with a
3.4 GHz Core i7-3770 CPU, 16 GB RAM, and Windows 10.

Table 1: Parameter search range of SDM and DDM (LB: lower bound; UB: upper bound). Parameter | RTC France | Photowatt-PWP201
---|---|---
LB | UB | LB | UB
$I_{ph}$ (A) | 0 | 1 | 0 | 2
$I_{0},I_{01},I_{02}$ ($\mu$A) | 0 | 1 | 0 | 50
$n,n_{1},n_{2}$ | 1 | 2 | 1 | 50
$R_{s}$ ($\Omega$) | 0 | 0.5 | 0 | 2
$R_{p}$ ($\Omega$) | 0 | 100 | 0 | 2000

###  5.2 Global optimality analysis via interval B&B

In this part, we investigate the global minimum of the nonlinear least squares
problem (10) for both the SDM and the DDM using the interval B&B method
described in Section 4.1. As mentioned above, lots of metaheuristics have been
employed in the literature to extract parameter values for solar PV models.
For example, the results of 69 different methods, including the SDM’s five
parameters and the resultant RMSE values, for the RTC France solar cell
dataset are collected in Table 1 of a latest article (Ćalasan et al., 2020).
It is reasonable that the primary objective of most studies is to reduce
further the fitting error evaluated conventionally by the RMSE. Nevertheless,
none of these 69 methods can certify the the global minimum or quantify the
optimality gap even for the simple SDM case. From a practical standpoint and
to be compatible with the vast majority of studies, we restrict our discussion
of RMSE values to five significant digits. A notable observation of existing
methods, e.g., those summarized in (Ćalasan et al., 2020, Table 1), (Li et
al., 2019, Table 3), (Gnetchejo et al., 2019, Table 11), (Yousri et al., 2019,
Table 1), and (Yu et al., 2019, Table 4), is that the minimum RMSE attained so
far remains 9.8602E-4 when the SDM is utilized on the RTC France cell dataset.
In line with the prior observation, a natural speculation is whether 9.8602E-4
is the global minimum RMSE value in the SDM case. Similar observations exist
regarding the DDM and the Photowatt-PWP201 dataset. Such a conjecture
motivates our work present in this section.

Before proceeding to global optimization, we note that several studies have
reported RMSE values even smaller than 9.8602E-4 when fitting the SDM to the
RTC France solar cell, e.g., some of the results collected in (Gnetchejo et
al., 2019, Table 2). There are primarily two reasons for these unexpectedly
small values, either due to accidental errors in calculation as pointed out by
(Gnetchejo et al., 2019) and (Ćalasan et al., 2020) or because of an
evaluation metric different from the one defined by (8) and (9) (Yousri et
al., 2019). Regarding the latter specifically, instead of Eq. (8), some
studies decide to calculate the estimated current more precisely by solving a
nonlinear equation despite the increased computational cost (see Section 3),
which often contributes to smaller RMSE values (Jordehi, 2016b, Table 1).
Nevertheless, such precise but expensive computation is adopted in only a
small portion of existing studies, which is thus not considered here. Another
subtle yet influential reason for such discrepancy is possibly that the RMSE
value is quite sensitive to roundoff errors (Jordehi, 2016b), especially to
those appearing in the exponential term of (4) and (2.2). To be as precise as
possible and consistent with the mainstream studies, the values of the two
constants $q$ and $k$ refer to exactly the numbers listed below Eq. (3) (see,
e.g., (Chen et al., 2016), (Yu et al., 2019), (Long et al., 2020), and (Liang
et al., 2020b)). In our discussion below, we have excluded those erroneous and
incompatible results. Finally, note that the interval B&B method in Algorithm
1 is a deterministic one with no control parameters to be tuned (except the
timeout and desired precision), which means multiple runs with the same inputs
will always produce the same result (Chenouard and El-Sehiemy, 2020).

####  5.2.1 SDM results

We first applied the interval B&B algorithm to solve problem (10) specialized
for the SDM and the RTC France PV cell data. The bound constraints, i.e., the
initial box in Algorithm 1, are specified by the search range in Table 1.
Following (Chenouard and El-Sehiemy, 2020), we limited the maximum runtime of
ibexopt to 20000 seconds. On the other hand, unlike (Chenouard and El-Sehiemy,
2020), we did not use the default absolute and relative precision but set them
to smaller values in order to obtain tighter bounds of the global minimum
(that is, $[\underline{f},\bar{f}]$ in Algorithm 1). In our setting, the
absolute precision was 1E-13, and the relative precision was 1E-9 (in contrast
to the default values 1E-7 and 1E-3, respectively; see also Remark 3). Note
that our objective value in optimization denotes the SSE (10), and the
corresponding RMSE is obtained afterwards by (9) once the optimization
finishes. The results obtained with the SDM and the RTC France dataset are
reported in the first column of Table 2. Recall that, unlike metaheuristics,
the interval B&B algorithm yields a tight interval that contains the global
minimum, while the displayed values of the five parameters correspond to the
upper bound of the produced interval. On the other hand, despite its
theoretical preciseness, the main drawback of the interval B&B algorithm is
the exceptionally long execution time, that is, around 4 hours versus a few
seconds of most metaheuristics.

Table 2: Optimization results for SDM using interval B&B. The bounds of the global minimum RMSE value are reported, whose upper bound is attained by the parameter values above it. Variable | RTC France | Photowatt-PWP201
---|---|---
$I_{ph}$ (A) | 0.760779120136 | 1.03052020484
$I_{0}$ ($\mu$A) | 0.322873926858 | 3.48287904343
$n$ | 1.48113747635 | 48.6435574734
$R_{s}$ ($\Omega$) | 0.0363792207867 | 1.20123680201
$R_{p}$ ($\Omega$) | 53.7009537057 | 981.263690780
RMSE |  | [9.860250397955652E-4,
---
9.860250417458982E-4]
| [2.425076598320144E-3,
---
2.425076599532477E-3]
Gap | 1.950333050615427E-12 | 1.2123329007351913E-12
Time (s) | 13547 | 38924

We notice from the “RTC France” column of Table 2 that the optimality gap
between the lower and the upper bound of the RMSE is extremely small. By
focusing on five significant digits only, we can certify safely that 9.8602E-4
is indeed the global minimum value of RMSE by fitting the SDM to the RTC
France dataset, which is in agreement with our observation of the literature.
The five parameters’ values are also very close to those acquired with various
metaheuristics in the literature (see, e.g., (Ćalasan et al., 2020, Table 1),
(Yu et al., 2019, Table 4), (Long et al., 2020, Table 5), and Table 6 in the
next section), though these metaheuristic algorithms cannot recognize the
global minimum even if they have actually found it. Such an agreement also
confirms the correctness of our implementation.

Next, we applied the same methodology to the Photowatt-PWP201 PV module
dataset. The major difference lies in the considerably widened parameter
search range in Table 1, which may pose a big challenge to the interval B&B
algorithm. The results of fitting the SDM to the Photowatt-PWP201 dataset are
reported also in Table 2. As expected, it took ibexopt much more time to solve
the minimization problem (10) in the Photowatt-PWP201 case probably due to the
enlarged search space. As far as we know, the best minimum value of the SDM’s
RMSE acquired in the literature for this PV module is 2.4250E-3 (truncated to
five significant digits), as summarized in (Gnetchejo et al., 2019, Table 15),
(Long et al., 2020, Table 13), (Li et al., 2019, Table 6), (Chen et al., 2016,
Table 8), (Yu et al., 2019, Table 3), and (Liang et al., 2020b, Table 7),
which shows exact agreement with the global minimum value’s bounds reported in
Table 2. Furthermore, the metaheuristics that have succeeded in acquiring the
global minimum all yield parameter values extremely close to those in Table 2,
such as results presented in (Gong and Cai, 2013), (Chen et al., 2016, Table
8), (Yu et al., 2019, Table 6), (Liang et al., 2020a, Table 7), and (Long et
al., 2020, Table 7). Again, such a consensus indicates the correctness of each
other.

###### Remark 3

The ibexopt toolkit has also been utilized in a recent study (Chenouard and
El-Sehiemy, 2020). In contrast to our setting, the default precision
tolerances are used in (Chenouard and El-Sehiemy, 2020), which leads to a
fairly larger optimality gap roughly equal to 1.95E-6. Consequently, the
result of (Chenouard and El-Sehiemy, 2020) appears insufficient to certify the
global minimum value 9.8602E-4 in the RTC France SDM case. More interestingly,
it is surprising that the extra running time required for the higher precision
in our setting was essentially negligible, less than one second. Besides, the
prior study (Chenouard and El-Sehiemy, 2020) does not report results for the
Photowatt-PWP201 PV module.

####  5.2.2 DDM results

Following the success of interval B&B in the above SDM case, we proceed to the
DDM case, which is more challenging due to the additional two parameters (see
(2.2)). In the previous study (Chenouard and El-Sehiemy, 2020), ibexopt failed
to generate a tight optimality bound subject to a timeout of 20000 s. More
specifically, when optimizing the seven parameters of the DDM for the RTC
France dataset, the lower bound of the global minimum obtained by ibexopt in
20000 s was 0. However, the zero lower bound is trivial and useless since we
are minimizing a summation of squared terms in (10a). From this particular
standpoint, we decided to allow ibexopt more running time attempting to get
tighter optimality bounds. The overall workflow is identical to the one
presented in the above SDM scenario except that the SDM is replaced with the
DDM.

We first examined the DDM for the RTC France dataset with the search space
given by Table 1. The timeout was set to 86400 s, that is, 24 hours. The
obtained results are listed in the left side of Table 3. We see that,
unfortunately, the significant increase of the allowed running time did not
help much. In particular, the lower bound still remains zero. In a more
detailed investigation, we found that the upper bound was reduced to 9.8657E-4
at the end of the first hour. In other words, ibexopt made little progress in
the subsequent 23 hours, only decreasing the upper bound from 9.8657E-4 to
9.8358E-4, while the lower bound stayed at zero. It thus seems hopeless that
the optimality bound can be tightened substantially by allowing even longer
execution time777We have actually run ibexopt for another 5 hours but observed
no perceivable change in both lower and upper bounds.. The considerable
slowdown of the solver at a later stage of the search process is probably
caused by the so-called cluster effect, a pain that B&B methods frequently
suffer from (Montanher et al., 2018; Floudas, 2013). This effect is roughly
characterized by the excessive splitting of the search space in the
neighborhood of the global optimum. The resultant boxes accumulate as a
cluster around the globally optimal solution, and, what is worse, the pattern
repeats itself recursively at increasingly smaller scales if we want high
accuracy. Upon closer inspection, we notice that the number of boxes that have
been explored at the end of 24 hours is already more than 25 million. In
addition, the multiple occurrences of each parameter incidental to the SSE
objective function (10a), known as the dependence problem, can also compromise
the efficiency of interval calculation and lead to over-estimation of the
value range (Hansen and Walster, 2003; Chenouard and El-Sehiemy, 2020).

Table 3: Optimization results for DDM using interval B&B. The bounds of the global minimum RMSE value are reported, whose upper bound is attained by the parameter values above it. Variable | RTC France | Photowatt-PWP201
---|---|---
$I_{ph}$(A) | 0.760815738919 | 1.0339286971
$I_{01}$($\mu$A) | 0.217867184041 | 1.86575472010E-23
$I_{02}$($\mu$A) | 0.781454995330 | 0.535399234849
$n_{1}$ | 1.44827388213 | 9.58860778809
$n_{2}$ | 1.98183166760 | 42.6724488388
$R_{s}$($\Omega$) | 0.0367359827333 | 1.63619822583
$R_{p}$($\Omega$) | 55.8931982861 | 607.690281231
RMSE |  | [0,
---
9.83581875679E-4]
| [0,
---
1.61865668151E-3]
Gap | 9.83581875679E-4 | 1.61865668151E-3
Time (s) | 86400 | 43200

Though we did not get an extraordinarily tight optimality bound like that in
Table 2, which is mainly attributed to the zero lower bound, the upper bound
is still informative: it is very close to the best known result in the
literature. For example, the best known minimum RMSE so far for the RTC France
benchmark with the DDM is 9.8248E-4 (see, e.g, results listed in (Chen et al.,
2016, Table 6), (Zhang et al., 2020, Table 11), (Zhang et al., 2020, Table
10), (Li et al., 2019, Table 4), (Yu et al., 2019, Table 3), (Yousri et al.,
2019, Table 1), and (Gnetchejo et al., 2019, Table 11)). In view of these
facts, the acquired upper bound 9.8596E-4 in Table 3 implies that 9.8248E-4 is
likely to be the global minimum value despite the lack of theoretical
guarantee. In addition, just like the above SDM case, the seven parameter
values reported in previous tables we have just cited all reside in a small
neighborhood of the result located by ibexopt in Table 3.

As for the Photowatt-PWP201 benchmark, most existing studies only employ the
SDM probably due to the increased difficulty of the DDM. Nevertheless, there
is no technical reason that prevents the application of the DDM to a PV
module. However, to our best knowledge, the limited existing studies that
consider the DDM and the Photowatt-PWP201 dataset, such as those reported in
(Nunes et al., 2018, Table 12) and (Kler et al., 2019, Table 14), adopt a
procedure distinct from the prevailing one: they use the Newton-Raphson method
rather than (8) to compute current values iteratively (refer to Section 3 for
more details). Consequently, unlike previous cases, we did not find any
compatible results that can be used for DDM’s perforamnce comparison with
respect to the Photowatt-PWP201 PV module here. In spite of the lack of
references, the result we have obtained in Table 3 looks reasonable if we
compare it with its SDM counterpart in Table 2. Besides, unlike the RTC France
case, the upper bound calculated by ibexopt did not reduce any longer after 12
hours in the PV module case.

###### Remark 4

Apart from the disparate current calculation procedures as emphasized above,
the parameter ranges given in (Nunes et al., 2018, Table 2) and (Kler et al.,
2019, Table 2), especially ranges of $I_{01}$ and $I_{02}$, are totally
different from the most widely used ones in Table 1, which thus blocks
reasonable comparison between these parameter estimation results.

####  5.2.3 Summary of global optimality analysis

Table 4: Overview of the optimization results of interval B&B in four
benchmark cases.

| SDM | DDM
---|---|---
| RTC France | Photowatt-PWP201 | RTC France | Photowatt-PWP201
RMSE |  9.8602E-4* |  2.4250E-3* |  9.8358E-4† |  1.6186E-3†
Time (s) | 13547 | 38924 | 86400 | 43200

  * *

Certified global minimum value with respect to five significant digits

  * †

Certified upper bound of the global minimum value in five significant digits

(a) RTC France

(b) Photowatt-PWP201

Figure 3: Measured data and calculated data using both the SDM and the DDM
with parameters optimized by interval B&B.

The results of the interval B&B method for the four representative cases,
reported previously in Table 2 and Table 3, are collected here in Table 4 to
give a concise overview. Notably, the interval B&B algorithm (via ibexopt) is
capable of locating the $\epsilon$-global minimum rigorously for both
benchmark datasets using the five-parameter SDM, though its running time is
exceedingly long compared with metaheuristic methods. For the more complicated
seven-parameter DDM, the interval B&B method failed to obtain an adequately
tight optimality bound in 24 hours, but the acquired upper bound can still
play the role of a good reference value. To demonstrate visually the accuracy
of the extracted parameters listed in Table 2 and Table 3, we insert them into
(8) to reconstruct the I-V curves, which are shown in Fig. 3. Note that the
negative current and voltage values therein simply imply a reverse direction
(Easwarakhanthan et al., 1986). Obviously, there is excellent consistency
between the measured data and the simulated data, even using the DDM whose
parameter values correspond to the RMSE upper bound in Table 3. Comparing the
RMSE values in Table 2 and Table 3 (summarized in Table 4), we notice that the
RMSE (or its upper bound) of the DDM is consistently smaller than that of the
SDM. This situation is expected and makes complete sense because the several-
parameter DDM is more flexible than the SDM and, in particular, the SDM (4) is
essentially a special instance of the DDM by fixing $I_{02}=0$ in (2.2).

###### Remark 5

We emphasize that, to avoid possible misunderstanding, the global minimum
value $J^{*}$ is defined in terms of the objective function value (10a), and
it does not imply that there is necessarily a unique global minimum point
$\bm{\theta}^{*}$ in (10b) corresponding to $J^{*}$. Moreover, since our
discussion is conventionally limited to five significant digits, it is not
surprising that many close but not identical parameter values can lead to the
same objective value. For example, slightly different values of $R_{p}$ are
identified by three methods in (Li et al., 2019, Table 3) though they all
obtain the same global minimum 9.8602E-4 for the first case in Table 4. This
fact is also noticeable by comparing parameter values in Table 2 and Table 6.

###  5.3 Optimization via simple differential evolution (DE)

As we have mentioned above, despite its capability to guarantee global
optimality (or optimality bounds), the interval B&B algorithm is not suitable
for practical estimation of PV parameters owing to its excessively long
execution time. In addition, a guaranteed global minimum may not be necessary
for real-world industrial applications. By contrast, a variety of
metaheuristic methods with distinct methodology and complexity can obtain a
reasonably small objective value in a far shorter time. Nevertheless, in view
of the abundance of metaheuristic algorithms, an interesting but rarely
studied question, as we have asked in Section 1, is whether such PV parameter
estimation really demands the increasingly complicated metaheuristics in the
recent literature. The answer to this puzzle is of both theoretical
implication and practical importance. Unfortunately, it is extremely
difficult, if ever possible, to perform theoretical analysis of stochastic
search based metaheuristics (Yang, 2020). In this section, we try to get some
insights from another perspective by examining whether the intentionally
simple differential evolution (DE) method in Algorithm 2 can achieve
comparable performance for PV parameter estimation using both the SDM and the
DDM. The two widely used benchmark datasets introduced in Section 5.1 are
considered again for comparison purposes. In addition to the informative
reference values acquired by interval B&B in Section 5.2, the accuracy and
efficiency of DE will also be evaluated against other popular metaheuristic
methods.

Starting with the recommendation in (Das and Suganthan, 2011, Section III), we
quickly determined appropriate control parameter values for DE in Algorithm 2
after several preliminary trials, which are listed in Table 5. Note that these
values are kind of canonical in the literature, and little time was spent in
parameter tuning thanks to the simplicity of Algorithm 2. Remarkably, both
models share almost the same control parameters on the two datasets: the only
difference lies in the number of generations. Intuitively, the more
complicated DDM requires more generations of DE than the simpler SDM.

Table 5: Control parameters of simple DE in Algorithm 2. | $N_{p}$ | $C_{r}$ | $F$ | $G$
---|---|---|---|---
SDM | 50 | 0.6 | 0.9 | 800
DDM | 50 | 0.6 | 0.9 | 1600

####  5.3.1 SDM and DDM results

Since DE in Algorithm 2 is a stochastic search algorithm, we followed the
convention in the literature, e.g., (Liang et al., 2020b; Jordehi, 2016b;
Yousri et al., 2020; Yu et al., 2019), and executed DE 30 times for each case
in order to obtain a trustworthy judgment of its accuracy and consistency. The
parameter values obtained by DE in a typical run are listed in Table 6. The
corresponding I-V curves reconstructed with these parameters are illustrated
in Fig. 4 to facilitate visual inspection. Since the RMSE values are fairly
small in all cases (see Table 4 and Table 6), there is, unsurprisingly, almost
no visual difference between Fig. 3 and Fig. 4.

Table 6: Parameter values obtained by DE in a typical run. (RT: RTC France; PW: Photowatt-PWP201) | SDM | DDM
---|---|---
| RT | PW | RT | PW
$I_{ph}$(A) | 0.760775 | 1.03051 | 0.760781 | 1.03051
$I_{0}/I_{01}$$(\mu\text{A})$ | 0.323021 | 3.48226 | 0.225974 | 9.77791E-3
$I_{02}$($\mu$A) | — | — | 0.749344 | 3.47248
$n/n_{1}$ | 1.481184 | 48.6428 | 1.45101 | 48.64284
$n_{2}$ | — | — | 1.99999 | 48.64283
$R_{s}$($\Omega$) | 0.036377 | 1.20127 | 0.0367404 | 1.20127
$R_{p}$($\Omega$) | 53.71852 | 981.982 | 55.4854 | 981.982
RMSE | 9.8602E-4 | 2.4250E-3 | 9.8248E-4 | 2.4250E-3

(a) RTC France

(b) Photowatt-PWP201

Figure 4: Measured data and calculated data by both the SDM and the DDM with
parameter values identified by DE.

The convergence curves of this simple DE in a typical run on the RTC France
dataset using both models are shown in Fig. 5. As we have observed in many
runs, just like Fig. 5 here, the DE often takes far fewer generations to
converge than the number $G$ specified in Table 5, which is a conservative
choice to ensure better convergence. The convergence curves on the other
dataset share a similar character and are presented in our online materials.

(a) SDM

(b) DDM

Figure 5: Convergence curves of DE on the RTC France dataset for both models
in a typical run.

The statistics of RMSE in the 30 runs are analyzed and reported in Table 7.
Overall, we see from Table 7 that the performance of our DE algorithm is
remarkably stable despite its stochasticity in nature. When applying to the
SDM on the two datasets and to the DDM on the Photowatt-PWP201 dataset, the DE
algorithm always yields the same minimal RMSE (in terms of five significant
digits) in all 30 trials. Even in the worst case (DDM and RT), the maximum and
minimum RMSE values obtained in 30 runs are still very close to each other,
which is also implied by the small standard deviation in Table 7.

Table 7: Statistics of RMSE values yielded by DE in 30 runs. (RT: RTC France; PW: Photowatt-PWP201)  | SDM | DDM
---|---|---
| RT | PW | RT | PW
Min | 9.8602E-4 | 2.4250E-3 | 9.8248E-4 | 2.4250E-3
Mean | 9.8602E-4 | 2.4250E-3 | 9.8267E-4 | 2.4250E-3
Max | 9.8602E-4 | 2.4250E-3 | 9.8602E-4 | 2.4250E-3
Std | 4.3929E-17 | 2.9525E-17 | 7.1027E-7 | 2.3955E-17

The simple DE algorithm managed to find the global minimum for the first two
cases in Table 7, i.e., for SDM on the two datasets. As for the more
challenging DDM cases, DE achieved the best known result in the literature for
the RTC France PV cell, i.e., 9.8248E-4. Nevertheless, comparing the last
case’s RMSE in Table 7 with its upper bound in Table 3, we see that DE failed
to accomplish an RMSE value below the upper bound (2.4250E-3 vs 1.6186E-3).
The obtained minimum RMSE thus cannot be the global minimum value in the last
case. Even so, we should not depreciate the simple DE’s performance hastily in
this particular case due to lack of usable comparison in the current
literature (recall the last paragraph of Section 5.2.2 for elucidation): no
comparable results, to our best knowledge, have been reported towards modeling
the Photowatt-PWP201 PV module with the DDM.

More specifically, by focusing on the last column of Table 6 (i.e., the DDM
fitted to the Photowatt-PWP201 module), we notice that the obtained RMSE is
equal to its SDM counterpart. This result is not surprising since the SDM is a
specialization of the DDM as we have mentioned in Remark 5. The main reason
that DE failed to locate an RMSE value within the optimality bounds revealed
by interval B&B in Table 3 is presumably attributed to the extremely small
value of $I_{01}$ therein, whose order of magnitude is -23. Such a tiny number
is commonly treated as zero in general computation, though it still has
perceivable impact on the RMSE in this specific case. Unfortunately, it is
really challenging for DE to pinpoint such a small value within a relatively
large search range $[0,50]$ (see Table 1). In fact, in the existing yet
incompatible studies (see Section 5.2.2) that consider the last case (“DDM +
PW”), a greatly narrowed search range is usually set for $I_{01}$ and $I_{02}$
(Nunes et al., 2018; Kler et al., 2019; Jordehi, 2016b), e.g., [1E-15, 1E-3]
in (Kler et al., 2019), which facilitates the discovery of good parameter
values. More interestingly, in the two PV module (“PW”) cases in Table 6, the
parameter values acquired by DE essentially degrade the DDM to the SDM by
$n_{1}=n_{2}\approx n$ and $I_{01}+I_{02}\approx I_{0}$, which can be figured
out easily by comparing Eq. (2.2) to (4). This argument is also justified by
the indistinguishable overlap between the SDM and DDM curves in Fig. 4(b) (as
opposed to the small but still observable gap between the two curves in Fig.
3(b)).

####  5.3.2 Comparison with existing algorithms

The parameter values obtained by DE in Table 6 are quite close to those
reported in existing studies and a detailed comparison of the values
themselves is omitted here to save space. Interested readers may refer to
tables in some latest articles like (Liang et al., 2020b; Long et al., 2020)
to examine parameter values obtained in various studies. In general, almost
identical parameter values are acquired as long as the resultant RMSE values
are similar, a fact that verifies the correctness of each other. We devote the
subsequent comparison to different algorithms’ performance statistics and time
efficiency, both being important factors in real-world applications.

For fair comparison, we mainly copied the results directly reported in recent
papers that have published state-of-the-art algorithms. No existing
investigations, however, target the Photowatt-PWP201 module using the DDM that
are compatible with our study here (see Section 5.2.2 for an explanation).
Besides, the relevant metaheuristic studies in the literature rarely provide
their source code to the public, which renders it cumbersome and inconvenient
to test their algorithms on new cases (Chenouard and El-Sehiemy, 2020). To our
best efforts, we only got the source code of two recent methods, namely
performance-guided JAYA (PGJAYA) (Yu et al., 2019) and self-adaptive ensemble-
based differential evolution (SEDE)888Available at
https://github.com/cilabzzu/Publications/tree/master/PGJAYA and
https://github.com/cilabzzu/Publications/tree/master/SEDE, respectively.
(Liang et al., 2020a). We ran the two algorithms to deal with the last case,
i.e., to estimate DDM parameters for the PV module, using their original
control parameter settings and report the results at the bottom of Table 8.

Table 8: Comparison of statistical results of various algorithms in four
benchmark case studies. (RT: RTC France; PW: Photowatt-PWP201)

Model | Data | Algorithm | RMSE | Time (s)
---|---|---|---|---
Min | Mean | Max | Std
SDM | RT |  MLBSA (Yu et al., 2018) | 9.8602E-4 | 9.8602E-4 | 9.8602E-4 | 7.0800E-11 | 28
TLABC (Chen et al., 2018) | 9.8602E-4 | 9.9417E-4 | 1.0308E-3 | 1.1896E-5 | 43
IJAYA (Yu et al., 2017c) | 9.8602E-4 | 9.8605E-4 | 9.8684E-4 | 1.4931E-7 | 25
PGJAYA (Yu et al., 2019) | 9.8602E-4 | 9.8602E-4 | 9.8603E-4 | 2.8029E-9 | 25
SATLBO (Yu et al., 2017a) | 9.8602E-4 | 9.8879E-4 | 1.0067E-3 | 4.8133E-6 | 35
SGDE (Liang et al., 2020b) | 9.8602E-4 | 9.8602E-4 | 9.8603E-4 | 2.4746E-9 |  —*
SEDE (Liang et al., 2020a) | 9.8602E-4 | 9.8602E-4 | 9.8603E-4 | 4.2000E-17 | 28
CoDE (Wang et al., 2011) | 9.8602E-4 | 9.8602E-4 | 9.8602E-4 | 2.3100E-17 | 25
Simple DE | 9.8602E-4 | 9.8602E-4 | 9.8602E-4 | 4.3929E-17 | 0.34
PW | MLBSA | 2.4250E-3 | 2.4253E-3 | 2.4336E-3 | 1.5600E-6 | 28
TLABC | 2.4250E-3 | 2.4254E-3 | 2.4287E-3 | 8.7464E-7 | 42
IJAYA | 2.4250E-3 | 2.4251E-3 | 2.4253E-3 | 5.0766E-8 | 27
PGJAYA | 2.4250E-3 | 2.4251E-3 | 2.4260E-3 | 1.7859E-7 | 24
SATLBO | 2.4250E-3 | 2.4254E-3 | 2.4315E-3 | 1.1622E-6 | 38
SGDE | 2.4250E-3 | 2.4250E-3 | 2.4250E-3 | 4.1697E-10 |  —*
SEDE | 2.4250E-3 | 2.4250E-3 | 2.4250E-3 | 3.1400E-17 | 28
CoDE | 2.4250E-3 | 2.4250E-3 | 2.4250E-3 | 2.1700E-17 | 27
Simple DE | 2.4250E-3 | 2.4250E-3 | 2.4250E-3 | 2.9525E-17 | 0.31
DDM | RT | MLBSA | 9.8248E-4 | 9.8506E-4 | 9.8613E-4 | 1.2400E-6 | 33
TLABC | 1.0012E-3 | 1.2116E-3 | 1.9826E-3 | 2.1100E-4 | 44
IJAYA | 9.8249E-4 | 9.8686E-4 | 9.9941E-4 | 3.2211E-6 | 26
PGJAYA | 9.8260E-4 | 9.8603E-4 | 9.9599E-4 | 2.3666E-6 | 28
SATLBO | 9.8282E-4 | 1.0054E-3 | 1.2306E-3 | 5.0271E-5 | 37
SGDE | 9.8441E-4 | 9.8577E-4 | 9.8602E-4 | 4.0150E-7 |  —*
SEDE | 9.8248E-4 | 9.8289E-4 | 9.8602E-4 | 9.1700E-7 | 28
CoDE | 9.8249E-4 | 1.0036E-3 | 1.5496E-3 | 1.0300E-4 | 25
Simple DE | 9.8248E-4 | 9.8267E-4 | 9.8602E-4 | 7.1027E-7 | 0.80
PW | PGJAYA | 2.4250E-3 | 2.4272E-3 | 2.4485E-3 | 5.4346E-6 |  373†
SEDE | 2.4250E-3 | 2.4250E-3 | 2.4250E-3 | 6.6661E-17 |  324†
Simple DE | 2.4250E-3 | 2.4250E-3 | 2.4250E-3 | 2.3955E-17 | 0.84

  * *

Since SGDE extracts three model parameters simultaneously, its runtime is not
directly comparable here and thus not reported by (Liang et al., 2020a).

  * †

Time measured on our PC using MATLAB R2020a and the source code provided by
the original studies. All the other runtime results of competing methods are
extracted from the literature.

Since metaheuristic algorithms are stochastic in nature, a common practice is
to evaluate their statistical characteristics in 30 independent runs (Liang et
al., 2020b; Jordehi, 2016b; Liang et al., 2020a; Li et al., 2019). In light of
the abundance of studies on this topic, we extracted the statistical results
of selected methods, especially the state-of-the-art ones, from those reported
in four latest articles, (Li et al., 2019, Table 12), (Yu et al., 2019, Table
3), (Liang et al., 2020b, Table 9), and (Liang et al., 2020a, Table 9). The
selection was guided by the following key criteria: first, choose state-of-
the-art ones according to the RMSE they have attained; second, select
deliberately algorithms of distinct methodology and particularly DE variants
for a comprehensive comparison. Interested readers may refer to the cited
tables for other well-established methods. The average CPU time of these
methods were mainly read from (Liang et al., 2020a, Fig. 7) (except TLABC) and
(Yu et al., 2019, Fig. 5) (for only TLABC). The results are collected and
listed in Table 8. We have the following interesting observations and
discussions regarding Table 8.

  * 1.

Overall, it is clear that our simple DE in Algorithm 2 and several strong
competitors like SEDE can attain the best RMSE values. Moreover, if we take
the complete statistics into consideration, DE exhibits superior performance,
e.g., with regard to the mean RMSE and the worst RMSE values. This argument is
in agreement with the fact that RMSE values acquired by DE exhibit tiny
standard deviations, among the top three smallest ones in all cases.

  * 2.

We compare the simple DE intentionally with three more complicated DE
variants, namely, SGDE, SEDE, and CoDE. In the two SDM cases, the four DE
algorithms exhibit almost identical accuracy in terms of RMSE values, while
CoDE demonstrates the highest stability measured by the standard deviation and
SGDE is the least stable one. By contrast, in the more challenging seven-
parameter DDM cases, our simple DE outperforms both CoDE and SEDE with its
enhanced performance stability. The favorable consistency of DE’s performance
across multiple runs is probably attributed to its simplistic mechanism that
introduces randomness in (13) and (15).

  * 3.

Compared with the other algorithms, the biggest advantage of the simple DE
highlighted in Table 8 is definitely its substantially reduced running time.
DE stands out from these competing techniques by taking virtually only one
percent of the CPU time consumed by others. This impressive speedup is mainly
brought by its extreme simplicity that includes only four computationally
cheap equations in the main loop of Algorithm 2. Moreover, a common metric to
evaluate the efficiency of metaheuristics is the number of objective function
evaluations termed $N_{FE}$. There is no standard value of $N_{FE}$ in the PV
parameter estimation literature though a typical value is 50000 (Liang et al.,
2020a; Yu et al., 2019; Li et al., 2019; Chen et al., 2016). We acknowledge
that an inevitable tradeoff of the simplicity of our DE is the (moderately)
increased number of function evaluations for comparable performance. To be
specific, the maximum number of function evaluations in DE is 40000 for the
SDM and 80000 for the DDM, respectively (see Table 5). Nonetheless, note that
the calculation of the objective (10a) with a few dozens of data points is
inexpensive. Consequently, it is the algorithm’s internal computation burden
rather than the objective function evaluation that tends to dominate the
overall time consumption, at least in the four cases here. This fact is
implied particularly by the significantly longer runtime of the three more
complex DE variants in Table 8.

  * 4.

Focusing on the last case (i.e., “DDM+PW”) in Table 8, we notice that the best
RMSE value attained by the three algorithms all turns to be 2.4250E-3, though
this value is certainly not the global minimum, since an upper bound of the
global minimum has been identified as 1.6186E-3 in Table 4. A possible reason
for this common failure, as we have speculated in Section 5.3.1, is that the
ideal value of $I_{01}$ is unreasonably small, which poses a tremendous
challenge to these metaheuristic methods since extraordinarily high precision
is required to get such a value for $I_{01}$ inside a wide range.

###### Remark 6

It should be pointed out that, among the runtime listed in Table 8, only the
bottom two of PGJAYA and SEDE in the “DDM+PW” case were measured on our
desktop PC. We understand that different computing environments and certain
implementation details may affect more or less the execution time. Note,
however, that the computing power of the server that produces other reference
runtime in Table 8 from (Yu et al., 2019) and (Liang et al., 2020a) is at
least as strong as our PC. Overall, the time reduction of our simple DE is
around two orders of magnitude, whose efficiency advantage is indeed
compelling regardless of these possible external disturbances.

###### Remark 7

Usually, far fewer generations are enough for the simple DE to converge in
reality (recall Fig. 5). Techniques like early stopping may help reduce
further the execution time, though they are not considered in this study. As a
side note, another “DE” algorithm (Ishaque and Salam, 2011) is also mentioned
in (Liang et al., 2020b, Table 9). However, that DE (Ishaque and Salam, 2011)
depends exclusively on the manufacturer’s datasheet information and is
distinct from our methodology here.

##  6 Conclusion

In this paper, we conducted a critical study of solar PV parameter estimation
from a perspective distinct from the majority of existing work. We tried to
address two essential issues that have seldom been attempted in the current
literature. In regard to the two most widely used benchmark datasets, the
globally minimum RMSE for the SDM and a reasonably tight upper bound for the
DDM were certified rigorously with a deterministic interval analysis based
branch and bound (B&B) algorithm. However, the running time of this interval
B&B algorithm is overly long for practical usages despite its theoretical
guarantee. Next, we showed through extensive examination that, for the first
time and somewhat surprisingly, a simple DE algorithm was able to locate the
global minimum or at least attain the best known result. More importantly, the
simple DE algorithm is indeed striking due to its favorable performance
stability and its unmatched efficiency. Therefore, we suggest that a
practitioner should start with the simple DE outlined in Algorithm 2 as the
first off-the-shelf tool, especially for real-time parameter estimation
applications, where the environment keeps changing (Chen et al., 2016).

According to the results obtained above in this study, we present the
following comments and suggestions regarding the two questions in Section 1
and for future work.

  1. 1.

The globally minimal RMSE values related to the SDM for both the RTC France PV
cell and the Photowatt-PWP201 PV module have been certified in this study. It
has already been reached by various metaheuristics. It is thereupon futile for
researchers to develop more sophisticated metaheuristics to pursue even
smaller RMSE. On the one hand, this certified global minimum may be viewed as
a golden reference when examining a metaheuristic’s fundamental capability. On
the other hand, since the global minimum can be attained by many recent
metaheuristic methods (see Table 8), the benchmark role of the SDM with the
two datasets in measuring the performance of various methods, still a de facto
standard in the literature, seems to be undermined and not fully convincing.

  2. 2.

The global minimum in the two DDM cases has not been identified, since the
interval B&B method can only yield a not so tight optimality bound (see Table
3). Even so, no metaheuristics can reach the pinpointed upper bound yet in the
“DDM+PW” case (recall Table 8 and Table 4). Overall, the more challenging DDM
cases may serve as better benchmarks than the SDM to assess various
metaheuristics. Besides, the certification of global optimality demands more
advanced global optimization techniques (Floudas and Gounaris, 2008) and
deserves more investigations.

  3. 3.

In all the four test cases, the simple DE depicted by Algorithm 2 performs
surprisingly equally well when compared with a variety of sophisticated
metaheuristic methods apart from its unparalleled efficiency advantage. We
surely admire the power of advanced metaheuristics like those in Table 8, but
it seems unfortunate that the competence of simpler algorithms has been
largely underestimated in the current literature for PV parameter estimation.
We suggest that future researchers justify their methods by comparing against
the simple DE here, which is an appropriate baseline for forthcoming studies.
Besides, further examination of other classic metaheuristics, like the
canonical particle swarm optimization, is also meaningful.

  4. 4.

It should be emphasized that we are definitely not underrating the value of
sophisticated metaheuristics in a general sense. Such complexity is necessary
for many difficult real-world tasks (Yang, 2020). However, note that the
normal PV parameter estimation is essentially a straightforward nonlinear
least squares problem with only five or seven variables. A reasonable
speculation is that the present problem is not tough enough to highlight the
real power of such advanced metaheuristics. Instead, the simple DE in
Algorithm 2 appears to be a better match for this specific problem from a
pragmatic view. We advocate the examination of more challenging scenarios in
future investigations, where there is still vast room for improvement, e.g.,
the DDM with the Photowatt-PWP201 PV module dataset, datasets of larger sizes,
datasets in different temperatures, and the less commonly studied but more
demanding three diode model (Yang et al., 2020; Yousri et al., 2020; Qais et
al., 2019).

  5. 5.

Given the abundance of metaheuristics in the literature, we believe that more
efforts may be better put into other aspects apart from fitting accuracy,
which may include rigorous certification of global optimality, satisfaction of
practical industry needs, and construction of new high-quality datasets, among
others.

## Acknowledgements

We would like to thank Dr. Gilles Chabert for development of ibex and helpful
discussions. This work was supported by the Exploiting Distributed Generation
(EDGE) programme grant (No: EDGE-GC-2018-002).

## Data Availability

The two benchmark datasets and our code are publicly accessible at
https://github.com/ShuhuaGao/rePVest.

## References

  * Abdulrazzaq et al. (2020) Abdulrazzaq, A.K., Bognár, G., Plesz, B., 2020.  Evaluation of different methods for solar cells/modules parameters extraction.  Solar Energy 196, 183–195.
  * Almonacid et al. (2010) Almonacid, F., Rus, C., Hontoria, L., Muñoz, F.J., 2010\.  Characterisation of PV CIS module by artificial neural networks. A comparative study with other methods.  Renewable Energy 35, 973–980.
  * Askarzadeh and Rezazadeh (2013) Askarzadeh, A., Rezazadeh, A., 2013\.  Artificial bee swarm optimization algorithm for parameters identification of solar cell models.  Applied Energy 102, 943–949.
  * Bilal et al. (2020) Bilal, Pant, M., Zaheer, H., Garcia-Hernandez, L., Abraham, A., 2020\.  Differential Evolution: A review of more than two decades of research.  Engineering Applications of Artificial Intelligence 90, 103479.
  * Burer and Letchford (2012) Burer, S., Letchford, A.N., 2012\.  Non-convex mixed-integer nonlinear programming: A survey.  Surveys in Operations Research and Management Science 17, 97–106.
  * Ćalasan et al. (2020) Ćalasan, M., Abdel Aleem, S.H.E., Zobaa, A.F., 2020.  On the root mean square error (RMSE) calculation for parameter estimation of photovoltaic models: A novel exact analytical solution based on Lambert W function.  Energy Conversion and Management 210, 112716.
  * Chen et al. (2018) Chen, X., Xu, B., Mei, C., Ding, Y., Li, K., 2018\.  Teaching–learning–based artificial bee colony for solar photovoltaic parameter estimation.  Applied Energy 212, 1578–1588.
  * Chen et al. (2016) Chen, Z., Wu, L., Lin, P., Wu, Y., Cheng, S., 2016\.  Parameters identification of photovoltaic models using hybrid adaptive Nelder-Mead simplex algorithm based on eagle strategy.  Applied Energy 182, 47–57.
  * Chenouard and El-Sehiemy (2020) Chenouard, R., El-Sehiemy, R.A., 2020\.  An interval branch and bound global optimization algorithm for parameter estimation of three photovoltaic models.  Energy Conversion and Management 205, 112400.
  * Chin et al. (2015) Chin, V.J., Salam, Z., Ishaque, K., 2015.  Cell modelling and model parameters estimation techniques for photovoltaic simulator application: A review.  Applied Energy 154, 500–519.
  * Cotfas et al. (2013) Cotfas, D.T., Cotfas, P.A., Kaplanis, S., 2013.  Methods to determine the dc parameters of solar cells: A critical review.  Renewable and Sustainable Energy Reviews 28, 588–596.
  * Das and Suganthan (2011) Das, S., Suganthan, P.N., 2011\.  Differential Evolution: A Survey of the State-of-the-Art.  IEEE Transactions on Evolutionary Computation 15, 4–31.
  * Dizqah et al. (2014) Dizqah, A.M., Maheri, A., Busawon, K., 2014.  An accurate method for the PV model identification based on a genetic algorithm and the interior-point method.  Renewable Energy 72, 212–222.
  * Easwarakhanthan et al. (1986) Easwarakhanthan, T., Bottin, J., Bouhouch, I., Boutrit, C., 1986\.  Nonlinear minimization algorithm for determining the solar cell parameters with microcomputers.  International journal of solar energy 4, 1–12.
  * El-Naggar et al. (2012) El-Naggar, K.M., AlRashidi, M., AlHajri, M., Al-Othman, A., 2012\.  Simulated annealing algorithm for photovoltaic parameters identification.  Solar Energy 86, 266–274.
  * Floudas (2013) Floudas, C.A., 2013.  Deterministic global optimization: theory, methods and applications. volume 37.  Springer Science & Business Media.
  * Floudas and Gounaris (2008) Floudas, C.A., Gounaris, C.E., 2008\.  A review of recent advances in global optimization.  Journal of Global Optimization 45, 3.
  * Gao et al. (2018) Gao, S., Xiang, C., Qin, K., Sun, C., 2018.  Mathematical modeling reveals the role of hypoxia in the promotion of human mesenchymal stem cell long-term expansion.  Stem cells international 2018.
  * Gnetchejo et al. (2019) Gnetchejo, P.J., Ndjakomo Essiane, S., Ele, P., Wamkeue, R., Mbadjoun Wapet, D., Perabi Ngoffe, S., 2019\.  Important notes on parameter estimation of solar photovoltaic cell.  Energy Conversion and Management 197, 111870.
  * Gong and Cai (2013) Gong, W., Cai, Z., 2013.  Parameter extraction of solar cell models using repaired adaptive differential evolution.  Solar Energy 94, 209–220.
  * Guo et al. (2016) Guo, L., Meng, Z., Sun, Y., Wang, L., 2016.  Parameter identification and sensitivity analysis of solar cell models with cat swarm optimization algorithm.  Energy Conversion and Management 108, 520–528.
  * Hansen and Walster (2003) Hansen, E., Walster, G.W., 2003\.  Global optimization using interval analysis: revised and expanded. volume 264.  CRC Press.
  * Ishaque and Salam (2011) Ishaque, K., Salam, Z., 2011\.  An improved modeling method to determine the model parameters of photovoltaic (PV) modules using differential evolution (DE).  Solar Energy 85, 2349–2359.
  * Jordehi (2016a) Jordehi, A.R., 2016a.  Parameter estimation of solar photovoltaic (PV) cells: A review.  Renewable and Sustainable Energy Reviews 61, 354–371.
  * Jordehi (2016b) Jordehi, A.R., 2016b.  Time varying acceleration coefficients particle swarm optimisation (TVACPSO): A new optimisation algorithm for estimating parameters of PV cells and modules.  Energy Conversion and Management 129, 262–274.
  * Kler et al. (2019) Kler, D., Goswami, Y., Rana, K.P.S., Kumar, V., 2019\.  A novel approach to parameter estimation of photovoltaic systems using hybridized optimizer.  Energy Conversion and Management 187, 486–511.
  * Li et al. (2019) Li, S., Gong, W., Yan, X., Hu, C., Bai, D., Wang, L., Gao, L., 2019.  Parameter extraction of photovoltaic models using an improved teaching-learning-based optimization.  Energy Conversion and Management 186, 293–305.
  * Li and Xu (2018) Li, Z., Xu, Y., 2018.  Optimal coordinated energy dispatch of a multi-energy microgrid in grid-connected and islanded modes.  Applied Energy 210, 974–986.
  * Liang et al. (2020a) Liang, J., Qiao, K., Yu, K., Ge, S., Qu, B., Xu, R., Li, K., 2020a.  Parameters estimation of solar photovoltaic models via a self-adaptive ensemble-based differential evolution.  Solar Energy 207, 336–346.
  * Liang et al. (2020b) Liang, J., Qiao, K., Yuan, M., Yu, K., Qu, B., Ge, S., Li, Y., Chen, G., 2020b.  Evolutionary multi-task optimization for parameters extraction of photovoltaic models.  Energy Conversion and Management 207, 112509.
  * Long et al. (2020) Long, W., Cai, S., Jiao, J., Xu, M., Wu, T., 2020\.  A new hybrid algorithm based on grey wolf optimizer and cuckoo search for parameter extraction of solar photovoltaic models.  Energy Conversion and Management 203, 112243.
  * McCormick (1976) McCormick, G.P., 1976.  Computability of global solutions to factorable nonconvex programs: Part i—convex underestimating problems.  Mathematical programming 10, 147–175.
  * Montanher et al. (2018) Montanher, T., Neumaier, A., Domes, F., 2018.  A computational study of global optimization solvers on two trust region subproblems.  Journal of Global Optimization 71, 915–934.
  * Navabi et al. (2015) Navabi, R., Abedi, S., Hosseinian, S.H., Pal, R., 2015\.  On the fast convergence modeling and accurate calculation of PV output energy for operation and planning studies.  Energy Conversion and Management 89, 497–506.
  * Nunes et al. (2018) Nunes, H.G.G., Pombo, J.A.N., Mariano, S.J.P.S., Calado, M.R.A., Felippe de Souza, J.A.M., 2018.  A new high performance method for determining the parameters of PV cells and modules based on guaranteed convergence particle swarm optimization.  Applied Energy 211, 774–791.
  * Oliva et al. (2017) Oliva, D., Abd El Aziz, M., Hassanien, A.E., 2017.  Parameter estimation of photovoltaic cells using an improved chaotic whale optimization algorithm.  Applied Energy 200, 141–154.
  * Price et al. (2006) Price, K., Storn, R.M., Lampinen, J.A., 2006.  Differential evolution: a practical approach to global optimization.  Springer Science & Business Media.
  * Qais et al. (2019) Qais, M.H., Hasanien, H.M., Alghuwainem, S., 2019.  Identification of electrical parameters for three-diode photovoltaic model using analytical and sunflower optimization algorithm.  Applied Energy 250, 109–117.
  * Sheng et al. (2019) Sheng, H., Li, C., Wang, H., Yan, Z., Xiong, Y., Cao, Z., Kuang, Q., 2019.  Parameters extraction of photovoltaic models using an improved moth-flame optimization.  Energies 12, 3527\.
  * Tong and Pora (2016) Tong, N.T., Pora, W., 2016.  A parameter extraction technique exploiting intrinsic properties of solar cells.  Applied Energy 176, 104–115.
  * Villalva et al. (2009) Villalva, M.G., Gazoli, J.R., Filho, E.R., 2009.  Comprehensive Approach to Modeling and Simulation of Photovoltaic Arrays.  IEEE Transactions on Power Electronics 24, 1198–1208.
  * Wang et al. (2011) Wang, Y., Cai, Z., Zhang, Q., 2011.  Differential evolution with composite trial vector generation strategies and control parameters.  IEEE transactions on evolutionary computation 15, 55–66.
  * Yang et al. (2020) Yang, B., Wang, J., Zhang, X., Yu, T., Yao, W., Shu, H., Zeng, F., Sun, L., 2020.  Comprehensive overview of meta-heuristic algorithm applications on PV cell parameter identification.  Energy Conversion and Management 208, 112595.
  * Yang (2020) Yang, X.S., 2020.  Nature-inspired optimization algorithms: Challenges and open problems.  Journal of Computational Science , 101104\.
  * Yousri et al. (2020) Yousri, D., Abd Elaziz, M., Oliva, D., Abualigah, L., Al-qaness, M.A.A., Ewees, A.A., 2020\.  Reliable applied objective for identifying simple and detailed photovoltaic models using modern metaheuristics: Comparative study.  Energy Conversion and Management 223, 113279.
  * Yousri et al. (2019) Yousri, D., Elaziz, M.A., Razaee, A., Merchaoui, M., Rana, K.P.S., Babu, T.S., Oliva, D., Ram, P., Rajasekar, N., Alam, D.F., Eteiba, M.B., Kler, D., Goswami, Y., Kumar, V., 2019\.  Comment on “Important notes on parameter estimation of solar photovoltaic cell”, by Gnetchejo et al. [Energy Conversion and Management, https://doi.org/10.1016/j.enconman.2019.111870].  Energy Conversion and Management 201, 112131.
  * Yu et al. (2017a) Yu, K., Chen, X., Wang, X., Wang, Z., 2017a.  Parameters identification of photovoltaic models using self-adaptive teaching-learning-based optimization.  Energy Conversion and Management 145, 233–246.
  * Yu et al. (2017b) Yu, K., Liang, J., Qu, B., Chen, X., Wang, H., 2017b.  Parameters identification of photovoltaic models using an improved jaya optimization algorithm.  Energy Conversion and Management 150, 742–753.
  * Yu et al. (2017c) Yu, K., Liang, J.J., Qu, B.Y., Chen, X., Wang, H., 2017c.  Parameters identification of photovoltaic models using an improved JAYA optimization algorithm.  Energy Conversion and Management 150, 742–753.
  * Yu et al. (2018) Yu, K., Liang, J.J., Qu, B.Y., Cheng, Z., Wang, H., 2018\.  Multiple learning backtracking search algorithm for estimating parameters of photovoltaic models.  Applied Energy 226, 408–422.
  * Yu et al. (2019) Yu, K., Qu, B., Yue, C., Ge, S., Chen, X., Liang, J., 2019.  A performance-guided JAYA algorithm for parameters identification of photovoltaic cell and module.  Applied Energy 237, 241–257.
  * Zaimi et al. (2019) Zaimi, M., El Achouby, H., Ibral, A., Assaid, E.M., 2019\.  Determining combined effects of solar radiation and panel junction temperature on all model-parameters to forecast peak power and photovoltaic yield of solar panel under non-standard conditions.  Solar Energy 191, 341–359.
  * Zhang et al. (2020) Zhang, H., Heidari, A.A., Wang, M., Zhang, L., Chen, H., Li, C., 2020.  Orthogonal Nelder-Mead moth flame method for parameters identification of photovoltaic modules.  Energy Conversion and Management 211, 112764.