Lyapunov Exponent as Physics-Informed Dense Reward: RL Discovery of Stabilization Beyond the Kapitza Pendulum Slava Andrejev We suggest using the Lyapunov characteristic exponent (LCE) as a dense reward signal for the reinforcement learning problem of stabilizing the inverted pendulum with vertical motion. With LCE, the agent not only successfully found the oscillatory motion known as the Kapitza pendulum but also damped the pendulum’s pivoting, leaving it in a strictly upright position.
1 Introduction
Whenever we write LCE, we will assume the maximal LCE. The LCE is an important parameter of a dynamiOne of the classic problems for testing reinforcement cal system. It measures how fast infinitesimally close learning methods is the inverted pendulum problem. trajectories in the phase space diverge. For the purpose In the textbook formulation, a mass on a light, rigid of stabilization, a negative LCE indicates the system is stick is attached to a cart on a rail, with the mass point- stable. If the LCE is positive, infinitesimally small perturing upward. An RL agent must find a strategy for cart bations of the initial state grow and diverge chaotically. movement to keep the pendulum upright. For the inverted pendulum, the LCE is positive. As a Besides this straightforward way of balancing the consequence, we cannot balance the inverted pendupendulum with a horizontal motion, we can also set it lum; we have to apply some form of control to keep the upright with a vertical motion. If we oscillate the pivot system upright. point up and down fast enough, the pendulum will have If we compute an average LCE over a window, it will an equilibrium point at the top. This type of pendulum, give the agent a direct signal about whether the system with an oscillating pivot, is called the Kapitza pendulum. is stable over time. The squared angle with vertical does To restore justice, Pyotr Kapitza didn’t discover it; he not carry this integral information. If the pendulum falls found a brilliant way to describe its motion [1]. The fact in the middle of an episode, we don’t know how well we that this pendulum has a stable equilibrium upright has performed before the fall. We just get a large negative been known before [2]. reward with no indication that oscillatory movement A fair question to ask is: if we apply the reinforcement would have led to a stable system. learning methods we use for the inverted pendulum on a This work will explore whether an agent, armed with cart to a pendulum with the pivot moving up and down, LCE as a reward signal, can discover inverted pendulum will the AI agent find the Kapitza pendulum? In one stabilization with vertical movement of the pivot point. attempt [3], the authors used the Deep Deterministic Policy Gradient (DDPG) method [4]. The agent failed to find oscillatory movement. It either let the pendulum 2 Motion Equations fall most of the time at the end of an episode, or started accelerating downward without bound. The authors Figure 1 shows the physical model of our pendulum. It hypothesized that the reward sparsity was the reason is not difficult to show that the Lagrangian of the system for the failure. They used the squared angle with vertical is as a reward signal. They didn’t verify their hypothesis. 𝑚 To improve their result, we suggest using the maxi𝐿 = (𝑥𝑐̇ 2 + 𝑦𝑐̇ 2 + 𝑙 2 𝜑̇2 + 2𝑙 𝜑̇ (𝑥𝑐̇ cos 𝜑 + 𝑦𝑐̇ sin 𝜑)) − 2 mal Lyapunov Characteristic Exponent (LCE) [5] as a − 𝑚𝑔𝑙 cos 𝜑, reward signal. From here on, we will omit “maximal”.
1
Lyapunov Exponent as Dense Reward 𝑥
has at most 𝑛 = dim 𝑇𝑥 𝑀 different values. These values are 𝜈1 > 𝜈2 > … > 𝜈𝑠 , with 1 ⩽ 𝑠 ⩽ 𝑛. All 𝑔 that produce 𝜈𝑠 in (3) form a linear space 𝐿𝑠 . Those 𝑔 that produce 𝜈𝑠 and 𝜈𝑠−1 form a linear space 𝐿𝑠−1 ⊃ 𝐿𝑠 . We can continue this process until we reach 𝐿1 . Therefore, we have
𝑚 𝑙
𝜑 𝑦
𝐿𝑠 ⊂ 𝐿𝑠−1 ⊂ … ⊂ 𝐿1 = 𝑇𝑥 𝑀,
Figure 1: Pendulum model.
and for a vector 𝜈 ∈ 𝐿𝑖 ∖ 𝐿𝑖+1 where 𝑥𝑐 and 𝑦𝑐 are the coordinates of the pivot point. Using the Euler–Lagrange equation for 𝜑, we will get 𝜑̈ +
𝑦𝑐̈ − 𝑔 𝑥𝑐̈ cos 𝜑 + sin 𝜑 = 0. 𝑙 𝑙
Since we have only vertical motion, 𝑥𝑐̈ = 0, and 𝑦𝑐̈ will be our action (RL, not physical action), which we will denote as 𝑎. Finally, 𝜑̈ +
𝑎−𝑔 sin 𝜑 = 0. 𝑙
(1)
For simplicity, we will measure acceleration in units of 𝑔 and distance in units of 𝑙. With this choice, the natural time scale becomes √𝑙/𝑔 seconds. (1) will take the following form 𝜑 ̈ = (1 − 𝑎) sin 𝜑,
𝜈𝑖 = lim
1
𝑡→∞ 𝑡
𝛼𝑘 = ‖𝑑Φ𝑠𝑥(𝑘−1)𝑠 𝑤𝑘−1 ‖ ,
(2)
𝑑Φ𝑠𝑥(𝑘−1)𝑠 𝑤𝑘−1 𝛼𝑘 𝑘
We will follow Benettin et al. [6, 7] in briefly describing the LCE theory and the numerical method for computing the LCE. We have a compact connected Riemannian manifold 𝑀 and a diffeomorphism Φ of 𝑀 onto itself. At each point 𝑥 ∈ 𝑀, 𝑇𝑥 𝑀 is the tangent space of 𝑀 at point 𝑥. A flow {Φ𝑡 } is a one-parameter family of diffeomorphisms satisfying Φ𝑠+𝑡 = Φ𝑠 ∘ Φ𝑡 , and Φ0 = 𝟏. 𝑑Φ𝑡𝑥 ∶ 𝑇𝑥 𝑀 → 𝑇Φ𝑡 𝑥 𝑀 is a linear map telling us how infinitesimally nearby points get stretched and rotated when we apply Φ. Finally, {𝑑Φ𝑡𝑥 } is a family of linear maps. For any nonzero vector 𝑔 ∈ 𝑇𝑥 𝑀 one can define (3)
1 ∑ ln 𝛼𝑖 . 𝑘→∞ 𝑘𝑠 𝑖=1
𝜈1 = lim
(5)
An ODE for 𝑑Φ𝑡𝑥 can be derived using the Jacobian matrix. If we have a dynamic system 𝑥 ̇ = 𝑓 (𝑥), then Φ𝑡𝑥 is a solution that maps the initial condition 𝑥(0) = 𝑥 to 𝑥(𝑡). We can write this fact down as 𝑑 𝑡 Φ = 𝑓 (Φ𝑡𝑥 ) . 𝑑𝑡 𝑥
(6)
𝑑Φ𝑡𝑥 by definition is the spatial derivative of Φ𝑡𝑥 with respect to the initial condition. If we apply the spatial derivative to (6) and use the chain rule, we will get an ODE for 𝑑Φ𝑡𝑥 𝑑 [𝑑Φ𝑡𝑥 ] = 𝐷𝑓 (𝑥) ⋅ 𝑑Φ𝑡𝑥 , 𝑑𝑡
This number is called the LCE of vector 𝑔 with respect to the family 𝑑Φ𝑡𝑥 . It can be shown [6] that {𝜒 (𝑔)}𝑔∈𝑇𝑥 𝑀 where 𝐷𝑓 is the Jacobian matrix.
2
.
Finally,
3 The Maximal LCE
𝑡→∞
(4)
𝑤0 = 𝜈 ∈ 𝑇𝑥 𝑀,
where 𝑎(𝑡) = 𝑦𝑐̈ (𝑡).
1 ln ‖𝑑Φ𝑡𝑥 𝑔‖. 𝑡
1 ⩽ 𝑖 ⩽ 𝑠.
If we randomly choose vector 𝑔 ∈ 𝑇𝑥 𝑀 and plug it to (4), we will obtain the maximal LCE 𝜈1 because 𝑔 has zero probability to end up in 𝐿2 , ..., 𝐿𝑠 . This is the basis of computing the maximal LCE numerically as it is shown in [7]. For numeric computation, we have to take care of the exponential growth of the evolving tangent vector 𝑑Φ𝑡𝑥 𝜈. To avoid floating point overflow, Benettin et al. [7] suggest periodically renormalizing the tangent vector after some fixed time 𝑠. That is, we choose a random 𝜈 ∈ 𝑇𝑥 𝑀 with ‖𝜈‖ = 1 and compute recursively
𝑤𝑘 =
𝜒 (𝑔) = lim sup
ln ‖𝑑Φ𝑡𝑥 𝜈‖,
Slava Andrejev For the case of the inverted pendulum, 𝑀 consists of and 𝑦𝑐 (𝑡) are fast-oscillating periodic functions of period ̇ From (1), the Jacobian matrix is pairs 𝑥 = (𝜑, 𝜑). 𝑇 . Following [8], we further require that the changes in 𝜑 are small, and we separate 𝜑 into slow and fast parts, 0 1 𝐷𝑓 = [ 1 ]. (𝑔 − 𝑎 (𝑡)) cos 𝜑 0 𝜑 = 𝜑 𝑠 + 𝜑𝑓 . (9) 𝑙 ̇ Using the Angle 𝜑𝑓 is a periodic function and its value over one Similarly, 𝑇𝑥 𝑀 is a linear space of pairs (𝛿, 𝛿). period of oscillation remains small. Angle 𝜑𝑠 can have Jacobian above, we get any value but during the same time 𝑇 it doesn’t change 𝑔 − 𝑎 (𝑡) ̈ 𝛿= 𝛿 cos 𝜑. (7) much. If we average these values over the period 𝑇 and 𝑙 mark this average with a bar, we get With the units we introduced earlier, this equation be𝜑 ≅ 𝜑𝑠 , 𝜑 𝑓 ≅ 0. (10) comes 𝛿 ̈ = (1 − 𝑎) 𝛿 cos 𝜑. (8) If we substitute (9) into (2), drop the slow 𝜑 ̈ , and expand 𝑠 Instead of using a random initial vector 𝑤0 = (𝛿0 , 𝛿0̇ ) sin 𝜑, we get for an LCE estimate in (5), we can directly derive an 𝜑𝑓̈ = (1 − 𝑎) (sin 𝜑𝑠 + 𝜑𝑓 cos 𝜑𝑠 ) . approximation of the direction of the maximum ‖𝑤0 ‖ change. We set our pendulum at some small initial angle Since 𝜑𝑓 is small, we can leave out 𝜑𝑓 cos 𝜑𝑠 , and the term 𝜑0 and let it go. The solution of (2) for small angles is 1 ⋅ sin 𝜑𝑠 is a slow drift, which belongs to 𝜑𝑠 . Therefore, 𝜑 ≈ 𝜑0 cosh(𝑡) ≈ 𝜑0 . Since ‖𝑤0 ‖ = 1, we can write a solution of (8) as sin 𝛼0 𝛿 (𝑡) = cos 𝛼0 cosh(𝑘𝑡) + sinh(𝑘𝑡), 𝑘 where 𝑘 = √cos 𝜑0 , and 𝛼0 is the initial direction of 𝑤0 . We need to find 𝛼0 ̇ that maximizes ‖(𝛿, 𝛿)‖. 4 𝑑 ̇ 2 = 1 − 𝑘 sin(2𝛼) sinh2 (𝑘𝑡)+ ‖(𝛿, 𝛿)‖ 𝑑𝛼0 𝑘2 1 + 𝑘2 sinh(2𝑘𝑡) cos(2𝛼). 𝑘
From which it follows that at 𝛼0 =
2𝑘 cosh(𝑘𝑡) 𝜋 1 − atan( ) 2 2 (1 − 𝑘 2 ) sinh(𝑘𝑡)
𝜑𝑓̈ = −𝑎 sin 𝜑𝑠 , and 𝜑𝑓 = −𝜂 sin 𝜑𝑠 ,
(11)
where 𝜂(𝑡) is the zero-mean shift of 𝑦𝑐 , 𝜂 = 𝑦𝑐 − 𝑦𝑐 , due to the requirement 𝜑 𝑓 = 0 from (10). The next step is to use (11) in (2) and (9), 𝜑𝑠̈ + 𝜑𝑓̈ = (1 − 𝑎) (sin 𝜑𝑠 + 𝜑𝑓 cos 𝜑𝑠 ) , and average the result over one period 𝑇 , 𝜑𝑠̈ = sin 𝜑𝑠 + 𝑎𝜂 cos 𝜑𝑠 sin 𝜑𝑠 . After integration by parts, 𝜑𝑠̈ = sin 𝜑𝑠 − 𝜂2̇ cos 𝜑𝑠 sin 𝜑𝑠 . This gives a necessary condition for stabilizing the pendulum at a nonzero constant angle 𝜑𝑠 ,
we will observe the maximum increase of ‖𝑤0 ‖. 𝑡 in the 𝜂2̇ = cos−1 𝜑𝑠 . (12) above will be the time period over which we estimate the LCE. We can follow a similar technique for 𝛿 in (8). We separate it into slow and fast parts,
4 LCE Solution For a Pendulum With an Oscillating Pivot Let’s verify that a solution for (2) and (8) with a constant angle 𝜑 and negative LCE exists. We assume that 𝑎(𝑡)
𝛿 = 𝛿 𝑠 + 𝛿𝑓 . Finding 𝛿𝑓 is very similar to finding 𝜑𝑓 , 𝛿𝑓 = −𝜂𝛿𝑠 cos 𝜑𝑠 .
3
Lyapunov Exponent as Dense Reward The slow part is more convoluted, 𝛿𝑠̈ = 𝛿𝑠 cos 𝜑𝑠 + 𝛿𝑓 cos 𝜑𝑠 − 𝑎𝛿𝑓 cos 𝜑𝑠 − 𝜑𝑓 𝛿𝑠 sin 𝜑𝑠 − 𝜑𝑓 𝛿𝑓 sin 𝜑𝑠 + 𝑎𝛿𝑠 𝜑𝑓 sin 𝜑𝑠 + 𝑎𝛿𝑓 𝜑𝑓 sin 𝜑𝑠 .
angle, and the pivot motion won’t be a single cycle periodic motion. In the stationary mode, fast oscillation will be accompanied by slow periodic movement (periodic because of the time symmetry), and in the transition phase, the agent may try to damp the oscillations of 𝜑𝑠 by tracking 𝛿𝑠̇ .
Dropping the terms that vanish after averaging over 𝑇 ,
5 Learning Method
𝛿𝑠̈ = 𝛿𝑠 cos 𝜑𝑠 − 𝑎𝛿𝑓 cos 𝜑𝑠 − 𝜑𝑓 𝛿𝑓 sin 𝜑𝑠 + 𝑎𝛿𝑠 𝜑𝑓 sin 𝜑𝑠 + 𝑎𝛿𝑓 𝜑𝑓 sin 𝜑𝑠 . Finally, let’s remove the higher order terms and average, 𝛿𝑠̈ = 𝛿𝑠 (cos 𝜑𝑠 − 𝜂2̇ cos 2𝜑𝑠 + 𝜂𝜂̈ 2 sin2 𝜑𝑠 cos 𝜑𝑠 ) . If we assume that 𝜑𝑠 = const and take into account (12), 𝛿𝑠̈ = 𝛿𝑠
sin2 𝜑𝑠 (1 + 𝜂𝜂̈ 2 cos2 𝜑𝑠 ) . cos 𝜑𝑠
(13)
A generic solution is 𝛿𝑠 = 𝐴𝑒 −𝑘𝑡 + 𝐵𝑒 𝑘𝑡 , where 𝑘 = 𝜈1 =
1/2 |sin 𝜑𝑠 | (1 + 𝜂𝜂̈ 2 cos2 𝜑𝑠 ) √cos 𝜑𝑠
(14)
is the maximal LCE. We’ve made several important discoveries. First, from (14), it follows that for 𝜑𝑠 = 0, the minimum LCE is zero. This means that truly upright stabilization is impossible. The upright equilibrium cannot damp oscillations. It makes intuitive sense: we cannot create any torque by vertical movement if 𝜑 = 0. We can also see it in (11); if 𝜑𝑠 = 0, then the fast term 𝜑𝑓 vanishes. Second, we can see from (13) that the only way to stabilize the pendulum is to have 𝜂𝜂̈ 2 ≠ 0. This can only happen if the pivot motion is asymmetric, that is, the upward motion has a different acceleration profile than the downward motion. Finally, with strictly periodic motion, when 𝜂𝜂̈ 2 = const, true stabilization is also impossible. To achieve stabilization and have negative LCE, we must have a damping term in (13), 𝜂𝜂̈ 2 ∼ −𝛿𝑠̇ . Therefore, we can expect that if we use the LCE as a reward signal, the agent will essentially go beyond the Kapitza pendulum. To make the LCE negative, it probably will stabilize the pendulum at some small, nonzero
4
Our environment consisted of the inverted pendulum described by (1). The length of the pendulum rod was 30 cm. ODE (1) was combined with (7) into one system of equations to solve them simultaneously. Note that the immediate output of (7) at time 𝑡 is useless; we have to see how it evolves through time and use (5) to derive the LCE. To do that, we calculated the maximal LCE according to (5) over a short time window of 0.25 s. For training, we used Soft Action-Critic (SAC) algorithm [9] from Stable Baselines 3 [10]. We ran 1 s episodes, and each sample in an episode increased time by 5 ms. Action in each sample was acceleration 𝑎 from (1). Since the reward in our environment was calculated from the LCE over a time window, the state had to include all values of the pendulum speed and location over the LCE window. That is, our window had 50 time steps, and at each step we recorded 𝑦𝑐 , 𝑦𝑐̇ , sin 𝜑, cos 𝜑, and 𝜑.̇ Therefore, the state was a set of 250 values: {(𝑦𝑐,𝑖 , 𝑦𝑐,𝑖 ̇ , sin 𝜑𝑖 , cos 𝜑𝑖 , 𝜑𝑖̇ ), 𝑖 = 1 … 50}, where 𝑖 is the index of a time step inside the window. To reset the environment, we set the pendulum at a small random angle within ±5.7° and let it fall for the duration of the LCE time window, filling in all 𝜑 and 𝜑̇ in the initial state. The SAC implementation was used with default parameters, except the learning rate 𝛼 = 1 × 10−4 . The default parameters were: discount 𝛾 = 0.99, soft update rate for target networks 𝜏 = 5 × 10−3 , replay buffer size 1 × 106 , batch size 256, train frequency was one gradient step per environment step, target entropy −1, and 100 random steps before training began. The implementation used fully connected neural networks with two hidden layers, 256 neurons each, for all function approximations: actor, critics, and target critics. We couldn’t use pure LCE as a reward signal. The agent would have found a trivial solution: spatially unbounded acceleration down with 𝑎 > 𝑔. To encourage the agent to stay in the center, we added a smooth negative reward that was almost constant when 𝑦𝑐 stayed
Slava Andrejev
Reward
0 −500 −1,000 −1,500
0
200
400
600
800
1,000
Evaluation number
and 𝑐 = 25 was selected to keep the reward roughly the same order of magnitude as the LCE. The result is shown in fig. 4. We can see that the agent managed to keep the pendulum upright, but the pendulum moves chaotically around the equilibrium with a far larger amplitude than the pendulum driven by an agent that was trained on the LCE signal.
7 Conclusion
Figure 2: Training progress.
We suggested a physics-informed reward signal based on the LCE in an attempt to let the agent discover the within ±𝑦max and rapidly increased when the agent went Kapitza pendulum. The trained model outperformed outside these bounds: it. Not only did it find an equilibrium point at the top, it also managed to damp the oscillations and bring the 20 (𝑦𝑐 /𝑦max ) 𝑟𝑡 = −𝜈1 − 100 . (15) pendulum to almost a halt with residual oscillations of 1 + (𝑦𝑐 /𝑦max )20 ±0.01°. The remaining question is why the stable mode keeps If the agent went outside the ±𝑦max bounds, we terminated the episode with 𝑟𝑡 = −1000. We had 𝑦max = the pendulum at 0.25° from vertical. Is it a fundamental constant that minimizes the LCE? Can we make it zero? 10 cm in our experiments. Can we stabilize the pendulum at an arbitrary angle from vertical?
6 Results
Figure 2 shows the training progress. The quick rise indicates that the agent rapidly learned how to stay within the ±𝑦max boundary. Then it slowly learned how to stabilize the pendulum. The work of the trained agent is shown in fig. 3a. The agent’s performance is remarkable. Not only did it find a mode that provides an equilibrium point at the top, it damped the angular oscillations, significantly improving upon the Kapitza pendulum. For comparison, the same picture shows how the Kapitza pendulum behaves when we oscillate the pivot point with the same frequency as the trained model. The Kapitza pendulum started from 1° and continued oscillating. The agent started from 10° with the pendulum falling. It stopped the fall at 15° and then brought the pendulum to oscillations as small as ±0.01° (see fig. 3d). Note that training episodes were 2 s. In this time span, the pendulum has not yet stabilized. At the end of the episode the agent was still damping the oscillations. However, the learned policy was general enough to let the agent keep the pendulum upright outside the training window. Since in the previous work [3] the authors used the DDPG learning algorithm, we verified whether SAC would stabilize the inverted pendulum with the same reward signal as they used. We changed 𝜈1 in (15) to 𝑐𝜑 2 ,
References [1]
Пётр Леонидович Капица. «Маятник с вибрирующим подвесом». В: Успехи физических наук 44 (1951), с. 7—20.
[2]
Harold Jeffreys and Bertha Swirles Jeffreys. Methods of Mathematical Physics. 2nd ed. Cambridge University Press, 1950. Chap. Solution of Linear Differential Equations of the Second Order, pp. 488–489.
[3]
Slava Andrejev and Nate Ruble. “RL Discovery of Kapitza Pendulum.” Course project, CS-541: Artificial Intelligence, Portland State University. 2023.
[4]
Timothy P. Lillicrap et al. “Continuous Control with Deep Reinforcement Learning.” In: Computing Research Repository abs/1509.02971 (2015). doi: 10 . 48550/arXiv.1509.02971. url: https://arxiv. org/abs/1509.02971.
[5]
A. M. Lyapunov. “The general problem of the stability of motion.” In: International Journal of Control 55.3 (1992), pp. 531–534. doi: 10 . 1080 / 00207179208934253.
[6]
Giancarlo Benettin et al. “Lyapunov Characteristic Exponents for smooth dynamical systems and for Hamiltonian systems; a method for computing all of them. Part 1: Theory.” In: Meccanica 15.1 (1980), pp. 9–20. issn: 1572-9648. doi: 10.1007/BF02128236.
5
Lyapunov Exponent as Dense Reward ⋅10−2 5
LCE pendulum Kapitza pendulum
𝑦𝑐 , m
10 𝜑, °
5 0
0 −5
−5 0
2
4 6 Time, s
8
10
0
0.5
1 Time, s
2
1.5
(a) Vertical angle for the LCE pendulum and the Kapitza (b) The agent is damping oscillations at the beginning of an pendulum for comparison. The initial position of the latter episode. was 1°.
⋅10−2
1
0.27 𝜑, °
𝑦𝑐 , m
0 −1 −2
0.26 0.25 0.24
9
9.2
9.4 9.6 Time, s
9.8
10
9
(c) Pivot point vertical motion in stable mode.
9.2
9.4 9.6 Time, s
9.8
10
(d) Residual 𝜑 oscillations in stable mode.
Figure 3: Inverted pendulum stabilization with a trained agent. 𝑟𝑡 = 𝜈1 𝑟𝑡 = 𝜑 2
𝜑, °
10
[9]
Tuomas Haarnoja et al. “Soft Actor-Critic: Off-Policy Maximum Entropy Deep Reinforcement Learning with a Stochastic Actor.” In: Proceedings of the 35th International Conference on Machine Learning. 2018, pp. 1861–1870.
[10]
Antonin Raffin et al. “Stable-Baselines3: Reliable Reinforcement Learning Implementations.” In: Journal of Machine Learning Research 22.268 (2021), pp. 1–8. url: https://jmlr.org/papers/v22/20-1364.html.
5 0 −5 0
2
4 6 Time, s
8
10
Figure 4: Comparison of models trained with the LCE and 𝜑 2 reward signals. [7]
Giancarlo Benettin et al. “Lyapunov Characteristic Exponents for smooth dynamical systems and for Hamiltonian systems; A method for computing all of them. Part 2: Numerical application.” In: Meccanica 15.1 (1980), pp. 21–30. issn: 1572-9648. doi: 10.1007/ BF02128237.
[8]
L. D. Landau and E. M. Lifshitz. “§30. Motion in a rapidly oscillating field.” In: Mechanics. 3rd ed. Vol. 1. Course of Theoretical Physics. Elsevier, 1976, pp. 93– 95. isbn: 978-0-7506-2896-9.
6