Fabien B. Analytical System Dynamics: Modeling and Simulation

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

268 6 Dynamic System Analysis and Simulation

6.1.3 Lyapunov direct method

In the linearization method described above the stability of the dynamic sys-

tem is determined by examining the eigenvalues of the linearized equations of

motion. In this section the stability of the equilibrium is determined by exam-

ining a positive deﬁnite scalar ‘energy function’ for the system. To motivate

why such an analysis may be fruitful, consider an unforced dynamic system

consisting of ideal inductors, capacitors and resistors. In such a system any

initial energy (due to perturbation of the system) will be dissipated by the

resistors, and the system will eventually return to an equilibrium point. The

Lyapunov direct method involves the construction and analysis of suitable

scalar energy-like function for the system.

Before stating the Lyapunov stability theorems we need some deﬁnitions.

Let x ∈ R

be some vector, and Ω be region in R

that contains the origin.

Then;

• The scalar function V(x) is said to be locally positive deﬁnite if V(0) = 0

and V(x) > 0 for all x 6= 0 in Ω.

• The scalar function V(x) is said to be locally negative deﬁnite if V(0) = 0

and V(x) < 0 for all x 6= 0 in Ω.

• The scalar function V(x) is said to be locally positive semi-deﬁnite if

V(0) = 0 and V(x) ≥ 0 for all x 6= 0 in Ω.

• The scalar function V(x) is said to be locally negative semi-deﬁnite if

V(0) = 0 and V(x) ≤ 0 for all x 6= in Ω.

• If V(x) ≥ 0 for some values of x ∈ Ω, and V(x) ≤ 0 for other values of

x ∈ Ω then, V(x) is said to be indeﬁnite.

If Ω is the entire space R

then these deﬁnitions apply globally instead of

locally.

Lyapunov Stability Theorems

Consider the autonomous dynamic system

˙x(t) = f(x(t)), (6.11)

where the state x(t) ∈ R

, and x

∗

= 0 is an equilibrium. Let Ω be region in

that contains the origin, x

∗

= 0. Then, the following theorems give suﬃ-

cient conditions for x

∗

= 0 to be stable, or asymptotically stable equilibrium.

I Lyapunov stability theorem.

If for x ∈ Ω there is a continuous, scalar, positive-deﬁnite function V(x)

such that

dV(x)

V(x) ≤ 0.

Then, the equilibrium x

∗

= 0 is stable in the sense of Lyapunov.

6.1 System Analysis 269

II Lyapunov asymptotic stability theorem.

If for x ∈ Ω there is a continuous, scalar, positive-deﬁnite function V(x)

such that

V(x) < 0.

Then, the equilibrium x

∗

= 0 is asymptotically stable in the sense of Lya-

punov.

III Krasovsky asymptotic stability theorem.

If for x ∈ Ω there is a continuous, scalar, positive-deﬁnite function V(x)

such that

V(x) < 0, x /∈ Ω

V(x) = 0, x ∈ Ω

where Ω

⊂ Ω, (excluding the origin) is deﬁned by the scalar function

F (x) = 0. Then, the equilibrium x

∗

= 0 is asymptotically stable in the

sense of Lyapunov if

U(x) =



dF (x)



f(x) 6= 0, x ∈ Ω

If the conditions in Theorem I or Theorem II or Theorem III are satisﬁed

then, V(x) is called a Lyapunov function.

Theorem I indicates that if we can ﬁnd a scalar function V(x) that is

positive deﬁnite and

V(x) is negative semi-deﬁnite then, x

∗

= 0 is a stable

equilibrium point.

Theorem II indicates that if we can ﬁnd a scalar function V(x) that is pos-

itive deﬁnite and

V(x) is negative deﬁnite then, x

∗

= 0 is an asymptotically

stable equilibrium point.

Theorem III indicates that x

∗

= 0 is an asymptotically stable equilibrium

point if the following conditions hold;

• The scalar function V(x) is positive deﬁnite.

•

V(x) = 0 in some subset of Ω, say Ω

which is deﬁned by the surface

F (x) = 0.

•

V(x) is negative deﬁnite outside of Ω

• U(x) 6= 0, x ∈ Ω

. This condition implies that the system trajectory does

not stay in Ω

forever. Note that Ω

excludes the origin.

Since the system does not remain in Ω

for all time, and

V(x) is negative

deﬁnite outside of Ω

, it can be inferred that the systems ‘energy’, V(x), will

decrease until it reaches the equilibrium.

In Theorem II asymptotic stability requires a monotonic decrease in the

function V(x), while Theorem III allows V(x) to decrease piecewise fashion.

See Merkin (1997) for the proofs of these theorems.

Note that these stability theorems are all local and apply only for x ∈ Ω.

However, similar statements hold globally if in addition to the conditions

270 6 Dynamic System Analysis and Simulation

stated above, we have V(x) → ∞ as kxk → ∞. The utility of these theorems

will be demonstrated via the examples below.

Example 6.5.

The equations of motion for the unforced linear

mass-spring system shown here are

˙x

= x

˙x

= −

where x

is the position, x

the velocity, m is the mass, and k is the spring

stiﬀness. The equilibrium for this system is x

∗

= 0, x

∗

= 0.

Now, select

V(x) =

(kx

+ mx

)

as a Lyapunov function candidate. (The function V(x) is called a Lyapunov

function only if its satisﬁes one of the stability theorems given above.)

First, note that V(x) is a scalar positive deﬁnite function, and it is the

sum of the kinetic coenergy and the potential energy of the system.

Next,

V(x) = kx

˙x

+ mx

˙x

= kx

+ mx

(−

)

= 0.

Since, V(x) is positive deﬁnite and

V(x) = 0, the Lyapunov stability Theorem

I, indicates that the equilibrium is stable. The same result can be obtained

using the linearization method.

Example 6.6.

For the dynamic system

˙x = λx,

where λ < 0, the point x

∗

= 0 is the equilibrium.

Let the scalar positive deﬁnite function

V(x) =

be a Lyapunov function candidate. Then,

V(x) = x ˙x

= λx

< 0,

6.1 System Analysis 271

for all x 6= 0. Therefore, according to the stability Theorem II we can conclude

that the equilibrium, x

∗

= 0, is asymptotically stable.

Example 6.7.

The equations of motion for the unforced linear

mass-spring-damper system shown here are

˙x

= x

˙x

= −

−

where x

is the position, x

the velocity, m is the mass, k is the spring

stiﬀness, and b is the damping coeﬃcient. The equilibrium for this system is

∗

= 0, x

∗

= 0.

Let

V(x) =

(kx

+ mx

)

be the Lyapunov function candidate. Then,

V(x) = kx

˙x

+ mx

˙x

= kx

+ mx



−



= −bx

Since,

V is negative semi-deﬁnite we can conclude from the stability Theorem

I that the equilibrium is stable. However, using the stability Theorem III it

can be shown that, in fact, the equilibrium is asymptotically stable. To do

so let Ω

= {x ∈ R

| x

= 0}, excluding the origin. That is, Ω

is all points

in R

where x

= 0, except the origin. Noting this exception, it can be seen

that

V(x) = 0, for x ∈ Ω

, and we can deﬁne Ω

using the function

F (x) = x

= 0.

Therefore,

U(x) =



dF (x)



f(x)

= [ 0 1 ]

−

= −

−

6= 0, x ∈ Ω

272 6 Dynamic System Analysis and Simulation

That is, U(x) is not identically zero for all points in Ω

. Thus, the system

trajectory is such that x does not remain in Ω

for all time, and since

V(x)

is negative deﬁnite outside of Ω

the system eventually approaches the equi-

librium. Hence, from the stability Theorem III we can conclude that the

equilibrium is asymptotically stable.

6.1.4 Lagrange’s stability theorem

Let q(t) ∈ R

denote the generalized displacements of a dynamic system, and

let ˙q(t) be the corresponding ﬂows. The system is called a natural dynamic

system if the kinetic coenergy is of the form T

∗

(q, ˙q) =

˙q(t)

M(q(t)) ˙q(t),

where M(q(t)) ∈ R

n×n

is a symmetric positive deﬁnite matrix.

If the natural dynamic system is conservative then, the equilibrium and

stability of the system can be determined by evaluating the potential energy

function. In this case all the applied eﬀorts, including gravity forces, must be

included in the potential energy function.

Equilibrium

Let

= −

∂V (q)

∂q

, i = 1, 2, ···, n,

be the eﬀorts due to the capacitors and gravity. Here, V (q(t)) is the potential

energy function. Let

= −

∂D( ˙q)

∂ ˙q

, i = 1, 2, ···, n,

be the eﬀorts due to ideal resistors, where D( ˙q) is the dissipation function. Let

, i = 1, 2, ···, n be the eﬀorts due to sources. If the system is conservative

then the eﬀorts due to the resistors and sources are zero, i.e., e

= 0, and

= 0, i = 1, 2, ···, n.

In the equilibrium position the the virtual work done by all the eﬀorts

applied to the system must sum to zero, i.e.,

δW(q) =

i=1

+ e

)δq

i=1

−

∂V

∂q

δq

= 0.

6.1 System Analysis 273

Since, the displacements are independent this expression will hold if and only

∂V/∂q

= 0, i = 1, 2, ···, n.

Thus, a necessary condition for q

∗

to be an equilibrium of a conservative

natural system is that the potential energy function must be stationary at

∗

Stability

Writing the potential energy function as a Taylor series centered at q

∗

gives

V (q) = V (q

∗

) +

∂V (q

∗

)

∂q

(q − q

∗

) +

(q − q

∗

)

H(q

∗

)(q − q

∗

)

+ higher order terms,

where

H(q) =

∂

V (q)

∂q

is the Hessian of the potential energy function. At the equilibrium, q

∗

, we have

∂V (q

∗

)/∂q = 0. Also, we can assume that V (q

∗

) = 0, since a constant can be

added to the potential energy function without aﬀecting the eﬀorts applied to

the system. Then, in a neighborhood of the equilibrium, the potential energy

function can be approximated by

V (q) =

(q − q

∗

)

H(q

∗

)(q − q

∗

Now, suppose that the the Hessian matrix H(q

∗

) is positive deﬁnite at the

equilibrium q = q

∗

then, the potential energy function is a minimum at q

∗

Moreover, V (q) is positive deﬁnite in a neighborhood of the equilibrium.

To evaluate the stability of the equilibrium, consider the Lyapunov func-

tion candidate

V(q, ˙q) =

i=1

∂T

∗

∂ ˙q

˙q

− T

∗

+ V

˙q

M(q) ˙q + V (q).

Since, M (q) is positive deﬁnite, and V (q

∗

) is a minimum then, there is a

neighborhood of the point q = q

∗

, ˙q = 0 where V(q, ˙q) is positive deﬁnite.

Also,

V =

i=1

∂T

∗

∂ ˙q

˙q

− T

∗

+ V

274 6 Dynamic System Analysis and Simulation

i=1



∂T

∗

∂ ˙q



˙q

i=1

∂T

∗

∂ ˙q

¨q

−

i=1

∂T

∗

∂q

˙q

−

i=1

∂T

∗

∂ ˙q

¨q

i=1

∂V

∂q

˙q

i=1



∂T

∗

∂ ˙q

−

∂T

∗

∂q

∂V

∂q



˙q

= 0.

Thus, according to the Lyapunov stability Theorem I, the equilibrium q = q

∗

˙q = 0 is stable. Note that to obtain the result

V = 0, we have used the fact

that

∂T

∗

∂ ˙q

−

∂T

∗

∂q

∂V

∂q

= −

∂D

∂ ˙q

+ e

= 0, i = 1, 2, ···, n.

The result developed above can be summarized as follows.

Stability Theorem of Lagrange

If the potential energy function of a conservative, natural system, has a min-

imum at an equilibrium, then the equilibrium is stable in the sense of Lya-

punov.

The stability theorem of Lagrange provides suﬃcient conditions for the

equilibrium to be stable. The inverse of this result is given by two theorems

attributed to Lyapunov (Meirovitch (1970)). Speciﬁcally, for a conservative,

natural system, suﬃcient conditions for the equilibrium to be unstable are

given by the following results.

• Lyapunov’s theorem.

The equilibrium q

∗

is unstable if the potential energy function is a maxi-

mum at q

∗

• Extended Lyapunov’s theorem.

If the potential energy function has no minimum at the equilibrium point,

the equilibrium is unstable.

6.1 System Analysis 275

Example 6.8.

In the system shown on the right the mass m

is attached to a spring with stiﬀness k > 0. The

displacement of the spring, as measured from its

free length is x. The mass m

is attached to m

via a inertialess rod with length l. The angular

displacement of the rod is given by the angle θ,

and g is the acceleration due to gravity.

Using x and θ as the generalized coordinates the

potential energy for the system is

V =

+ m

gl cos θ.

The ﬁrst term is the potential energy due to the

spring force, and the second term is the potential

energy due to the force of gravity on m

lθ

At the equilibrium the potential energy function is stationary, hence we

require

∂V

∂x

= kx = 0,

∂V

∂θ

= −m

gl sin θ = 0.

These equations are satisﬁed at the points



∗





±nπ



, n = 0, 1, 2, ···.

We will examine the two unique equilibrium positions, i.e., (a) x

∗

= 0,

∗

= 0, and (b) x

∗

= 0, θ

∗

= π. In either case the Hessian of the potential

energy function is

H(x, θ) =







∂

∂x

∂

∂x∂θ

∂

∂θ∂x

∂

∂θ









k 0

0 −m

gl cos θ



Therefore, at the ﬁrst equilibrium point the Hessian is

H(0, 0) =



k 0

0 −m



Since k > 0, m

> 0, g > 0 and l > 0, the Hessian matrix, H(0, 0), is

indeﬁnite, i.e., x

∗

= 0, θ

∗

= 0 is a saddle point of the potential energy

function. By the extended Lyapunov theorem it can be concluded that this

equilibrium is unstable. Intuitively this can be shown to be true since, at

276 6 Dynamic System Analysis and Simulation

θ = 0 any small perturbation in the angle θ will cause the rod to tip over,

and thus move ‘far’ away from the equilibrium.

At the second equilibrium point the Hessian is

H(0, π) =



k 0

0 m



which is a positive deﬁnite matrix. In this case we can use the Lagrange

stability theorem to argue that the equilibrium point x

∗

= 0, θ

∗

= π is

stable.

6.2 System Simulation

One of the attractive features of the modeling technique described in this

book is that the process can be readily automated on computer systems, so as

to obtain the equations of motion and simulate the behavior of the model. In

this section we describe two computer programs that facilitate the numerical

solution of the Lagrangian DAEs. We also present several simulation examples

using these programs.

The program ldaetrans, which is described in section 6.2.1, takes a system

described by q, f, e

, e, s, T

∗

, V , D, φ, ψ, Γ , and Σ, and generates the LDAEs

(4.8) via symbolic diﬀerentiation of the appropriate terms. This program also

produces a ﬁle that can be used to simulate the equations of motion using

MATLAB or Octave.

The function ride, which is described in section 6.2.2, is an implementa-

tion of an implicit Runge-Kutta method for the solution of implicit diﬀeren-

tial equations (IDEs). The solution technique is based on the 3-stage Radau

IIA implicit Runge-Kutta method as presented in Chapter 5. The function

is written in the MATLAB/Octave programming language.

Both computer codes can be found on the web site:

http://abs-5.me.washington.edu/pub/fabien/asd

6.2.1 The translator ldaetrans

The program ldaetrans translates an ASCII input ﬁle description of the

model into the LDAEs (4.8), and it also produces ﬁles suitable for the numer-

ical integration of the model using MATLAB or Octave. (Here, we used the

function ride to integrate the LDAEs.) The program ldaetrans is written

6.2 System Simulation 277

in ANSI C, and has been compiled and executed on several UNIX based op-

erating systems including Linux, Solaris, Mac OSX, OpenBSD and FreeBSD.

The input to this translator takes advantage of the fact that the modeling

approach developed in this book leads to a set of highly structured equations.

The ASCII input ﬁle uses keywords to enter the system variables, scalar

energy functions, constraints, etc. A list of the keywords recognized by the

program are shown in Table 6.1. These keywords are used as follows.

DisplacementVariables FlowVariables EffortVariables

DynamicVariables KineticCoEnergy PotentialEnergy

DissipationFunction GeneralizedEfforts DisplacementConstraints

FlowConstraints EffortConstraints DynamicConstraints

InitialDisplacements InitialFlows InitialEfforts

InitialDynamic InitialTime FinalTime

OutputPoints

Table 6.1 ldaetrans keywords

• DisplacementVariables: This keyword is used to input the displace-

ment variables in the model. For example, the displacement variables

q = [q

, q

]

are entered using the statement

DisplacementVariables = [q1,q2,q3];

Note that, all vector inputs to the program are enclosed in square brackets,

and the elements of the vector are separated by commas. In addition, all

input lines to the program end with a semicolon.

Similar input statements are used to enter the ﬂow variables f, the eﬀort

variables e, and the dynamic variables s, via the keywords FlowVariables,

EffortVariables and DynamicVariables, respectively.

• KineticCoEnergy: This keyword is used to enter the kinetic coenergy for

the system. For example, T

∗

can be entered using the statement

KineticCoEnergy = (m*f1*f1)/2;

Similar statements are used to enter the potential energy V , and the dis-

sipation function D via the keywords PotentialEnergy and

DissipationFunction, respectively.

• GeneralizedEfforts: The eﬀorts e

= [e

, e

, ···, e

]

are entered

using this keyword. For example, consider a system with displacements

q = [q

, q

]

, and the virtual work due to the applied eﬀorts is δW =

τ δq

+ F δq

. In this case e

= [τ, 0, F ]

, and can be described using the

statement

GeneralizedEfforts = [tau,0,F];

Note that the eﬀorts e

are distinct from the eﬀorts associated with the

eﬀort constraints Γ , i.e., the eﬀort variables e. Also, in a strict sense these