Modified Interactive Algorithm Based on Runge Kutta Optimizer for Photovoltaic Modeling: Justification under Partial Shading and Varied Temperature Conditions

The accuracy of characteristic the PV cell/module/array under several operating conditions of radiation and temperature mainly relies on their equivalent circuits sequentially; it is based on identified parameters of the circuits. Therefore, this paper proposes a modified interactive variant of the recent optimization algorithm of the rung-kutta method (MRUN) to determine the reliable parameters of single and double diode models parameters for different PV cells/modules. The results of the MRUN optimizer are validated via series of statistical analyses compared with five new meta-heuristic algorithms including aquila optimizer (AO), electric fish optimizer (EFO), barnacles mating optimizer (BMO), capuchin search algorithm (CapSA), and red fox optimization algorithm (RFSO) moreover, twenty-five state-of the art techniques from literature. Furthermore, the identified parameters certainty is evaluated in implementing the characteristics of an entire system consists of series (S), and series-parallel (S-P) PV arrays with numerous dimensions. The considered arrays dimensions are three series (3S), six series (6S), and nine series (9S) PV modules. For the investigated arrays, three-dimensional arrays are recognized. The first array comprises 3S-2P PV modules where two parallel strings (2P) have three series modules in each string (3S). The second array consists of six series-three parallel (6S-3P) PV modules, and the third one has nine series-nine parallel (9S-9P) PV modules. The results prove that the proposed algorithm precisely and reliably defines the parameters of different PV models with root mean square error and standard deviation of 7.7301e-4 ± 4.9299e-6, and 7.4653e-4±7.2905e-5 for 1D, and 2D models, respectively meanwhile the RUN have 7.7438e-4 ± 3.5798e-4, and 7.5861e-4 ± 4.1096e-4, respectively. Furthermore MRUN provided extremely competing results compared to other well-known PV parameters extraction methods statistically as it has.


I. INTRODUCTION
Because of the rising and unstable prices of fossil fuels and pollution and solid wastes, renewable energy sources have emerged as a viable option for sustaining energy while reducing pollution [1]. Solar energy can be considered the  Number of parallel strings in the array P max maximum power at SOC most efficient alternative to fossil fuels and coal when power system operators meet energy needs. Solar energy is the most potential renewable energy source because of its global distribution, low maintenance, noise-free operation, and near-conventional production processes. Over the last three decades, efforts have been made to transition from small to large-scale solar cells or from cell size to module size [2]. This trend reflects advancements in solar cell technology, which have increased utilization levels from a small scale (few cells) to a big scale (many modules). Many limitations that previously limited its use in largescale power system applications have now been removed.
Photovoltaic (PV) power is one of the preeminent universal renewable electricity methods behind wind power transforming methods worldwide [3]. Lately, active investments target the PV production outstandingly, and many revenues are advised to research in PV energy. Usually, such regulation arises from the severe cost reduction in the PV operation parts, which designates the enormous development of the PV industry [4], [5]. The cost of PV machines pointedly fell within the recent decade, including the PV operation itself with its design and adapters [6]. Thus, the PV price adjustment increases its defiance with other renewable energy operations. According to the statistics estimation budget of the universal requirement for PV operations, a fifty percent achievement rate of PV investments was returned with 306.5 GW in 2016 [7].
In this context, the PV system has gotten a lot of press in recent years [8], [9]. However, to examine the dynamic conversion behavior of a PV system, one must first understand how to model the PV cell, which is the system core component [10]. To represent PV cells, a variety of methodologies have been devised, the most prevalent of which is the use of comparable circuit models. The single diode and double diode types are the most often utilized circuit types [11]. The correctness of the parameters associated with the model structure is critical for estimating, sizing, performance assessment, management, efficiency computations, and reactive power control of solar PV systems after finding a suitable model structure [12].
Because of the nonlinear nature of photovoltaic models and the rising number of parameters that must be evaluated, metaheuristic approaches inspired by numerous natural phenomena have been widely employed to find parameters of PV models as a possible alternative to deterministic methodologies [13], [14]. Metaheuristic methods are employed to extract parameters with high precision. They do not put any constraints on the problem characteristics; therefore, they are simple to apply to a variety of real-world challenges [15]. The Grasshopper Optimization Algorithm (GOA) is used in [16] to extract parameters from a three-diode PV model. The simulation results are run at different temperatures and levels of irradiation. The GOA-based PV model efficiency is assessed by comparing its numerical findings to other optimization method-based PV models. A new optimization approach termed turbulent flow of water-based optimization is proposed in [17] to extract the characteristics of three photovoltaic (PV) cell models (TFWO). Compared to existing approaches, the suggested TFWO achieves a high degree of similarity between the estimated voltage-power (V-P) and current-voltage (I-V) curves. A gradient-based optimizer (GBO) is used in [18] to estimate the parameters of solar cells and PV modules in an efficient and precise manner. To show the GBO ability to estimate the parameters of solar cells, three popular solar cell models are used. Compared to the experimental, the suggested GBO achieves a high degree of closeness between the simulated V-P and V-I curves.
The work in [19] provided an improved chaotic JAYA algorithm for precisely and reliably classifying the parameters of various photovoltaic models, such as singlediode and double-diode models. Furthermore, on multiple phases of the search space, the proposed method includes a self-adaptive weight to govern the trend and obtain the ideal solution while avoiding the worst result. Compared to previous algorithms in the literature, comprehensive analysis and practical results show that the proposed algorithm achieved highly competitive efficiency in terms of accuracy and dependability. In [20], a numerical approach is proposed for finding the five-parameter model of photovoltaic cells. To investigate the relationships between parameters, explicit equations are used, which are then solved by an optimization method. The suggested approach produces accurate results, and it can be used with a variety of solar cells. This paper proposed a classified perturbation mutation-based particle swarm optimization algorithm [21]. The performance of each updated personal best position is evaluated and quantified as high-quality or low-quality during each generation of the suggested algorithm. The presented technique outperformed other well-known parameter extraction methods in terms of accuracy, stability, and speed, as demonstrated by the results of the experiments. This paper proposed an improved wind-driven optimization approach for identifying the nine unknown parameters [22]. The suggested technique is a hybrid of the differential evolution algorithm mutation strategy and the wind-driven optimization algorithm covariance matrix adaption evolution strategy. The findings showed that in terms of accuracy, convergence speed, and practicality, the revised wind-driven optimization model surpassed the other models.
Based on the literature, most researchers applied their algorithms to determine the cell/module parameters without testing those parameters in emulating the behavior of a complete array. Therefore, the certainty of the identified parameters has not been validated with simulating an entire PV string or array under different environmental conditions. Further, proposing the reliable and efficient optimizer to determine the PV models parameters with high certainty for implementing its physical behaviors is a challenging task because the actual values of the model parameters are unknown, so any improvement in fitting accuracy is regarded as highly beneficial from a modeling standpoint [12]. Furthermore, the No Free Lunch theorem in [23] states that no one method can handle all optimization issues. As a result, the search for an alternate optimization technique to consistently identify the PV model parameters is ongoing. Therefore, this paper introduced a new metaheuristic search method for parameter extraction of PV models using a modified variant of the Runge Kutta (RUN) method. RUN is, without a doubt, one of the most influential and versatile search methods used today. It has the benefit of being simple to execute and effective. However, although RUN has recently been used to tackle various real-world problems, it has a significant issue in discovering the search space efficiently with nonlinear optimization problems and exploiting the solutions. Therefore, in this work, a modified variant of RUN (MRUN) is proposed to tackle these drawbacks and handle the nonlinear optimization problem of PV modeling with high accuracy and consistency. The following lines sum up the analyses through the paper:-• A novel variant of the RUN optimizer is proposed for identifying the 1D and 2D models parameters utilizing a set of several experimental data for PV cells and modules under numerous operating conditions of radiation and temperature. • The proposed method results are validated via Lambert equations for the two models; 1D and 2D. • The proposed optimizer results are compared with recent meta-heretic algorithms are implemented on the same settings and well-known PV parameters extraction methods in the literature. • The quality of the extracted model parameters is justified in implementing the physical behavior of entire systems consisting of a set of strings (S), and seriesparallel (SP) arrays PV arrays subject to multiple levels of radiation and temperature.
The results show that the proposed algorithm precisely and reliably defines the parameters of different PV models besides providing extremely competing results compared to other well-known PV parameters extraction methods statistically. Furthermore, the results divulge that the proposed modification enhanced the core abilities of the optimizer. For example, MRUN is good at searching the search space and discovering the optimum global region. Moreover, it has an excellent ability at exploiting the solutions, so even when the population has not converged to a local optimum, the algorithm can find any better answers to advance the evolution. Finally, the results of the entire S and SP systems reveal the certainty of the identifying parameters in acting the PV arrays characteristics with high efficiency.
The rest sections of this paper is given as follows. Section II presents a description of the equivalent circuits. Problem definitions are given in Section III. The details of the proposed optimizer is presented in section IV. The results and discussions are established in section V. The VOLUME 4, 2016 justification process of the proposed optimizer results id presented in section 6, finally, Conclusion and future open direction are given in Section VII. In this section, the basic model of Photovoltaic (PV) equivalent circuits is introduced. In general, the most popular PV models are the single diode model (1D) and the double diode model (2D) (see Fig. 1). Each of one of these models has its own characteristics; for example, 1D is considered as the simplest equivalent circuit. Whereas 2D is more efficient than 1D since 2D simulates the physical quality of PV models at low level of irradiation conditions [24]. In this model (i.e., 2D), the first and second diode stand for the diffusion current and recombination effects, respectively. The SD consists of one diode and photo that produced shunt resistance (R sh ) and current. After that, the integration become in series form with a resistance (R s ). Then the output PV current (I) is calculated based on the law of kirchhoff current as given in Eqs. (1)-(2) [7], [10], [25].

II. PHOTOVOLTAIC EQUIVALENT CIRCUITS
In Eq. (1), a 1 is ideality factor of the diode. the I p , I o1 , and I d1 stand for leakage shunt, saturation diode and the diode currents, respectively. Also, V t represents the thermal voltage and it is measured at temperature (T in Kelvin) as KT q , where k = 1.35 × 10 −23 J/K is Boltzmann constant and q = 1.6 × 10 −19 C denotes the electron charge.
Besides Eq. (2), there are five parameters (i.e., I ph , I o1 , a 1 , R S , and R sh ) are required to be determined. From Fig. 1(b) it can be observed that the 2D is an extension of 1D by integrating a parallel second diode with the first diode in SD. This simulates the physical influences at the P-N junction and similar to 1D, the output current of PV is computed as [7], [10], [25].: In Eq. (3), a 1 (a 2 ) stands for the ideality factor of first(second) diode. I o2 stands for the saturation current. Besides Eq. (3), there are seven parameters (i.e., I ph , I o1 , I o2 , a 1 , a 2 , R s , and R sh ) are required to be identified.
In addition, the produced photocurrent is computed based on incident radiation value (G) at T as Eq. (4b). As well as, the reverse saturation currents of I o1,2 are computed as in Eq. (4c). While, R sh depends on G and it is defined in Eq. (4e) and Eq.(4f) represents the process of computing the open circuit voltage (V oc (T ) ) at temperature T [26], [27].
where V oc (s) , R ps , I o1,2 s , and I phs stand for the open circuit voltage, shunt resistance, reverse saturation currents, and photo current, respectively. The Cf v (%/ • C) and Cf i (%/ • C) denote the temperature coefficient of voltage and current, respectively. In addition, E g stand for the band-gap energy that defined as [26], [27]. : where E g (s) represents the band-gap energy at standard operating conditions (SOC).

III. IMPLEMENTED OBJECTIVE FUNCTION
The process of identify the parameters of the 1D and 2D is considered as a nonlinear optimization problem. The most popular objective function that used to achieve this process is the root mean square error (RMSE) that computed using the value of the estimated (I est ) and the measured (I meas ) current. To make the objective function more accurate and suitable for real-world application, the newton-raphson is used to solve the system of nonlinear equations and this can be formulated as: where Z stand for the vector that contains the set of estimated parameters. Whereas, M represents the length of the measured values. To compute the value of the estimated current (I estt+1 ), the extracted parameters are used with find the solution of Eqs. 2, and 3 using Newton-raphson method as formulated in Eq. ((7)): In Eq. (7), dI and dI represent the a difference function of I and its first derivative. For clarity, the definition of dI and dI for 1D is formulated as: According to Eqs. 8-9 and substituting them in Eq. 7, then the estimated current is computed as: Rs a1V th exp V +Iest t Rs a1V th − Rs R sh − 1 (10) Similar to computed I estt+1 for 1D, we can compute the I estt+1 for 2D where five iterations are enough to find the solution using Newton-raphson method.

A. VERIFYING THE RESULTS WITH LAMBERT FORM
To evaluate the performance of the identified parameters using the developed method, we used the Lambert W function (LWF) to compute the currents of 1D and 2D. The Lambert W function in mathematics is a set of functions that are the branches of the inverse relation of the β function as shown below y = f exp y (11) where exp y represents the exponential function, and y is any complex number. The y of Eq.(12) can be written based on Lambert W form as follows: where W is the solution of the Lambert equation. By using the same concept several engineering applications can be formulated. One of those applications is solving the PV characteristic equations [28]- [30]. In general, the LWF of 1D PV model (Eq. (2)) is defined as: In addition, the LWF of 2D (Eq. (2)) is defined as: where I Lambert denotes the output current obtained using LWF.
So, the RMSE has been recomputed for I est using LWF according to the estimated parameters and I meas . In case there is large difference (Diff RM SE ) between the RMSE (as in Eq.6) and RMSE based on LWF (i.e., RMSE Lambert ), so the process of estimated parameters is inefficient [28]- [30]. The mathematical form of RM SE Lambert is computed as below:

IV. MODIFIED ALGORITHM BASED ON RUNGE KUTTA METHOD
In this section, the proposed enhanced Algorithm Based on Runge Kutta Method (MRUN) is proposed to modify the classical RUn optimizer performance. The details of the proposed are described below in the following subsections.

A. OVERVIEW OF THE BASIC ALGORITHM BASED ON RUNGE KUTTA (RUN) METHOD
The RUN optimization algorithm is inspired by the Runge Kutta method (RKM) [31] that was applied to find the solution for the ordinary differential equations. In general, RKM generates a high-precision numerical value based VOLUME 4, 2016 on the functions only and doesn't require any gradient information (Zheng & Zhang, 2017). Therefore, the RUN optimization algorithm depends on computing the slope in RKM, that used as a searching strategy to emulate the exploration ability in swarm-based optimization. The mathematical formulation of the RUN algorithm contains a set of stages that are discussed as follows: • Initialization stage: In this stage, the initial solutions of N agents are constructed based on the boundaries of the search landscape [LB, U B]. This step is conducted using the following formula: where D represents the dimension of the given problem. The LB j and U B j are the lower and upper boundaries of j th variable in the solution set Z i,j ; i=1,2,...,N , the symbol N refers to the total number of search agents. • Updating solutions stage : The RUN algorithm uses a search mechanism (SM) based on the Runge Kutta method to update the position of current solution at each iteration, which is defined as: is an integer number used to change the search process direction. The symbols of g ∈ [0, 2] and µ ∈ [0, 1] are random numbers. The SF represents adaptive factor that defined as: where M axt denotes the total number of iterations. Tthe Z c and Z m given in Eq. (18) are defined as below: In Eq. (21), the ϕ ∈ [0, 1] denotes a random number. The symbols of Z b and Z pb stand for the best-so-far agent and the best one at each iteration, respectively. For the SM parameter that defined in Eq. (18) is updated using the following formula: where rand 1 and rand 2 stand for random numbers. The value of ∆Z is computed as: where the value of Z b and Z w are updated as: • Enhanced solution quality stage: In this stage, the quality of solutions is enhanced using different operators to improve the convergence rate with skipping the local optima. This process is defined as: (25), β ∈ [0, 1] represents random number and r ∈ {1, 0, −1} is an integer number.
Following [31], in case the fitness value of Z new2 not better than the fitness value of Z i then there is another chance to modify the the value of Z i . The solutions can be updated using the following formula: where v = 2 × r 3 is a random number stands for the interval of 2 × [0 1]. r 1 , r 2 and r 3 are random values.
The pseudo-code of the classical RUN is exhibited in algorithm 1 to summarize the main steps of the algorithm structure. Compute fitness for each solution with determine the best (Z b ), worst Z w and Z pb solutions. 5: for (i=1 to N) do 6: for (j=1 to D) do 7: Using Eq. (18) to update Z i . 8: if rand < 0.5 then 9: Using Eq. (24) to obtain Z new2 10: if rand < w then

B. STEPS OF THE PROPOSED MODIFIED RUN (MRUN)
As described in the previous section of classical RUN, the main two stages are updating the solutions and enhancing them. In updating the solution stage, the transition process between the exploration and exploitation stage based on a random process where if rand ≤ 0.5, the exploration phase is performed else the exploitation phase is implemented. Such of this approach may be cause inconsistency in the attained solutions by the RUN. For modifying this process, in MRUN, the transition process between the exploration and exploitation phases is divide in the search agents where in the first third of the search agents group the exploration is performed and for the rest of the agents the exploitation is executed as reported in Algorithm 2. Furthermore, in the MRUN, the main controller variable in the enhanced solution stage ω is modified via using pareto heavy-tied distribution of Eq.27 rather than uniform distribution rand(0,2) to enhance the new solutions Z new2 .
Algorithm 2 Pseudo-code of enhancing the updated solution stage of the classical RUN 1: for (i=1 to N) do 2: if i< N/3 then 3: Use the following formula Use the following formula Pareto distribution (PF): a random variable has been shown to follow Pareto distribution if it has the following tail function PF: where a, and b are the scale and shape parameters and have values of 0.0001 and 2, respectively. the x is random vector with dimension (1,D). Then the updated ω based on PF is written as follows: The Fig.2 depicts the flowchart of the proposed MRUN. The algorithm starts with set of random solution while these solutions are modified using Algorithm 2 in the updated solution stage and employing pareto front in the phase of enhancing the solutions of Eqs. 24, 28, 26.

V. SIMULATIONS AND DISCUSSIONS
The proposed ERIN is assessed throughout several sets of experimental series. Firstly, the proposed optimizer is applied to identify the parameters of 1D and 2D models using RTC France datasets. Then, in the second series of experiments, the parameters of 2D models of three modules including Kyocera Solar KC200GT PV module, Sharp NU-(Q250W2) panel and Pythagoras Solar Large PVGU Window under different operating conditions of temperature and irradiation. The manufacture properties of the realized panels at SOC are reported in Table 1. The MRUN technique results are compared with set of optimization techniques including aquila optimizer (AO) [32], electric fish optimizer (EFO) [33], barnacles mating optimizer (BMO) [34], capuchin search algorithm (CapSA) [35], and red fox optimization algorithm (RFSO) [36] to investigate the proposed performance. The setting values of the population size, number of iterations and number of independent runs are 50, 500 and 25, respectively. These   values are applied for all the implemented algorithms for achieving unbiased comparison. The upper (B max ) and lower (B min ) boundaries for the realized models unknown parameters are tabulated in Table 2.

A. MODELING OF PV SOLAR CELL
In this section a set of 26 (V-I) measured pairs for a commercial RTC France silicon solar cell at 1000W/m 2 and 33 • C is used while evaluating the performance of the proposed optimizer. The identified parameters of 1D and 2D models by MRUN, other peers and set of recent literature are reported in Table 3 with the established fitness function value (RMSE) as a primary metric for the accuracy of the results. Moreover, the accuracy of the estimated parameters is evaluated though implementing Lambert forms for 1D and 2D models of section III-A to compute the RMSE lambert then calculating the deviation (Diff RM SE ) between the RMSE lambert and the attained RMSE of Eq. (6), the large deviation refers to an inefficient identified parameters. Furthermore, the absolute error at MPP (AE M P P ) is reported in the Table 3 for providing an extensive analysis.
From the results given in Table 3, it can be noticed that RMSE and Diff RM SE by MRUN is smaller than other implemented techniques (RUN, AO, EFO, BMO, CaSA, RFSO) in the two test cases (i.e., 1D and 2D). This indicates the high performance of the developed MRUN approach and affirm its superiority over the other approaches. The AE M P P values by the MRUN show its reliability in modeling the PV solar cell with high accuracy.
For further analysis of the performance of MRUN, different statistical measures are used. These measures are the best, worst, average, median, variance, and standard deviation (StD). The Wilcoxon signed-rank test is also applied to if the difference between the MRUN and other methods is significant or not, and this is evaluated at a significant level of 0.05. From these tabulated results, it can be observed that MRUN allocates the first rank in terms of best, worst, average, and median. In addition, its stability is better than other methods. The traditional RUN can provide better results than the different algorithms; however, BMO is competitive. Moreover, the p-value obtained using the Wilcoxon signed-rank test indicates a significant difference between MRUN and other methods.
For affirm the certainty of the identified parameters, Fig.  3 and Fig. 4 are reported to depict the current-voltage (V-I) characteristic, Power-voltage (V-P) characteristic, absolute error curve between measured and estimated datasets, Mean convergence curve and RMSE throughout 25 independent times using 1D and 2D model of RTC France solar cell, respectively. From these graphs it can be reached to the following observations; the identified parameters by MRUN provides a closely matching between the measured and estimated that sets. This observation is confirmed from the absolute error curves of Figs.3(c)-4(c) for 1D and 2D models, respectively. The mean convergence rate of MRUN of Figs. 3(d)-4(d) converges to the highest quality solutions compared to the other techniques. In addition, the boxplots of RMSE of Figs.3(e)-4(e) indicate the high stability of MRUN in 1D and 2D models. Followed by EFO and CapSA, however, AO is the worst algorithm in both models.
The results of Table 6 reveal that the developed MRUN provides the most powerful results in comparison to most of the recent state-of-the-art techniques. Meanwhile, the MRUN has the same RMSE similar to TVACPSO and CPSO when applied to determine the parameters of the 1D model. For the 2D, the MRUN statistical metrics values prove the algorithm superiority in providing the most accurate results with comparable execution time.

B. MODELING OF PV SOLAR COMMERCIAL MODULES
Within this section, the developed method has been applied to determine the parameters of 2D model of KC200GT solar panel, Sharp NU-(Q250W2) panel, and Pythagoras Solar Large PVGU Window under three different operating conditions for each module. Table 9 shows the results of the developed MRUN and other methods using KC200GT, Sharp NU-(Q250W2), and Pythagoras Large PVGU Window PV module for 2D model. It can be noticed from the reported data that the high efficiency of the MRUN algorithm at the three different cases (i.e., at different temperatures). Since the Diff RM SE of MRUN is zeros among all the tested cases. However, it has been observed that the performance of other models are decreased at 1000W/m 2 45 • C in the three PV module. In some cases, the quality of other compared algorithms at 1000W/m 2 60 • C is better than at 1000W/m 2 45 • C and 1000W/m 2 25 • C.
Furthermore, the V-I, V-P characteristics, and AE curves are plotted in Figs. 5-6, and 7 respectively. The figures show the high matching between the measured and estimated datasets of the three studied models, manning the identified parameters high accuracy. The mean convergence curves of MRUN is better than the competitive algorithms (i.e., RUN, AO, EFO, BMO, CapSA, and RFSO).
For the statistical evaluation part, the difference between the developed MRUN and RUN is significant in all cases, as noticed from the P-value obtained by Wilcoxon signedrank test. In addition, the MRUN can get results better than RUN in terms of Best, Worst, Average, Median. Also, the developed MRUN is more stable than traditional RUN among the tested cases.
The accuracy of the extracted parameters is tested for further assessment while emulating the PV modules charac-terises under different operating irradiation, and temperature conditions as depicted in Figs. 8-9-10 for the studied three modules. The curves reveal the high qualified identified parameters by the MRUN as the most point of the AE curves of Figs. 8(c)-9(c)-10(c) are less than 0.02. Therefore the developed MRUN is recommended to provide efficient equivalent circuit parameters of the PV solar modules under different operating conditions.

VI. JUSTIFICATION UNDER PARTIAL SHADING AND VARIED TEMPERATURE CONDITIONS
The simulation implemented above and debates divulge the efficiency and accuracy of the extracted parameters for considered PV cells/panels under various operations cases. Thus, in this section, we are motivated to justify the recognized model parameters certainty in modeling connected strings with (N × 1, N is the number of series panels) and series-parallel arrays with (N × M ) panels under uniform and partially shaded operating conditions with considering temperature impact. For a brief, the data of the KC200GT PV panel is the established case in this section. The considered PV strings consist of three series (3S), six series (6S), and nine series (9S) PV modules. For the investigated arrays, three-dimensional arrays are recognized. The first array comprises 3S-2P PV modules where two parallel strings (2P) have three series modules in each string (3S). The second array consists of 6 series-3 parallel (6S-3P) PV modules, and the third one has 9 series-9 parallel (9S-9P) PV panels. The described series of experiments are listed in Table 11.

A. MATHEMATICAL FORMULATION FOR THE CURRENT, VOLTAGE AND POWER OF THE PV ARRAY
In this part, the strategies that followed in computing the string and SP interconnected arrays current and voltage under PS conditions are represented as follows: • In the first step, the parameters of the used PV equivalent circuit (1D or 2D) are identified at SOC. Under the different operating conditions, those parameters are normalized using equations of Section II. By this way, the V-I curves of the PV module are computed for the corresponding operating conditions. • For the string current and voltage values, suppose having a string with connected three modules receive three different levels of radiations that are 900W/m 2 , 400W/m 2 , 100W/m 2 as shown in Fig. 11(a). The V-I and V-P characteristics of each module in the string are plotted in Fig. 11(b), it's obvious that the modules currents levels are varied because of the PS. The series connection between these modules in a string as in Fig. 11(a) with bypass diode and blocking diode generates V-I, V-P characteristics with three levels (ladder shape) as in Fig. 11(c). To understand the strategy that followed to conducting these laddered characteristics, the following steps are reported. CWOA [12] PSO-WOA [12] SATLBO [45] HFAPS [46] MLBSA [38] TVACPSO [47] CPSO [47] GSA [51] CSA [40] ICSA [40] GA [48] CS [37] 8.0621e − 04 ± 4.3109e − 04 9.8602e − 04 ± 1.4485e − 09 9.98678e − 04 ± 6.33831e − 03 1.07101e − 03 ± 1.17001e − 03 9.86022e − 04 ± 2.3002e − 06 9.8602e − 04 ± − 9.8602e − 04 ± 9.1461e − 12    where Diff RM SE is equalled to ( RMSE lambert -RMSE of Eq. 6). where Diff RM SE is equalled to ( RMSE lambert -RMSE of Eq. 6).    -For zone 1 of Fig. 12(a),the current of M2 and M3 is lower than that of M1 as illustrated in Fig.  11(b). Then, the M1 current is the dominant in this zone while the other modules are bypassed. The current of the bypassed modules flow through their parallel diodes (conducting links are appeared by red line), as illustrated in Fig.12(a). Thus the current of the string in this zone is controlled by the condition of I SC M 2 ≤ I str Zone1 < I SC M 1 . The corresponding voltage of the string in this zone is calculated by Eq.(ps1). 29) where V str Zone1 is the voltage of the string of Zone 1, V M 1 is the first module voltage, and the V Bypassdiode2 and V Bypassdiode3 are the voltage values of the parallel-connected bypass diodes of the second and third modules, respectively. Finally, the V Blockingdiode refers to the voltage of the blocking diode. The voltage drop V D across the diode can be computed using Eq.30.
where V D is an implementation to the V Bypassdiode2 , V Bypassdiode3 and V Blockingdiode . The V f wd symbolizes to the diode forward voltage, the R refers to the diode equivalent resistance in the conducing mode and the I D is the current flows through the diode. For the ideal diode the V f wd and R are equaled 0. -For Zone 2 of Fig.12(b), M3 is generating a lower current than that of the other modules (M1 and M2 ) as illustrated in Fig.11(b). Thus, the currents of the M1 and M2 are the dominant ones in this zone; in contrast, the M3 is bypassed. Thus, the current of this zone is controlled by the condition of The corresponding voltage of the string in this zone is expressed as follows: -For Zone 3 of Fig. 12(c), the three modules of the string (M1,2,3) are generating currents (the conducting elements are colored in red). Thus the current of the string in this zone is controlled by the rule of 0 ≤ I str Zone3 ≤ I SC M 3 accordingly, the string voltage can be modeled through Eq.33.
The total harvested string power (P str ) is computed by multiplying the string voltage and current as reported in Eq. 33.
• To compute the voltage and current of the whole array composed of M parallel strings, the previously described steps are conducted for each string. Then, the current and voltage of the array are calculated as follows.  This part presents a case study for implementing the V-I and V-P characteristics of Kyocera Solar (KC200GT) based on the module identified parameters of section V-B by MRUN. Fig. 13 shows the V-I and V-P characteristics of the KC200GT PV interconnected strings and arrays under uniform irradiation levels.
The highlighted values of current, voltage, and power at short circuit point and maximum power one are illustrated in Fig.13(a) illustrate the high certainty of the identified model parameters that conducted to detect the short circuit and maximum power current, voltage, and power with nearly 99.9% accuracy as the expected and detected values of the short circuit and maximum current are (8.21 (A), 7.61 (A)), and (8.207 (A), 7.609 (A)) respectively. The expected value of the string voltage at the maximum power point is nearly equal to the summation of the voltage of the three modules that is nearly equal to 78.9 (3*26.3). The detected value of the voltage is 78.75 that is about 99.8% from the expected value. The accuracy of the identified parameters by MRUN is also affirmed from the figures of the 6x1 string, 9x1 string, and the considered arrays as shown in Fig. 13. For the non-uniform incident irradiation condition of Fig.14, the V-I and V-P characteristics are highly matched with the calculated current, voltage, and power based on equations of Section VI-A. For fig.15(a), a 3x1 string is subjected to three different radiations levels of 1000W/m 2 , 600W/m 2 , 200W/m 2 , accordingly the PV string characteristics has three ladders. The first ladder of the characteristics is controlled by the current of M1 that is subjected to 1000W/m 2 until the flowing current becomes equal to the generates current of M2 (600/1000 * 8.21 = 4.9260); in this step, the M2 current is the dominant one, and finally, M1 generated current (100/1000 * 8.21 = 1.6420) becomes the dominant one. The output voltage is the summation of the three modules voltage. The displayed curves of Fig. 14 are closely matched with these calculations that approve the accuracy of the identified parameters. For considering the change in the temperature and irradiation, the characteristics of PV string and arrays are exhibited in Fig. 15. Based on the previous discussions and observations, the MRUN is a recommended optimizer for identifying the 1D and 2D PV models parameters.

VII. CONCLUSION
Providing efficient PV solar cell/module modeling is the first essential step before the system installation. Therefore, this paper proposed an interactive modified version of an algorithm based on runge kutta optimizer (MRUN). The proposed algorithm performance has been enhanced to identify the PV single (1D) and double (2D) diode models parameters using several solar cells/modules datasets. The MRUN results have been validated through intensive comparisons with recent meta-heuristic optimization algorithms, including basic RUN, aquila optimizer (AO), electric fish optimizer (EFO), barnacles mating optimizer (BMO), capuchin search algorithm (CapSA), red fox optimization algorithm (RFSO) moreover recent twenty-five state-ofthe-art-techniques from literature. The power-voltage (V-P), voltage-current (V-I) characteristics, the absolute error between the estimated and measured datasets, the mean convergence curves, and the fitness function (RMSE throughout a set of separate runs are utilized for the comparison purpose. Furthermore, a group of statistical metrics and nonparametric analysis of the Wilcoxon signed-rank test have been computed. The comparisons and analyses reveal the quality enhancements achieved via applying the proposed RUN optimizer in identifying the PV modules parameters compared to the basic Run optimizer and other recent stateof-the-art techniques. The accuracy, consistency, and convergence properties of the proposed optimizer are remarkable points in the results. The efficiency of the extracted parameters has been justified in modeling large interconnected solar systems in series, parallel or series-parallel schemes; several sizes of strings and PV arrays have been considered under several operating conditions of temperature and radiation. The attained V-P and V-I characteristics of the studied strings and arrays divulge the efficient emulating for the physical performance of the considered systems.