#### Abstract

The distributed order concept, which is a parallel connection of fractional order integrals and derivatives taken to the infinitesimal limit in delta order, has been the main focus in many engineering areas recently. On the other hand, there are few numerical methods available for analyzing distributed order systems, particularly under stochastic forcing. This paper proposes a novel numerical scheme for analyzing the behavior of a distributed order linear single input single output control system under random forcing. The method is based on the operational matrix technique to handle stochastic distributed order systems. The existing Monte Carlo, polynomial chaos, and frequency methods were first adapted to the stochastic distributed order system for comparison. Numerical examples were used to illustrate the accuracy and computational efficiency of the proposed method for the analysis of stochastic distributed order systems. The stability of the systems under stochastic perturbations can also be inferred easily from the moment of random output obtained using the proposed method. Based on the hybrid spectral framework, the optimal design was elaborated on by minimizing the suitably defined constrained-optimization problem.

#### 1. Introduction

Fractional/distributed order calculus is applied widely across a range of disciplines, such as physics, biology, chemistry, finance, physiology, and control engineering [1–6]. The memory property of fractional order calculus provides a novel tool to model real-world plants better than integer order ones such as diffusion plants [5]. Fractional calculus has been used for modeling of turbulence in [2]. In [3], the concept of fractional calculus is used for interpreting the underlying mechanism of dielectric relaxation. A method for design fractional order controllers for deterministic systems is proposed in [6].

The distributed order (DO) equation, which is a generalized concept fractional order, was first proposed by Caputo in 1969 [7] and solved by him in 1995 [8]. The general solution of linear DO was then discussed systematically [9]. Later, the DO concept was used to examine the diffusion equation [10], the rheological properties of composite materials, and other real complex physical phenomena [11–14]. Several different methods for the time domain analysis of DO systems have been reported [15–18]. On the other hand, a numerical method for the analysis of a DO operator is still immature and requires further development. In particular, there are few methods to analyze DO systems under the excitation of random processes. This motivated the theme of this study: the development of a computational scheme for the analysis basic of a DO system with stochastic settings. The operational matrix (OP) has attracted considerable attention for the analysis of a range of dynamic systems [19–21]. The main characteristic of this technique is that different analysis problems can be reduced to a system of algebraic equations using different types of orthogonal functions, which greatly simplifies the problem [19]. On the other hand, to the best of the author’s knowledge, there are no reports on the analysis of stochastic DO systems using an OP. Many natural systems often suffer from stochastic noise that causes fluctuations in their behavior, making them deviate from deterministic models. Therefore, it is important to examine the statistical characteristics of states (mean, variance) for those stochastic systems. This problem is often called statistical analysis (or uncertainty quantification) of a system [22–24]. This paper proposes a numerical scheme based on the OP technique for the statistical analysis of DO systems.

The Monte Carlo (MC) method is commonly used to simulate a stochastic model [25, 26]. The method relies on the sampling of independent realizations of random inputs according to their prescribed probability distribution. The data is fixed for each realization and the problem becomes deterministic. Solving the multiple deterministic realizations builds an ensemble of solutions, that is, the realization of random solutions, from which statistical information can be extracted, for example, the mean and variance. Nevertheless, this method typically reveals slow convergence and has a large computational demand. For example, the mean values typically converge as , where is the number of samples.

Generalized polynomial chaos (gPC) [27–32] represents a more recent tool for quantifying the uncertainty within system models. The approach involves expressing stochastic quantities as the orthogonal polynomials of random input parameters. This method is actually a spectral representation in random space and converges rapidly when the expanded function depends smoothly on random parameters. On the other hand, the stochastic inputs of many systems involve random processes parameterized by truncated Karhunen-Loeve (KL) expansions, and the dimensionality of the KL expansions depends on the correlation lengths of these processes. For input processes with low correlation lengths, the number of dimensions required for an accurate representation can be extremely large.

The OP method [29], where a system is described by a stochastic operator (operational matrix), is an alternative approach for the simulation of stochastic integer order systems. This method involves the inverse of the stochastic operators as Neumann series and is most effective for systems with inputs with low correlation lengths. On the other hand, it is restricted to small random parametric uncertainty.

In a recent study [33], the authors introduced a hybrid spectral method, which combines the advantages of both the OP and polynomial chaos (PC), to simulate single input single output (SISO) stochastic fractional order systems. In the present study, the method reported in [33] was extended to the statistical analysis of DO systems affected by stochastic fluctuations. Here, the stochastic operator was approximated using PC instead of a Neumann series. This method provides the algebraic relationships between the first- and second-order stochastic moments of the input and output of a system, hence bypassing the KL expansions that can require large dimensions for accurate results. In contrast to the traditional OP method, the proposed method is not limited by the magnitude of the uncertainty.

Section 2 briefly introduces a DO system and the OP technique for uncertainty quantification in this system, leading to computation of the moments of random matrices. Section 3 summarizes the process of calculating the moments of the random matrices using a stochastic collocation. Section 4 defines the suitable performance objectives coupled with the spectral method for the design of a stochastic linear DO system. Section 5 provides examples to demonstrate the use of the proposed method. The results of the proposed deterministic system with a DO were compared with those of other existing numerical and analytical methods. To assess a stochastic DO system, the MC, gPC, and frequency methods were first adopted to the stochastic DO system for comparison because the analytical results were unavailable. The results from the proposed method were then compared with the numerical results from the MC, gPC, and frequency methods.

#### 2. Preliminary of Fractional and Distributed Order System

In this section, we give some necessary definitions and preliminaries of the fractional calculus theory which will be used in this paper.

##### 2.1. Governing Equation for System Dynamics with Fractional Order Dynamics

Fractional calculus considers the generalization of the integration and differentiation operator to a noninteger order [34, 35]:where is the order of the operator.

Among many formulations of the generalized derivative, the Riemann-Liouville (RL) definition is used most often:where denotes the gamma function and is an integer satisfying .

The RL fractional integral of a function is defined as follows:Another popular definition of a fractional order derivative is the Caputo () definition [36],The Laplace transform for a fractional order derivative under zero initial conditions can be defined as .

Note that, under a zero initial condition, the two Riemann-Liouville and Caputo definitions are equivalent.

Therefore, a fractional order single input single output (SISO) system can be described by a fractional order differential equation as or by a transfer function:where and are the arbitrary real positive numbers and and are the input and output of the system, respectively.

##### 2.2. Distributed Order Systems

The DO differential operation is defined as follows [17]:where denotes the distribution function of order .

Therefore, the general form of the DO differential equation isFor time domain analysis of the DO system, the integral in (7) is discretized using the quadrature formula as follows [16, 17]:where , are the node and weight from the quadrature formula, respectively. In other words, the DO equation is approximated as a multiterm fractional order equation and can be rearranged as (5).

##### 2.3. Operational Matrices of Block Pulse Function for the Analysis of Distributed Order Systems

Block pulse functions (BPFs) are a complete set of orthogonal functions that are defined over the time interval, ,where is the number of block pulse functions.

Therefore, any function that can be absolutely integrated on the time interval can be expanded to a series from the block pulse basis as follows:where constitutes the block pulse basis. From here, the subscript of is dropped out for the convenience of notation.

The expansion coefficients (or spectral characteristics) can be evaluated as follows:Furthermore, any function that is absolutely integrable on the time interval can be expanded as follows:with expansion coefficients (or spectral characteristics) ofEquation (3) can be expressed in terms of the OP [19],where the generalized OP integration of the block pulse function, , isThe elements of the generalized OP integration can be given byThe generalized OP of a derivative of order iswhere is the identity matrix.

The generalized OP of derivative can be used to approximate (2) as follows:Therefore, using the OP, the discretization of DO can be expressed asThe DO system in (7) can be rewritten in terms of the OP, , as follows:The input and output are related by the following equation:

##### 2.4. Stochastic Analysis of Distributed Order Systems

Consider the system described by (7), which has the spectral characteristics of input and output linked by (21). Assume that the system is excited by random forcing with a given mean and covariance function as follows:where denotes the expectation operator and the spectral characteristics of the mean and covariance function of the input are calculated in (11) and (13).

Using the OP, the spectral characteristics of the mean and covariance of the output are given by the following [33] (the details are available in Appendices):Therefore,The random parameters , result in the random OP in (23) and (24), and its (OP ) moment can be estimated using a stochastic collocation method, which is described in the next section.

When the parameters, , , are deterministic, (23) and (24) become

*Remarks*. The relationship in (25) is invariant with respect to the orthogonal polynomial used to construct the OP of the fractional order integral and derivative. The relationships in (24) and (25) are only available for a linear system.

#### 3. Stochastic Collocation for the Operational Matrix

A stochastic collocation method, which is described briefly below, is based on the gPC and can easily estimate the means and variances of complex dynamics. Therefore, it has been used to estimate the moment of the random matrix in (24).(i)Assume that a random OP has the form, , where is a vector of independent random parameters with the probability density functions . Vector has the joint probability density function, , with the support, .(ii)Choose a suitable quadrature set for each random parameter according to the probability density so that a one-dimensional integration can be approximated as accurately as possible by , where is the th node and is the corresponding weight.(iii)Construct a multidimensional cubature set by tensorization of the one-dimensional quadrature set over all the combined multi-index . Because manipulation of the multi-index is cumbersome in practice, a single index is preferable for manipulating these equations. The multi-index is often replaced by a graded lexicographic order index, [27]. Because the weighting functions of the cubature are the same as the probability density functions, the moment of the random matrix can be approximated byThe Matlab suite, OPQ, can be used to obtain the one-dimensional quadrature sets and their corresponding orthogonal polynomials (polynomial chaos) with respect to the different density weights [36].

The algorithm of the proposed method for the analysis of a stochastic system can be summarized as follows:(a)Calculate the coefficients , of the expansions of the mean and covariance of the input as shown in (11) and (13).(b)Rewrite the DO differential equation in (7) in terms of OP as (20).(c)The coefficients of expansions of the mean and covariance function of the output are obtained from (23). In (23), the moments of several random matrices need to be calculated. The moment of a random matrix is calculated by the stochastic collocation method as (26).(d)Finally, the mean and covariance of the output are obtained as (24).For a clearer understanding of the algorithm, a similar algorithm is depicted graphically in [33] for the analysis of stochastic linear fractional order systems.

#### 4. Optimal D Controller Design

Assume that the system is described by (7), where coefficients , are independent random variables with given distributions. The set point input is a random process with a given mean and covariance function as follows:The system is in the closed loop configuration, as shown in Figure 1, with a fractional order controller [35]

The OP for this controller iswhere , , and are the identity matrix, the integration OP of fractional order, , and the OP of the fractional order derivative, , respectively.

Denote the OP for the system as . Using block algebra for OP operator, the OP for a closed loop system can be obtained as follows:This closed loop OP can be used to obtain the first- and second-order moment of the random output from (24).

The parameters of the controller can be obtained by optimizing the cost function defined as [37]where and are obtained from (24).

#### 5. Examples

Before going to detailed examples, let us give some information about the existing methods. Ito calculus can be used for the statistical analysis of integer order linear/nonlinear systems only with ideal white noise (noise with direct delta covariance function and infinite bandwidth). On the other hand, the PC method can be used for a system with low bandwidth noise. However, in the PC method the computational load increases significantly as the bandwidth of noise increases. The MC and Quasi-MC methods can be used for arbitrary cases (i.e., with arbitrary type of noise). However, they require a large computational effort for obtaining accurate results. To overcome these limitations in each existing method, a hybrid spectral method [33] was proposed for the statistical analysis of fractional order linear SISO systems with arbitrary type of random input. In this paper, the methodology in [33] was extended to the DO case. Several different case studies were considered to show the efficiency of the proposed method handling different kind of random inputs in a unified frame work: band-limited white noise (noise with low bandwidth), ideal white noise, and fractional Brownian motion with Hurst parameter . It should be noted that when , fractional Brownian motion becomes Brownian motion, whose derivative is ideal white noise.

##### 5.1. Examples 1(a) and 1(b): Band-Limited White Noise Input

Because an integer order system can be considered as a special case of DO systems, this example considers a simple linear integer order system, from [38] with a band-limited white noise as the input.

Let the input have a zero mean and covariance function of , where the sinc function is defined as The power spectral density of the input iswhere is known as bandwidth of the noise. As approaches to infinity, the process will become ideal white noise.

Therefore, the power spectral density function of stochastic output is given byThis is a linear system; the means of the input and output are zero.

Using the frequency method, the exact steady state variance of output can be expressed as [38]The OP for this linear system iswhere is the identity matrix; is the OP of integration (of order one), which was calculated using (15) and (16).

This system lacks random parameters. The covariance function of the system output is approximated by (25). From (25), it can be seen that the mean of the output by the proposed method is zero.

Two numerical cases are considered.

(a) Consider ; . Figure 2 shows the output variance obtained by the proposed method for . For comparison, the results by the gPC, MC, and frequency methods are also given. Note that the frequency method in (34) can provide the exact (analytical) steady state variance. The random process input, , was parameterized using a noncanonical decomposition [33, 39] (the details are available in Appendices). The results from the proposed method were quite satisfactory. Table 1 lists the simulation parameters and computational times required for each method.

(b) Consider ; . Figure 3 presents the variance obtained by the proposed method for . If the number of cubature nodes is kept as in case (a), the gPC method cannot obtain an accurate result in the steady state. The result indicates that in the gPC method the number of cubature nodes needs to increase with increasing bandwidth of the noise, and the computational load increases for obtaining the same accurate result accordingly.

Table 1 presents the computational times and simulation parameters for all methods. From this table and Figures 2 and 3, the proposed method provides better performance in terms of accuracy and computational load.

*Remark*. The gPC approaches can be divided into two subcategories: intrusive Galerkin [27, 31] approaches and nonintrusive projection approaches [27–30, 32]. The advantage of nonintrusive approach is ease of implementation. For this reason, nonintrusive (collocation) methods have become very popular. The intrusive Galerkin method offers the most accurate solutions involving the least number of equations in multidimensional random spaces, but it is more cumbersome to implement. Thus, in this paper, the nonintrusive method is referred to as the gPC method.

##### 5.2. Examples 2(a), 2(b), 2(c), and 2(d): Ideal White Noise

*(a) Example 2(a): Double Delta Function Distributed Order System*. This example considers the statistical analysis of a special case of DO integrator taken from the literature [40] as follows:where and is the Dirac delta function. Therefore, (36) is actually a double fractional integrator,The case where the input is an ideal white noise with a zero mean and covariance function was considered:The exact variance of the output is given by [40]where is the Mittag-Leffler function, which can be calculated using the Matlab mlf.m function [41].

The OP for system (37) is given bywhere is the OP of derivative of order .

This system does not have random parameters. Therefore, the covariance function of system output can be approximated by (25). The regularization technique [33] is used to approximate the Dirac delta covariance function. Figure 4 presents the variance obtained using the proposed method for , , and . The relative error of the proposed method with respect to (39) is also shown in Figure 4. The result by the proposed method is quite satisfactory. Figure 5 compares the absolute error for the output variance by the proposed and MC methods (with respect to the exact variance in (39)). The simulation times are listed in Table 1 for both methods. For MC simulations, the Matlab code, fode_sol.m, from [42] was used. Again, Table 1 and Figure 5 show that the proposed method has better accuracy with less computational burden than the MC method.

*Remarks*. The fact that the gPC method becomes computationally intractable for ideal white noise input makes the proposed method more attractive.

*(b) Example 2(b): Uniform Distributed Order Integrator*. This example considers the statistical analysis of a DO integrator taken from [40]:Again, this study considered the case where the input is an ideal white noise with a zero mean and covariance function as shown in (38).

The variance of the output is given by the following [40]:The OP for system (41) is given bywhere is the OP of derivative of order and are the nodes and weights from the Legendre quadrature.

Figure 6 shows the variance of the output by the proposed method. The absolute error with respect to the exact variance (42) is also shown.

*(c) Example 2(c): ** Controller Design for Uniform Distributed Order Integrator with Stochastic Input*. From the examples above, it can be seen that the proposed method for predicting the mean and variance of the system output provides better accuracy and lower computational load than the other methods such as the MC and gPC. Therefore, it is more suitable for direct optimal design by minimization of the objective function in (30).

Consider a controller as a closed loop configuration with the uniform DO integrator above. The set point is a random process with a unit mean and covariance function . This input can be viewed as a combination of the deterministic set point and zero mean measurement noise [37].

The control objective is to track the deterministic unit step input. This can be achieved by minimizing objective function defined in Section 4. The search space for the optimal parameters of the controller was limited to ; ; ; ; for simplicity, as in other studies on the probabilistic approach [29, 37]. The resulting controller can be expressed as .

Figure 7 shows mean and variance of system output with this controller. Figure 8 shows 500 possible responses of the uncertain system with the proposed controller. From the finiteness of the output variance, the stability of system can be determined.

*(d) Example 2(d): Improved Mean Tracking Control with Iterative Learning Control*. Since the proposed method allows lower computational time for prediction of the mean of system output under random forcing, it can be used with iterative learning control in which input sequence is refined from one trial to next trial [43].

Consider a problem where the mean of closed loop system in Example 2(c) needs to track desired mean . The following iterative learning control scheme can be used for refining the mean of set point input:where and is the closed loop OP.

Figure 9 shows the simulation results by the proposed method. It can be seen from the figure that the iterative learning algorithm in (44) improves the tracking error in the mean as increases. The MC simulations with 512 sample responses are shown in Figure 10. It can be seen that the designed control algorithm can track the desired mean despite the random forcing.

*Remarks*. Note that the low computational cost of the proposed method enables using the iterative learning control algorithm.

##### 5.3. Example 3: Linear Integer Order with Fractional Brownian Motion Input

For each , there is a real-value Gaussian process such that and the covariance function of , . This process is called standard fractional Brownian motion with the Hurst parameter . If , then the corresponding standard fractional Brownian motion is the well-known standard Brownian motion.

Recently, there has been growing interest in linear systems with fractional Brownian motion input [44–46]. On the other hand, in contrast to standard Brownian motion, where the moment of output can be obtained by Ito calculus, there are very few methods available for obtaining the moment of the stochastic output of a general linear system with a fractional Brownian motion input.

Consider a random process that satisfies the following:where is a fractional Brownian motion with the Hurst parameter . The mean of is . The variance of satisfies a differential equation [44]:where is given byFigure 11 shows the evolution of variance of for obtained from (46) and (47).

Equation (45) can be rewritten in terms of OP as follows:where is the OP of the derivative of order one. Therefore, one can easily obtain the covariance of the random process, , by utilizing the OP for this system, , and covariance function of fractional Brownian motion and (25). Figure 11 shows the variance of obtained using the proposed method. Finally, the variance of by a MC estimation is also given in the same figure.

##### 5.4. Example 4: Linear Distributed Order with Stochastic Parametric and Additive Uncertainties

*(a) Example 4(a)*. This case considers a DO system with both random parameters and random forcing. The system can be expressed aswhere and are uniform random variables in . The system is in a closed loop configuration with a fractional PI controller, from [46]. Note that the controller, , was designed for a nominal system, that is, for . The input is a band-limited white noise process with a unit mean and covariance function of the following: , where .

Figure 12 compares the mean and variance of the output calculated by the proposed method using (24) with the results from the gPC and MC methods. In this study, the gPC and MC methods were first applied to the DO systems under stochastic forcing for comparison. In the MC and gPC methods, the DO term was discretized first and the routine fode_sol.m [42] was then used to integrate the multiterm fractional order version of the DO system. Again, the random process input, , was parameterized using a noncanonical decomposition. Table 1 lists the simulation parameters and computational times needed for each method. Figure 12 and Table 1 show that the proposed method can provide similar accuracy with much less computational effort than the other methods. The advantage of the proposed method lies in its use of operational matrices: the mean and covariance of the output can be obtained directly from those of the input without parameterization of the input.

*(b) Example 4(b)*. The proposed method was applied to the design of fractional order PID controller for this system. The covariance of the output can be obtained directly from those of the input without parameterization of the input. The control objective is to track the deterministic unit step input. The search space for the parameters of the controller was limited to ; ; ; ; . The result is as follows: .

Figure 13 shows the mean and variance of the system output by the proposed controller and the controller, from [46]. Figure 14 shows a bounded region for 1000 possible responses of the uncertain system with the proposed and initial controllers. The proposed controller outperformed the initial controller.

#### 6. Conclusions

A hybrid spectral method was proposed to analyze DO systems in a stochastic setting with arbitrary random forcing and parametric uncertainties. To analyze the system with stochastic parameter perturbation, the stochastic collocation was used to estimate the random operator. This combines the advantages of both the OP technique and PC method. The use of operational matrices explicitly provides the relationship between the first- and second-order moment for the input and output of a system, bypassing parameterization of the random input when predicting the statistical characteristics and reducing the dimensions of the random space. This can also effectively handle a system with a low correlation length input (i.e., ideal white noise) by regularization. The numerical examples show that the proposed method provides superior accuracy and computational efficiency for analyzing stochastic DO systems over other existing methods, such as the gPC, MC, and frequency methods: the frequency method can give the only result at the steady state; the accuracy and efficiency of the gPC method are degraded for a wideband process. Although the MC method is straightforward, its accuracy and computational burden are problematic. On the other hand, the explicit relationship in (24) is only available for a linear system; the applicability of the proposed method is restricted to linear systems.

#### Appendices

#### A. Derivation of (24)

Consider a system with its input and output linked bywhere is the OP of the system. The input and parameters , are random.

Therefore, the mean of the input and the output in (A.1) was calculated aswhere denotes the expectation operator; ; .

The statistical independence of and leads toTherefore, the spectral characteristics (or expansion coefficients) of the mathematical expectations of input and output are related byIntroducing the system’s signal in the spectral form leads to an equation defining the correlation function of the output to be written as follows:Therefore, (A.5) becomeswhere is the square matrix of expansion coefficients (spectral characteristics) of the input’s correlation function, which is given byThe covariance function of the system’s input is defined as Expanding (A.8) in terms of the orthogonal functions gives the following:The spectral characteristics of the input signal’s moments are given bySubstituting (A.10) into (A.6) givesThe covariance function of the system’s output is then given byor in spectral form,Equation (A.13) gives the relationship between the spectral characteristics of the moments of the system’s output and input. Equations (A.4) and (A.13) are then combined to form (24) in the paper

#### B. Noncanonical Decomposition of Stationary Random Processes [39]

Consider a stationary random process with mean , covariance , and variance . This process can be represented as withwhere and are independent, is Gaussian, and is a random variable with a probability density function (pdf) given in (B.1).

*Proof. * has a mean ofand a covariance function ofwhere is the central component of the random process . In (B.3), the properties of and and the independence of are used to simplify the equation.

The covariance function also can be calculated as the inverse Fourier transform of the power spectral densityComparing (B.3) and (B.4) gives the pdf of in (B.1). Because , is a proper pdf.(A)A first-order Markov process with a mean and exponential covariance can be parameterized as , where is Gaussian, as in (B.1), and , (B)Band-limited white noise with a mean and covariance can be parameterized as , where is Gaussian, as in (B.1), and , .

#### Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.

#### Acknowledgments

This research was supported by Basic Science Research Program through the National Research Foundation of Korea (NRF) funded by the Ministry of Education (2015R1D1A3A01015621). This work was also supported by Priority Research Centers Program through the National Research Foundation of Korea (NRF) funded by the Ministry of Education (2014R1A6A1031189).