Schulz M. Control Theory in Physics and Other Fields of Science: Concepts, Tools, and Applications

Подождите немного. Документ загружается.

8.3 L´evy Processes 227

If the initial probability distribution functions have diﬀerent tails, C

=

−

but equal exponents, p

(ξ) converges to the asymmetric L´evy distribution

with the exponent γ, the asymmetry parameter (8.59)anda ∼ R(C

−

)/2

If the asymptotic exponents γ

of the elementary probability density p (ξ)are

diﬀerent but min (γ

,γ

−

) < 2, the convergence is to a completely asymmetric

L´evy distribution with an exponent γ = min (γ

,γ

−

)andb =1forγ

−

<γ

or b = −1forγ

−

>γ

Finally, upon a suﬃciently large number of convolutions, the Gaussian

distribution attracts also all the probability distribution functions decaying

as or faster than |ξ|

−3

at large |ξ|. Therefore, L´evy laws with γ<2are

sometimes denoted as true L´evy laws.

Unfortunately, all L´evy distributions with γ<2 have inﬁnite variances.

That limits its physical, but not its mathematical, meaning. Physically, L´evy

distributions are meaningless with respect to ﬁnite systems. But in complex

systems with an almost unlimited reservoir of hidden irrelevant degree of

freedom, such probability distribution functions are quite possible at least

over a wide range of the stochastic variables. Well-known examples of such

wild distributions [13, 41] have been found to quantify the velocity-length

distribution of the fully developed turbulence (Kolmogorov law) [14, 20, 21],

the size–frequency distribution of earthquakes (Gutenberg–Richter law) [25,

26], or the destruction losses due to storms [22]. Further examples related

to social and economic problems are the distribution of wealth [23, 24] also

known as Pareto law, the distribution of losses due to business interruption

resulting from accidents [15, 16] in the insurance business, or the distribution

of losses caused by ﬂoods worldwide [17] or the famous classical St. Petersburg

paradox discussed by Bernoulli [18, 19]

8.3.3 Truncated L´evy Distributions

As we have seen, L´evy laws obey scaling relations but have an inﬁnite variance.

A real L´evy distribution is not observed in ﬁnite physical systems. However, a

stochastic process with ﬁnite variance and characterized by scaling relations in

a large but ﬁnite region close to the center is the truncated L´evy distribution

[43]. For many realistic problems, we have to ask for a distribution which in

the tails is a power law multiplied by an exponential

p (ξ) ∼

|ξ|

γ+1

exp



−

|ξ|



. (8.62)

The characteristic function of L´evy laws truncated by an exponential as in

(8.62) can be written explicitly as [42, 43]

ln ˆp (k)=a



1+k



γ/2

cos (γ arctan (kξ

)) − 1

cos (πγ/2)



1+ib tan (γ arctan (|k|ξ

))

|k|



. (8.63)

228 8 Filters and Predictors

After R convolutions we get the characteristic distribution function

ln ˆp

(k)=−Ra



1+k



γ/2

cos (γ arctan (kξ

)) − 1

cos (πγ/2)



1+ib tan (γ arctan (|k|ξ

))

|k|



. (8.64)

It can be checked that (8.63) recovers (8.53)forξ

→∞. The behavior of

(ξ) can be obtained from an inverse Fourier transform (6.51). In order

to determine the characteristic scale of the probability distribution p

(ξ), we

have to consider the main contributions to the inverse Fourier transform. This

condition requires that the characteristic wave-number k

char

is of an order of

magnitude satisfying ln ˆp

char

)  1. This relation is equivalent to





+ ξ

−2



γ/2

− ξ

−γ



 1 . (8.65)

For R  ξ

,(8.65) is satisﬁed if k

char

 1. Thus we obtain immediately

char

∼ (Ra)

−1/γ

and therefore the characteristic scale ξ

char

∼ (Ra)

1/γ

,which

characterizes an ideal L´evy distribution.

When, on the contrary, R  ξ

, the characteristic value of k

char

becomes

much smaller than ξ

−1

, and we ﬁnd now the relation k

char

∼ (Ra)

−1/2

γ/2−1

The characteristic scale ξ

char

∼ (Ra)

1/2

1−γ/2

corresponding to what we ex-

pect from the Gaussian behavior.

Hence, as expected, a truncated L´evy distribution is not stable. It ﬂows

to an ideal L´evy probability distribution function for small R and then to

the Gaussian distribution for large R. The crossover from the initial L´evy-like

regime to the ﬁnal Gaussian regime occurs if the characteristic scale of the

L´evy distribution reaches the truncation scale ξ

char

∼ ξ

, i.e., if Ra ∼ ξ

8.4 Rare Events

8.4.1 The Cram´er Theorem

The central limit theorem states that the Gaussian law is a good description

of the center of the probability distribution function p

(ξ) for suﬃciently

large R. We have demonstrated that the range of the center increases with

increasing R but it is always limited for ﬁnite R. A similar statement is valid

for the generalized version of the central limit theorem regarding the conver-

gence behavior of L´evy laws. Fluctuations exceeding the range of the center

are denoted as large ﬂuctuations.

Of course, large ﬂuctuations are rare events. The behavior of possible large

ﬂuctuations is not, or is only partially, aﬀected by the predictions of the central

limit theorem so that we should ask for an alternative description. We start

our investigation from the general formulae (8.12) for a one-component event.

8.4 Rare Events 229

The characteristic function can also be calculated for an imaginary k → iz so

that the Fourier transform becomes a Laplace transform

ˆp (z)=



dξp(ξ)exp{−zξ} , (8.66)

which holds under the assumption that the probability distribution function

decays faster than an exponential for |ξ|→∞. We obtain again an algebraic

relation for R convolution of the elementary probability distribution function

p (ξ),

ˆp

(z)=[ˆp (z)]

. (8.67)

On the other hand, we assume that for suﬃciently large R the probability

density p

(ξ) may be written as

(ξ) = exp



−RC





, (8.68)

where C (x) is the Cram´er function [44, 45]. We will check by a construction

principle, whether such a function exists for the limit R →∞.Tothisaimwe

calculate the corresponding Laplace transform

ˆp

(z)=R



dx exp {−R [C (x)+zx]} (8.69)

by using the method of steepest descent. This method approximates the in-

tegral by the value of the integrand in a small neighborhood around its max-

imum ˜x. The value of ˜x depends not on R and is a solution of

∂

∂˜x

C (˜x)+z =0. (8.70)

With the knowledge of ˜x we can expand the Cram´er function in powers of x

around ˜x

C (x)+zx = C (˜x)+z˜x +

∂

∂˜x

C (˜x)[x − ˜x]

+ ··· . (8.71)

Note that the ﬁrst-order term vanishes because of (8.70). Substituting (8.71)

into (8.69), we obtain the integral

ˆp

(z)=R exp {−R [C (˜x)+z˜x]}



dy exp



−R



∂

C (˜x)

∂˜x

+ ···



(8.72)

with y = x − ˜x. The leading term in the remaining integral is a Gaussian

law of width δy ∼ R

−1/2

. With respect to this width all other contributions

of the series expansion can be neglected for R →∞. Therefore, we focus

here in the second-order term. The corresponding Gaussian integral exists if

∂

C/∂x

> 0. In this case we obtain

ˆp

(z) ∼

R/ (∂

C (˜x) /∂ ˜x

)exp{−R [C (˜x)+z ˜x]} . (8.73)

230 8 Filters and Predictors

For R →∞, the leading term of the characteristic function is given by

ˆp

(z) ∼ exp {−R [C (˜x)+z˜x]} . (8.74)

Combining (8.67), (8.74), and (8.70), we obtain the equations

C (˜x)+z˜x +lnˆp (z)=0 and

∂

∂˜x

C (˜x)+z = 0 (8.75)

which allow the determination of

♥

C (x). These two equations indicate that

the Cram´er function is the Legendre transform of ln ˆp (z). Hence, in order to

determine C (˜x)wemustﬁndthevalueofz which corresponds to a given ˜x.

The diﬀerentiation of (8.75) with respect to ˜x leads to

∂

∂˜x

C (˜x)+z +˜x

∂z

∂˜x

∂ ln ˆp (z)

∂z

∂˜x



˜x +

∂ ln ˆp (z)

∂z



∂z

∂˜x

=0. (8.76)

Because of ∂z/∂˜x = −∂

C (˜x) /∂ ˜x

< 0 (see above), we ﬁnd the relation

˜x = −

∂ ln ˆp (z)

∂z

(8.77)

from where we can calculate z = z(˜x). Having C (˜x), the Cram´er theorem

reads

(ξ) = exp



−RC





for R →∞. (8.78)

This theorem describes large ﬂuctuations outside the central region of ˆp

(ξ).

The central region is deﬁned by the central limit theorem, which requires

ξ ∼ R

with α<1(see(8.36) and (8.37)). Thus, the central region collapses

to the origin in the Cram´er theorem.

But outside of the center we have |ξ|/R > 0. Obviously, the scaling of the

variables diﬀers between the central limit theorem and the Cram´er theorem.

While the rescaling ξ → ξ/

√

R leads to the form-stable Gaussian behavior

of p

(ξ) in the limit R →∞, the rescaling ξ/R yields another kind of form

stability concerning the expression R

−1

ln p

(ξ).

Furthermore, the properties of the initial elementary probability distribu-

tion disappear close to the center for R →∞. Therefore, the central limit the-

orem describes a universal phenomenon. The Cram´er function conserves the

properties of the elementary probability distribution functions due to (8.75)

so that the large ﬂuctuations show no universal behavior.

8.4.2 Extreme Fluctuations

The Cram´er theorem provides a concept for the treatment of large ﬂuctuation

as a sum of an inﬁnite number of successive events. This limit R →∞cor-

responds to the fact that the rescaled accumulated ﬂuctuations ξ/R remains

ﬁnite.

Another important regime is the extreme ﬂuctuation regime [46]. Here we

have to deal with ﬁnite R but ξ/R →∞.

8.4 Rare Events 231

In order to quantify this class of ﬂuctuations, we start again from (8.12)

and consider one-component events.

We use the representation p (ξ) = exp {−f (ξ)} and obtain

(ξ)=





j=1

dξ





ξ −



j=1





exp







−



j=1

f (ξ

)







. (8.79)

In order to simplify, we restrict ourselves on the case of an extreme positive

ﬂuctuation ξ → +∞. We have now two possibilities. On the one hand, the

asymptotic behavior of the function f (ξ) can be concave. Then we have f(x)+

f(y) >f(x + y) so that the dominant contributions to (8.79) are obtained

from conﬁgurations with all ﬂuctuations are very small except of one extreme

ﬂuctuation being almost equal to ξ. Therefore, we get

ln p

(ξ) ∼ ln p (ξ) ∼−f (ξ) . (8.80)

On the other hand, if the asymptotic behavior of f (ξ)isconvex,f (x)+

f(y) <f(x + y), the minimum of the exponentials is given by the symmetric

conﬁguration ξ

= ξ/R for all j =1,...,R. The convexity condition requires

a global minimum of the sum of all exponentials in (8.79) so that



j=1

f (ξ

) ≥ Rf





. (8.81)

We apply again the method of the steepest descent. To this aim we introduce

the deviations δξ

= ξ

− ξ/R and expand the sum in (8.81) around its

minimum



j=1

f (ξ

)=Rf













j=1

(δξ

)

+ o



|δξ|



, (8.82)

where we have used the constraint δξ

+δξ

+···+δξ

= 0. We substitute this

expression into (8.79). Then, with the assumption of convexity, f



(ξ/R) > 0,

the integral (8.79) can be estimated. We get the leading term

(ξ) ∼ exp



−Rf





. (8.83)

This approximate result approaches the true value for ξ/R →∞. Apparently,

(8.83) and (8.68) are identical expressions. But we should reminder that (8.68)

holds for R →∞but ﬁnite ξ/R, while (8.83) requires ξ/R →∞. However,

the Cram´er function C (x) becomes equal to f (x)forx →∞.

In summary, the knowledge of the tails of an elementary probability distri-

bution p (ξ) allows the determination of the tails of the probability distribution

function p

(ξ)via

(ξ) ∼







(8.84)

232 8 Filters and Predictors

if ln p

−1

(ξ) is a convex function in ξ. On the other hand, if ln p

−1

(ξ) is concave,

we get p

(ξ) ∼ p (ξ)forξ/R →∞.

8.5 Kalman Filter

8.5.1 Linear Quadratic Problems with Gaussian Noise

Let us now study a stochastic system under control which is described by the

linear evolution equation of the type (8.3).

X(t)=A(t)X(t)+B(t)u(t)+ξ(t) , (8.85)

where ξ(t)istheN-component noise vector modeling the uncertainty of the

system. Because of the central limit theorem (see Sect. 8.2), the probability

distribution functions of the components of ξ(t) are assumed to be those of a

Gaussian stochastic process. Furthermore, we introduce a p-component output

Y (t)=C(t)X(t)+η(t) , (8.86)

where C(t) is a matrix of type p × N and η(t) represents the p-component

Gaussian random observation error. Both noise vectors have zero mean

ξ(t)=0 and η(t) = 0 (8.87)

while the correlation functions are given by

(t)ξ



)=Ω

αβ

(t)δ(t − t



) η

(t)η



)=Θ

αβ

(t)δ(t − t



) (8.88)

and

(t)η



)=0. (8.89)

The initial value of the state vector, X

= X(0), may have the mean X

while

the covariance matrix is given by

− X

)

− X

)

= σ

αβ

. (8.90)

Obviously, we have a double problem. The ﬁrst part must be the reconstruc-

tion of the state X(t) from the knowledge of the observations Y (t) while the

second problem is the control of the system.

8.5.2 Estimation of the System State

The problem of the optimal estimate of the state of a system from the available

observations is also called a ﬁltering procedure. To solve this problem, we split

the state and the observation variable

X(t)=

X(t)+X

(t)andY (t)=

Y (t)+Y

(t) (8.91)

with

8.5 Kalman Filter 233

X(t)=A(t)

X(t)+ξ(t)and

(t)=A(t)X

(t)+B(t)u(t) (8.92)

and

Y (t)=C(t)

X(t)+η(t)andY

(t)=C(t)X

(t) (8.93)

while the initial conditions are

X(0) = δX

and X

(0) = X

. (8.94)

Note that the initial ﬂuctuations are given by δX

= X

− X

. Now we

consider the evolution of

X(t) and try to reconstruct this state at the current

time t from the knowledge of

Y (t



) with t



<t. To this aim we deﬁne a certain

basic {e

,...,e

} spanning the phase space P. Then the projections of the

current state onto this basic

(t)=

X(t)e

(8.95)

(k =1,...,N) represent the dynamics of the system completely. On the other

hand, we may introduce the scalar quantities

(t)=





)

Y (t



) (8.96)

(k =1,...,N) and ask for certain p-component vector functions Λ

(t) satis-

fying the N minimum problems

(t)=

(t) − θ

(t))

→ min (8.97)

at the current time

t>0. It means that we have decomposed the ﬁlter-

ing problem into N separate minimum problems leading to N optimal p-

component ﬁlter functions Λ

. In order to solve these minimum problems, we

consider the diﬀerential equations (k =1,...,N)



)=−A



)+C



)Λ



) (8.98)

for t



∈ [0,t] with matrices A(t



)andC(t



)from(8.85) and (8.86), respectively,

and the ﬁnal condition

Z(t)=e

. (8.99)

We transform this equation by the application of (8.85) and (8.86):

+ ξZ

YΛ

− ηΛ

+ ξZ

. (8.100)

By integrating both sides of this equation between 0 and t, we get with (8.99),

(8.95), and (8.96)

Note that the initial time is t =0.

234 8 Filters and Predictors

(t) − θ

(t)=δX

(0) +





[ξ(t



) − η(t



)Λ



)] . (8.101)

Hence, by squaring both sides, performing the average and considering (8.88),

we obtain

(t)=

(0)δX

δX

(0)











)ξ(t



)ξ



)











)η(t



)η



)Λ



)

(0)σZ

(0)





)ΩZ



)+Λ



)ΘΛ



)] . (8.102)

Thus, the ﬁltering problem is reduced to a deterministic linear quadratic con-

trol problem with performance J

(t), the constraints (8.98), and the ﬁnal

conditions (8.99). However, the roles of the ﬁnal and initial times have been

interchanged in the performance functional. This fact can be managed by the

reﬂection of the time direction. Then the comparison of the problem with the

results of Sect. 3.1.4 requires now the solution

= Θ

−1

CGZ

, (8.103)

where the symmetric N × N matrix G is a solution of the Ricatti equation

G − GA

− AG + GC

−1

CG = Ω (8.104)

with the initial condition

G (0) = σ (8.105)

while the function Z

is the solution of



−1

CG − A

(8.106)

with the ﬁnal condition

(t)=e

. (8.107)

Thus the wanted estimation #x

(t) of the state vector with respect to the basic

vector e

is given by

The changed sign is also a consequence of the mentioned time reﬂection.

8.5 Kalman Filter 235

(t)=θ

(t)=





)

Y (t



)





)G(t



)Θ

−1

Y (t



)

= e





,t)G(t



)Θ

−1

Y (t



) (8.108)

and therefore

X(t)=





,t)G(t



)Θ

−1

Y (t



) . (8.109)

Here, Γ (t



,t) is the Green’s function solving the diﬀerential equation (8.106).

This solution may be formally written as

Γ (t



,t) = exp













(τ) Θ

−1

C (τ) G (τ ) − A

(τ)

dτ







. (8.110)

Thus we get



,t) = exp







−







G (τ) C

(τ) Θ

−1

C (τ) − A (τ)

dτ







Γ (t, t



) , (8.111)

where

Γ (t, t



) is Green’s function associated to A−GC

−1

C. In other words,

the optimal estimation

X(t) satisﬁes the diﬀerential equation

X(t)=



A − GC

−1

X(t)+G(t)C

(t)Θ

−1

Y (t) (8.112)

with the initial condition

X(0) = 0. We remark that the optimal estimation

depends essentially on the strength of the noise concerning the state evolution

of the system, the observation error, and the uncertainties of the initial state

(see (8.112) and (8.104)).

The estimation

X(t) is taken for u = 0. In order to obtain the estimation

(t) for the presence of a ﬁnite control, we must add the deterministic so-

lution X

(t), obtained from the solution of the second group of equations of

(8.92) and (8.93),

(t)=

X(t)+X

(t). Because of (8.92) and (8.93), the

complete estimation fulﬁlls the diﬀerential equations



A − GC

−1

+ Bu + GC

−1

Y, (8.113)

236 8 Filters and Predictors

where we have taken into account (8.91) for the elimination of

Y (t). In order

to complete the estimation equation (8.113), we have to consider the initial

condition

(0) = X

. (8.114)

Equation (8.113) is the famous Kalman ﬁlter equation [90, 91]. Since the

estimation of a state from the knowledge of continuous observations is a very

important problem for a large class of physical experiments, we will illustrate

the algorithm with a simple example which belongs to the standard classical

measurement. Let us assume that a one-component quantity X(t) follows a

well-deﬁned deterministic law

X(t)=u(t), where u(t) is the expected time-

dependent trajectory while the observations, Y (t)=X(t)+η(t), are ﬂawed

with an error η(t) of the variance Θ.Thus,wehaveA =0,B =1,C =1,

and Ω = 0. The Ricatti equation (8.104) becomes now

G + Θ

−1

= 0 with G(0) = σ, (8.115)

where σ is the variance of the initial data. The solution of this equation is

simply G = σΘ/(σt + Θ). Thus, we get the estimation equation

σt + Θ

= u +

σt + Θ

Y (8.116)

and therefore the optimal estimated state

(t)=



[σt



u(t



)+Θu(t



)+σY (t



)] + X

σt + Θ

. (8.117)

This formula gives the optimal estimation of state on the basis of the obser-

vations Y (t). In principle, the Kalman ﬁlter can be interpreted as a special

regression model. In order to obtain the deviations of the estimated state from

the real state, we substitute Y = X + η and u =

X in (8.117) and obtain after

integration by parts

(t)=X(t)+

σt + Θ







− X







ση





. (8.118)

Obviously, the estimation error,



(t) − X(t)



, is dominated by the integral

over the noise process. Hence, we get for a Gaussian noise the relative error



(t) − X(t)

X(t)



∼ t

−1/2

|X(t)|

−1/2

. (8.119)

Thus, the Kalman ﬁlter estimation leads to an asymptotically convergence to

the true behavior if X(t) decays not faster than t

−1

for t →∞.

Finally, we will deﬁne the meaning of the function G(t). To this aim we

make use of (8.106) and (8.104) and evaluate the derivative