To get familiar with the methodology of control systems simulations in Mathematica, methods for results visualization and analysis.

Scientific & Engineering Programming

II Year Electronics and Computer Engineering, FoE, WUST

Laboratory Class 7 – Control systems in Mathematica

The scope

To get familiar with the methodology of control systems simulations in Mathematica, methods for results visualization and analysis.

Prerequisites

Before the classes you should know, how to:

• model simple physical systems,

• represent and define differential equations,

• solve numerically differential equations,

• visualize functions being the differential equations solutions,

• visualize the results on animations.

Tasks

Read carefully the tasks preliminaries. Analyze the to dos. Choose the task more suitable for you.

1. Controlling the Inverted Pendulum with Moving Vertically Base

Preliminaries

Differential equations, describing the dynamical behavior of a depicted below inverted pen- dulum of the length l and the mass m, with a moving vertically base,

x y

θ y

l

x

m

y

has a form

l ¨ θ = g sin θ − c ˙ θ + sin θ ¨ y, (1) where θ denotes the angle, the pendulum link makes with the vertical, y is the base offset, c viscous damping coefficient, and g stands for the gravitational acceleration. The dot operator means a time derivative ˙  =

.

1

Scientific & Engineering Programming, LC7, II Year ECE, FoE, WUST 2

y

CONT PLANT y

u

Figure 1: Open loop control scheme

Let’s assume, the pendulum base is moving in a oscillatory way, described by the equation

y(t) = A sin(ωt), (2)

where A is the oscillations amplitude and ω their frequency. Introduction of the oscillatory base (2) into the system (1) allows to control its behavior exploiting the idea of an open loop (feedforward) control scheme, illustrated in Figure 1. In the figure PLANT stands for the object to be controlled, y is its output, and u its input. CONT is a controller, which on the base of desired system behavior tries to generate the control input u, causing the system to behave as desired. In the case of the analyzed inverted pendulum with a moving vertically base, the role of the input u plays the second derivative of the base position ¨ y, and the angle θ is taken as the output y. The oscillator takes the role of the controller. Changing its oscillations frequency and amplitude one can control the behavior of the pendulum and decide, which equilibrium the pendulum stabilizes at.

To dos

For the described system do the following:

(a) Prepare the Mathematica model of the inverted pendulum with a moving vertically base (1):

i. define a function which returns the pendulum base position, ii. define a function which returns the mass m position, iii. define the differential equation 1.

Analyze the model behavior for an initial θ value equal to 0

and 0.2, and the steady base (y=0). Plot the joint angle as a function of time. Animate the pendulum movement.

Investigate the influence of the dumping coefficient on the pendulum behavior.

(b) Equip the pendulum with the oscillator (2):

i. define the rule which introduces the controller (2), ii. apply the defined rule to the system.

Examine the correctness of the provided solution.

(c) Simulate the system behavior for both the initial values with different values of oscil- lations amplitude and frequency. Try

A = 0.2, ω = 4; A = 0.2, ω = 40; A = 0.2, ω = 400; A = 0.02, ω = 400; A = 0.002, ω = 400; A = 0.002, ω = 4000. Try to derive some conclusions about the influence of the oscillator amplitude and its frequency on the system behavior

.

∗An inverted pendulum with such a base is know as Kapitza’s pendulum.

†For the practical purposes you can observe also the pendulum behavior for the initial value very close to 0 (i.e. 10⁻⁴.)

‡With the pendulum length l = 1 and the dumping coefficient c = 0.2.

§Do not forget to change the animation speed if necessary, to observe the actual system behavior

Scientific & Engineering Programming, LC7, II Year ECE, FoE, WUST 3

2. Controlling the Double Pendulum

Preliminaries

Differential equations, being the dynamical model of a presented during the lecture and depicted below double pendulum,

x y

θ

y

l

x

x

y

l

θ

m

m

consisting of two masses m

and m

, attached by rigid massless links of lengths l

and l

, can be simplified to a form

( m

l

θ ¨

+ m

l

l

θ ¨

cos(θ

− θ

) + m

l

l

θ ˙

sin(θ

− θ

) + gm

l

sin θ

= u

m

l

θ ¨

+ m

l

l

θ ¨

cos(θ

− θ

) − m

l

l

θ ˙

sin(θ

− θ

) + gm

l

sin θ

= u

, (3)

where θ

, θ

denote the angles, the links make with the vertical, u

, u

describe the external, control torques applied at the pendulum joints, m

= m

+ m

, and g stands for the gravitational acceleration. The dot operator stands for a time derivative ˙  =

. The equations (3) can be rewritten in a general, vector form

Q(q)¨ q + B(q, ˙ q) = u, (4)

where the vectors q =

2

, u = (

), and the matrix Q(q) and the vector B(q, ˙ q) are expressed as follows:

Q(q) =

 m

l

m

l

l

cos(θ

− θ

) m

l

l

cos(θ

− θ

) m

l

 ,

B(q, ˙ q) =  m

l

l

θ ˙

sin(θ

− θ

) + gm

l

sin θ

−m

l

l

θ ˙

sin(θ

− θ

) + gm

l

sin θ

 .

To control the position of the double pendulum end, which is a function of the link angles θ

and θ

x

y



=

 l

sin θ

+ l

sin θ

−l

cos θ

− l

sin θ



and has a form of the vector output

y = k(q)

with y = (

), one can apply the idea of a closed loop feedback linearization, illustrated in Figure 2. In the figure PLANT stands for the object to be controlled, y is its output, and u its input. LIN is a linearizing block, which makes the compound system with the input v and the output y to be linear. PD is just a linear PD controller, providing a control loop feedback,

¶The system equations in this form, together with all the accompanying matrices and vectors are defined in the Mathematica notebook file, provided on the course web page in the laboratory classes table.

Scientific & Engineering Programming, LC7, II Year ECE, FoE, WUST 4

v y

LIN PLANT

u

PD _y

d

Figure 2: Linearizing closed loop control scheme

allowing the output y to follow a desired value y

. In the case of the double pendulum the LIN module has to realize the relation

u = −D

(q)P (q, ˙ q) + D

(q)v, (5) with D(q) = J (q)Q

(q), P (q, ˙ q) = ˙ J (q) ˙ q − J (q)Q

(q)B(q, ˙ q), J (q) =

. The PD controller follows the well known rule of control, described by

v = ¨ y

− K

( ˙ y − ˙ y

) − K

(y − y

), (6) with diagonal gain matrices K

, K

, containing nonnegative elements.

To dos

For the described system do the following:

(a) Compare the double pendulum vector equations (4) and their matrix Q and vector B with these introduced in the provided file – put the attention to the equations structure.

(b) Simulate the system behavior for the initial condition θ

(0) = P i/2, θ

(0) = ˙ θ

(0) = θ ˙

(0) = 0, and the controls u = (

) equal to zero. Plot the joint angles as a function of time. Animate the pendulum movement.

(c) Set the pendulum parameters l

= l

= m

= m

= 1. Try the controls (

) = (

), (

) = (

), (

) = (

), (

) = (

).

(d) Analyze the way, the linearizing module equations (5) are introduced in the provided file.

(e) Simulate the linearized system behavior (with the input v) for the initial condition θ

(0) = P i/2, θ

(0) = ˙ θ

(0) = ˙ θ

(0) = 0, and the controls

(

) = (

), (

) =

, (

) = (

), (

) = (

), (

) =

, (

) = (

). Try to interpret the meaning of the controls v

, v

.

(f) Repeat the simulation for the last value of controls v with the nonzero initial joint speed θ ˙

(0) = −0.1. With this initial speed try the control (

) =

.

(g) Analyze the way, the control loop (6) is defined in the provided file. Apply the defined control rules to the linearized system. Analyze the control system behavior for an exemplary desired trajectory y

(t) =

0.5+0.3 sin t

(the parametric equation of a circle).

Add to the animation the desired and actual trajectories (y

and y.)

(h) Analyze the control system behavior for different values of the gains, different initial conditions. Change the desired trajectory. What happens when the desired trajectory (part of it) is located outside the region, which the pendulum end can reach?

kThe formula together with the matrices is included in the Mathematica notebook file, provided on the course web page in the laboratory classes table.

∗∗Also included in the Mathematica notebook file, provided on the course web page in the laboratory classes table.

††In the cases when the NDSolve function doesn’t give the solution for the whole time domain you’ve specified, visualize just the part which has been obtained. Try to explain, why no the whole solution was obtained.