Parameter and Order Identiﬁcation of Fractional Systems with Application to a Lithium-Ion Battery

: This paper deals with a method for the parameter and order identiﬁcation of a fractional model. In contrast to existing approaches that can either handle noisy observations of the output signal or systems that are not at rest, the proposed method does not have to compromise between these two characteristics. To handle systems that are not at rest, the parameter, as well as the order identiﬁcation, are based on the modulating function method. The novelty of the proposed method is that an optimization-based approach is used for the order identiﬁcation. Thus, even if only noisy observations of the output signal are available, an approximate identiﬁcation can be performed. The proposed identiﬁcation method is, then, applied to identify the parameters and orders of a lithium-ion battery model. The experimental results illustrate the practical usefulness and verify the validity of our approach.


Introduction
In recent years, fractional models have been increasingly used to describe a broad range of real world problems [1]. Fractional analysis provides models that accurately describe the long-term phenomena of systems as the fractional operators take into account the entire system's past. In contrast to the integer order operators, the information of the system's past is not concentrated in the initial conditions; however, the system's past is represented by a non-constant function [2]. In the context of parameter and order identification, this non-constant function has to be considered, otherwise it leads to incorrect parameter and order identification [3].
To overcome this drawback, the majority of approaches for parameter and order identification assume the system to be at rest (see e.g., [4][5][6][7][8]). This means that the input and output signal are assumed to be zero before the start of the identification. However, fractional models are primarily used to describe long-term phenomena, such as the diffusion process in a lithium-ion battery. For a lithium-ion battery to be considered at rest, this diffusion process must have completely decayed, which can take several hours [9]. Therefore, the assumption that the system is at rest is restrictive in terms of applicability.
In [3,10], two approaches were proposed that dropped the assumption of a system at rest. The authors of [3] applied the Caputo definition for the fractional derivatives. Integrating the underlying fractional differential equation with the highest derivative order leads to the initial values in the system description. These initial values are interpreted as additional parameters of the system and collected in an extended parameter vector. Under the assumption that the fractional orders are known, a method for the bias-free identification of the extended parameter vector under noisy observations of the input and output signal is presented.
A shortcoming of the method from [3] is that the initial values of the Caputo derivatives do not initialize the system properly [11]. Instead, as with the Riemann-Liouville definition, the entire past of the input and output signal has to be taken into account to initialize the Caputo derivative properly [12]. This is considered in the approach of [10], which is based on the findings from [5]. In opposite to [3], not only the parameters but also the orders of the fractional system are identified in the method from [10]. The approach is derived from the modulating function method, which allows elimination of the influence of the entire system past. Even though the approach is demonstrated with a noisy observation of the output signal, the orders are determined by solving a root-finding problem. Thus, the noisy observation of the output signal leads to an ill-posed problem.
To determine the orders in the presence of noise, an optimization-based approach is proposed in this paper. The approach is also based on the modulating function method to keep the possibility of eliminating the influence of the entire system's past. For the sake of practical implementation, the Grünwald-Letnikov definition is used instead of the Riemann-Liouville definition.
The remainder of this paper is structured as follows: First, Section 2 outlines the preliminaries. In the first part of Section 3, the parameter identification is briefly introduced. The second part of Section 3 presents the order identification under noisy observations of the output signal, which is the first main contribution of the paper. The application to a lithium-ion battery and the identification of the parameters and orders are the second main contribution, which is described in Section 4. In Section 5, we draw the main conclusions from this paper.

Preliminaries
In the literature, many definitions for the fractional operator exist. With regard to practical implementation, the presented approach is based on the definition of Grünwald-Letnikov. Further, the fractional model and the identification equation on which the approach is based on are given in this section.
The definition of the left-sided Grünwald-Letnikov derivative is where T is the sampling time.
The right-sided Grünwald-Letnikov derivative is defined as: where T is the sampling time.
Following the notation of [16], the fractional derivative operators are given by d α t 0 t . The derivative order α is given as the right upper index d α . The lower left index d t 0 for the left-sided case or the lower right index d t 3 for the right-sided case marks the time before t 0 or after t 3 the function is zero. The third parameter is the actual time t and is given as lower right respectively left index. The derivative of the right-sided Grünwald-Letnikov definition with respect to the order is given in the following lemma, which was originally stated in [10].

Lemma 1.
Derivative of the Right-Sided Grünwald-Letnikov Definition. The derivative of the right-sided Grünwald-Letnikov definition is The proof can be found in [10].

Remark 1.
Assuming a continuous and bounded function h and a finite sampling time T, then the derivative of the right-sided Grünwald-Letnikov definition (3) does not have singularities and is bounded.

Fractional Model
Consider the general non-commensurable fractional model where u(t) is the time-continuous input signal andỹ(t) is a noisy observatioñ of the noise-free output signal y, which is superimposed by a white noise with a zero mean. The derivative operator D indicates that the system is not at rest. RL marks that the system is interpreted in terms of Riemann-Liouville [13][14][15]. The unknown parameters are collected in the vector p = [a n , . . . , a 0 , b m , . . .
where a i , b k ∈ R, and the fractional orders in the vector where α i , β k ∈ R, α i > 0 and β k > 0. Without loss of generality, the fractional orders are assumed to be ordered 0 ≤ α 0 < · · · < α n , 0 ≤ β 0 < · · · < β m with α n ≥ β m . The number of unknown parameters n, m ∈ N are assumed to be known. Regarding the noise , the input signal u, the output signalỹ, and the parameters p, the following assumptions are made.

Assumption 2. Bounded Signals.
The input u and output signalỹ are assumed to be bounded.

Assumption 3. Parameter Normalization.
Without loss of generality, the normalization a n = 1 holds in (4).

Simulation of a Fractional Model
In the parameter identification method, a simulation y sim of the system output signal y is needed. To simulate a fractional model (4), a closed-form solution with zero initial conditions is given in [18] (p. 113). In general, if the system is not at rest and, hence, no zero initial conditions are present, an error between the simulated output signal with and without zero initial conditions occurs [19]. To reduce this error, an extension based on a closed form solution that considers the latest past of the system in sense of the short-memory principle [20] is proposed in [19].

Definition 3. Closed-Form Solution Considering the Short-Memory Principle.
A closed-form solution that considers the short-memory principle is defined as: where T is the sampling time, t 1 is the time where the measurement starts, t 2 = t 1 + LT is the time where the simulation starts, and L is the memory length.

Remark 2.
For the definition 3, it is assumed that the system is at rest before t 0 .
The input signal u and the output signal of the system y in t ∈ [t 1 , t 2 ) are used to initialize the closed-form solution (10). From t = t 2 onward, the simulation of the system output signal y sim is calculated recursively.

Modulating Function Method and Identification Equation
The modulating function method was first proposed in [21]. The basic idea is to transfer the derivatives of the initial signals to an arbitrary and continuously differentiable function by means of integration by parts. Therefore, the fractional model (4) is multiplied with a function that is continuously differentiable. The resulting products of the function and the input as well as output signal are integrated over t ∈ [t 1 , t 2 ]. Applying integration by parts leads to the result that the derivatives are swapped but at the expense of arising boundary terms.

Definition 4.
Modulating Function [21]. With regard to a fractional model (4), a modulating function is a function satisfying The property (P1) ensures that all necessary derivatives of the modulating function exist. The property (P2) eliminates the boundary terms which arise due to the integration by parts. Considering a third property, viz.
which is stated in [22]. The identification of the modulating function method for a fractional model with initialized derivative operators results in The derivatives of the model (4) are considered to be described by the fractional operator of Riemann-Liouville. Applying the modulating function method, the derivative type changes to the Caputo definition for the modulating function [22]. Therefore, the approach can be formulated using the Grünwald-Letnikov derivative, and the Caputo definition has to be replaced by the Grünwald-Letnikov definition. If the modulating function is α n + 1-times continuously differentiable, the derivative operator of Riemann-Liouville and Grünwald-Letnikov are interchangeable [23]. The Riemann-Liouville derivative is linked with the Caputo derivative through ( [24], p. 91) which is the identification equation for the approach of this paper.

Combined Iterative Parameter and Order Identification
In this section, a method for parameter and order identification is introduced. The novelty is that the method can handle noisy observations of the output signal. The identification is performed iteratively whereby the orders are updated in a first step, and afterwards the parameters are identified for the orders of the actual iteration. This procedure is shown in Figure 1. For a consistent parameter estimation under noisy measurements, we apply the approach from [19]. This approach is outlined in Section 3.1. The focus of the paper is on the order identification under noisy observations of the output signal. Instead of solving a root-finding problem, as in [10], an optimization-based approach is used to identify the derivative orders of a fractional model (4). This main contribution of the paper is given in Section 3.2.

Parameter Identification
The method introduced in [19] is an extension of the instrumental variable method for integer order models to fractional order models. To apply the instrumental variable method, we require N p ≥ N identification equations, where N = n + m + 1 is the number of model parameters. The N p identification equations are generated by evaluating the integrals in (13) for different time intervals. In general, the lower bound of the integral can be expressed by t 1,h and the upper bound by t 2,h where h ∈ 0, 1, . . . , N p − 1 . These time points can be derived, e.g., by shifting the integrals by a fixed time T ∆ . Assuming that the identification process is started at t = t 1 and the fixed integral width is T I , the lower and the upper bound are given by t 1,h = t 1 + hT ∆ and t 2,h = t 1,h + T I , respectively. initial order θ 0 parameter identification p θ q order identification θ q convergence? yes no parameter identification p (θ 0 ) q = q + 1 is set up for the parameter identification. The parameters depend on the latest estimates of θ q . The vector Y := y 0 , y 1 , . . . , y N p −1 describes the α n -th fractional derivative of the modulated output signal The measurement matrix M := m 0 , m 1 , . . . , m N p −1 consists of the vectors m h ∈ R N×1 , which are the remaining modulated derivatives of the output and input signal: where h ∈ 0, 1, . . . , N p − 1 . Due to the noisy observations of the output signal, the leastsquares method leads to biased estimates of the parameters [25]. Therefore, a matrix of instrumental variables W := w 0 , w 1 , . . . , w N p −1 where is used to solve the linear system (14): Instead of the noisy observations of the system output signalỹ, the simulation of the output signal y sim is used to set up (17). Thus, the instrumental variables fulfill the requirements of the instrumental variable method in that they are uncorrelated to the noise but strongly correlated to the undisturbed input signal and output signal of the system. The approach can be summarized as follows: 1.
The parameters are estimated using the least-squares method after N iterations where N is the number of parameters for the first time.

2.
These parameters are used to simulate the output signal y sim with the current fractional orders θ q . In general, the entire past of the system is not available. Thus, the closed-form solution based on the short-memory principle (10) is used to calculate y sim . 3.
The calculated output signal y sim is used to set up the instrumental variable matrix W and to refine the parameter estimates iteratively via (18) until a maximum number of iterations N p is reached or the difference between the two estimates is almost zero.

Order Identification
For the parameter identification, we assumed that the derivative orders are available in each iteration. In this subsection, an order identification method is proposed. Therefore, the identification equation of the modulating function (13) is interpreted as a nonlinear function that depends on the orders of the system For a noisy output signal, the nonlinear equation for the order identification results in In the sequel, we first investigate the influence of the noise on the order identification. Afterwards, an optimization-based method for the order identification is proposed, and the combination with the parameter identification is outlined.

Ill-Posed Root-Finding Problem
For the original parameters a i and b k and the orders α i and β k of the model, f o is zero. In [10], this fact is used to transfer the order determination into a root-finding problem.
In contrast to f o , in general, f is not zero for the original parameters and orders but takes the value It is possible that, for small deviations of the original parameters and orders, f is zero. However, the possibility also exists that f is always different from zero. In this case, an identification of the orders using root-finding algorithms is not possible. This is exemplified in Figure 2. A system described by where a 0 = 10, b 0 = 1 and α = 0.8 is excited with a pseudo random binary signal. The identification Equation (20) is evaluated in a brute-force manner for α = {0.5, 0.01, . . . , 0.9}, while the parameters are kept at a 0 = 10 and b 0 = 1. The first evaluation is performed with the noise-free output signal, which is represented by the blue line in Figure 2. The red line, on the other hand, is obtained from an output signal, which is superposed by a time-variant noise ε. The root disappears in this case, and thus the order can not be derived by the approach proposed in [10].

Optimization-Based Approach
In order to identify at least the derivative orders in the presence of noisy observations of the output signal approximately, an objective function where N θ ≥ n + m + 2 is defined. As can be seen, the objective function is formulated in the sense of best fitting. In (24), the function To identify the orders, the objective function has to be minimized with respect to the orders: min The optimization problem (26) is solved using the Gauß-Newton method ( [26], S. 214 ff.).
where ∇J is the gradient: and H J is the approximated Hessian of the Gauß-Newton method:

Convergence Analysis
In this subsection, the derivative of the order identification Equation (20) is calculated first. This derivative is needed for gradient (28) as well as Hessian (29). Secondly, the convergence of the Gauß-Newton method for the stated problem (26) with (24) is analyzed. (25). Suppose N o = n + m + 2, J = { j ∈ N|j ≤ N o }, j ∈ J and an independent equation f h (25) of f in (20). Further, let {θ} j describe the j-th element of the order vector θ in (7).

Lemma 2. Partial Derivative of Order Identification Equation
The partial derivative of f h with respect to {θ} j is (30) where In (31), the partial derivative of the Grünwald-Letnikov derivative is given in Lemma 1.
Proof. If (25) is derived with respect to {θ} j , it must be considered that the parameters as well as the derivative of Grünwald-Letnikov depend on the orders. To calculate the derivative, the product rule is used. While the parameters depend on all orders and are stated in (30), the Grünwald-Letnikov derivatives only depend on one specific order. Therefore, a case distinction depending on the derivative order {θ} j has to be performed, which is given in (31). The needed derivative of the Grünwald-Letnikov derivative with respect to the order is given in Lemma 1.
The derivatives of the parameter ∂a i (θ) /∂{θ} j and ∂b k (θ) /∂{θ} j in (30) depend on the specific method for the parameter identification. In this paper, the instrumental variable method from [10] is used to achieve a consistent parameter estimation. Based on the parameter identification Equation (18), the corresponding partial derivative of the parameters is calculated. Lemma 3. Partial Derivative of the Parameter. Suppose N = n + m + 1, N p ∈ N, N p ≥ N and the parameter identification Equation (18): where M(θ) = m 0 (θ), m 1 (θ), . . . , m N p −1 (θ) is the measurement matrix with Further, let {θ} j describe the j-th element of the order vector θ (7). The partial derivative of p with respect to {θ} j can be calculated as Proof. Starting from the linear system (14), the first step is the multiplication of the linear system (14) with the transposition of the instrumental variable matrix W W (θ)M(θ)p(θ) = W (θ)Y(θ).
Deriving (37) with respect to {θ} j yields To obtain this equation, the product rule is used. Equation (38) is rearranged such that all expressions except for the expression depending on the partial derivative of the parameter are on one side Under the assumption of independent equations for the parameter identification problem, the inverse of W (θ)M(θ) must exist. Otherwise, the requirements for the parameter identification are not met. The inversion of W (θ)M(θ) is the last step to calculate the derivative of the parameters.
Due to the computational effort of the non-recursive instrumental variable method, often the recursive realization is used in practical applications ( [17], p. 269). In this case, (36) cannot be calculated, and the difference quotient is used to approximate the derivative of the parameters. As an example, the difference quotient for the κ-th element of the parameter vector p and the j-th element of the order vector θ result in where κ = 1, 2, . . . , n + m + 1, j = 1, 2, . . . , n + m + 2. In (40), ∆ ∈ R >0 is a sufficiently small constant that represents the difference interval. The vector e j is the j-th unit vector of R n+m+2 and is used to select the element of the vector θ. Under the assumption that the Hessian H J is non-singular, the Hessian of the Gauß-Newton method is a positive definite matrix ( [26], S. 215). Hence, if the derivatives ∂ f h (θ q )/∂ {θ} j exist and ∇J has no singularity, the Gauß-Newton method converges.

Lemma 4. Convergence Analysis.
The derivatives ∂ f h (θ q )/∂ {θ} j (30) exist and ∇J (28) has no singularity if (a) the sampling time T is finite, and (b) the modulating function γ, (c) the input signal u, the output signalỹ, the parameters a i and b k , and the derivatives of the parameters ∂a i (θ) /∂{θ} j and ∂b k (θ) /∂{θ} j are bounded.
Hence, if these requirements are fulfilled, the Gauß-Newton method converges against a minimum for the stated problem (26) with (24).
Proof. If the sampling time is equal to zero, the partial derivative of the Grünwald-Letnikov definition (3) has a singularity. Due to this, the sampling time must be restricted to finite values, which is ensured by requirement (a). If (30) is not bounded, the inverse of the Hessian does not exist. For (30) to be bounded, all terms of (30) must be bounded, which is covered by the requirements (b) and (c).
With regard to practical applications, the requirements (a)-(c) are not restrictive. The requirement (a) is always fulfilled for implementation on processing units. The fulfillment of requirement (b) can be ensured by the choice of the modulating function without any restriction by the properties (P1) and (P2). The boundedness of the input signal, output signal, parameters, and the derivatives of the parameters regardless of the chosen method is guaranteed by the practical Assumption 2, which fulfills requirement (c) [19]. However, the Gauß-Newton method generally does not ensure that the global minimum is reached. This results from the fact that the Gauß-Newton method is based on a linear approximation of the original problem. The minimum the Gauß-Newton method converges against depends on the chosen initial solution and on the specific nonlinearity ( [26], S. 214 ff.).

Parameter and Order Identification for a Model of a Lithium-Ion Battery
In this section, the proposed parameter and order identification approach is applied to a lithium-ion battery model. To this end, it is crucial that the modulating function satisfies property (P3) as, otherwise, the system is required to be at rest. Considering the lithiumion battery, the diffusion process has to decay, which can take several hours. The choice of such a modulating function is stated in Section 4.1. In Section 4.2, the fractional model of a lithium-ion battery is given, and the experimental setup is explained in Section 4.3. The ground truth for the identification is determined in Section 4.4. The identification is carried out in Section 4.5.

Choice of Modulating Function
One advantage of the approach from Section 3 is that the system need not to be at rest at the beginning of the identification process. This advantage is at the cost of the modulating function fulfilling property (P3) in addition to (P1) and (P2). In [22], they proved that the spline-type modulating function also fulfilled property (P3).

Definition 5.
Spline-Type Modulating Function [27]. The spline-type modulating function is the integration over a weighted sequence of impulses where s ∈ N is the number of splines, and ν ∈ N represents the order of the modulating function.

Model of a Lithium-Ion Battery
Several models to describe a lithium-ion battery exist in the literature. Depending on the modeling goal, the range is from detailed electrochemical models to behavior models, see [28,29]. The electrochemical models are based on particle models, and partial differential equations are used, whereas the behavior models are described by RC circuits. Fractional models are a compromise between the electrochemical models and behavior models, see [30] (p. 141 ff.) and [16,29,31,32]. In general, the fractional impedance model of a lithium-ion battery is used to describe the overvoltages of the different loss processes and the overvoltage of the diffusion process in a lithium-ion battery [31].
The model where s depicts the complex number frequency parameter and α diff ∈ (0, 1) is the fractional order used in this paper is based on [29]. Due to the minimal sampling time T min = 0.07 s is realizable with the experimental setup, which is explained in the next subsection, and only the low frequency range of the model proposed in [29] can be correctly reconstructed. The loss processes of the high and mid frequency range are collected inR 0 . The differential capacity C diff describes the amount of charge that has to be removed or inserted to change the open circuit voltage (OCV) by a specified quantity [33]: For the parameter and order identification, the model (42) has to be rearranged so that it appears in the form of (4). First, the fractions are extended Second, using the Laplace transform ( [13], p. 79) yields Comparing (46) with (4), the parameters and order to identify are

Experimental Setup
Due to safety reasons, the lithium-ion battery was operated in a climate chamber VC 3 4018 from Vötsch. The lithium-ion battery was excited with a BOP 20-20M from Kepco, and the excitation was realized as a current excitation. The current was measured with a digital multimeter 34410A from Agilent, which had a standard deviation of σ i = 1.1 mA. The voltage was measured with a DS2004 High-Speed A/D Board from dSpace. The DS2004 High-Speed A/D Board was integrated in a real-time system and had a standard deviation of σ u = 0.14 mV.
The data transfer of the multimeter to the host computer took T t = 0.05 − 0.07 s. While the data were transferred, no new data could be recorded from the multimeter. To avoid an incorrect parameter and order identification due to measurement losses, the sampling time was chosen to be T = 0.1 s.
The measurement was recorded in a voltage-correct circuit, which is shown in Figure 3. With this, it is taken into account that the parameters and orders of the fractional model of the lithium-ion battery (44) depended on the OCV and were only valid for small-signal excitation ( [9], p. 63 f.). To achieve a consistent parameter identification, the measured voltage had to be corrected by the instantaneous OCV ( [17], p. 255). A current-correct circuit leads to systematic errors in voltage measurements, such that the correction of the measured voltage is also incorrect. Even if the voltage across the ammeter is very low, the voltage drop suggests a great change in the state of charge (SOC), especially in the range of SOC ∈ [40%, 60%]. This is due to the flat course of the SOC-OCV curve, see Figure 4.

Experiment and Ground Truth
The fractional battery model (44) is only valid in a single operating point, which depends on the temperature and state of charge. To prevent any changes from the operating point due to temperature changes, the temperature was kept constant at 20°C inside the climate chamber. Thus, the state of charge was not varied during the measurement. A pseudo random binary signal with a pulse duration of T P = 10 s with the following two properties was used. The first property concerns the peak-to-peak amplitude, which was chosen asî = 400 mA. Next to this, the pseudo random binary signal was designed such that it was zero mean over the entire measuring range.
lithium-ion battery  The measured input signal is shown in Figure 5, and the measured output signal of the lithium-ion battery is shown in Figure 6. In both figures, a dashed red line indicates the point in time from which the data were used for the combined parameter and order identification. The data before T B = 80 s were only available for the determination of the reference parameters.
The lithium-ion battery was not excited for the first T R = 20 s to check if the lithiumion battery was at rest and to determine the instantaneous OCV. As the standard deviation of the measurement σ u,mea = 0.146 mV is comparable to the standard deviation of the manual σ u = 0.14 mV, we assumed that the lithium-ion battery was at rest. The OCV was calculated as mean value of the first T R = 20 s and was In combination with Figure 4, the SOC was determined to be SOC = 34.9%. Considering (43), the corresponding differential capacity can be derived from the SOC-OCV curve. The differential capacity is equivalent to the reciprocal of the derivative of the SOC-OCV curve. With the given SOC-OCV curve (cf. Figure 4) and SOC value, the differential capacity was determined as The parameter to be identified is The parameter b 1 =R 0 = 0.039 Ω and the order α = 0.395 are estimated by a complex nonlinear least squares regression. As an objective function, the quadratic error between a simulated output signal and the measured output signal is used. For the simulation, the complete input signal is used. Thus, we ensured that the system was at rest and no error occurred in the parameter and order due to neglecting a part of the system's history.

Results of the Combined Parameter and Order Identification
The instrumental variable method as described in Section 3.1 was applied for the parameter identification. To set up the instrumental variables, a simulated output signal y sim is needed. The closed-form solution (10) considering the short-memory principle was used so that a strong correlation was achieved with the undisturbed output signal, but the whole past of the system did not have to be taken into account. The memory length was chosen to be L = 200. The less computationally expensive recursive instrumental variable method instead of the non-recursive method was applied. Hence, an initial covariance matrix must be provided. As two parameters have to be identified, the covariance matrix is a 2 × 2 matrix and is chosen to be The values in (50) are motivated by the very small reference parameters. Before the identification can be performed, the spline-type modulating function (41) has to be parameterized. The number of splines was chosen to be s = 5, and the order was ν = 3. The identification horizon was T I = 40 s. Further, to generate independent Equation (25), the start of the integral is shifted by T ∆ = 4 s for each equation.
For the first parameter identification q = 0 and, hence, before the iterative parameter and order identification is started (see Figure 1), the first two iterations within the parameter identification are performed in sense of the least-squares method as described at the end of Section 3.1. After the first two iterations, the output signal is simulated using the actual estimates of the parameters and the instrumental variables are calculated for the parameter identification. Within the iterative parameter and order identification (q > 0), the parameters of the iteration before (q − 1) are used to set up the instrumental variables for the first iteration of the parameter identification. Afterwards, the updated estimates are always used.
As in the case of parameter identification, the identification Equation (25) of the order identification is also based on the modulating function method. The spline-type modulating function with the same parameterization as for the parameter identification is used. The initial value of the order is specified to be α diff,0 = 0.8. Since the identification is performed in the non-recursive form, the derivative of the parameters can not be calculated analytically. Instead, the numerical approximation (40) is used. The deviation from the derivation order of the actual iteration, which is required to evaluate the difference quotient, is set to ∆ = 0.001 (see (40)).
The iterative character of the parameter identification and of the combined parameter and order identification approach makes stopping criteria necessary. In this paper, two stopping criteria were formulated for each iteration. One stopping criterion is the maximum number of iterations q max . The maximum number of iterations is limited by the measurement duration T D = 120 s (cf. Figures 5 and 6), the memory length L in combination with the sampling time T, and the chosen identification horizon T I and shifting time T ∆ of the modulating function method. Using the chosen data, the maximum number of iterations is A second stopping criterion considers the change of the estimated parameter and orders between two consecutive iterations. If the change of the parameters is smaller than ε p = 10 −6 , the instrumental variable method terminates, and the next step of the order identification is performed. If the change of the order is smaller than ε o = 10 −6 , the combined parameter and order identification process terminates.
To distinguish the identified parameters from the reference parameters, these are denoted asp. The notationα is used to mark the identified order. The identification results are presented in Figure 7. The reference parameters, which were determined in the section before, are illustrated as dotted blue lines, while the identified parameters and orders are given as a solid red line. Even though only noisy observations of the input and output signal are available, all parameters and the order converge against the reference values. Nonetheless, the exact values are not reached within the considered time. This is in line with the discussion about the influence of the noise on the nonlinear equation for the order identification in Section 3.2. The relative error after termination, i.e., the relative error between the data of the eleventh iteration and reference values, for the parameters are e r,1 =b 1 (α diff,11 ) − b 1 b 1 × 100 % = 0.6 %,

Conclusions
In this paper, a combined parameter and order identification method for fractional systems under noisy observations of the output signal was proposed.
The order identification was performed by an optimization-based approach that relied on the modulating function method. The optimization-based approach ensured the order identification to be robust against measurement noise. The parameter identification was based on the instrumental variable approach from [19]. In contrast to the methods from the literature, our method did not require the system to be at rest at the beginning of the identification process. Furthermore, we provided practical implementable expressions of the necessary equations and derivatives, which were derived in analytical form.
We applied the combined identification method to a fractional model of a lithium-ion battery. The experimental results revealed that the proposed method yielded the correct parameters in the presence of noisy observations of the output and input signal. These identification results demonstrated that the method was practically applicable.