A single diode model can describe the static behavior of a PV module and is proposed in many kinds of literature such as in24,25, the electric circuit for the single diode model is shown in Fig. 2.

From Fig. 2, the mathematical equation which describes the static PV model can be driven as follow:

Where \({I}_{ph}\) is the photocurrent, \({I}_{sat}\) is the reverse saturation current of the diode, q is the electron charge, A is the diode ideality factor, K is the Boltzmann constant, T is the PV module temperature, \({R}_{s}\) is the series resistor, \({R}_{sh}\) is the shunt resistor, and \({N}_{s}\) is the number of connected PV cells in series, the photocurrent (\({I}_{ph}\)) which is produced by the PV module is influenced by the solar irradiance falling on the PV module plane and its influenced by the ambient temperature of the PV module as shown below26,27:

Where \({I}_{phs}\) is the photocurrent at standard conditions (STC), \({I}_{rrs}\) is the irradiation at STC and is 1000 \(\frac{W}{{m}^{2}}\),\({ K}_{i}\) is the temperature coefficient, \(\varDelta T=T-{T}_{s}\) is the difference between the PV temperature (\(T\)) and STD temperature (\({T}_{s}=25℃\)), from equations (4) and (5) the output current (\({I}_{PV}\) ) depends on the amount of sunlight and the PV temperature, the characteristics curves(P-V), and(I-V) for different irradiance and temperature are shown in Fig. 2,

**PV System Modeling**

The modeling process for dynamical systems is used to drive mathematical representation for the system, it depends on two techniques: mathematical representation is derived from physics or measured input-output data, in general, the modeling process can be organized into three techniques: white box, grey box, and black box modeling28, in white box modeling the physics of the system, are completely known; and the mathematical model can be constructed entirely from prior knowledge and physical insight29, in grey box modeling, some of the physical parameters are known, the other parameters are determined from measured input-output data, but the black box modeling is used when the physical properties of the model are unknown, and the model is constructed from measured input-output data, but the chosen model structure belongs to the candidate models that are known to have good flexibility and have been successful in the past 30. The back box modeling is chosen for the system identification methodology in this paper, in the back box modeling method, many models are compared with the real data, and the best is determined by an optimization criterion. Four basic stages characterize the system identification technique (I) The measured input-output data, (II) the proposed model, (III) a criterion of fit and parameter estimation, and (IV) Validate the model which is obtained.

**4.1. The measured input-output data.** The data used in system identification is obtained from the PV system, shown in figure .1, where the input data is the PV current and output data is the PV system output power.

**4.2. The proposed model structure**. a fractional order model is considered in this paper, A generalized fractional-order dynamic system can be represented by the next Eq. 31:

The dynamic system in (6) can be written in the next form:

$$\sum _{k=0}^{n}{a}_{k}{\mathcal{D}}^{{\alpha }_{k}}\mathcal{Y}\left(t\right)= \sum _{k=0}^{m}{b}_{k}{\mathcal{D}}^{{\beta }_{k}}\mathfrak{u}\left(t\right) \left(7\right)$$

Where( \({a}_{k}\),\({ b}_{k})\in \mathcal{ }\mathcal{R}\), and (\({\alpha }_{k},{\beta }_{k}\)) \(\in {\mathcal{R}}_{+}^{2}\). in (6) if all the derivation orders are integer multiples of a base order q such that \({\alpha }_{k}\), \({\beta }_{k}=qk\), \(q\in {\mathcal{R}}_{ }^{+}\), the system is said to be of commensurate order. If \(q=\frac{1}{\mathcal{g}}\), \(\mathcal{g}\in {\mathbb{Z}}_{+}\), the system is said to be of rational order.

The Laplace transform of (7) can be written in the form:

$$\sum _{k=0}^{n}{a}_{k}{s}^{{\alpha }_{k}}Y\left(s\right)= \sum _{k=0}^{m}{b}_{k}{s}^{{\beta }_{k}}U\left(s\right), \left(8\right)$$

Then the fractional order transfer function can be written in the form:

$$G\left(s\right)=\frac{Y\left(s\right)}{U\left(s\right)}=\sum _{k=0}^{m}{b}_{k}{s}^{{\beta }_{k}}/\sum _{k=0}^{n}{a}_{k}{s}^{{\alpha }_{k}}, \left(9\right)$$

By considering the system is commensurate order *q*, We can write \(\text{ɦ}={s}^{q}\), then the continuous time rational transfer function is:

$$H\left(\lambda \right)= \frac{\sum _{k=0}^{m}{b}_{k}{ɦ}^{k}}{\sum _{k=0}^{n}{a}_{k}{ɦ}^{k}} . \left(10\right)$$

In fractional calculus the differential operator has many definitions, The three most frequently used definitions for the general fractional differ integral are the Grunwald-Letnikov (GL) definition, the Riemann-Liouville (RL), and the Caputo definitions32, Grünwald-Letnikov definition is considered in this work for the drive numerical solutions to fractional-order differential equations, which is given by the next equation:

$${}_{\text{a}}{}^{\text{ }}{\text{D}}_{\text{t}}^{ \alpha }=\underset{h\to 0}{\text{lim}}\frac{1}{{h}^{\alpha }} \sum _{j}^{\left[\frac{t-a}{h}\right]}{\left(-1\right)}^{k}\left(\begin{array}{c}\alpha \\ k\end{array}\right)f\left(t-kh\right), \left(11\right)$$

Where \(h\) is the step size, and the binomial coefficient can be expressed with Gamma function \(\varGamma \left(.\right)\)in the following expression:

$$\left(\begin{array}{c}\alpha \\ k\end{array}\right)=\frac{\alpha !}{k!\left(\propto -k\right)!}=\frac{\varGamma \left(\alpha +1\right)}{\varGamma \left(k+1\right)\varGamma \left(\alpha -k+1\right)} ,$$

12

Laplace transform of the Grünwald-Letnikov fractional operator is:

$$\mathcal{L}\left[{\mathcal{D}}^{\alpha }f\left(t\right)\right]={S}^{\alpha }F\left(S\right)$$

13

.

Oustaloup’s recursive filter is an approximation for fractional differentiator and integrator, in a specified frequency range (\({\omega }_{b},{\omega }_{h}\)) and is capable of producing a very good fitting to the fractional-order elements 33,34, the algorithm for Oustaloup’s recursive filter is shown in the following form:

$${S}^{\alpha }\approx \prod _{K=1}^{N}\frac{S+{\omega }_{k}^{\text{'}}}{S+{\omega }_{k}} \left(14\right)$$

Where

\({\omega }_{k}^{\text{'}}= {\omega }_{b}* {\omega }_{u}^{\frac{2k-1-\alpha }{N}}\) , \({\omega }_{k}= {\omega }_{b}* {\omega }_{u}^{\frac{2k-1+\alpha }{N}}\), \(K={\omega }_{h}^{\alpha }\), \({\omega }_{u}= \sqrt{\frac{{\omega }_{h}}{{\omega }_{b}}}\)

And *N* is the approximation order in the frequency range (\({\omega }_{b},{\omega }_{h}\))

**4.3. Criterion of fitting and parameters estimation**

The system identification method depends on the error minimization between model output and system output, the least-square method is used in this work, and the problem is expressed as follows:

$$\mathcal{F}\left(x\right)=\sum _{p=1}^{n}{℮}_{p}$$

15

Where \({℮}_{p}={y}_{p}-{\widehat{y}}_{p}\) and \({\widehat{y}}_{p}=\widehat{\mathcal{H}}\left({\mu }_{p},\theta \right)\) is the output of the estimated model\(\widehat{\mathcal{H}}\) for the input signal (\({\mu }_{p}\)), \(p\) is the iteration, and \(\theta\) is the estimated parameters of the model \(\widehat{\mathcal{H}}\), from Eq. (6) the model parameter\(\theta\) is formed by:\({a}_{r}=\left[{a}_{0} {a}_{1}\dots { a}_{m}\right]\), \({\alpha }_{r}=\left[{\alpha }_{0} {\alpha }_{1}\dots { \alpha }_{m}\right]{b}_{\mathcal{F}}=\left[{b}_{0} {b}_{1} \dots {b}_{n}\right]\), \({\beta }_{\mathcal{F}}=\left[{\beta }_{0} {\beta }_{1} \dots { \beta }_{n}\right]\).

In this work the FOMCON toolbox for MATLAB/Simulink is used, the toolbox is based on fractional calculus for system modeling and control design. The FOMCON toolbox uses the Levenberg-Marquardt algorithm for optimization of the model parameters θ under the nonlinear least squares technique. The Levenberg-Marquardt algorithm is described by the Eq. (16) According to the references21,35,36 :

$${-J}_{p}^{T}{℮}_{p}=\varDelta {\theta }_{p} \left({J}_{p}^{T}{J}_{p}+\lambda I\right)$$

16

Where \(\lambda\) is a positive scaler, called a damping factor, \(I\) is the identity matrix, and \(J\) is the Jacobian matrix of the function \(\mathcal{F}\) and is defined as

$$J=\left[\begin{array}{ccc}\frac{\partial {℮}_{\text{1,1}}}{\partial {\theta }_{1}}& \frac{\partial {℮}_{\text{1,1}}}{\partial {\theta }_{2}} \dots & \frac{\partial {℮}_{\text{1,1}}}{\partial {\theta }_{N}}\\ \begin{array}{c}\frac{\partial {℮}_{\text{1,2}}}{\partial {\theta }_{1}}\\ \dots \end{array}& \begin{array}{c}\frac{\partial {℮}_{\text{1,2}}}{\partial {\theta }_{1}}\\ \dots \end{array} \dots & \begin{array}{c}\frac{\partial {℮}_{\text{1,2}}}{\partial {\theta }_{N}}\\ \dots \end{array}\\ \frac{\partial {℮}_{1,m}}{\partial {\theta }_{1}}& \frac{\partial {℮}_{1,m}}{\partial {\theta }_{1}} \dots & \frac{\partial {℮}_{1,m}}{\partial {\theta }_{N}}\end{array}\right] \left(17\right)$$

the update parameters are defined in Eq. (18).

$${\theta }_{p+1}={\theta }_{p}+\varDelta {\theta }_{p}$$

18

From equations (16) and (18) the update parameters can be determined by Eq. (19):

$${\theta }_{p+1}={\theta }_{p} - \frac{{-J}_{p}^{T}{℮}_{p}}{\left({J}_{p}^{T}{J}_{p}+\lambda I\right)}. \left(19\right)$$

**4.4. The model quality**

The quality of the identified model will be discussed in this section. let \({y}_{s}\) denote The output data of the system, and \({y}_{d}\) denote The output data of the identified model, for single input single output (SISO)system, both \({y}_{s}\) and \({y}_{d}\) are vectors of dimension *N X 1*, and the difference between \({y}_{s}\) and \({y}_{d}\) is the model output error since \({y}_{s}\) and \({y}_{d}\) are vectors, then the error is a vector of size *N X 1* which is called residuals

$$e= {y}_{s}-{y}_{s}$$

20

.

The percentage fit can be described as

$$fit=\left(1-\frac{\parallel e\parallel }{{y}_{s}-\stackrel{-}{{y}_{s}}}\right)*100\% , \left(21\right)$$

where \(\parallel .\parallel\) is the Euclidean norm, and \(\stackrel{-}{{y}_{s}}\) is the mean value of the system output data\({y}_{s}\)

**4.5 stability analysis**

To study the fractional system stability, we consider Matignon’s stability theorem37 :

**Theorem 1**

(*Matignon’s stability theorem*) The fractional transfer function G(s) = Z(s) = P (s) is stable if and only if the following condition is satisfied in the σ-plane:

$$\left|\text{arg}\left(\sigma \right)\right|>q\frac{\pi }{2} , \forall \sigma \in \complement , P\left(\sigma \right)=0$$

22

,

*where* \(0 < q < 2\)*and* \(\sigma ={s}^{q}\). *When* \(\sigma =0\) *is a single root of* \(P \left(s\right),\) *the system cannot be stable. For q* = 1, *this is the classical theorem of pole location in the complex plane: no pole is in the closed right plane of the first Riemann sheet.*

The algorithm for determining the system stability in (10) is shown in figure (4), according to Matignon’s stability theorem, the fractional order system stability in the complex plane is shown in figure.5.