*Article* **A Fast-Tracking Hybrid MPPT Based on Surface-Based Polynomial Fitting and P&O Methods for Solar PV under Partial Shaded Conditions**

**Catalina González-Castaño 1,† , Carlos Restrepo 2,\*,† , Javier Revelo-Fuelagán 3,† , Leandro L. Lorente-Leyva <sup>4</sup> and Diego H. Peluffo-Ordóñez 5,6,†**


**Abstract:** The efficiency of photovoltaic (PV) systems depends directly on solar irradiation, so drastic variations in solar exposure will undoubtedly move its maximum power point (MPP). Furthermore, the presence of partial shading conditions (PSCs) generates local maximum power points (LMPPs) and one global maximum power point (GMPP) in the P-V characteristic curve. Therefore, a proper maximum power point tracking (MPPT) technique is crucial to increase PV system efficiency. There are classical, intelligent, optimal, and hybrid MPPT techniques; this paper presents a novel hybrid MPPT technique that combines Surface-Based Polynomial Fitting (SPF) and Perturbation and Observation (P&O) for solar PV generation under PSCs. The development of the experimental PV system has two stages: (i) Modeling the PV array with the DC-DC boost converter using a real-time and high-speed simulator (PLECS RT Box), (ii) and implementing the proposed GMPPT algorithm with the double-loop controller of the DC-DC boost converter in a commercial low-priced digital signal controller (DSC). According to the simulation and the experimental results, the suggested hybrid algorithm is effective at tracking the GMPP under both uniform and nonuniform irradiance conditions in six scenarios: (i) system start-up, (ii) uniform irradiance variations, (iii) sharp change of the (PSCs), (iv) multiple peaks in the P-V characteristic, (v) dark cloud passing, and (vi) light cloud passing. Finally, the experimental results—through the standard errors and the mean power tracked and tracking factor scores—proved that the proposed hybrid SPF-P&O MPPT technique reaches the convergence to GMPP faster than benchmark approaches when dealing with PSCs.

**Keywords:** maximum power point tracking; photovoltaic system; partial shading conditions; surfacebased polynomial fitting

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

#### **1. Introduction**

Global energy demand has increased substantially in recent decades, and along with a greater need for electric power comes the expansion of alternative energy solutions that promote environmental protection and sustainability. For example, the solar energy industry has rapidly grown resulting in decreased manufacturing costs and the proliferation of affordable photovoltaic (PV) systems [1]. However, PV systems face significant challenges, mainly the irregular solar irradiation due to partial shading that affects their efficiency.

PV systems are prone to fluctuations in their efficiency related to the operating environment, so they need to work at their maximum power point (MPP) all the time regardless

**Citation:** González-Castaño, C.; Restrepo, C.; Revelo-Fuelagán, J.; Lorente-Leyva, L.L.; Peluffo-Ordóñez, D.H. A Fast-Tracking Hybrid MPPT Based on Surface-Based Polynomial Fitting and P&O Methods for Solar PV under Partial Shaded Conditions. *Mathematics* **2021**, *9*, 2732. https:// doi.org/10.3390/math9212732

Academic Editor: Nicu Bizon

Received: 29 September 2021 Accepted: 26 October 2021 Published: 28 October 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

of atmospheric conditions. PV-generated power varies significantly due to rapid changes in irradiance caused by the shadows of passing clouds. The PV curve presents multiple peaks, several local maximum power points (LMPPs) and a global maximum power point (GMPP), under these partial shading conditions (PSCs).

Several works use the maximum power point tracking (MPPT) technique in their practical implementation to optimize PV output under PSCs [2–7]. In this vein, an experimental study about MPP characteristics of partially shaded strings is presented in [2]. Based on over 26,000 measured current-voltage curves, this work tests six and seventeen series-connected PV modules. The results proved that following the MPP closest to the nominal MPP voltage can significantly reduce the wide operating voltage range at the expense of small energy losses. To circumvent these shortcomings, a novel MPPT control for PV systems based on the search and rescue (SRA) optimization algorithm is presented in [3]. The improvements exhibited by the proposed technique are enhanced PV system performance, very low oscillations at global maximum (GM) and quick and efficient tracking of GM). Additionally, robustness, up to 99.93% power tracking efficiency in steady-state and implementation simplicity are other outstanding features of the proposed SRA control strategy.

In [8], a new framework based on a sliding-mode controller (SMC) is developed for the MPPT algorithm. This technique delivers precise tracking under changing weather conditions and performs better than conventional methods. Similarly, a new MPPT approach based on the adaptive fuzzy logic controller (AFLC) is introduced in [9]. In the AFLC method, with the goal to produce the optimal duty cycle for MPPT, the membership functions (MFs) are optimized through Grey Wolf Optimization (GWO). Testing under four shading patterns proves that the AFLC approach can track the global MPP for all conditions with enhanced speed, efficiency, and reduced oscillations.

As seen in [10,11], bio-inspired optimization has also been a source for developers to draw inspiration for designing MPPT algorithms. Improving the squirrel search algorithm (ISSA), Ref. [10] introduces a novel MPPT technique that reduces tracking time by half compared to the conventional SSA algorithm. Moreover, the results showed faster convergence and fewer power oscillations when tracking the GMPP.

An alternative method to the classical Marine Predator Algorithm (MPA) is proposed in [11] to cope with its implicit weaknesses. The MPAOBL-GWO method integrates the Opposition Based Learning (OBL) strategy with the Gray Wolf Optimizer (GWO), hence the name. This combination enhances the efficiency of the MPA and prevents it from descending into local points, as can be observed by the results.

Research similar to that of this study can also be found in the literature. For example, an MPPT optimal design based on a surface-based polynomial fit (MPPT-SPF) for a PV system is presented in [12]. The hardware-in-the-loop system is implemented using a high-speed, real-time simulator (PLECS RT Box 1) and a digital signal controller (DSC). In addition, this work applies an optimized version of the SPF technique for partial shading conditions. Likewise, a novel two-stage MPPT method is presented in [13]. In the first stage, the presence of (PSCs) is detected, and then, in the second stage, the global maximum power point (MPP) is reached using a new algorithm based on ramp change of the duty cycle and continuous sampling from the P-V characteristic of the array. Finally, the "Perturb and Observe" algorithm traces small changes of the new MPP.

The comprehensive review of online, offline, and hybrid optimization MPPT algorithms conducted in [14,15] indicates that most conventional MPPT algorithms track the GMPP correctly under conditions of uniform solar irradiance. Otherwise, under conditions of partial shading, rapidly switching, and conditions of nonuniform irradiance, it fails to obtain an accurate GMPP. However, under conditions of rapid solar irradiance change and PSCs, hybrid optimization algorithms are quick and precise in GMPP tracking. Unfortunately, these models are complex and, therefore, difficult to implement using integrated technologies. Therefore, this paper proposes a fast-tracking hybrid MPPT based on Surface-Based Polynomial Fitting and P&O for solar PV under PSCs.

This paper introduces the SPF-P&O GMPPT hybrid approach that exploit simultaneously the benefit of the classical P&O approach and a more data-driven curve-fitting method, as its name implies. This synergistic and complementary strategy improves the P&O algorithm efficiency using a polynomial approximation for global data fitting; the algorithm outputs the polynomial coefficients that best fit the PV panel characteristic curves. Specifically, the proposed approach is a curve-based type, which is considered to be data-driven since its parameters are optimized over the structure of the data. Following from this premise, curve-based fitting approaches have been successfully applied in different application, such as: variation of cosmic ray intensity with atmospheric pressure using a straight line fitting, calibration of a prism spectrometer using a polynomial curve, variation of viscosity of water with temperature using polynomial with equally spaced observations, variation of vapour pressure of ethyl alcohol with temperature using a generalized linear function, and the counting rate of a type I counter using non-linear functions, among others [16]. Finally, as for the MPPT stage, the resulting coefficients are used as the basis aimed at achieving a more accurate estimation and simplifying digital implementations into low-cost digital signal controllers.

The experimental results that validate the GMPPT algorithm come from a simulation in six testing scenarios and different transients from the 133 different cases available in the P-V characteristic data (shown in Figure 1). The analyzed case studies are: (i) system startup, (ii) uniform irradiance variations, (iii) sharp change of the (PSCs), (iv) multiple peaks in the P-V characteristic, (v) dark cloud passing, and (vi) light cloud passing. The results for the six scenarios are evaluated using the standard errors and the mean power tracked and tracking factor scores.

The main contributions of this paper are as follows:


The rest of this paper is structured as follows. First, Section 2 shows the PV system and its controllers. Then, the proposed SPF-P&O GMPPT algorithm is explained in Section 3. Next, the numerical simulations and Real-Time HIL Results are detailed in Section 4. Lastly, conclusions and future work are presented in Section 5.

**Figure 1.** P-V characteristics of four PV module KC200GT under uniform and nonuniform irradiance conditions and connected in series.

#### **2. PV Modeling**

Through a setup of four PV modules, with its respective parallel bypass diodes and all connected in series, the essence of PV string characteristics is studied as shown in Figure 1. This PV string employs a DC-DC boost switching converter to charge a battery while a nested control loop tracks the MPP of the PV system.

Figure 1 depicts the MPPT algorithm and the double-loop controls. Using PLECS, the PV units are modeled as an array of four series-connected KC200GT solar modules whose electrical characteristics are presented in Table 1. Likewise, Figure 1 details the nonlinear P-V characteristic for 133 different cases of shading patterns. There are 24 possible permutations for each of the 133 cases, totaling 3192 possible PSCs. Thus, the PV system can be studied under uniform irradiance levels (cases 1 to 10) with a unique peak in the P-V

characteristic corresponding to the GMPP. On the other hand, different combinations of PSCs generate multiple peaks in the P-V characteristic, from two to three. Figure 1 shows the nonuniform irradiance levels on the PV modules that produce these PSCs. The GMPP moves with respect to the LMPPs at the multiple P-V characteristic peaks to test different situations with the MPPT algorithms.

The following are the differential equations for the boost converter [17–19]:

$$\frac{di\_L(t)}{dt} = \frac{v\_\mathcal{g} - (1 - \mu)\,v\_o}{L} \tag{1}$$

$$\frac{dv\_o(t)}{dt} = \frac{-v\_o}{R\_L \mathcal{C}} + \frac{(1-u)i\_L}{\mathcal{C}},\tag{2}$$

where *vo* is the output voltage, *u* represents the control variable ∈ {0, 1} and *iL* is the inductor current. For the boost converter, the duty cycle is [17–19]:

$$
\mathfrak{u} = 1 - \frac{\mathfrak{v}\_{\mathfrak{g}}}{\mathfrak{v}\_{\mathfrak{o}}}.\tag{3}
$$

**Table 1.** Electrical characteristics of PV module KC200GT.


#### *2.1. Discrete-Time Sliding-Mode Current Control*

Designing the inner-loop control for the DC-DC boost converter is complex due to the inherent non-linearity of the converter. Thus, a robust sliding-mode controller is used in this study. This section outlines the discrete-time sliding-mode current control (DSMCC) strategy with a fixed frequency. In order to ensure that the control surface (4) is reached in the next sampling period (*fsamp* = *fs*), this approach computes the variable control *u*[*n*] in the *n*-th time sample period. This control has been implemented for switching systems in [20–22].

$$s[n] = i\_{Lref}[n-1] - i\_L[n].\tag{4}$$

Equation (2) is used to calculate the inductor current slopes of the boost converter, presented in Table 2. Assuming the averaged model of the inductor current slope of the converter *diL dt* <sup>≈</sup> *iL*[*n*+1]−*iL*[*n*] *<sup>T</sup>* , the Euler approximation heads to the next discrete-time inductor current expression:

$$i\_L[n+1] = i\_L[n] + T(m\_1 + m\_2)u[n] - m\_2T \tag{5}$$

being *T* the sampling or switching period. Consequently, the resulting duty cycle expression is:

$$
\mu[n] = \frac{1}{(m\_1 + m\_2)T} e[n] + \frac{m\_2}{m\_1 + m\_2} \,\prime \tag{6}
$$

where *e*[*n*] = *iLre f* [*n*] − *iL*[*n*], being *iLre f* [*n*] = *iL*[*n* + 1] (see Figure 1). Using the expressions for *m*<sup>1</sup> and −*m*<sup>2</sup> for the output current slopes from Table 2 in (6), the boost converter's control law is given by:

$$u[n] = \frac{L}{v\_o[n]T}e[n] + 1 - \frac{v\_\mathcal{g}[n]}{v\_o[n]}.\tag{7}$$

**Table 2.** Slope of the inductor current waveform.


#### *2.2. Discrete-Time PI Voltage Control*

The input voltage of the boost converter *vg* is regulated for the external loop using a proportional-integrator controller. The forward Euler method is used to express the controller transfer function in the z domain as follows:

$$G\_{vpi}(z) = K\_{pv} + \frac{K\_{iv}T\_{sample}}{z - 1},\tag{8}$$

where *Tsamp* = 1/ *fsamp*,

and

$$K\_{pv} = 2\pi \,\mathrm{C}\_{\mathrm{in}} f\_{\mathcal{L}'} \tag{9}$$

$$K\_{iv} = \frac{K\_{pv}}{T\_i} \,\prime \tag{10}$$

with *Cin* representing the input capacitor. For the voltage loop (*fc*), the crossover frequency (CF) value should be lower than the current loop CF. Furthermore, the PI zero should be below *fc* (1/(2π*Ti*) < *fc*).

#### **3. SPF-P&O MPPT**

#### *3.1. Maximum Power Point Tracking (MPPT) Algorithm*

The Maximum Power Point Tracking (MPPT) control keeps the power transfer at the highest efficiency, optimizing the PV system performance in any radiation and temperature conditions that the solar panels undergo. The MPPT approach allows the PV module to function at its maximum power point by controlling the switching converter [23]. This section briefly explains the "Perturb and Observe" benchmark method before introducing the offered MPPT method.

#### *3.2. Conventional "Perturb and Observe" Method*

The so-called "Perturb and Observe" (P&O) method is widely used due to its simplicity and low cost [24–26]. This algorithm provokes perturbations by either decreasing or increasing the reference voltage according to the output power PV module. If the current measured power *P*[*n*] is higher than its previous sampled value *P*[*n* − 1], the voltage change continues in the same direction. Otherwise, it is reversed. Next, the PV module voltage is compared to the maximum voltage, to predict the MPP. Finally, a power step of the PV module [24] is produced through a small step of reference voltage. The P&O-based MPPT is abbreviated as MPPT-P&O.

#### *3.3. Proposed MPPT Method*

This work proposes a polynomial curve-fitting approach in favor of a more accurate MPP estimation. This approach, known as Surface-based Polynomial Fitting (MPPT-SPF), works as explained below:

#### 3.3.1. Curve-Based Fitting

The fundamentals of a polynomial model (*y* = *f*(*x*)) for any curve can be expressed in mathematical terms as follows:

$$\{y\{\rho\_N\} = \sum\_{i=1}^{N+1} \mathbb{C}\_i \mathfrak{x}^{N+1-i}\} \tag{11}$$

where *x* is the input times series, *y*{*ρN*} is the output time series, *n* is the degree of the polynomial, such that 1 ≤ *N* ≤ 9, and *N* + 1 is the order of the polynomial. In this case, the order is the number of coefficients to be adjusted, while the degree represents the highest power of the predictor variable.

Below, and based on their degree, polynomials are presented. For example, a fourdegree polynomial is given by:

$$y\{\rho\_4\} = f(\mathbf{x}) = \mathbb{C}\_1 \mathbf{x}^4 + \mathbb{C}\_2 \mathbf{x}^3 + \mathbb{C}\_3 \mathbf{x}^2 + \mathbb{C}\_4 \mathbf{x} + \mathbb{C}\_5. \tag{12}$$

Polynomials are helpful when a simple empirical model is needed. Hence, a polynomial model is suitable for interpolation and extrapolation processes or characterizing data using a global fit.

Polynomial fitting is reasonably flexible when dealing with simple data structures. However, fitting can become unstable for high-degree polynomials. Moreover, although all polynomials fit correctly within a predefined data range, they diverge significantly outside of it.

Curve-fitting with high-degree polynomials can result in scale affectations because it uses predictor values as the basis for a matrix with high values. This problem can be addressed by preprocessing input data using z-score normalization, i.e., centering to zero mean and scaling to unit standard deviation [16].

#### 3.3.2. Surface-Based Fitting

The output time series is represented as *z* = *f*(*x*, *y*) when the fitting *f*(·) involves two input time series. For MPPT purposes, the variables are defined as follows:


The following notation is used for polynomial surfaces: *ρij* is the fitting type, where *j* represents the degree of *y* and, on the other hand, *i* the degree of *x*. Additionally, the maximum value for *i* and for *j* [27] is 5. The maximum between *i* and *j* is the overall degree of the polynomial. The degree of *x* is going to be less than or equal to *i* in each term. Likewise, in each term, the degree of *y* is going to be less than or equal to *j*. Accordingly, a surface with *i* and *j* degrees is denoted as *z*{*ρij*} = *f*(*x*, *y*). Some examples are mentioned in Table 3.

**Table 3.** Examples of polynomial models for surfaces.


Table 4 shows the degrees that make up the model terms. For example, for an *x* degree of 1 and a *y* degree of 3, the name of the model will be *ρ*13. The mathematical foundation of the numerical curve-fitting methods is further detailed in [16].


**Table 4.** Polynomial model terms.

A function is obtained from the *ρij* approach that precisely fits the behavior and general trend of the analyzed data. Next, the most accurate fitting is found, considering both the criteria to quantify the adjusting procedure suitability and the polynomial to be tuned to represent the input data. The relationship between the obtained polynomial degree, the curve adjustment, and the values to interpolate is obtained from this adjustment. Finally, a five-degree polynomial is generated from a given data sequence in the form (*x*[*n*], *y*[*n*]), as shown below:

$$\begin{aligned} f(\mathbf{x}, \mathbf{y}) &= \begin{aligned} \mathbb{C}\_{00} + \mathbb{C}\_{10}\mathbf{x} + \mathbb{C}\_{01}\mathbf{y} + \mathbb{C}\_{20}\mathbf{x}^2 + \mathbb{C}\_{11}\mathbf{x}\mathbf{y} \\ &+ \mathbb{C}\_{02}\mathbf{y}^2 + \mathbb{C}\_{30}\mathbf{x}^3 + \mathbb{C}\_{21}\mathbf{x}^2\mathbf{y} + \mathbb{C}\_{12}\mathbf{x}\mathbf{y}^2 + \mathbb{C}\_{03}\mathbf{y}^3 \\ &+ \mathbb{C}\_{40}\mathbf{x}^4 + \mathbb{C}\_{31}\mathbf{x}^3\mathbf{y} + \mathbb{C}\_{22}\mathbf{x}^2\mathbf{y}^2 + \mathbb{C}\_{13}\mathbf{x}\mathbf{y}^3 + \mathbb{C}\_{04}\mathbf{y}^4 \\ &+ \mathbb{C}\_{50}\mathbf{x}^5 + \mathbb{C}\_{41}\mathbf{x}^4\mathbf{y} + \mathbb{C}\_{32}\mathbf{x}^3\mathbf{y}^2 + \mathbb{C}\_{23}\mathbf{x}^2\mathbf{y}^3 + \mathbb{C}\_{14}\mathbf{x}\mathbf{y}^4. \end{aligned} \tag{13}$$

The afore-developed analysis allows assessing how well the curve fits the data. Goodness-of-fit is measured by determining the accuracy of the coefficients based on the 95% confidence limits. The polynomial coefficients determined from the robust fitting of the data are:

*C*<sup>00</sup> = 186.8(184.7, 188.9), *C*<sup>10</sup> = −14.14(−14.55, −13.73), *C*<sup>01</sup> = 6.969(6.883, 7.056), *C*<sup>20</sup> = 0.3872(0.3721, 0.4022), *C*<sup>11</sup> = −0.3113(−0.316, −0.3065), *C*<sup>02</sup> = 0.009194(0.008359, 0.01003), *C*<sup>30</sup> = −0.004332(−0.004533, −0.00413), *C*<sup>21</sup> = 0.006992(0.006861, 0.007124), *C*<sup>12</sup> = −0.0004896(−0.0005158, −0.0004635), *C*<sup>03</sup> = 7.545e − 06(4.966e − 06, 1.012e − 05), *C*<sup>40</sup> = 2.297e − 05(2.184e − 05, 2.409e − 05), *C*<sup>31</sup> = −7.728e − 05(−7.875e − 05, −7.581e − 05), *C*<sup>22</sup> = 1.026e − 05(9.913e − 06, 1.061e − 05), *C*<sup>13</sup> = −6.762e − 07(−7.171e − 07, −6.353e − 07), *C*<sup>04</sup> = 2.86e − 08(2.565e − 08, 3.155e − 08), *C*<sup>50</sup> = −4.633e − 08(−4.858e − 08, −4.407e − 08), *C*<sup>41</sup> = 3.104e − 07(3.049e − 07, 3.158e − 07), *C*<sup>32</sup> = −5.768e − 08(−5.93e − 08, −5.607e − 08), *C*<sup>23</sup> = 5.3e − 09(5.014e − 09, 5.586e − 09), and *C*<sup>14</sup> = −2.523e − 10(−2.799e − 10, −2.248e − 10).

generating the following polynomial:

$$\begin{aligned} f(x,y) &= -186.8 - 14.14x + 6.969y + 0.3872x^2 - 0.3113xy \\ &+ 0.009194y^2 - 0.004332x^3 + 0.006992x^2y - 0.0004896xy^2 \\ &+ 7.545e - 06y^3 + 2.297e - 05x^4 - 7.728e - 05x^3y \\ &+ 1.026e - 05x^2y^2 - 6.762e - 07xy^3 + 2.86e - 08y^4 \\ &- 4.633e - 08x^5 + 3.104e - 07x^4y - 5.768e - 08x^3y^2 \\ &+ 5.3e - 09x^2y^3 - 2.523e - 10xy^4. \end{aligned}$$

As a result, the obtained fitting reaches a remarkable 0.9507 R-square, indicating a significant data trend and a good model fitting. Moreover, a 38.74 root mean square error (RMSE) and a 9.477e + 07 sum of square error (SSE) estimation are reached. Figure 2 depicts the previous adjustment as a MPP surface in terms of *Ppv* and *vg*.

**Figure 2.** Plotting of the MPP surface in terms of *Ppv* and *vg*.

Validation and implementation of the MPPT-SPF aim to maximize the available energy of the connected solar modules at any given moment during operation. Hence, MPPT continually samples the output of the PV cells and adjusts the voltage and current so that the PV system generates maximum power regardless of environmental conditions.

Additionally, by following the approach described in [28], uncertainty values for all the coefficients can be calculated. To do so, consider the vector **C** holding the estimates of the coefficients given by:

$$\tilde{\mathbf{C}} = \left(\boldsymbol{\Phi}^{\top}\boldsymbol{\Phi}\right)^{-1}\boldsymbol{\Phi}^{\top}\mathbf{E},\tag{15}$$

where

$$\Phi = \left[ \mathbf{1}\_{\mathfrak{M}}, \mathbf{x}, \mathbf{y}, \mathbf{x}^2, \mathbf{x}\mathbf{y}, \mathbf{y}^2, \mathbf{x}^3, \mathbf{x}^2\mathbf{y}, \mathbf{x}\mathbf{y}^2, \mathbf{y}^3, \mathbf{x}^4, \mathbf{x}^3\mathbf{y}, \mathbf{x}^2\mathbf{y}^2, \mathbf{y}^3, \mathbf{y}^4, \mathbf{x}^5\mathbf{y}, \mathbf{x}^3\mathbf{y}^2, \mathbf{x}^2\mathbf{y}^3, \mathbf{x}\mathbf{y}^4 \right], \tag{16}$$

and

$$\mathbf{E} = [\varepsilon\_1, \dots, \varepsilon\_m]\_\prime \tag{17}$$

with *m* denoting the number of samples, *varepsiloni* being the error of approximation and *i* ∈ {1, . . . , *m*}. Then, the uncertainty *u*(*Cj*) for the coefficient *j* can be estimated as:

$$
\mu(C\_j) = \sqrt{\frac{\left(\boldsymbol{\Phi}\tilde{\mathbf{C}} - \mathbf{E}\right)^\top \left(\boldsymbol{\Phi}\tilde{\mathbf{C}} - \mathbf{E}\right)}{m - 3}} \sqrt{\theta\_{j\mid j}}\tag{18}
$$

where

$$
\Theta = (\Phi^\top \Phi)^{-1} \,, \tag{19}
$$

*θjj* is the entry *jj* of matrix **Θ** and *j* ∈ {1, ... , 20}. In this case, all the uncertainty values are admissible as can be noted in the following:

*<sup>u</sup>*(*C*1) = *<sup>u</sup>*(*C*00) = <sup>−</sup>8.92 <sup>×</sup> <sup>10</sup>11, *<sup>u</sup>*(*C*2) = *<sup>u</sup>*(*C*10) = <sup>−</sup>4.09 <sup>×</sup> <sup>10</sup>9, *<sup>u</sup>*(*C*3) = *<sup>u</sup>*(*C*01) = 2.49 <sup>×</sup> 105, *<sup>u</sup>*(*C*4) = *<sup>u</sup>*(*C*20) = <sup>−</sup>2.79 <sup>×</sup> <sup>10</sup>6, *<sup>u</sup>*(*C*5) = *<sup>u</sup>*(*C*11) = 7.36 <sup>×</sup> 103, *<sup>u</sup>*(*C*6) = *<sup>u</sup>*(*C*02) = 3.61 <sup>×</sup> 101, *<sup>u</sup>*(*C*7) = *<sup>u</sup>*(*C*30) = <sup>−</sup>4.53 <sup>×</sup> <sup>10</sup>2, *u*(*C*8) = *u*(*C*21) = 4.33, *<sup>u</sup>*(*C*9) = *<sup>u</sup>*(*C*12) = 2.5 <sup>×</sup> <sup>10</sup><sup>−</sup>2, *<sup>u</sup>*(*C*10) = *<sup>u</sup>*(*C*03) = 1.09 <sup>×</sup> <sup>10</sup><sup>−</sup>4, *<sup>u</sup>*(*C*11) = *<sup>u</sup>*(*C*40) = <sup>−</sup>1.75 <sup>×</sup> <sup>10</sup><sup>−</sup>2, *<sup>u</sup>*(*C*12) = *<sup>u</sup>*(*C*31) = 3.77 <sup>×</sup> <sup>10</sup><sup>−</sup>4, *<sup>u</sup>*(*C*13) = *<sup>u</sup>*(*C*22) = 2.13 <sup>×</sup> <sup>10</sup><sup>−</sup>06, *<sup>u</sup>*(*C*14) = *<sup>u</sup>*(*C*13) = 5.14 <sup>×</sup> <sup>10</sup><sup>−</sup>08, *<sup>u</sup>*(*C*15) = *<sup>u</sup>*(*C*04) = 4.55 <sup>×</sup> <sup>10</sup><sup>−</sup>12, *<sup>u</sup>*(*C*16) = *<sup>u</sup>*(*C*50) = <sup>−</sup>1.05 <sup>×</sup> <sup>10</sup><sup>−</sup>07, *<sup>u</sup>*(*C*17) = *<sup>u</sup>*(*C*41) = 4.06 <sup>×</sup> <sup>10</sup><sup>−</sup>09, *<sup>u</sup>*(*C*18) = *<sup>u</sup>*(*C*32) = 2.55 <sup>×</sup> <sup>10</sup><sup>−</sup>11, *<sup>u</sup>*(*C*19) = *<sup>u</sup>*(*C*23) = 1.52 <sup>×</sup> <sup>10</sup><sup>−</sup>12, and *<sup>u</sup>*(*C*20) = *<sup>u</sup>*(*C*14) = 5.12 <sup>×</sup> <sup>10</sup><sup>−</sup>16.

#### *3.4. SPF-P&O Algorithm*

The SPF-P&O GMPPT algorithm is detailed in Algorithm 1, which works as follows: Its objective is to obtain the input voltage reference (*vgre f*) for the boost converter from the measurement of the output voltage and output current of the PV module array. In this way, the maximum power characteristic (*Pmax*) associated with the measured of current and voltage is obtaining from the expression (14). Once *Pmax* is estimated, the voltage reference for the P&O is selected by searching from a register of possible solutions for power point maximum power (*Pmpp*).

When power changes, greater than or equal to reference Δ*Ppv* [%] occur, the MPPT-SPF algorithm is executed, so:

$$\frac{\left|P\_{pvnew} - P\_{pvlast}\right|}{P\_{pvlast}} \ge \Delta P\_{pv} \text{ [\%]},\tag{20}$$

where *Ppvnew* is the actual power measurement and *Ppvlast* is the previous power measurement.

The GMPP search will again execute if the condition (20) is satisfied. This condition ensures that the expression (14) is evaluated to detect the optimal voltage solution for the GMPP even if there is any change in solar irradiance.

**Algorithm 1:** SPF-P&O GMPPT algorithm running at the microcontroller (see Figure 1).

```
Input: vg[n], iL[n]
Output: vgre f [n]
Main Function Main:
   Calculate Ppvnew[n] = vg[n] · iL[n] if Inequation (20) = true then
       Calculate Pmax for the solution vgre f by Equation (14)
        gop = ∞
        while s <number of solutions do
           g(s) = abs(Pmpp(s) − Pmax);
           if (g(s) < gop ) then
               gop = g(s);
                sop = s;
           end
       end
       vgre f = vmpp(sop);
   else
       Calculate Perror = Ppvnew[n] − Ppvlast[n] Calculate verror = vg[n] − vglast[n] vgre f [n]= P&O (Perror, verror, vgre f [n])
   end
   Update Ppvlast[n + 1] = Ppvnew[n] Update vglast[n + 1] = vg[n]
return
Subfunction P&O(Perror, verror, vref )
   // Dynamic step size
   if abs (Perror) > ΔP then
       Δv = Δvbig
   else
       Δv = Δvsmall
   end
   if (Perror > 0) then
       if (verror > 0) then
           vref = vref + Δv
       else
           vref = vref - Δv
       end
   else
       if (verror > 0) then
           vref = vref - Δv
       else
           vref = vref + Δv
       end
   end
   return vref
```
#### **4. Simulation and Real-Time HIL Results**

In this section, the SPF-P&O MPPT algorithm is validated using a DC-DC boost converter through hardware in the loop (HIL). PLECS RT Box 1 implements the stage power converter and the PV array, and 6.6 μs is the sampled time to model the converter. And the values of the boost converter components are: *Cin* = 200 μF, *L* = 1 mH, *Vo* = 160 V and *fs* = 25 kHz. Through the TI 28069M LaunchPad, which is a low-cost Texas Instrument microcontroller, different controls that integrate the PV global system control scheme are implemented as presented Figure 1. Figure 3 shows the setup for HIL experiments.

**Figure 3.** Hardware in-the-loop experimental setup: (a) oscilloscope, (b) PLECS RT box, (c) RT box launchPad interface, (d) Texas Instruments LAUNCHXL-F28069M, (e) Laptop.

The proposed SPF-P&O MPPT method is compared with a GMPPT P&O algorithm studied in [13], which realizes a voltage swept in partial shading conditions across the entire voltage range (from 0 to the open voltage value) for the PV system. The voltage reference that produces the maximum power is used to initialize after the voltage swept the P&O algorithm.

#### *4.1. Inner-Loop Current Control Results*

Figure 4 illustrates the transient responses for the inductor current of the boost converter in reaction to variations of the current reference. Where *vg*, *vo* and *iL* are the signals sampled for the control, and 500 μs is the sampling time. The current reference is changed from 4 A to 8 A and back to 4 A, as can be observed in Figure 4a,b. to ensure a boost operation, the output voltage is *Vo* = 160 V and the input voltage is set in 100 V. Smooth transitions during reference changes are found, settling times near 250 μs and no overshoot. As shown, the inductor current is well regulated. The inductor currents followed the change in the current reference perfectly. During the current step reference change, the current control performance is validated.

#### *4.2. Double-Loop Results*

To demonstrate the effectiveness of the external loop control, voltage reference variations from 100 V to 110 V with a step between variations of 5 V for external loop validation, are presented in Figure 4c. The crossover frequency (CF) corresponds to *fc* = 500 Hz, allowing the proportional gain calculation according to (9). Since the location of the PI zero in Equation (10) is lower than *fc* (1/(2π*Ti*) < *f <sup>c</sup>*), *Ti* = 3.18e−3 s was chosen.

**Figure 4.** Experimental (**a**,**b**), responses of the sliding digital current input control when the reference *iref* : (**a**) changes from 4 A to 8 A, and (**b**) from 8 A to 4 A. The converter operates with *Vg* = 100 V and *Vo* = 160 V of input and output voltage, respectively. CH1: *Vg* (60 V/div), CH2: *Vo* (60 V/div), CH3: *iL* (2 A/div) and a time base of 500 μs. (**c**) responses of the double-loop using sliding digital current control when the reference *vref* changes with steps of 5 V between 100 V to 110 V while the output voltage (*Vo* = 160 V) ensures a boost operation. CH1: *vg* (20 V/div), CH2: *Vo* 20 V/div), CH3: *iL* (20 A/div) and a time base of 156 ms.

Every 400 μs, the voltage regulator (*Gvpi*(*z*)) calculates the inductor current reference as shown in Figure 1. And current transitions, caused by the voltage changes, are smooth and the voltage reference is concisely tracked as can be appreciated from Figure 4c.

#### *4.3. GMPPT P&O and Proposed SPF-P&O Method Comparison*

The proposed method is compared with the GMPPT P&O studied in [13], every 100 ms the GMPPT algorithms are executed to provide a new voltage reference for the voltage loop control, as shown in Figure 1 for the GMPPT algorithm block. For the SPF-P&O, Δ*Ppv* [%] from expression (20) is set to 8%. The inner current loop is updating at 25 kHz, the outer voltage loop is calculating at 2.5 kHz, and the MPPT strategy is computing at 10 Hz, these strategies are implemented in the DSC.

These scale time differences between the systems are challenging but at the same time advantageous to implement an MPPT based on a finite-state machine.

For the validation of the proposed algorithm, six scenarios with different transient of the cases presented in Figure 1 are described below.

#### 4.3.1. Scenario 1: System Start-Up

This scenario presents the start-up with the panels' irradiance levels corresponding to case 56 of Figure 1.

Figure 5 shows, corresponding to the maximum power, the transient behavior from zero current until an equilibrium point, where three modules have an irradiance of 1000 W/m<sup>2</sup> and a module of 400 W/m2.

In Figure 5, both GMPPT algorithms reach the steady-state close to 1.3 s, with the proposed MPPT a uniform step voltage reference to tracking the MPP (599.9 W at 79.1 V) while the system starts up.

The proposed SPF-P&O algorithm works at the optimum point, without oscillation, as is shown in the input voltage signal after it has been tracked. Due to the GMPPT P&O realized a swept of voltage, the reference voltage is increased until 120 V to find the GMPPT at 79.1 V.

In Figure 6, a quantitative analysis can be observed of the proposed GMPPT method and GMPPT P&O method, and the results can be seen in Figure 5. These results confirm a similar performance for both methods through the start-up, where the GMPPT P&O algorithm has a higher tracking factor because the mean power tracked value is closer to the global power maximum.

**Figure 5.** Simulated and experimental dynamic behavior of the MPPT algorithms for Scenario 1 and an output voltage *Vo* = 160 V. The proposed MPPT algorithm (**left**) is compared with GMPPT P&O (**right**). CH1: *vg* (50 V/div), CH2: *iL* (10 A/div), CH3: Maximum power (200 W/div), CH4: Measured power (200 W/div) and a time base of 220 ms.

**Figure 6.** Comparative analysis of the MPPT methods for Scenario 1 shown in Figure 5.

4.3.2. Scenario 2: Uniform Irradiance Variations

Scenario 2 studies the results for the MPPT techniques under uniform irradiance variations as shown in Figure 7. The irradiance sequence corresponds to cases 1, 3, 5, 7, and 9 shown in Figure 1.

**Figure 7.** Simulated and experimental dynamic behavior of the MPPT algorithms for Scenario 2 with an output voltage *Vo* = 160 V. The proposed MPPT algorithm (right) is compared with the GMPPT P&O algorithm (right). CH1: *vg* (50 V/div), CH2: *iL* (10 A/div), CH3: Maximum power (200 W/div), CH4: Measured power (200 W/div) and a time base of 960 ms.

The proposed algorithm outputs the optimal voltage reference for the case and tracks the MPP faster while the GMPPT P&O has a slow convergence.

The proposed model facilitates the calculation and estimation of the variable of interest with high levels of accuracy; however, its proper evaluation through the results obtained to ensure its accuracy is crucial. Certainly, this evaluation can be performed by analyzing the errors, such as the mean absolute error (MAE), the relative error (RE) and the root mean square error (RMSE) and evaluate the performance of the results obtained [29,30], whose equations can be written as follows:

$$\text{RE} = \frac{\sum\_{i=1}^{m} (P\_{pvi} - P\_{mpp})}{P\_{mpp}} \mathbf{100\ [\%]},\tag{21}$$

$$\text{MAE} = \frac{\sum\_{i=1}^{m} |P\_{pvi} - P\_{mpp}|}{m}, \text{ and} \tag{22}$$

$$\text{RMSE} = \sqrt{\frac{\sum\_{i=1}^{m} (P\_{pvi} - P\_{mpp})^2}{m}},\tag{23}$$

where *Ppvi* represents the measured power of the PV module, *Pmpp* is the available MPP power of the solar module, and *m* the total number of sampling data. Figure 8 shows the sensitivity of the MPPT algorithms through MAE, RE, and RMSE for the results shown in Figure 7.

**Figure 8.** Comparative analysis of the MPPT methods under different uniform irradiance conditions for Scenario 2 shown in Figure 7.

Standard error values indicate that the performance of the proposed MPPT algorithm has higher effectiveness in tracking the maximum power point. This statistical analysis shows that the SPF-P&O method reaches the lowest value for the error compared to the GMPPT P&O algorithm. The proposed SPF-P&O algorithm has a tracking factor of 98.07%, while for the GMPPT P&O method, the tracking factor is 80.39%.

#### 4.3.3. Scenario 3: Sharp Change of the PSCs

Scenario 3 presents a sharp change between case 103 and case 68, producing high PV power variations from 588 W to 355 W. Figure 9 shows simulated and HIL results of the GMPP tracking performance. The MPPT tracking efficiency for the GMPPT P&O method is 75.65%, while the 90.48% is achieved by the SPF-P&O proposed method. Figure 10 show the measure of the standard error for each GMPPT method, where the SPF-P&O GMPPT method has minimum error. As can be seen in the input voltage of the converter in Figure 9 and by the inductor current, the PV array always operates in an oscillating mode for the GMPPT P&O method.

**Figure 9.** Simulated and experimental dynamic behavior of the MPPT algorithms for Scenario 3 with an output voltage *Vo* = 160 V. The proposed MPPT algorithm (right) is compared with the GMPPT P&O algorithm (right). CH1: *vg* (50 V/div), CH2: *iL* (10 A/div), CH3: Maximum power (200 W/div), CH4: Measured power (200 W/div) and a time base of 450 ms.

**Figure 10.** Comparative analysis of the MPPT methods for Scenario 3 shown in Figure 9.

Therefore, the proposed GMPPT method achieves a superior performance during abrupt irradiation variations than the GMPPT P&O method. Please note that during the change from case 68 to 103, the experimental result of the SPF-P&O algorithm did not estimate an optimal reference voltage. Nonetheless, when P&O operates, GMPPT is achieved, demonstrating the robustness of the SPF-P&O method.

4.3.4. Scenario 4: Multiple Peaks in the P-V Characteristic

Scenario 4 has the sequence of case 3 with only one peak (GMPP), case 52 with two peaks (One GMPP and one LMPP), and case 85 with three peaks (One GMPP and two LMPP) in the P-V characteristic. The results for each MPPT algorithms are shown in Figure 11.

**Figure 11.** Simulated and experimental dynamic behavior of the MPPT algorithms for Scenario 4 with an output voltage *Vo* = 160 V. The proposed MPPT algorithm (right) is compared with the GMPPT P&O algorithm (right). CH1: *vg* (50 V/div), CH2: *iL* (10 A/div), CH3: Maximum power (200 W/div), CH4: Measured power (200 W/div) and a time base of 540 ms.

For case 3, the irradiance is uniform with a unique MPP at 104.28 V, which is tracked faster by the proposed algorithm. When the irradiance changes in case 52, there are two peaks, the proposed algorithm misidentifies the optimum voltage but archived the maximum power.

Finally, when the irradiation is changed in case 85, there are three peaks, where the algorithm has quickly identified and tracked the second peak (55.3 V, 81.3 W) as the GMPP. The SPF-P&O presents an overall GMPPT tracking efficiency of 90.86%, while the GMPPT P&O method presents a tracking efficiency of 79.85%.

The results for the cases' sequence are evaluated through standard errors (21)–(23), and the scores of mean power tracked and tracking factor as shown in Figure 12.

**Figure 12.** Comparative analysis of the MPPT methods for Scenario 4 shown in Figure 11.

From this figure it can inferred that the proposed GMPPT method has high MPP tracking capability regarding to GMPPT P&O method under multiple peaks in the P-V characteristic due to different PSCs.

#### 4.3.5. Scenario 5: Dark Cloud Passing

A dark cloud passes in this scenario, obscuring each one of the PV modules. The sequence of the cases included in this scenario correspond to the cases 7, 17, 26, and 35 as presented in Figure 13.

**Figure 13.** Simulated and experimental dynamic behavior of the MPPT algorithms for Scenario 5 with an output voltage *Vo* = 160 V. The proposed MPPT algorithm (right) is compared with the GMPPT P&O algorithm (right). CH1: *vg* (50 V/div), CH2: *iL* (10 A/div), CH3: Maximum power (200 W/div), CH4: Measured power (200 W/div) and a time base of 780 ms.

Cases 17, 26 and 35 correspond to (PSCs), which has two peaks, and the proposed GMPPT correctly identifies the optimal voltage for the maximum power point for the transition between cases 7, 17, and 26. Nevertheless, the reference voltage value for case 35 is incorrectly identified by the SPF-P&O algorithm.

This error is because the power between the GMPP (77.1 W at 26.1 V) and the LMPP located in the second peak (72.26 W at 100 V) are similar. The proposed MPPT tracked the GMPP much faster than the GMPPT P&O algorithm. Figure 14 shows the standard error and comparison indicators for the results presented in Figure 13.

The SPF-P&O algorithm exhibits smaller error values in comparison to the GMMP P&O method. The SPF-P&O algorithm presents a tracking factor of 93.53% while the GMPPT P&O method a tracking factor of 72.29%.

**Figure 14.** Comparative analysis of the MPPT methods under different nonuniform irradiance conditions for Scenario 5 shown in Figure 13.

4.3.6. Scenario 6: Light Cloud Passing

In this scenario a light cloud passes, partially obscuring one by one the PV modules. The results of scenario 6 has the sequence of cases: 8, 63, 95, and 127, and the results

are presented in in Figure 15. Figure 16 illustrates that the proposed method provides the lowest value error for the results seen in Figure 15.

When the irradiance level is reduced by the transition of the cloud, the SPF-P&O algorithm identifies the new GMPPT power for the different cases without oscillations while the GMPPT P&O method has a big oscillation around the GMPPT, and this is reflected in the current and voltage waveforms.

The SPF-P&O algorithm presents a tracking factor of 97.1% while for the GMPPT P&O method a value of 76.27%.

**Figure 15.** Simulated and experimental dynamic behavior of the MPPT algorithms for Scenario 6 with an output voltage *Vo* = 160 V. The proposed MPPT algorithm (right) is compared with the GMPPT P&O algorithm (right). CH1: *vg* (50 V/div), CH2: *iL* (10 A/div), CH3: Maximum power (200 W/div), CH4: Measured power (200 W/div) and a time base of 780 ms.

**Figure 16.** Comparative analysis of the MPPT methods for Scenario 6 shown in Figure 15.

#### **5. Conclusions**

A fast-tracking hybrid MPPT technique based on Surface-Based Polynomial Fitting and P&O has been presented for solar PV under PSCs. The SPF-P&O MPPT uses a polynomial model from the characterization of the PV module data, which is evaluated during irradiance variations. Meanwhile, the conventional P&O tracks the MPPT under uniform irradiance. The power circuit model was implemented using an RT BOX 1 tool. A low-cost commercial DSC was used to implement the proposed MPPT algorithm and the double-loop strategies in a C programming software.

As experimentally demonstrated in this work, the introduced SPF-P&O GMPPT approach can tie together the estimation capability of the conventional P&O with a curvefitting-based approach (namely surface-based polynomial fitting). Thus, it can be said that the hybrid approach represents a synergistic and complementary strategy to leverage the effectiveness of a P&O-type method with global data fitting using a polynomial approach. Different profile tests proved that the proposed hybrid method is robust and performs a faster and more effective MPP tracking with no steady-state oscillations for partial shading conditions.

In the future, a more extensive comparison between different MPPT techniques under PSCs using a low-cost commercial microcontroller will be conducted to accomplish deeper insights about the efficiency of PV systems, furthering the growth of photovoltaics as an affordable and sustainable energy source.

**Author Contributions:** Conceptualization, D.H.P.-O., C.R, L.L.L.-L., J.R.-F. and C.G.-C.; methodology, D.H.P.-O., C.R., L.L.L.-L., J.R.-F. and C.G.-C.; software, C.R. and C.G.-C.; validation, C.R. and C.G.-C.; formal analysis, D.H.P.-O., L.L.L.-L. and J.R.-F.; investigation, D.H.P.-O., C.R., L.L.L.-L. and C.G.-C.; resources, D.H.P.-O., C.R., L.L.L.-L., J.R.-F. and C.G.-C.; data curation, D.H.P.-O., C.R. and C.G.-C.; writing—original draft preparation, D.H.P.-O., C.R. and C.G.-C.; writing—review and editing, D.H.P.-O., C.R., L.L.L.-L., J.R.-F. and C.G.-C.; visualization, D.H.P.-O., C.R. and C.G.-C.; supervision, D.H.P.-O., C.R. and L.L.L.-L.; project administration, D.H.P.-O., L.L.L.-L. and C.R.; funding acquisition, D.H.P.-O., C.R. and J.R.-F. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work was supported in part by the Chilean Government under projects ANID/ FONDECYT/ 1191680 and SDAS Research Group, www.sdas-group.com (accessed on 29 September 2021.)

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Not applicable.

**Conflicts of Interest:** The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

#### **References**


**Stefan Jovˇci´c \* and Petr Pr ˚uša**

Department of Transport Management, Marketing and Logistics, Faculty of Transport Engineering, University of Pardubice, Studentská 95, 532 10 Pardubice, Czech Republic; petr.prusa@upce.cz **\*** Correspondence: stefan.jovcic@upce.cz; Tel.: +420-466-036-385

**Abstract:** Third-party logistics (3PL) is becoming more and more popular because of globalization, e-commerce development, and increasing customer demand. More and more companies are trying to move away from their own account transportation to third-party accounts. One reason for using 3PLs is that the company can focus more on its core activities, while the 3PL service provider can provide distribution activities in a more professional way, save costs and time, and increase the level of customer satisfaction. An emerging issue for companies in the logistics industry is how they can decide on the 3PL evaluation and selection process for outsourcing activities. For the first time, the entropy and the criteria importance through intercriteria correlation (CRITIC) methods were coupled in order to obtain hybrid criteria weights that are of huge importance to decide on the 3PL provider evaluation and selection process. The obtained criteria weights were further utilized within the additive ratio assessment (ARAS) method to rank the alternatives from the best to the worst. The introduced hybrid–ARAS approach can be highly beneficial, since combining two methods gives more robust solutions on one hand, while on the other hand eliminating subjectivity. Comparative and sensitivity analyses showed the high reliability of the proposed hybrid–ARAS method. A hypothetical case study is presented to illustrate the potentials and applicability of the hybrid–ARAS method. The results showed that 3PL-2 was the best possible solution for our case.

**Keywords:** 3PL logistics; decision making; ARAS; entropy; CRITIC

#### **1. Introduction**

Third-party logistics (3PL) selection is becoming more and more popular because of globalization, e-commerce development, and increasing customer demand. Today's society, through the needs of different entities, represents a source of numerous new requirements and expectations for companies in the postal and logistics industries [1]. Systems for the distribution of goods, both at the national and international levels, are very important for appropriate business functioning and for the normal life of citizens [2]. Because of e-commerce development, there is increasing pressure on 3PL service providers all over the world. Wang et al. [3] emphasized that more trustworthy delivery, high inventory turnover, and inventory staged in forwarding locations near consumers are all effects of the e-commerce trend. According to Wang [4], e-commerce was sped up by the COVID-19 outbreak. According to Wang et al. [5], today's business globalization, customer satisfaction, and strong competition have forced many companies to work closely with external business partners. Third-party logistics service providers have had a significant impact on society on a global scale. They are not responsible not only for moving goods from the point of origin to the point of the destination, but in some cases for packaging, storage, etc. Third-party logistics services depend on the contracts they sign with collaborating companies.

Based on the aforementioned facts, the selection and evaluation process of 3PL service providers, is not an easy task for logisticians, since multiple factors affect the decisionmaking process. A poor choice of business partner can greatly negatively affect a company. The resulting losses can be financial, material, loss of reputation, loss of users, and many

**Citation:** Jovˇci´c, S.; Pr ˚uša, P. A Hybrid MCDM Approach in Third-Party Logistics (3PL) Provider Selection. *Mathematics* **2021**, *9*, 2729. https://doi.org/10.3390/math9212729

Academic Editor: Nicu Bizon

Received: 16 September 2021 Accepted: 26 October 2021 Published: 27 October 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

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

others. Furthermore, Soh [6] stated that if an appropriate 3PL provider was not selected, serious problems could occur, such as low-quality logistics services and contract nonfulfillment. He also emphasized that the decision-making problem for selecting the best 3PL provider had been receiving much attention recently among scholars as well as business practitioners. Nevertheless, Hsu et al. [7] stated that engaging 3PL could reduce fixed costs and increase flexibility, allowing organizations to focus on their core competencies and thereby enhancing their efficiency. Since 3PL evaluation and selection is a multidisciplinary field, there are many research questions that have been addressed by various authors in the field.

In this article, two research questions were dealt with: research question 1 (RQ1) refers to the combination of objective methods that allow the evaluation criteria for the 3PL provider selection process to obtain the best possible result; research question 2 (RQ2) is directed to the 3PL provider selection process.

To answer the aforementioned research questions, knowledge from various multidisciplinary fields was applied. Given the fact that no complete numerical data were available to create a real-life case study, given the time constraints of this research, some hypothetical data were used. The 3PL evaluation and selection process, in this paper, starts with a discussion with the experts in the field of logistics. The experts, according to their knowledge and experience, helped the authors define some of the criteria that should be of vital importance for companies that deal with the decision-making process. The number of alternatives (possible 3PL providers) depends on the company, but in this hypothetical case, there were five possible alternatives. Not all criteria were equally important, so it was necessary to assign them degrees of importance. For this part of the paper, the entropy and the CRITIC methods were coupled, and the new hybrid criteria importance is obtained. To obtain the final rank of the alternatives considered, the ARAS MCDM method is used.

The main contribution of this article lies in the proposed entropy–CRITIC (hybrid)– ARAS methodology, which, according to the authors' knowledge and the reviewed literature in the field, has not previously applied to this problem. The proposed methodology is general and can be applied to any other MCDM problem dealing with the selection of collaboration partners as well as any interrelated criteria. In addition, this paper contributes a combination of two objective methods to obtain hybrid criteria weights, since combining two methods gives more robust solutions on one hand while eliminating subjectivity on the other hand.

In this article, the methodology was applied only to a hypothetical example, which should be emphasized as a limitation. However, the authors intend to test the methodology for the real-life case study. Another weak point of the paper is that the considered criteria for 3PL selection were mostly part of the economic pillar. Nevertheless, future directions include addressing the other aspects such as environmental, social, and technical ones.

Since there is an increasing number of logistics companies that provide various logistics services as third parties, it is not easy to evaluate and select the best collaborative partner. The authors of this paper had the main motivation to propose methodology that on the one hand should be easy to implement and on the other hand should help decision makers select the best 3PL. The methodology upgraded the ARAS method by combining objective methods in order to eliminate subjectivity in the assessment of criterion weights and obtain more robust solutions using the ARAS method.

This paper is organized as follows: After the introductory section, Section 2 provides the current state of the methods used in the 3PL logistics field. The methodology proposed to solve the 3PL evaluation and selection problem is elaborated in Section 3. The application of the hybrid MCDM method to a hypothetical example is elaborated in Section 4. Section 5 gives some managerial insights into the 3PL logistics field. Section 6 is the conclusion and gives some future research directions.

#### **2. Literature Review**

As pointed out above, 3PL service provider evaluation and selection are not easy tasks for decision makers, given the fact that multiple criteria and many existing methods ought to be taken into consideration. Researchers have created many methods to solve the 3PL evaluation and selection problem. Most of these methods have belonged to the class of multicriteria decision-making methods. In addition to multicriteria analysis methods, many other methods, such as statistical or mathematical programming methods and integrated approaches, have been used. This section provides a review of the literature in the field of 3PL service providers based on the methods that other authors used to solve the 3PL evaluation and selection problem. Jovˇci´c et al. [8] provided an extensive review of the literature regarding the most commonly used methods for 3PL provider evaluation and selection.

One of the most often used multicriteria analysis methods is the analytic hierarchy process (AHP). Saaty [9] originally developed this method. After its introduction, the AHP was widely used in many fields to solve multicriteria decision-making problems; one of those fields was logistics. Korpela and Touminen [10] applied the AHP to find the best possible solution of 3PL warehousing in the processing industry. Yahya and Kingsman [11] used the AHP to determine priorities in selecting suppliers. Akarte et al. [12] proposed a web-based AHP system to evaluate casting suppliers. Muralidharan et al. [13] developed a five-step AHP method to rank suppliers. Liu and Hai [14] applied the AHP to evaluate and select suppliers. So et al. [15] used the AHP to assess the quality of service of suppliers in Korea, while Göl and Çatay [16] applied the AHP to select the best 3PL service provider in a Turkish automotive company. Chan et al. [17] considered the supplier selection issue in the airline industry by using the AHP. Hou and Su [18] assessed and selected suppliers in the mass-customization environment by utilizing the AHP. Gomez et al. [19] proposed a model to evaluate the performance of suppliers by using the AHP. Hudymácov ˇ á [20] applied the AHP in supplier selection. Asamoah [21] applied the AHP in a pharmaceutical manufacturing company in Ghana. Hruška et al. [22] solved a 3PL selection problem by AHP in the production company in the Czech Republic. Jayant and Singh [23] applied an AHP–VIKOR hybrid MCDM approach for 3PL selection. Tuljak-Suban and Bajec [24] upgraded the AHP with the graph theory and matrix approach (GTMA). Aguezzoul and Pache [25] combined the AHP with the ELECTRE I methodology to solve the 3PL selection problem.

The analytic network process (ANP) is also a frequently used method for 3PL evaluation and selection problems. Meade and Sarkis [26] proposed a conceptual model to evaluate and select a third-party reverse logistics provider (3PRLP). Sarkis and Talluri [27] applied the ANP to evaluate and select the best supplier, considering seven evaluating criteria. Bayazit [28] applied the ANP to tackle the supplier selection problem. Jkharkharia and Shankar [29] applied the ANP to select the best logistics service provider. Further research regarding 3PRLP evaluation and selection was proposed by Zareinejad and Javanmard [30]. They applied ANP, intuitionistic fuzzy sets (IFS), and grey relation analysis (GRA). In their study, the ANP was used to identify the most important attributes in the selection and evaluation of 3PRLP. The technique for order of preference by similarity to ideal solution (TOPSIS) is one of the most frequently applied approaches in 3PL. Mostly, this method is coupled with fuzzy logic, ANP, AHP, etc. There have been many studies in the literature that may confirm it. Chen and Yang [31] proposed restricted fuzzy-AHP and fuzzy-TOPSIS to assess and choose the best supplier. Zeydan et al. [32] used a combination of fuzzy-AHP, fuzzy-TOPSIS, and DEA methods in the automotive industry to evaluate and select suppliers. Singh et al. [33] applied the TOPSIS method for supplier selection in the automotive industry as well. Jayant et al. [34] evaluated and selected reverse third-party logistics service providers (R3PLs) in the mobile phone industry by coupling the AHP (to evaluate the criteria for R3PLs) and TOPSIS (to select the best one). Laptate [35] used fuzzy modified TOPSIS for supplier selection problems in the supply chain. ELECTRE (ELimination Et Choice Translating REality) is a family of multicriteria

decision analysis methods. Aguezzoul and Pires [36] used the ELECTRE method for 3PL performance evaluation and selection in a complex strategic decision process that involved various qualitative and quantitative criteria. Before that, Govindan et al. [37] used the fuzzy-ELECTRE method to rank 3PRL providers. The method was applied to a battery recycling case. When it comes to the combination of fuzzy logic with the multicriteria decision-making methods, there have been many studies in the scientific literature. According to Cheng [38] and Cheng et al. [39], a fuzzy-AHP approach handles issues that use the theory of fuzzy sets as well as hierarchical structure analysis. On the other hand, Ayhan [40] declared that the fuzzy-AHP approach was an extended AHP model into a fuzzy domain. This method has found application in various fields. For example, Kilincci and Onal [41] applied fuzzy-AHP to select suppliers for a washing machine company. When it comes to a low carbon supply chain, Shaw et al. [42] coupled the fuzzy-objective linear programming (LP) with the fuzzy-AHP in order to choose an optimal supplier solution. First, to determine the weights of the predetermined criteria, the Fuzzy-AHP was used. Second, the best supplier was determined by using fuzzy-objective LP. Zhang et al. [43], Zhang and Feng [44], Göl and Catay [16], and Soh [6] combined the AHP and fuzzy approaches to solve a 3PL service provider assessment issue. To compute the importance of individual parameters and subcriteria in fourth-party logistics (4PL), Cheng et al. [45] utilized the fuzzy-AHP approach. Arikan [46] dealt with the fuzzy-AHP method for multiple-objective supplier selection problems. Jagannath et al. [47] evaluated and selected 3PL providers from a sustainability perspective using the interval-valued fuzzy-rough approach. Rezaeisaray et al. [48] conducted a study on a pipe and fittings manufacturing company using a novel hybrid MCDM model for outsourcing supplier selection. They concluded that among the selective criteria for outsourcing, business development, focus on basic activities, and order delays were the three most important ones. They also ranked suppliers to facilitate decision making for selection. Sremac et al. [49] assessed logistics providers by combining the rough stepwise weight assessment ratio analysis and rough weighted aggregated sum product assessment approaches. Zarbakhshnia et al. [50] proposed a multiple attribute decision-making (MADM) model to rank and select 3PRLPs, using fuzzy stepwise weight assessment ratio analysis (SWARA) to weight the evaluation criteria. To rank and select sustainable 3PRLPs, the COPRAS (complex proportional assessment of alternatives) method was used. Özcan and Ahıskalı [51] solved the 3PL service provider selection problem by combining multicriteria decision-making methods with linear programming models. Some statistical methods that deal with the 3PL supplier selection problem can be found in the literature. For example, the correlation method was used by various authors [52–54]. Lai [55] conducted cluster analysis, which analyzed the service capability and performance of logistics service providers. Sinkovics and Roath [56] used descriptive statistics in 3PL relationships, considering six parameters: customer orientation, competitor orientation, operational flexibility, collaboration, logistics performance, and market performance. Knemeyer and Murphy [57] evaluated the performance of 3PL arrangements from a marketing perspective. Regarding mathematical programming methods (linear and nonlinear programming, dual and multiobjective programming, data envelopment analysis (DEA)), various research papers in the field of logistics service providers can be found. For example, Falsini et al. [58] carried out a study regarding logistic service provider evaluation and selection based on an integration of AHP, DEA, and linear programming methods. Zhou et al. [59] used the DEA method to evaluate the efficiency of Chinese 3PL. Hamdan and Rogers [60] evaluated the efficiency of 3PL operations using the DEA method. Kumar et al. [61] solved a multiobjective 3PL allocation problem for fish distribution. Tsai et al. [62] applied the new fuzzy DEA model to solving MCDM problems in supplier selection. Liu et al. [63] compared suppliers from a collaboration perspective for the new energy vehicle manufacturers in China. Hoseini et al. [3] used the combination of the fuzzy-best-worst method and the fuzzy inference system to solve the supplier selection issue in the construction industry. Kurniawan and Puspitasari [64] evaluated the criteria for supplier selection by applying fuzzy logic and the best-worst method. Whang et al. [65]

used the fuzzy-AHP and fuzzy-VIKOR Methods in sustainable supply chain third-party logistics. For better transparency, the aforementioned literature review on the methods for 3PL evaluation and selection is summarized in Table 1.

**Table 1.** Review based on the methods for 3PL evaluation and selection.


Based on an extensive review of the literature in the field of third-party logistics, the most often used methods were multicriteria decision-making methods in combination with fuzzy logic. At present, in order to support the 3PL evaluation and selection process, there are many new multicriteria methods, such as SWARA, EDAS, MACBAC, and WASPAS, that are based on group decision making in a fuzzy environment. The main advantage of the methods mentioned above is the fact that when coupled with fuzzy logic, they have the power to help decision makers decide in an uncertain environment. In other words, the methods can help decision makers decide when the input data are not defined as crisp values but are given descriptively through linguistic statements. In this study, a hybrid MCDM approach is proposed to evaluate and rank the best 3PL service provider when crisp input data are given. Future research based on this paper should address the hybrid

MCDM approach in regard to the fuzzy domain as well. No previous research has been conducted that applies the hybrid–ARAS method in the way that is proposed here. One of the advantages of the method we propose is that it combines two objective methods to obtain more robust solutions and at the same time eliminate subjectivity.

#### **3. Methodology**

This paper combined three possible methods to solve a hypothetical example of the 3PL selection problem. The data for 3PL providers is usually, as a rule, privately owned. Moreover, some data are not freely available to the general public or the scientific community, probably because of corporate policy to protect proprietary information. Furthermore, the COVID-19 crisis additionally hampered obtaining data. However, the input data for 3PL service providers were formulated based on interviews with experts. The methodology proposed is presented in Figure 1. The first phase was problem formulation. In this case, the 3PL service provider selection problem was considered. From the extensive literature review as well as the Experts' opinions, the authors of this paper identified the criteria for 3PL provider selection. The criteria that were taken into consideration were mostly part of the economic pillar. The CRITIC (criteria importance through intercriteria correlation) and entropy methods were used to find the criteria importance for the 3PL service provider selection. By combining those two methods, the hybrid criteria weights were determined. The main reason for coupling the CRITIC and entropy methods was that both are used to obtain objective criteria weights, and when they are coupled together, a more robust solution is obtained (the subjectivity is eliminated). To find the final rank of the best possible 3PL provider, the criteria importance obtained by the hybrid method were further used within the ARAS method. The ARAS method is a relatively new MCDM method used to rank alternatives. and, according to the authors' knowledge and the literature review, it has not been previously applied and coupled with the hybrid criteria weights as it was here.

**Figure 1.** A flowchart of the methodology used for 3PL service provider evaluation and selection (Source: Authors).

In the second part of the paper, the hybrid criteria importance was integrated into the ARAS method to obtain the best 3PL solution.

#### *3.1. CRITIC (Criteria Importance through Intercriteria Correlation) Method*

In a decision-making process, the importance of criteria plays an important role, since not all criteria are equally important [66]. In this article, we applied the CRITIC method, since it is very useful in obtaining objective criteria importance. The CRITIC method was originally developed by Diakoulaki et al. [67]. The criteria importance calculated by the CRITIC method considers both conflicts among criteria and contrast intensity of each criterion [66]. Ghorabaee et al. [68] declared that in this method, the correlation coefficient is used to observe the conflict between criteria, while the standard deviation is used to consider the contrast intensity of each criterion. According to Diakoulaki et al. [67], the CRITIC method should be described as follows:

Step 1 computes the transformations of performance values (*xij*) and obtains criteria vectors. It is calculated by following Equation (1):

$$\mathbf{x}\_{ij}^{-T} = \begin{cases} \frac{\mathbf{x}\_{ij} - \mathbf{x}\_{j}^{-}}{\mathbf{x}\_{j}^{\*} - \mathbf{x}\_{j}^{-}} \text{ if } j \in B; \\\frac{\mathbf{x}\_{j}^{-} - \mathbf{x}\_{ij}}{\mathbf{x}\_{j}^{-} - \mathbf{x}\_{j}^{\*}} \text{ if } j \in N; \end{cases} \tag{1}$$

where *xij <sup>T</sup>* explains the transformed value, *xj* represents the vector of *j*th criterion, and *x*<sup>∗</sup> *j* and *x*− *<sup>j</sup>* represent the ideal and anti-ideal values with respect to *j*th criterion. If *j B* then *x*∗ *<sup>j</sup>* = *maxixij* and *x*<sup>−</sup> *<sup>j</sup>* = *minixij*. If *j N* then *x*<sup>∗</sup> *<sup>j</sup>* = *minixij* and *x*<sup>−</sup> *<sup>j</sup>* = *maxixij*.

Step 2 calculates the standard deviation Ჿ*<sup>j</sup>* of each criterion utilizing the corresponding vector. Step 3 defines a square (*m* × *m*) matrix *R* with *rjk* elements, where *k* = 1, 2, . . . , *m*:

$$R = \left[ r\_{jk} \right]\_{m \times m} \tag{2}$$

The elements of this matrix are the linear correlation coefficients between the *xj* and *xk* vectors.

Step 4 computes the information measure of each criterion by using Equation (3):

$$H\_{\rangle} = \mathfrak{G}\_{\rangle} \sum\_{k=1}^{m} \left(1 - r\_{jk}\right) \tag{3}$$

Step 5 calculates the criteria importance by utilizing Equation (4):

$$\mathcal{W}\_{\dot{l}} = \frac{H\_{\dot{j}}}{\sum\_{k=1}^{m} H\_k} \tag{4}$$

#### *3.2. Entropy Method*

According to Zhang [69], the entropy weight method was originally a concept of thermodynamics, which was first added into the information theory by C.E. Shannon and is now applied widely in the fields of engineering technology, social economy, etc. When it comes to multicriteria, Randelovi´ ¯ c et al. [70] emphasized that entropy was mainly used to determine the priority of an alternative. According to Randelovi´ ¯ c et al. [70] the method of entropy is described as follows: let us assume that *cj* = (*a*1*j*, *a*2*j*, ... *amj*) describes a priority vector according to an exact criterion *j*, *j* = 1, . . . , *n*. The entropy value for this vector can be calculated by applying Equation (5):

$$H\_{Wj} = -\sum\_{i=1}^{m} a\_{ij} \ln \left( a\_{ij} \right), j = 1, \dots, n \tag{5}$$

In addition, Randelovi´ ¯ c et al. [70] emphasized that in the theory of information, the entropy value *HWj* could be defined as a unit of discrete random variable *X* uncertainty, which could have a value from the fixed set (*x*1, *x*2,... *xn*) in such a way that the feasibility that *X* is equal to *xj* is given by *wj* and may be presented by Equation (6):

$$P(X = x\_j) = w\_j \tag{6}$$

#### *3.3. Hybrid Criteria Weights*

To obtain the hybrid criteria weights, the authors of this paper combined the criteria weights obtained by the entropy and CRITIC methods. The combination of those two methods is demonstrated in the following equation:

$$\text{Hydrotid Weight (n\*)} = 0.5 \cdot \text{Entropy Weight (n\*)} + 0.5 \cdot \text{CRITIC Weight (n\*)} \tag{7}$$

where (*n\**) represents the hybrid weight of the *n*th criterion.

The hybrid criteria weights were chosen to test how those weights may affect the ranking of the final alternatives, as well as to notice whether there was any difference if the entropy and CRITIC methods were coupled with the ARAS method separately. In addition, the combination of two objective methods gave more robust and stable results.

#### *3.4. The Additive Ratio Assessment (ARAS) Method*

The additive ratio assessment (ARAS) method is an MCDM method originally developed by Zavadskas and Turskis [71]. Boškovi´c et al. [66] declared that the ARAS method was very efficient and easy to implement in situations where multiple criteria are considered. There are several steps of the ARAS method described by Zavadskas and Turskis [71]:

Step 1 formulates an initial decision-making matrix which consists of *m* alternatives (rows) compared on *n* criteria (columns). The initial decision-making matrix is presented below:

$$X = \begin{bmatrix} \mathbf{x\_{01}} & \cdots & \mathbf{x\_{0j}} & \cdots & \mathbf{x\_{0n}} \\ \vdots & \ddots & \vdots & \ddots & \vdots \\ \mathbf{x\_{i1}} & \cdots & \mathbf{x\_{ij}} & \cdots & \mathbf{x\_{in}} \\ \vdots & \ddots & \vdots & \ddots & \vdots \\ \mathbf{x\_{m1}} & \cdots & \mathbf{x\_{mj}} & \cdots & \mathbf{x\_{mn}} \end{bmatrix}; i = \overline{0, m}, j = \overline{1, n} \tag{8}$$

where *m* represents a number of alternatives, *n* represents a number of criteria describing each alternative, *xij* describes the performance value of the *i*th alternative in terms of the *j*th criterion, and *x*0*<sup>j</sup>* shows an optimal value of *j*th criterion.

If the optimal value of *j*th criterion is unknown, then:

$$\begin{array}{l} \mathbf{x}\_{0j} = \max \mathbf{x}\_{i} \mathbf{x}\_{ij} \text{ if } \max\_{i} \mathbf{x}\_{ij} \text{ is preferred};\\ \mathbf{x}\_{0j} = \min\_{i} \mathbf{x}\_{ij'}^{\*} \text{ if } \min\_{i} \mathbf{x}\_{ij}^{\*} \text{ is preferred}.\end{array} \tag{9}$$

Usually, the performance values *xij* and the criteria weights *Wj* are considered as the entries of a DMM. The system of criteria and the values and initial weights of criteria were determined by experts. The information can be corrected by the interested parties by considering their goals and opportunities.

Step 2 is the normalization of the input data of the initial decision-making matrix from the step 1. The normalization means that all the input data should be between an interval from 0 to 1. The normalized values *xij* in the normalized decision-making matrix *X* are obtained by applying Equations (11) and (12):

$$
\overline{X} = \begin{bmatrix}
\overline{x}\_{01} & \cdots & \overline{x}\_{0j} & \cdots & \overline{x}\_{0n} \\
\vdots & \ddots & \vdots & \ddots & \vdots \\
\overline{x}\_{i1} & \cdots & \overline{x}\_{ij} & \cdots & \overline{x}\_{in} \\
\vdots & \ddots & \vdots & \ddots & \vdots \\
\overline{x}\_{m1} & \cdots & \overline{x}\_{mj} & \cdots & \overline{x}\_{mn}
\end{bmatrix}; i = \overline{0, m}, j = \overline{1, n}; \tag{10}
$$

For the criteria with the highest preferable merits, the normalization is computed by utilizing Equation (11):

$$\overline{\mathfrak{X}\_{ij}} = \frac{\mathfrak{x}\_{ij}}{\sum\_{i=0}^{m} \mathfrak{x}\_{ij}};\tag{11}$$

For the criteria with the lowest preferable merits, the normalization is obtained by Equation (12):

$$\mathfrak{x}\_{\mathrm{ij}} = \frac{1}{\mathfrak{x}\_{\mathrm{ij}}^{\*}} ; \overline{\mathfrak{x}\_{\mathrm{ij}}} = \frac{\mathfrak{x}\_{\mathrm{ij}}}{\sum\_{i=0}^{m} \mathfrak{x}\_{\mathrm{ij}}} ; \tag{12}$$

Step 3 involves defining a normalized-weighted matrix *X*ˆ . It is possible to evaluate criteria with weights 0 < *Wj* < 1. Only well-founded weights should be used, because weights are always subjective and influence the solution. The values of weight *Wj* are usually determined by the expert evaluation method. The sum of the weights *Wj* is limited as follows:

$$\sum\_{j=1}^{n} w\_j = 1;\tag{13}$$

$$\hat{X} = \begin{bmatrix} \mathfrak{k}\_{01} & \cdots & \mathfrak{k}\_{0j} & \cdots & \mathfrak{k}\_{0n} \\ \vdots & \ddots & \vdots & \ddots & \vdots \\ \mathfrak{k}\_{i1} & \cdots & \mathfrak{k}\_{ij} & \cdots & \mathfrak{k}\_{in} \\ \vdots & \ddots & \vdots & \ddots & \vdots \\ \mathfrak{k}\_{m1} & \cdots & \mathfrak{k}\_{mj} & \cdots & \mathfrak{k}\_{mn} \end{bmatrix}; i = \overline{0, m}\_{\prime} \text{ } j = \overline{1, n}$$

Normalized-weighted values of all the criteria are calculated as follows:

$$\pounds\_{ij} = \overline{\pi}\_{ij} \cdot \mathcal{W}\_{\mathfrak{j}}; i = \overline{0, m}; \tag{15}$$

where *Wj* is the weight (importance) of the *j*th criterion and *xij* is the normalized rating of the *j*th criterion.

Step 4 determines the value of optimality function:

$$S\_i = \sum\_{j=1}^{n} \pounds\_{ij}; \ i = \overline{0, m}; \tag{16}$$

where the optimality function of *i*th alternative is marked with *Si*.

The maximum value of *Si* reflects the best option, while on the contrary, the minimum value reflects the worst option. In other words, the greater the value of the optimality function *Si*, the more effective the alternative. The preferences of alternatives can be observed according to the value *Si*.

Step 5 computes the level of the alternative utility. To do so, it is necessary to make a comparison of the solutions with is ideal solution (*S*0). The calculation of the utility level *Ki* of an alternative *ai* is calculated by Equation (17):

$$K\_{\bar{i}} = \frac{S\_{\bar{i}}}{S\_0}; i = \overline{0, m}; \tag{17}$$

where *Si* and *S*<sup>0</sup> are the optimality criterion values. The calculated values *Ki* are between 0 and 1.

#### **4. Application of the Hybrid–ARAS Method to 3PL Evaluation and Selection**

The previously described methodology was applied to a hypothetical example, and the results are described in this section. The main reason for the hypothetical example was, as aforementioned, because it was hard to obtain real data for 3PL selection given the time constraints of this research, etc. However, three experts helped the authors define the criteria that influence the decision-making process and agreed on their ranking by the entropy and *CRITIC* methods. In this case, the authors selected five 3PL providers as the possible alternatives. According to the experts' opinions, five criteria that should be taken into consideration when evaluating and selecting 3PL providers were selected. These criteria were price (C1), delivery service (C2), quality of service (QoS) from customer experience (C3), territorial coverage of the EU (C4), and flexibility (C5). The selected criteria are completely expressed in the following table (Table 2). It is important to point out that the considered criteria were mostly included in the economic aspect. When it comes to the consulted experts' knowledge and expertise, expert 1 held a managerial position in a tire manufacturing company with four years of experience and had a Ph.D. in the field of logistics; expert 2 held a managerial position in a multinational beverage company with five years of experience and had a Master's degree; and Expert 3 was a manager of a cold

chain company with more than eight years of experience and a Ph.D. in the field of logistics and supply chains. Because of the COVID-19 outbreak, the authors interviewed the experts by telephone. The experts gave some limited information about themselves but required complete anonymity because of the business policies of their firms. The hypothetical input data used in the ARAS decision-making matrix were formulated in collaboration with the experts as well; the experts agreed that the data should respond to real conditions in the 3PL logistics market.

**Table 2.** Criteria for 3PL evaluation and selection (Source: Authors).


After the description of the criteria, the entropy and CRITIC methods were used in order to obtain the criteria weights. First, the weights were obtained by each separate method. Second, hybrid weights were obtained by combining the results of the two methods by applying Equation (7). The entropy method was applied to find criteria importance. After applying the entropy method, the following criteria weights were obtained (Tables 3–5).


**Table 3.** Initial entropy decision-making matrix (Source: Authors).

**Table 4.** Normalization of the entropy decision-making matrix (Source: Authors).



**Table 5.** Computed entropy value (h) and final weights (Source: Authors).

According to the entropy method, the highest importance was assigned to price and QoS, while lesser importance was assigned to territorial coverage of the EU, delivery service, and flexibility. By applying the previously described CRITIC method, the following criteria weights were calculated (Tables 6–9).


**Table 6.** Initial CRITIC decision-making matrix (Source: Authors).

**Table 7.** Initial CRITIC decision-making matrix with standard deviations (Source: Authors).



**Table 8.** *m* × *m* matrix (Source: Authors).

**Table 9.** *m* × *m* matrix with the final weights (Source: Authors).


According to this method, the highest importance was assigned to price, followed by QoS from customer experience, territorial coverage of the EU, delivery service, and flexibility. The following table (Table 10) compares the criteria weights obtained from both methods and the hybrid weights obtained by applying Equation (7).

**Table 10.** Obtained hybrid criteria weights (Source: Authors).


The hybrid method ranked the criteria in a following way: the highest importance was assigned to price (0.2983), second place was related to QoS from the customer perspective (0.2486), third place was related to territorial coverage of the EU (0.1947), and flexibility (0.1755) and delivery service (0.0828) had less importance. Using the criteria weights obtained by the hybrid method, the final ranking of the 3PL providers was obtained by applying the ARAS method. For better clarity, the obtained criteria weights by the hybrid method are presented in Figure 2.

**Figure 2.** Obtained hybrid criteria weights (Source: Authors).

#### *4.1. Application of the Hybrid–ARAS Method to 3PL Evaluation and Selection Problem*

After the criteria weights were determined, the ARAS method was applied to obtain the final ranking of the 3PL providers. The input data used in the ARAS decision-making matrix are presented in Table 11.


**Table 11.** Initial ARAS decision-making matrix (Source: Authors).

The normalization of the input data is presented in Table 12.



The weighted decision-making matrix is described in Table 13.


#### **Table 13.** Weighted D–M matrix (Source: Authors).

As shown in Figure 3, by applying the hybrid–ARAS method, the best alternative was shown to be 3PL-2, with the preference value of 0.9873, followed by 3PL-1, with the preference of 0.9270; 3PL-3, with the preference of 0.9141; 3PL-4, with the preference of 0.8930; and in the last place 3PL-5, with the preference of 0.8747.

**Figure 3.** Final ranking of 3PL providers obtained by the hybrid–ARAS method (Source: Authors).

#### *4.2. Sensitivity Analysis*

To test the stability of the proposed hybrid–ARAS method, a sensitivity analysis was conducted. The main purpose of the analysis was to examine how the change in trade-off parameter ξ affected the final ranking of the alternatives. In this regard, the parameter ξ was changed within the interval of [0, 1] with an increment value of 0.1 (Figure 4).

**Figure 4.** The sensitivity analysis to changes in the trade-off parameter ξ (Source: Authors).

When ξ = 1, only the entropy ARAS method was applied to prioritize the 3PL providers. When ξ = 0, only the CRITIC ARAS method was used to evaluate the 3PL providers. Therefore, in the base case scenario, ξ was set to 0.5 to equally appraise both methods and generate hybrid criteria importance. According to Figure 4, 3PL-2 was the best alternative under all ξ values. In addition, there was no change in the ranks of any 3PL provider in all 10 new test cases; i.e., the 3PL service provider ranking order was 3PL-2 > 3PL-1 > 3PL-3 > 3PL-4 > 3PL-5. The sensitivity analysis revealed that the proposed model has a high level of stability.

#### *4.3. Comparative Analysis*

A comparative analysis was performed to check the reliability of the results obtained through the hybrid–ARAS method for 3PL selection. The 3PL selection process was solved with two state-of-the art MCDM approaches, WASPAS [72] and EDAS [73]. The result of the comparative analysis is presented in Figure 5.

**Figure 5.** The comparative analysis of the hybrid–ARAS approach with WASPAS and EDAS.

The proposed hybrid–ARAS and EDAS methods ranked 3PL-2 as the best solution and 3PL-5 as the worst one. When it comes to 3PL-4, it was ranked as the fourth best option according to hybrid–ARAS and EDAS, but the WASPAS method ranked it as the best solution. According to the results, the model is reliable.

#### **5. Managerial Insights**

Since third-party logistics play an important role in the logistics market, managers all around the world should carefully monitor and identify the most important parameters that may affect this field. Not all the parameters are equally important, so sorting out the most important ones may greatly influence the 3PL business. When identifying the criteria for 3PL selection, each company should consider the whole picture of its business environment and select the most important criteria to attract the best collaboration partner according to its expectations. In addition, collaboration between experts is of crucial importance for 3PL selection. The wrong choice of a 3PL provider can have negative consequences for a company such as loss of profit, business image, reputation, customer loyalty, etc. The recommendation for managers should be to eliminate subjectivity and combine one or more objective methods to identify the priority of the criteria with MCDM methods in order to obtain the best possible 3PL provider to collaborate with. The results of using more objective methods are more robust and confident, which leads to a more stable managerial solution.

#### **6. Discussion and Conclusions**

This paper aimed to propose a possible hybrid–ARAS approach to the 3PL evaluation and selection process. The main reason for the methodology proposed was to give a theoretical contribution in the 3PL logistics field that could help managers and scientists think about the possible approaches for evaluating and selecting 3PL service providers. Three methods were combined in order to obtain a ranking of the best alternatives. The first method was the entropy method, applied to obtain objective criteria importance. The second was the CRITIC method, utilized to obtain objective criteria importance. The

third was the ARAS method combined with the entropy and CRITIC methods in order to obtain a ranking of the possible alternatives. Furthermore, the combination of the entropy and CRITIC methods resulted in obtaining *hybrid* criteria weights, which were further coupled with the ARAS method to rank the best alternative. In each case, either separate (entropy–ARAS; CRITIC–ARAS) or coupled (hybrid–ARAS), the final ranking results were not changed. The main idea of coupling the former two methods was to obtain objective criteria weights, which should eliminate subjectivity and give more robust results.

To examine the stability of the proposed hybrid–ARAS method, a sensitivity analysis was conducted. The results of the analysis revealed that there was no changing in the ranking alternatives when 10 variations in criteria were taken. In addition, a comparative analysis was carried out to compare the obtained results. The results of the comparative analysis revealed that the model is reliable.

The major contributions of this paper are: (i) for the first time, a new combination of methods was introduced to rank 3PL service providers; (ii) a combination CRITIC–entropy (hybrid) method was employed to prioritize the criteria that influence the 3PL evaluation and selection process—a more robust solution was obtained by coupling two objective methods into a hybrid one, and subjectivity was eliminated; (iii) a new hybrid–ARAS method was developed to rank 3PL service providers; (iv) the methodology was very simply described and easy to implement and should be beneficial for the 3PL logistics field.

In our hypothetical case, the hybrid–ARAS MCDM approach for 3PL evaluation and selection process generated the following ranking order: *A2* (3PL-2) > *A1* (3PL-1) > *A3* (3PL-3) > *A4* (3PL-4) > *A5* (3PL-5). This approach identified 3PL-2 as the best possible alternative. On the other hand, the worst-ranked alternative was 3PL-5. As a result, it would be strongly recommended to select 3PL-2, since it was shown as the best alternative according to the methodology.

Limitations of this paper can indicate possible areas for its extension. The limitations included: (1) the methodology was not applied to a real-life case, but only to a hypothetical example; (2) the criteria influencing the decision-making process were filtered not only according to a review of the literature, but by the consideration of experts' opinions indeed, filtering was mostly performed according to discussion with experts in the field. However, the methodology is general and can be applied with any other criteria influencing the 3PL evaluation and selection process. The main point of the paper was to show the applicability of the methodology in order to respond to the aforementioned research questions; (3) there were only considered criteria from the economic pillar. However, there is space to include any other aspects that should matter in the decision-making process, such as environmental, technical, and social aspects; (4) the authors did not consider how 3PL logistics may fit with the industry 4.0.

This study can be seen as an important trigger for future papers in the field. The future directions inspired by this paper should be to (a) apply the methodology to real-life cases; (b) include the criteria most often used by other authors in the field; (c) apply the methodology in the fuzzy and picture fuzzy environment, since fuzzy logic deals with uncertainty, which often factors into the decision-making process; (d) examine how 3PL logistics may fit with the industry 4.0.

**Author Contributions:** Conceptualization, S.J. and P.P.; methodology: S.J.; validation: S.J. and P.P.; investigation: S.J. and P.P.; writing—original draft preparation, S.J. and P.P.; writing—review and editing: S.J. and P.P.; visualization, S.J. and P.P.; supervision, S.J. and P.P. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Not applicable.

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

#### **References**

