# Optimal Design of Stochastic Distributed Order Linear SISO Systems Using Hybrid Spectral Method

Authors: Pham Luu Trung Duong; Ezra Kwok; Moonyong Lee
Journal: Mathematical Problems in Engineering
 (2015)

---

## 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.


---

## Body

## 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 P
I
λ
D
μ 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 1
/
M, where M 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.Section2 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]:(1)
D
0
α
=
d
α
d
t
α
α
>
0
1
α
=
0
∫
0
t
d
τ
-
α
α
<
0
,where α
∈
R is the order of the operator.Among many formulations of the generalized derivative, the Riemann-Liouville (RL) definition is used most often:(2)
D
RL
0
α
f
t
=
1
Γ
m
-
α
d
d
t
m
∫
0
t
f
τ
t
-
τ
1
-
m
-
α
d
τ
,where Γ
(
x
) denotes the gamma function and m is an integer satisfying m
-
1
<
α
<
m.The RL fractional integral of a functionf
(
t
) is defined as follows:(3)
I
RL
0
α
f
t
=
1
Γ
α
∫
0
t
f
τ
t
-
τ
1
-
α
d
τ
.Another popular definition of a fractional order derivative is the Caputo (C) definition [36],(4)
D
C
t
α
=
1
Γ
m
-
α
∫
a
t
t
-
τ
m
-
α
-
1
f
m
τ
d
τ
.The Laplace transform for a fractional order derivative under zero initial conditions can be defined as L
{
D
0
α
f
(
t
)
}
=
s
α
F
(
s
).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 asa
0
D
0
α
0
y
t
+
a
1
D
0
α
1
y
t
+
⋯
+
a
l
D
0
α
l
y
t
=
b
0
D
0
β
0
u
t
+
⋯
+
b
m
D
0
β
m
u
(
t
) or by a transfer function:(5)
G
s
=
Y
s
U
s
=
b
m
s
β
m
+
⋯
+
b
0
s
β
0
a
l
s
α
l
+
⋯
+
a
0
s
α
0
,where α
i and β
i are the arbitrary real positive numbers and u
(
t
) and y
(
t
) are the input and output of the system, respectively.


### 2.2. Distributed Order Systems

The DO differential operation is defined as follows [17]:(6)
D
t
ρ
α
f
t
=
∫
γ
1
γ
2
ρ
α
D
t
α
f
t
d
α
,where ρ
(
α
) denotes the distribution function of order α.Therefore, the general form of the DO differential equation is(7)
∑
i
=
1
n
a
i
D
t
ρ
i
α
y
t
=
∑
j
=
1
m
b
j
D
t
ρ
j
α
u
t
.For time domain analysis of the DO system, the integral in (7) is discretized using the quadrature formula as follows [16, 17]:(8)
∫
γ
1
γ
2
ρ
α
D
t
α
f
t
d
α
≈
∑
l
=
1
Q
ρ
α
l
D
t
α
l
f
t
υ
l
,where α
l, υ
l 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,[
0
,
τ
],(9)
ψ
i
=
1
i
-
1
N
τ
≤
t
≤
i
N
τ
0
e
l
s
e
w
h
e
r
e
,where N is the number of block pulse functions.Therefore, any function that can be absolutely integrated on the time interval[
0
,
τ
] can be expanded to a series from the block pulse basis as follows:(10)
f
t
=
ψ
N
T
t
C
f
=
∑
i
=
1
N
c
f
i
ψ
i
t
,where ψ
N
T
(
t
)
=
[
ψ
1
(
t
)
,
…
,
ψ
N
(
t
)
] constitutes the block pulse basis. From here, the subscript N of ψ
N
T
(
t
) is dropped out for the convenience of notation.The expansion coefficients (or spectral characteristics) can be evaluated as follows:(11)
c
f
i
=
N
τ
∫
i
-
1
/
N
τ
i
/
N
τ
f
t
ψ
i
t
d
t
.Furthermore, any function g
(
t
1
,
t
2
) that is absolutely integrable on the time interval [
0
,
τ
]
×
[
0
,
τ
] can be expanded as follows:(12)
g
t
1
,
t
2
=
∑
i
=
1
N
∑
j
=
1
N
c
i
j
ψ
i
t
1
ψ
j
t
2
=
ψ
T
t
1
C
g
ψ
t
2
,with expansion coefficients (or spectral characteristics) of(13)
c
i
j
=
N
τ
2
∫
i
-
1
/
N
τ
i
/
N
τ
∫
i
-
1
/
N
τ
i
/
N
τ
g
t
1
,
t
2
ψ
i
t
1
ψ
j
t
2
d
t
1
d
t
2
.Equation (3) can be expressed in terms of the OP [19],(14)
I
0
α
f
t
=
ψ
t
T
A
α
C
f
,where the generalized OP integration of the block pulse function, A
α, is(15)
A
α
=
P
α
T
=
τ
N
α
1
Γ
α
+
2
f
1
f
2
f
3
⋯
f
N
0
f
1
f
2
⋯
f
N
-
1
⋮
⋱
⋱
⋱
⋮
0
⋯
⋯
⋯
f
1
T
.The elements of the generalized OP integration can be given by(16)
f
1
=
1
;
f
p
=
p
α
+
1
-
2
p
-
1
α
+
1
+
p
-
2
α
+
1
for
p
=
2,3
,
…
,
N
.The generalized OP of a derivative of order α is(17)
B
α
A
α
=
I
,where I is the identity matrix.The generalized OP of derivative can be used to approximate (2) as follows:(18)
D
0
α
f
t
=
ψ
t
T
B
α
C
f
.Therefore, using the OP, the discretization of DO can be expressed as(19)
D
t
ρ
α
f
t
=
∫
γ
1
γ
2
ρ
α
D
t
α
f
t
d
α
≈
∑
l
=
1
Q
ρ
α
l
D
t
α
l
f
t
υ
l
=
∑
l
=
1
Q
υ
l
ρ
α
l
ψ
t
T
B
α
l
C
f
.The DO system in (7) can be rewritten in terms of the OP, A
G, as follows:(20)
A
G
=
∑
i
=
1
n
a
i
∑
l
=
1
Q
υ
l
ρ
i
α
l
B
α
l
-
1
∑
j
=
1
m
b
j
∑
l
=
1
Q
υ
l
ρ
j
α
l
B
α
l
.The input and output are related by the following equation:(21)
C
Y
=
A
G
C
U
;
Y
t
=
C
Y
T
ψ
t
;
U
t
=
C
U
T
ψ
t
.


### 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:(22)
M
U
t
=
E
U
t
=
C
m
U
T
ψ
t
,
κ
U
U
=
E
U
t
1
-
M
U
t
1
U
t
2
-
M
U
t
2
=
∑
i
=
1
N
∑
j
=
1
N
ψ
i
t
1
ψ
j
t
2
c
i
j
=
ψ
t
1
T
C
K
U
U
ψ
t
2
,where E 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):(23)
C
m
Y
=
E
A
G
C
m
U
,
C
κ
Y
Y
=
E
A
G
C
κ
U
U
+
C
m
U
C
m
U
T
A
G
T
-
C
m
Y
C
m
Y
T
.Therefore,(24)
m
Y
t
=
ψ
t
1
T
C
m
Y
=
ψ
t
1
T
E
A
G
C
m
U
,
κ
Y
Y
t
1
,
t
2
=
ψ
t
1
T
C
κ
Y
Y
ψ
t
2
=
ψ
t
1
T
E
A
G
C
κ
U
U
+
C
m
U
C
m
U
T
A
G
T
ψ
t
2
-
ψ
t
1
T
C
m
Y
C
m
Y
T
ψ
t
2
.The random parameters a
i, b
j result in the random OP A
G in (23) and (24), and its (OP A
G) moment can be estimated using a stochastic collocation method, which is described in the next section.When the parameters,a
i, b
j, are deterministic, (23) and (24) become(25)
m
Y
t
=
ψ
t
1
T
A
G
C
m
U
,
κ
Y
Y
t
1
,
t
2
=
ψ
t
1
T
A
G
C
κ
U
U
+
C
m
U
C
m
U
T
A
G
T
ψ
t
2
-
ψ
t
1
T
C
m
Y
C
m
Y
T
ψ
t
2
.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.


## 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]:(1)
D
0
α
=
d
α
d
t
α
α
>
0
1
α
=
0
∫
0
t
d
τ
-
α
α
<
0
,where α
∈
R is the order of the operator.Among many formulations of the generalized derivative, the Riemann-Liouville (RL) definition is used most often:(2)
D
RL
0
α
f
t
=
1
Γ
m
-
α
d
d
t
m
∫
0
t
f
τ
t
-
τ
1
-
m
-
α
d
τ
,where Γ
(
x
) denotes the gamma function and m is an integer satisfying m
-
1
<
α
<
m.The RL fractional integral of a functionf
(
t
) is defined as follows:(3)
I
RL
0
α
f
t
=
1
Γ
α
∫
0
t
f
τ
t
-
τ
1
-
α
d
τ
.Another popular definition of a fractional order derivative is the Caputo (C) definition [36],(4)
D
C
t
α
=
1
Γ
m
-
α
∫
a
t
t
-
τ
m
-
α
-
1
f
m
τ
d
τ
.The Laplace transform for a fractional order derivative under zero initial conditions can be defined as L
{
D
0
α
f
(
t
)
}
=
s
α
F
(
s
).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 asa
0
D
0
α
0
y
t
+
a
1
D
0
α
1
y
t
+
⋯
+
a
l
D
0
α
l
y
t
=
b
0
D
0
β
0
u
t
+
⋯
+
b
m
D
0
β
m
u
(
t
) or by a transfer function:(5)
G
s
=
Y
s
U
s
=
b
m
s
β
m
+
⋯
+
b
0
s
β
0
a
l
s
α
l
+
⋯
+
a
0
s
α
0
,where α
i and β
i are the arbitrary real positive numbers and u
(
t
) and y
(
t
) are the input and output of the system, respectively.


## 2.2. Distributed Order Systems

The DO differential operation is defined as follows [17]:(6)
D
t
ρ
α
f
t
=
∫
γ
1
γ
2
ρ
α
D
t
α
f
t
d
α
,where ρ
(
α
) denotes the distribution function of order α.Therefore, the general form of the DO differential equation is(7)
∑
i
=
1
n
a
i
D
t
ρ
i
α
y
t
=
∑
j
=
1
m
b
j
D
t
ρ
j
α
u
t
.For time domain analysis of the DO system, the integral in (7) is discretized using the quadrature formula as follows [16, 17]:(8)
∫
γ
1
γ
2
ρ
α
D
t
α
f
t
d
α
≈
∑
l
=
1
Q
ρ
α
l
D
t
α
l
f
t
υ
l
,where α
l, υ
l 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,[
0
,
τ
],(9)
ψ
i
=
1
i
-
1
N
τ
≤
t
≤
i
N
τ
0
e
l
s
e
w
h
e
r
e
,where N is the number of block pulse functions.Therefore, any function that can be absolutely integrated on the time interval[
0
,
τ
] can be expanded to a series from the block pulse basis as follows:(10)
f
t
=
ψ
N
T
t
C
f
=
∑
i
=
1
N
c
f
i
ψ
i
t
,where ψ
N
T
(
t
)
=
[
ψ
1
(
t
)
,
…
,
ψ
N
(
t
)
] constitutes the block pulse basis. From here, the subscript N of ψ
N
T
(
t
) is dropped out for the convenience of notation.The expansion coefficients (or spectral characteristics) can be evaluated as follows:(11)
c
f
i
=
N
τ
∫
i
-
1
/
N
τ
i
/
N
τ
f
t
ψ
i
t
d
t
.Furthermore, any function g
(
t
1
,
t
2
) that is absolutely integrable on the time interval [
0
,
τ
]
×
[
0
,
τ
] can be expanded as follows:(12)
g
t
1
,
t
2
=
∑
i
=
1
N
∑
j
=
1
N
c
i
j
ψ
i
t
1
ψ
j
t
2
=
ψ
T
t
1
C
g
ψ
t
2
,with expansion coefficients (or spectral characteristics) of(13)
c
i
j
=
N
τ
2
∫
i
-
1
/
N
τ
i
/
N
τ
∫
i
-
1
/
N
τ
i
/
N
τ
g
t
1
,
t
2
ψ
i
t
1
ψ
j
t
2
d
t
1
d
t
2
.Equation (3) can be expressed in terms of the OP [19],(14)
I
0
α
f
t
=
ψ
t
T
A
α
C
f
,where the generalized OP integration of the block pulse function, A
α, is(15)
A
α
=
P
α
T
=
τ
N
α
1
Γ
α
+
2
f
1
f
2
f
3
⋯
f
N
0
f
1
f
2
⋯
f
N
-
1
⋮
⋱
⋱
⋱
⋮
0
⋯
⋯
⋯
f
1
T
.The elements of the generalized OP integration can be given by(16)
f
1
=
1
;
f
p
=
p
α
+
1
-
2
p
-
1
α
+
1
+
p
-
2
α
+
1
for
p
=
2,3
,
…
,
N
.The generalized OP of a derivative of order α is(17)
B
α
A
α
=
I
,where I is the identity matrix.The generalized OP of derivative can be used to approximate (2) as follows:(18)
D
0
α
f
t
=
ψ
t
T
B
α
C
f
.Therefore, using the OP, the discretization of DO can be expressed as(19)
D
t
ρ
α
f
t
=
∫
γ
1
γ
2
ρ
α
D
t
α
f
t
d
α
≈
∑
l
=
1
Q
ρ
α
l
D
t
α
l
f
t
υ
l
=
∑
l
=
1
Q
υ
l
ρ
α
l
ψ
t
T
B
α
l
C
f
.The DO system in (7) can be rewritten in terms of the OP, A
G, as follows:(20)
A
G
=
∑
i
=
1
n
a
i
∑
l
=
1
Q
υ
l
ρ
i
α
l
B
α
l
-
1
∑
j
=
1
m
b
j
∑
l
=
1
Q
υ
l
ρ
j
α
l
B
α
l
.The input and output are related by the following equation:(21)
C
Y
=
A
G
C
U
;
Y
t
=
C
Y
T
ψ
t
;
U
t
=
C
U
T
ψ
t
.


## 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:(22)
M
U
t
=
E
U
t
=
C
m
U
T
ψ
t
,
κ
U
U
=
E
U
t
1
-
M
U
t
1
U
t
2
-
M
U
t
2
=
∑
i
=
1
N
∑
j
=
1
N
ψ
i
t
1
ψ
j
t
2
c
i
j
=
ψ
t
1
T
C
K
U
U
ψ
t
2
,where E 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):(23)
C
m
Y
=
E
A
G
C
m
U
,
C
κ
Y
Y
=
E
A
G
C
κ
U
U
+
C
m
U
C
m
U
T
A
G
T
-
C
m
Y
C
m
Y
T
.Therefore,(24)
m
Y
t
=
ψ
t
1
T
C
m
Y
=
ψ
t
1
T
E
A
G
C
m
U
,
κ
Y
Y
t
1
,
t
2
=
ψ
t
1
T
C
κ
Y
Y
ψ
t
2
=
ψ
t
1
T
E
A
G
C
κ
U
U
+
C
m
U
C
m
U
T
A
G
T
ψ
t
2
-
ψ
t
1
T
C
m
Y
C
m
Y
T
ψ
t
2
.The random parameters a
i, b
j result in the random OP A
G in (23) and (24), and its (OP A
G) moment can be estimated using a stochastic collocation method, which is described in the next section.When the parameters,a
i, b
j, are deterministic, (23) and (24) become(25)
m
Y
t
=
ψ
t
1
T
A
G
C
m
U
,
κ
Y
Y
t
1
,
t
2
=
ψ
t
1
T
A
G
C
κ
U
U
+
C
m
U
C
m
U
T
A
G
T
ψ
t
2
-
ψ
t
1
T
C
m
Y
C
m
Y
T
ψ
t
2
.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,A
=
A
(
ξ
), where ξ
=
(
ξ
1
,
ξ
2
,
…
,
ξ
n
) is a vector of independent random parameters with the probability density functions ρ
i
(
ξ
i
)
:
Γ
i
→
R
+. Vector ξ has the joint probability density function, ρ
=
∏
i
=
1
n
ρ
i, with the support, Γ
≡
∏
i
=
1
n
Γ
i
∈
R
+
n.
(ii)
Choose a suitable quadrature setξ
i
m
,
w
m
m
=
1
q
i for each random parameter according to the probability density so that a one-dimensional integration can be approximated as accurately as possible by ∫
Γ
i
A
(
ξ
i
)
ρ
i
(
ξ
i
)
d
ξ
i
=
∑
i
=
1
q
i
A
(
ξ
i
(
m
)
)
w
i
(
m
), where ξ
i
(
m
) is the mth node and w
m is the corresponding weight.
(iii)
Construct a multidimensional cubature set by tensorization of the one-dimensional quadrature set over all the combined multi-index(
j
1
,
…
,
j
n
). Because manipulation of the multi-index (
j
1
,
…
,
j
n
) 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, j [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 by(26)
E
A
=
∫
Γ
A
ξ
ρ
ξ
d
ξ
=
∑
j
=
1
Q
A
ξ
j
w
j
=
∑
j
1
=
1
q
1
⋯
∑
j
n
=
1
q
n
A
ξ
1
j
1
,
…
,
ξ
n
j
n
w
1
j
1
,
…
,
w
n
j
n
.The 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 coefficientsC
m
R, C
κ
R
Rof 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. OptimalP
I
λDμ Controller Design

Assume that the system is described by (7), where coefficients a
i, b
j are independent random variables with given distributions. The set point input is a random process with a given mean and covariance function as follows:(27)
M
R
t
=
E
R
t
=
C
m
R
T
ψ
t
,
κ
R
R
=
E
R
t
1
-
M
R
t
1
R
t
2
-
M
R
t
2
=
∑
i
=
1
N
∑
j
=
1
N
ψ
i
t
1
ψ
j
t
2
c
i
j
R
=
ψ
t
1
T
C
K
R
R
ψ
t
2
.The system is in the closed loop configuration, as shown in Figure 1, with a fractional order P
I
λ
D
μ controller [35] C
s
=
K
p
+
(
K
i
/
s
λ
)
+
K
d
s
μ
.Figure 1
Closed loop control system.The OP for this controller is(28)
A
C
=
K
p
I
+
K
i
A
λ
+
K
d
B
μ
,where I, A
λ, and B
μ 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 asA
G. Using block algebra for OP operator, the OP for a closed loop system can be obtained as follows:(29)
A
=
I
+
A
G
A
C
-
1
A
G
A
C
.This closed loop OP can be used to obtain the first- and second-order moment of the random output from (24).The parameters of theP
I
λ
D
μ controller can be obtained by optimizing the cost function defined as [37](30)
min
Kp
,
Ki
,
Kd
,
λ
,
μ
⁡
J
=
min
Kp
,
Ki
,
Kd
,
λ
,
μ
⁡
E
∫
0
T
y
d
t
-
y
t
2
d
t
=
min
Kp
,
Ki
,
Kd
,
λ
,
μ
⁡
∫
0
T
E
y
2
d
t
+
E
y
t
2
-
2
E
y
d
t
E
y
t
d
t
,where E
y
(
t
)
2 and E
y
(
t
) 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 H. It should be noted that when H
=
1
/
2, 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,G
(
s
)
=
1
/
(
1
+
T
s
) from [38] with a band-limited white noise as the input.Let the input have a zero mean and covariance function ofκ
U
U
(
t
1
,
t
2
)
=
W
B
/
π
sinc
⁡
(
(
t
1
-
t
2
)
W
B
/
π
), where the sinc function is defined as (31)
sinc
⁡
x
=
sin
⁡
π
x
π
x
elsewhere
1
for
x
=
0
.The power spectral density of the input is(32)
S
U
ω
=
1
2
π
ω
≤
W
B
0
ω
>
W
B
,where W
B is known as bandwidth of the noise. As W
B approaches to infinity, the process will become ideal white noise.Therefore, the power spectral density function of stochastic output is given by(33)
S
Y
ω
=
1
2
π
1
1
+
T
ω
2
ω
≤
W
B
0
ω
>
W
B
.This 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](34)
D
y
ss
=
1
2
π
∫
-
W
B
W
B
1
1
+
T
ω
2
d
ω
=
1
π
T
arctan
⁡
W
B
T
.The OP for this linear system is(35)
A
G
=
T
I
+
A
1
-
1
A
1
,where I is the identity matrix; A
1 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) ConsiderW
B
=
π; T
=
2. Figure 2 shows the output variance obtained by the proposed method for W
B
=
π. 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, U
(
t
), 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.Table 1
Simulation parameters and time profiles for obtaining the statistical characteristics by the MC, gPC, and proposed methods in Examples 1(a),1(b), 2(a), 2(b), 3, and 4(a).
Example
Simulation parameters
Computational time (sec.)
MC (Halton sampling)
gPC
Proposed
MC
gPC
Proposed
Example1(a)
10000 samples
400 cubature nodes
512 basis functions
205.30
1.10
0.01
Example1(b)
10000 samples
1000 cubature nodes
512 basis functions
197.18
2.67
0.01
Example2(a)
10000 samples
N/A
512 basis functions
199.64
N/A
0.01
Example2(b)
10000 samples
N/A
2048 basis functions
827.76
N/A
3.22
Example3
10000 samples
N/A
512 basis functions
210.66
N/A
0.01
Example4(a)
10000 samples
625 cubature nodes
512 basis functions
4503.90
609.51
15.32Figure 2
Variances of the output in Example 1(a).(b) ConsiderW
B
=
4
π; T
=
2. Figure 3 presents the variance obtained by the proposed method for W
B
=
4
π. 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.Figure 3
Variances of the output in Example 1(b).Table1 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:(36)
D
t
ρ
α
y
t
=
u
,where ρ
(
α
)
=
a
1
δ
(
α
-
α
1
)
+
a
2
δ
(
α
-
α
2
) and δ is the Dirac delta function. Therefore, (36) is actually a double fractional integrator,(37)
a
1
D
0
α
1
y
t
+
a
2
D
0
α
2
y
t
=
u
t
.The case where the input u
(
t
) is an ideal white noise with a zero mean and covariance function was considered:(38)
κ
U
U
τ
=
δ
τ
=
δ
t
1
-
t
2
.The exact variance of the output is given by [40](39)
D
y
t
=
σ
2
t
=
1
a
2
2
∫
0
t
u
2
α
2
-
1
E
α
2
-
α
1
,
α
2
-
a
1
u
α
2
-
α
1
a
2
2
d
u
,where E
α
2
-
α
1
,
α
2 is the Mittag-Leffler function, which can be calculated using the Matlab mlf.m function [41].The OP for system (37) is given by(40)
A
G
=
∑
i
=
1
2
B
α
i
a
i
-
1
,where B
α
i is the OP of derivative of order α
i.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 a
1
=
a
2
=
1, α
1
=
3
/
4, and α
2
=
1. 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.Figure 4
Variances of the output in Example 2(a) witha
1
=
a
2
=
1, α
1
=
3
/
4, and α
2
=
1.Figure 5
Absolute errors by the proposed and MC methods in Example 2(a).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]:(41)
D
t
ρ
α
y
t
=
u
,
ρ
α
=
1
,
0
≤
α
≤
1
.Again, this study considered the case where the input u
(
t
) 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]:(42)
D
y
t
=
σ
2
t
=
∫
0
t
e
u
E
1
u
2
d
u
,
E
1
u
=
∫
u
∞
e
y
y
d
y
.The OP for system (41) is given by(43)
A
G
=
∑
i
=
1
5
B
α
i
w
i
-
1
,where B
α
i is the OP of derivative of order α
i and  α
i
,
w
i
i
=
1
5 are the nodes and weights from the Legendre quadrature.Figure6 shows the variance of the output by the proposed method. The absolute error with respect to the exact variance (42) is also shown.Figure 6
Output variance and absolute errors by the proposed and MC methods in Example 2(b).(c)  Example 2(c):
P
I
λ
D
μ
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 aP
I
λ
D
μ controller as a closed loop configuration with the uniform DO integrator above. The set point r
(
t
) is a random process with a unit mean and covariance function κ
R
R
(
τ
)
=
0.01
δ
(
τ
). This input r
(
t
) 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 Section4. The search space for the optimal parameters of the controller was limited to 0
≤
K
p
≤
5; 0
≤
K
i
≤
5; 0
≤
K
d
≤
5; 0.1
≤
λ
≤
1.9; 0.1
≤
μ
≤
1.9 for simplicity, as in other studies on the probabilistic approach [29, 37]. The resulting controller can be expressed as C
1
(
s
)
=
0.2805
+
1.1458
/
s
0.6422.Figure7 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.Figure 7
Mean and variances of the output of the system in Example 2(c) by the proposed controller.Figure 8
500 responses of stochastic system in Example 2(c) by the proposed controller.(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 meanm
Y
des
⁡
=
0.1
t
(
10
-
t
). The following iterative learning control scheme can be used for refining the mean of set point input:(44)
m
Y
k
t
=
ψ
t
T
C
m
Y
=
ψ
t
T
E
A
cl
⁡
C
m
R
k
,
e
k
t
=
m
Y
des
⁡
t
-
m
Y
k
t
=
ψ
t
T
C
m
e
k
,
m
R
k
+
1
t
=
m
R
k
t
+
e
k
t
=
ψ
t
T
C
m
R
k
+
0.5
ψ
t
T
C
m
e
k
,where k
∈
{
0,1
,
2
,
…
} and A
G is the closed loop OP.Figure9 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 k 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.Figure 9
Mean and variances of the output of the system in Example 2(d). Line with red o: desired mean.Figure 10
512 responses of stochastic system in Example 2(d) with iterative learning control algorithm.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 eachH
∈
(
0,1
), there is a real-value Gaussian process B
H
(
t
)
,
t
≥
0 such that M
B
H
(
t
)
=
E
B
H
(
t
)
=
0 and the covariance function of E
B
H
(
t
)
B
H
(
s
)
=
1
/
2
[
t
2
H
+
s
2
H
-
t
-
s
2
H
], s
,
t
∈
R
+. This process is called standard fractional Brownian motion with the Hurst parameter H. If H
=
1
/
2, 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 processX
(
t
) that satisfies the following:(45)
X
˙
t
=
-
X
t
+
B
H
t
,
X
0
=
0
,where B
H
(
t
) is a fractional Brownian motion with the Hurst parameter H. The mean of X
(
t
) is M
X
(
t
)
=
0. The variance of X
(
t
) satisfies a differential equation [44]:(46)
D
˙
X
t
=
-
2
D
X
t
+
2
d
t
,
t
,where d
(
t
,
t
) is given by(47)
d
s
,
t
=
1
2
t
2
H
1
-
e
-
s
+
j
1
s
+
j
2
s
,
t
,
j
1
s
=
e
-
s
∫
0
s
u
2
H
e
u
d
u
,
j
2
s
,
t
=
e
-
s
∫
0
s
t
-
u
2
H
e
u
d
u
.Figure 11 shows the evolution of variance of X
(
t
) for H
=
0.6 obtained from (46) and (47).Figure 11
Variances ofX
(
t
) in Example 3.Equation (45) can be rewritten in terms of OP as follows:(48)
C
X
=
I
+
B
-
1
C
B
=
A
G
C
B
,where B is the OP of the derivative of order one. Therefore, one can easily obtain the covariance of the random process, X
(
t
), by utilizing the OP for this system, A
G
=
(
I
+
B
)
-
1, and covariance function of fractional Brownian motion and (25). Figure 11 shows the variance of X
(
t
) obtained using the proposed method. Finally, the variance of X
(
t
) 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 as(49)
G
s
=
Y
s
U
s
=
k
τ
∫
0
1
s
α
d
α
+
1
,where k and τ are uniform random variables in [
0.5,1.5
]. The system is in a closed loop configuration with a fractional PI controller, C
1
(
s
)
=
0.187
+
36.35
/
s
1.216 from [46]. Note that the controller, C
1
(
s
), was designed for a nominal system, that is, for G
=
1
/
(
s
0.5
+
1
). The input R
(
t
) is a band-limited white noise process with a unit mean and covariance function of the following: κ
R
R
(
t
1
,
t
2
)
=
0.02
sinc
⁡
(
(
t
1
-
t
2
)
W
B
/
π
), where W
B
=
0.02
π.Figure12 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, U
(
t
), 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.Figure 12
Mean and variances of the output of the system in Example 4(a).(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 0
≤
K
p
≤
5; 0
≤
K
i
≤
50; 0
≤
K
d
≤
5; 0.5
≤
λ
≤
1.5; 0.1
≤
μ
≤
1.5. The result is as follows: C
2
(
s
)
=
5
+
50
/
(
s
)
+
5
s
0.1.Figure13 shows the mean and variance of the system output by the proposed controller and the controller, C
1
(
s
)
=
0.187
+
36.35
/
s
1.216 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.Figure 13
Mean and variance of system output for Example 4(b) by the proposed and initial controllers.Figure 14
Bounded regions for 1000 MC simulations of the stochastic system by the proposed and initial controllers in Example 4(b).


## 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,G
(
s
)
=
1
/
(
1
+
T
s
) from [38] with a band-limited white noise as the input.Let the input have a zero mean and covariance function ofκ
U
U
(
t
1
,
t
2
)
=
W
B
/
π
sinc
⁡
(
(
t
1
-
t
2
)
W
B
/
π
), where the sinc function is defined as (31)
sinc
⁡
x
=
sin
⁡
π
x
π
x
elsewhere
1
for
x
=
0
.The power spectral density of the input is(32)
S
U
ω
=
1
2
π
ω
≤
W
B
0
ω
>
W
B
,where W
B is known as bandwidth of the noise. As W
B approaches to infinity, the process will become ideal white noise.Therefore, the power spectral density function of stochastic output is given by(33)
S
Y
ω
=
1
2
π
1
1
+
T
ω
2
ω
≤
W
B
0
ω
>
W
B
.This 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](34)
D
y
ss
=
1
2
π
∫
-
W
B
W
B
1
1
+
T
ω
2
d
ω
=
1
π
T
arctan
⁡
W
B
T
.The OP for this linear system is(35)
A
G
=
T
I
+
A
1
-
1
A
1
,where I is the identity matrix; A
1 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) ConsiderW
B
=
π; T
=
2. Figure 2 shows the output variance obtained by the proposed method for W
B
=
π. 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, U
(
t
), 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.Table 1
Simulation parameters and time profiles for obtaining the statistical characteristics by the MC, gPC, and proposed methods in Examples 1(a),1(b), 2(a), 2(b), 3, and 4(a).
Example
Simulation parameters
Computational time (sec.)
MC (Halton sampling)
gPC
Proposed
MC
gPC
Proposed
Example1(a)
10000 samples
400 cubature nodes
512 basis functions
205.30
1.10
0.01
Example1(b)
10000 samples
1000 cubature nodes
512 basis functions
197.18
2.67
0.01
Example2(a)
10000 samples
N/A
512 basis functions
199.64
N/A
0.01
Example2(b)
10000 samples
N/A
2048 basis functions
827.76
N/A
3.22
Example3
10000 samples
N/A
512 basis functions
210.66
N/A
0.01
Example4(a)
10000 samples
625 cubature nodes
512 basis functions
4503.90
609.51
15.32Figure 2
Variances of the output in Example 1(a).(b) ConsiderW
B
=
4
π; T
=
2. Figure 3 presents the variance obtained by the proposed method for W
B
=
4
π. 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.Figure 3
Variances of the output in Example 1(b).Table1 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:(36)
D
t
ρ
α
y
t
=
u
,where ρ
(
α
)
=
a
1
δ
(
α
-
α
1
)
+
a
2
δ
(
α
-
α
2
) and δ is the Dirac delta function. Therefore, (36) is actually a double fractional integrator,(37)
a
1
D
0
α
1
y
t
+
a
2
D
0
α
2
y
t
=
u
t
.The case where the input u
(
t
) is an ideal white noise with a zero mean and covariance function was considered:(38)
κ
U
U
τ
=
δ
τ
=
δ
t
1
-
t
2
.The exact variance of the output is given by [40](39)
D
y
t
=
σ
2
t
=
1
a
2
2
∫
0
t
u
2
α
2
-
1
E
α
2
-
α
1
,
α
2
-
a
1
u
α
2
-
α
1
a
2
2
d
u
,where E
α
2
-
α
1
,
α
2 is the Mittag-Leffler function, which can be calculated using the Matlab mlf.m function [41].The OP for system (37) is given by(40)
A
G
=
∑
i
=
1
2
B
α
i
a
i
-
1
,where B
α
i is the OP of derivative of order α
i.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 a
1
=
a
2
=
1, α
1
=
3
/
4, and α
2
=
1. 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.Figure 4
Variances of the output in Example 2(a) witha
1
=
a
2
=
1, α
1
=
3
/
4, and α
2
=
1.Figure 5
Absolute errors by the proposed and MC methods in Example 2(a).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]:(41)
D
t
ρ
α
y
t
=
u
,
ρ
α
=
1
,
0
≤
α
≤
1
.Again, this study considered the case where the input u
(
t
) 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]:(42)
D
y
t
=
σ
2
t
=
∫
0
t
e
u
E
1
u
2
d
u
,
E
1
u
=
∫
u
∞
e
y
y
d
y
.The OP for system (41) is given by(43)
A
G
=
∑
i
=
1
5
B
α
i
w
i
-
1
,where B
α
i is the OP of derivative of order α
i and  α
i
,
w
i
i
=
1
5 are the nodes and weights from the Legendre quadrature.Figure6 shows the variance of the output by the proposed method. The absolute error with respect to the exact variance (42) is also shown.Figure 6
Output variance and absolute errors by the proposed and MC methods in Example 2(b).(c)  Example 2(c):
P
I
λ
D
μ
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 aP
I
λ
D
μ controller as a closed loop configuration with the uniform DO integrator above. The set point r
(
t
) is a random process with a unit mean and covariance function κ
R
R
(
τ
)
=
0.01
δ
(
τ
). This input r
(
t
) 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 Section4. The search space for the optimal parameters of the controller was limited to 0
≤
K
p
≤
5; 0
≤
K
i
≤
5; 0
≤
K
d
≤
5; 0.1
≤
λ
≤
1.9; 0.1
≤
μ
≤
1.9 for simplicity, as in other studies on the probabilistic approach [29, 37]. The resulting controller can be expressed as C
1
(
s
)
=
0.2805
+
1.1458
/
s
0.6422.Figure7 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.Figure 7
Mean and variances of the output of the system in Example 2(c) by the proposed controller.Figure 8
500 responses of stochastic system in Example 2(c) by the proposed controller.(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 meanm
Y
des
⁡
=
0.1
t
(
10
-
t
). The following iterative learning control scheme can be used for refining the mean of set point input:(44)
m
Y
k
t
=
ψ
t
T
C
m
Y
=
ψ
t
T
E
A
cl
⁡
C
m
R
k
,
e
k
t
=
m
Y
des
⁡
t
-
m
Y
k
t
=
ψ
t
T
C
m
e
k
,
m
R
k
+
1
t
=
m
R
k
t
+
e
k
t
=
ψ
t
T
C
m
R
k
+
0.5
ψ
t
T
C
m
e
k
,where k
∈
{
0,1
,
2
,
…
} and A
G is the closed loop OP.Figure9 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 k 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.Figure 9
Mean and variances of the output of the system in Example 2(d). Line with red o: desired mean.Figure 10
512 responses of stochastic system in Example 2(d) with iterative learning control algorithm.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 eachH
∈
(
0,1
), there is a real-value Gaussian process B
H
(
t
)
,
t
≥
0 such that M
B
H
(
t
)
=
E
B
H
(
t
)
=
0 and the covariance function of E
B
H
(
t
)
B
H
(
s
)
=
1
/
2
[
t
2
H
+
s
2
H
-
t
-
s
2
H
], s
,
t
∈
R
+. This process is called standard fractional Brownian motion with the Hurst parameter H. If H
=
1
/
2, 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 processX
(
t
) that satisfies the following:(45)
X
˙
t
=
-
X
t
+
B
H
t
,
X
0
=
0
,where B
H
(
t
) is a fractional Brownian motion with the Hurst parameter H. The mean of X
(
t
) is M
X
(
t
)
=
0. The variance of X
(
t
) satisfies a differential equation [44]:(46)
D
˙
X
t
=
-
2
D
X
t
+
2
d
t
,
t
,where d
(
t
,
t
) is given by(47)
d
s
,
t
=
1
2
t
2
H
1
-
e
-
s
+
j
1
s
+
j
2
s
,
t
,
j
1
s
=
e
-
s
∫
0
s
u
2
H
e
u
d
u
,
j
2
s
,
t
=
e
-
s
∫
0
s
t
-
u
2
H
e
u
d
u
.Figure 11 shows the evolution of variance of X
(
t
) for H
=
0.6 obtained from (46) and (47).Figure 11
Variances ofX
(
t
) in Example 3.Equation (45) can be rewritten in terms of OP as follows:(48)
C
X
=
I
+
B
-
1
C
B
=
A
G
C
B
,where B is the OP of the derivative of order one. Therefore, one can easily obtain the covariance of the random process, X
(
t
), by utilizing the OP for this system, A
G
=
(
I
+
B
)
-
1, and covariance function of fractional Brownian motion and (25). Figure 11 shows the variance of X
(
t
) obtained using the proposed method. Finally, the variance of X
(
t
) 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 as(49)
G
s
=
Y
s
U
s
=
k
τ
∫
0
1
s
α
d
α
+
1
,where k and τ are uniform random variables in [
0.5,1.5
]. The system is in a closed loop configuration with a fractional PI controller, C
1
(
s
)
=
0.187
+
36.35
/
s
1.216 from [46]. Note that the controller, C
1
(
s
), was designed for a nominal system, that is, for G
=
1
/
(
s
0.5
+
1
). The input R
(
t
) is a band-limited white noise process with a unit mean and covariance function of the following: κ
R
R
(
t
1
,
t
2
)
=
0.02
sinc
⁡
(
(
t
1
-
t
2
)
W
B
/
π
), where W
B
=
0.02
π.Figure12 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, U
(
t
), 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.Figure 12
Mean and variances of the output of the system in Example 4(a).(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 0
≤
K
p
≤
5; 0
≤
K
i
≤
50; 0
≤
K
d
≤
5; 0.5
≤
λ
≤
1.5; 0.1
≤
μ
≤
1.5. The result is as follows: C
2
(
s
)
=
5
+
50
/
(
s
)
+
5
s
0.1.Figure13 shows the mean and variance of the system output by the proposed controller and the controller, C
1
(
s
)
=
0.187
+
36.35
/
s
1.216 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.Figure 13
Mean and variance of system output for Example 4(b) by the proposed and initial controllers.Figure 14
Bounded regions for 1000 MC simulations of the stochastic system by the proposed and initial controllers in Example 4(b).


## 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.



---

*Source: 989542-2015-10-05.xml*