Week 1

Objective

The objectives of this laboratory experiment are as follows:

  1. Obtain the linear state-space representation of the rotary pendulum plant.

  2. Design a state-feedback control system that balances the pendulum in its upright position using Pole Placement.

  3. Simulate the closed-loop system to ensure the given specifications are met.

Equipment

Week 1

A. Modeling

This experiment involves modeling of the rotary inverted pendulum.

Rotary Inverted Pendulum Model

Figure 1 Rotary servo base with Rotary Pendulum module

Table 1 Rotary Pendulum components (Figure 1)

ID#
Component

1

Rotary Servo

2

Thumbscrews

3

Rotary Arm

4

Shaft Housing

5

Shaft

6

Pendulum T-Fitting

7

Pendulum Link

8

Pendulum Encoder Connector

9

Pendulum Encoder

Table 2 Main Parameters associated with the Rotary Pendulum module

Symbol
Description
Value
Unit

Mass of pendulum

0.127

kg\mathrm{kg}

Total length of pendulum

0.337

m\rm{m}

Distance from pivot to center of mass

0.156

m\rm{m}

Pendulum moment of intertia about center of mass

0.0012

kgm2\rm{kg \cdot m^2}

Pendulum viscous damping coeffi­cient as seen at the pivot axis

0.0024

Nms/rad\rm{N\cdot m\cdot s/rad}

Mass of rotary arm with two thumb screws

0.257

kg\rm{kg}

Rotary arm length from pivot to tip

0.216

m\rm{m}

Rotary arm length from pivot to center of mass

0.0619

m\rm{m}

Rotary arm moment of inertia about its center of mass

9.98 x 10^-4

kgm2\rm{kg \cdot m^2}

Rotary arm viscous damping coeffi­cient as seen at the pivot axis

0.0024

Nms/rad\rm{N\cdot m\cdot s/rad}

Rotary arm moment of inertia about pivot

0.0020

kgm2\rm{kg\cdot m^2}

Pendulum encoder resolution

4096

counts/rev\rm{counts/rev}

Model Convention

The rotary inverted pendulum model is shown in Figure 2. The rotary arm pivot is attached to the Rotary Servo system and is actuated. The arm has a length of , a moment of inertia of , and its angle, θ\bm\theta, increases positively when it rotates counterclockwise (CCW). The servo (and thus the arm) should turn in the CCW direction when the control voltage is positive, i.e. > 0.

The pendulum link is connected to the end of the rotary arm. It has a total length of LpL_p and it center of mass is . The moment of inertia about its center of mass is . The inverted pendulum angle, α\bm\alpha, is zero when it is perfectly upright in the vertical position and increases positively when rotated CCW.

Figure 2 Rotary inverted pendulum conventions

Nonlinear Equations of Motion

Instead of using classical (Newtonian) mechanics, the Lagrange method is used to find the equations of motion of the system. This systematic method is often used for more complicated systems such as robot manipulators with multiple joints.

Specifically, the equations that describe the motions of the rotary arm and the pendulum with respect to the servo motor voltage, i.e., the dynamics, will be obtained using the Euler-Lagrange equation:

The variables qi\bm q_i are called generalized coordinates. For this system let

where, as shown in Figure 2, θ(t)\theta(t) is the rotary arm angle and α(t)\alpha(t) is the inverted pendulum angle. The corresponding velocities are

With the generalized coordinates defined, the Euler-Lagrange equations for the rotary pendulum system are

The Lagrangian of the system is described by

where TT is the total kinetic energy of the system and VV is the total potential energy of the system. Thus the Lagrangian is the difference between a system’s kinetic and potential energies.

The generalized forces QiQ_i are used to describe the nonconservative forces (e.g. friction) applied to a system with respect to the generalized coordinates. In this case, the generalized force acting on the rotary arm is

and acting on the pendulum is

​See Table A in the Appendix for a description of the corresponding Rotary Servo parameters (e.g., such as the back-emf constant, ). Our control variable is the input servo motor voltage, . Opposing the applied torque is the viscous friction torque, or viscous damping, corresponding to the term . Since the pendulum is not actuated, the only force acting on the link is the damping. The viscous damping coefficient of the pendulum is denoted by .

Once the kinetic and potential energy are obtained and the Lagrangian is found, then the task is to compute various derivatives to get the EOMs. After going through this process, the nonlinear equations of motion for the Rotary Pendulum are:

The torque applied at the base of the rotary arm (i.e. at the load gear) is generated by the servo motor as described by the equation 7. Refer to Table A in the Appendix for the Rotary Servo parameters.

Linearization

Linearization of a nonlinear function about a selected point is obtained by retaining upto first order term in the Taylor Series expansion of the function about the selected point. For example, linearization of a two variable nonlinear function f(z) where

z=[z1z2]Tz = \begin {bmatrix} z_1 & z_2\end {bmatrix}^T

about the point

z0=[ab]Tz_0 = \begin {bmatrix} a & b\end {bmatrix}^T

can be written as

flin(z)=f(z0)+f(z)z1z=z0(z1a)+f(z)z2z=z0(z2b)f_{lin}(z) = f(z_0) + \frac{\partial f(z)}{\partial z_1}\Bigr|_{z = z_0} (z_1 - a)+\frac{\partial f(z)}{\partial z_2}\Bigr|_{z=z_0}(z_2 - b)

Linearization of Inverted Pendulum Equations

The nonlinear equations of the inverted pendulum system obtained as Equations (7) and (8) are linearized about the equilibrium point with the pendulum in the upright position, i.e.,

θo=0\theta_o = 0, αo=0\alpha_o = 0, θ˙o=0\dot{\theta}_o = 0, α˙o=0\dot{\alpha}_o = 0, θ¨o=0\ddot{\theta}_o = 0, α¨o=0\ddot{\alpha}_o = 0

Linearization of Equation (7) about the above equilibrium state gives

(Jr+mpLr2)θ¨12mpLpLrα¨=kVmbθ˙Brθ˙(10)(J_r + m_pL_r^2)\ddot{\theta} - \frac{1}{2}m_pL_pL_r\ddot{\alpha} = kV_m-b\dot{\theta} - B_r\dot{\theta} \quad\quad (10)

where, from Equation(9) the servo motor torque coefficient kk is

k=ηgKgηmktRmk = \frac{\eta_gK_g\eta_mk_t}{R_m}

and the back emf coefficient bb is

b=ηgKg2ηmktkmRmb = \frac{\eta_gK_g^2\eta_mk_tk_m}{R_m}

Likewise, linearization of Equation (8) about the equilibrium point with the pendulum in the upright position gives

12mpLpLrθ¨+(Jp+14mpLp2)α¨12mpLpgα=Bα˙(11)-\frac{1}{2}m_pL_pL_r\ddot{\theta} + (J_p + \frac{1}{4}m_pL_p^2)\ddot{\alpha} - \frac{1}{2}m_pL_pg\alpha= -B\dot{\alpha} \quad \quad (11)

where gg is the acceleration due to gravity. Note that the negative sign of the gravity torque (α\alpha term) in Equation 11 indicates negative stiffness with the pendulum in the vertically upright position.

Equation (10) and (11) can be arranged in the matrix form as

[M][θ¨α¨]+[D][θ˙α˙]+[K][θα]=[k0]Vm(12)[M]\begin{bmatrix}\ddot{\theta} \\ \ddot{\alpha} \end{bmatrix} +[D] \begin{bmatrix}\dot{\theta} \\ \dot{\alpha} \end{bmatrix} + [K] \begin{bmatrix}\theta \\ \alpha \end{bmatrix} = \begin{bmatrix}k \\ 0 \end{bmatrix}V_m \quad \quad (12)

Where mass matrix MM, damping matrix DD, and the stiffness matrix KK become

[M]=[(Jr+mpLr2)12mpLpLr12mpLpLr(Jp+14mpLp2)][M] = \begin{bmatrix} (J_r + m_pL_r^2) & -\frac{1}{2}m_pL_pL_r \\ -\frac{1}{2}m_pL_pL_r & (J_p + \frac{1}{4}m_pL_p^2) \end{bmatrix}
[D]=[b+Br00Bp][D] = \begin{bmatrix} b+B_r & 0 \\ 0 & B_p \end{bmatrix}
[K]=[00012mpLpg][K] = \begin{bmatrix} 0 & 0 \\ 0 & -\frac{1}{2}m_pL_pg \end{bmatrix}

Defining the state vector xx, output vector yy and the control input uu as

x=[θαθ˙α˙]T,y=[θα]T,u=Vmx = \begin{bmatrix} \theta & \alpha & \dot{\theta} & \dot{\alpha}\end{bmatrix}^T , \quad y = \begin{bmatrix} \theta & \alpha \end{bmatrix}^T, \quad u = V_m

Equation (12) can be rewritten in state space form as

x˙=Ax+Bu(13)\dot{x} = Ax + Bu \quad \quad (13)
y=Cx+Du(14)y = Cx + Du \quad \quad (14)

where the system state matrix AA, control matrix BB, output matrix CCand the feedthrough of input to the output matrix DD become

A=[O2×2I2×2[M]1[K][M]1[D]]A = \begin{bmatrix} O_{2\times2} & I_{2\times2} \\ -[M]^{-1}[K] & -[M]^{-1}[D] \end{bmatrix}
B=[O2×1M1[k0]]B = \begin{bmatrix} O_{2\times1} \\ M^{-1} \begin{bmatrix} k \\0 \end{bmatrix} \end{bmatrix}
C=[10000100]C = \begin{bmatrix} 1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \end{bmatrix}
D=[00]D = \begin{bmatrix} 0 \\ 0 \end{bmatrix}

In the equations above, note that O2×2O_{2\times2} is a 2×22 \times 2 matrix of zeros, O2×1O_{2\times1} is a 2×12 \times 1 matrix of zeros, and I2×2I_{2 \times 2} is a 2×22\times2 identity matrix. Using the parameters of the system listed in Table 2 and the appendix Table A, the linear model matricesAA and BB can be computed.

Analysis: Modeling

  1. Download the above zip file and extract it.

  2. Open rotpen_week1_student.mlx live script and run the Modelling section. It will automatically load the parameters required for the state-space representation, and subsequently generate the A, B, C and D matrices required for the upcoming analysis. Please refer to Table 2 and Table A for additional information regarding the parameters. Note: The representative C and D matrices have already been included. The actuator dynamics have been added to convert your state-space matrices to be in terms of voltage. Recall that the input of the state-space model is the torque acting at the servo load gear (i.e. the pivot of the pendulum). However, we control the servo input voltage instead of control torque directly. The script uses the voltage torque relationship given in Equation 9 to transform torque to voltage.

  3. Note down your state-space matrices for your report. Note: You may want to cross-check the state-space matrix with TAs before proceeding to balance control.

  4. Find the open-loop poles of the system. Hint: Use eig(A).

B. Balance Control

Specification

The control design and time-response requirements are: Specification 1: Damping ratio: 0.6 < ζ\zeta < 0.8 Specification 2: Natural frequency: 3.5 rad/s < < 4.5 rad/s Specification 3: Maximum pendulum angle deflection: α|\alpha| < 15 deg. Specification 4: Maximum control effort / voltage: < 10 V. The necessary closed-loop poles are found from specifications 1 and 2. The pendulum deflection and control effort requirements (i.e. specifications 3 and 4) are to be satisfied when the rotary arm is tracking a ±20\pm 20 degree angle square wave.

Stability

The stability of a system can be determined from its poles ([2]):

  • Stable systems have poles only in the left-hand plane.

  • Unstable systems have at least one pole in the righthand plane and/or poles of multiplicity greater than 1 on the imaginary axis.

  • Marginally stable systems have one pole on the imaginary axis and the other poles in the left-hand plane.

The poles are the roots of the system’s characteristic equation. From the state-space, the characteristic equation of the system can be found using where det()\mathrm{det}()is the determinant function, ss is the Laplace operator, and II the identity matrix. These are the eigenvalues of the state-space matrix AA.

Controllability

If the control input uu of a system can take each state variable, xix_i where i=1...ni = 1 ... n, from an initial state to a final state in finite time then the system is controllable, otherwise it is uncontrollable ([2]).

Rank Test The system is controllable if the rank of its controllability matrix

equals the number of states in the system,

Companion Matrix

If (A,B) are controllable and B is n×1, then A is similar to a companion matrix ([1]). Let the characteristic equation of A be Then the companion matrices of A and B are

and

Define where TT is the controllability matrix defined in Equation 16 and Then and

Pole Placement Theory

If (A,B)(A,B) are controllable, then pole placement can be used to design the controller. Given the control law u=Kxu = -Kx, the state-space model of Equation (13) becomes

We can generalize the procedure to design a gain KK for a controllable (A,B)(A,B) system as follows:

Step 1 Find the companion matrices A~\tilde{A} and B~\tilde{B}. Compute W=TT~1W = T\tilde{T}^{-1}. Step 2 Compute K~\tilde{K} to assign the poles of A~B~K~\tilde{A} - \tilde{B}\tilde{K} to the desired locations.

Step 3 Find K=K~W1K = \tilde{K}W^{-1} to get the feedback gain for the original system (A,B)(A,B). Remark-1: It is important to do the K~K\tilde{K} \rightarrow K conversion. Remember that (A,B)(A,B) represents the actual system while the companion matrices A~\tilde{A} and B~\tilde{B} do not.

Remark-2: The entire control design procedure using the pole placement method can be simply done in MATLAB using the function called 'place' or 'acker'. For a selected desired set of closed loop poles DP, the full state feedback gain matrix KK is obtained from

>> K = acker(A,B,DP);

Desired Poles

The rotary inverted pendulum system has four poles. As depicted in Figure 3, poles p1p_1 and p2p_2 are the complex conjugate dominant poles and are chosen to satisfy the natural frequency, ωn\omega_n, and damping ratio, ζ\zeta, as given in Specification. Let the conjugate poles be and where​ σ=ζωn\sigma = \zeta\omega_n and ωd=ωn1ζ2\omega_d = \omega_n \sqrt{1-\zeta^2} is the damped natural frequency. The remaining closed-loop poles, p3p_3 and p4p_4, are placed along the real-axis to the left of the dominant poles, as shown in Figure 3.

Figure 3 Desired closed-loop pole locations

Simulation Model with Feedback

The feedback control loop that balances the rotary pendulum is illustrated in Figure 4. The reference state is defined as where θd\theta_d is the desired rotary arm angle. The controller is Note that if xd=0x_d = 0 then u=Kxu = -Kx, which is the control used in the pole-placement algorithm.

Figure 4 State-feedback control loop

B.1 Experiment: Designing the Balance Control

  1. Select ζ\zeta and ωn\omega_n such that they satisfy the given design specifications.

  2. In the same rotpen_week1_student.mlx live script, go to the Balance Control section. Enter the chosen zeta and omega_n values.

  3. Determine the locations of the two dominant poles p1p_1 and p2p_2 based on the specifications and enter their values in the MATLAB live script. Ensure that the other poles are placed at p3 = -30 and p4 = -40. Hint: Use Equation 27.

  4. Find gain K using a predefined Compensator Design MATLAB command K = acker(A,B,DP), which is based on pole-placement design. Note: DP is a row vector of the desired poles found in Step 3.

B.2 Experiment: Simulating the Balance Control

The s_rotpen_bal SIMULINK diagram shown in Figure 5 is used to simulate the closed-loop response of the Rotary Pendulum using the state-feedback control described in Balance Control with the control gain K found above. The Signal Generator block generates a 0.1 Hz square wave (with an amplitude of 1). The Amplitude (deg) gain block is used to change the desired rotary arm position. The state-feedback gain K is set in the Control Gain gain block and is read from the MATLAB workspace. The SIMULINK State-Space block reads the A, B, C, and D state-space matrices that are loaded in the MATLAB workspace. The Find State X block contains high-pass filters to find the velocity of the rotary arm and pendulum.

Figure 5: s_rotpen_bal SIMULINK diagram used to simulate the state-feedback control
  1. Run rotpen_week1_student.mlx live script. Ensure the gain K you found is loaded in the workspace (type K matrix in the command window).

  2. Open and run the s_rotpen_bal.mdl for 10 seconds. The responses in the scopes shown in Figure 6 were generated using an arbitrary feedback control gain. Note: When the simulation stops, the last 10 seconds of data is automatically saved in the MATLAB workspace to the variables data_theta, data_alpha, and data_Vm.

    Figure 6: Balance control simulation
  3. Save the data corresponding to the simulated response of the rotary arm, pendulum, and motor input voltage obtained using your obtained gain K. Note: The time is stored in the data_theta(:,1) vector, the desired and measured rotary arm angles are saved in the data_theta(:,2) and data_theta(:,3) arrays. Similarly, the pendulum angle is stored in the data_alpha(:,2) vector, and the control input is in the data_Vm(:,2) structure.

  4. Measure the pendulum deflection and voltage used. Are the specifications given satisfied?

Results for Report

A) Modeling

  1. The linear state-space representation of the rotary inverted pendulum system, i.e., AA, BB, CC and DD matrices (numerical values).

  2. Open-loop poles of the system.

B.1) Balance Control Design

  1. Chosen ζ\zeta and ωn\omega_n based on design specifications.

  2. Corresponding locations of the two dominant poles p1p_1 and p2p_2.

  3. Gain vector KK.

B.2) Balance Control Simulation

  1. Plots of the commanded position of the rotary arm (θc\theta_c), simulated responses of the rotary arm (θ\theta), pendulum (α\alpha), and motor input voltage (VmV_m) obtained using your obtained gain K.

  2. Are Design Specifications 3 and 4 satisfied? Justify using the measured maximum pendulum deflection and motor input voltage values.

Questions for Report

B.1) Balance Control Design

  1. Based on your open-loop poles found in Result A.2, is the system stable, marginally stable, or unstable? Did you expect the stability of the inverted pendulum to be as what was determined?

  2. Determine the controllability matrix TT of the system. Is the inverted pendulum system controllable? Hint: Use Equation 17.

  3. Using the open-loop poles, find the characteristic equation of AA. Hint: The roots of the characteristic equation are the open-­loop poles.

    Instead of using det(sIA)=0\mathrm{det}(sI - A) = 0 , characteristic polynomials can also be found using MATLAB function poly.

  4. Find the corresponding companion matrices A~\tilde{A} and B~\tilde{B} . Hint: For A~\tilde{A} , use the characteristic equation of A found in Question B.3 and Equation 19. For B~\tilde{B}, use Equation 20.

  5. Determine the controllability matrix T~\tilde{T} of the companion system.

  6. Determine the transformation matrix WW.

  7. Check if A~=W1AW\tilde{A} = W^{-1}AW and B~=W1B\tilde{B} = W^{-1}B with the obtained matrices.

  8. Using the locations of the two dominant poles, p1p_1 and p2p_2, based on the specifications (Result B.1), and the other poles at p3=30p_3 = -30 and p4=40p_4 = -40, determine the desired closed-loop characteristic equation. Hint: The roots of the closed-loop characteristic equation are the closed-­loop poles.

  9. When applying the control u=K~xu = -\tilde{K}x to the companion form, it changes (A~,B~)(\tilde{A}, \tilde{B}) to (A~B~K~,B~)(\tilde{A}-\tilde{B}\tilde{K}, \tilde{B}). Find the gain K~\tilde{K} that assigns the poles to their new desired location. Hint: Use Equation 26 and find the corresponding characteristic equation. Compare this equation with the desired closed-loop characteristic equation found in Question B.5 to determine the gain vector K~\tilde{K}.

  10. Once you have found K~\tilde{K}, find KK using Step 3 in Pole Placement Theory.

  11. Compare the gain vector KKcalculated using Pole Placement Theory (Question B.1.10) with the gain vector KK obtained using MATLAB (Result B.1.3).

Appendix

Table A Main Rotary Servo Base Unit Specifications

VnomV_{\rm{nom}}

Motor nominal input voltage

6.0 V

Motor armature resistance

2.6 Ω ± 12%

Motor armature inductance

0.18 mH

Motor current-torque constant

7.68 × 10−3 N-m/A ± 12%

Motor back-emf constant

7.68 × 10−3 V/(rad/s) ± 12%

High-gear total gear ratio

70

Motor efficiency

0.69 ± 5%

Gearbox efficiency

0.90 ± 10%

rotor Rotor moment of inertia

3.90 × 10−7 kg-m2 ± 10%

High-gear equivalent moment of inertia without external load

1.823 × 10−3 kg-m2

High-gear Equivalent viscous damping coefficient

0.015 N-m/(rad/s)

Mass of bar load

0.038 kg

Length of bar load

0.1525 m

Mass of disc load

0.04 kg

Radius of disc load

0.05 m

Maximum load mass

5 kg

Maximum input voltage frequency

50 Hz

Maximum input current

1 A

Maximum motor speed

628.3 rad/s

Reference

[1] Bruce Francis. Ece1619 linear systems, 2001. [2] Norman S. Nise. Control Systems Engineering. John Wiley & Sons, Inc., 2008.

Last updated

Was this helpful?