1 of 3

Lab 3: Rotary Flexible Link

A. System Identification (Week 1)

Objective

The purpose of this experiment is to understand the existence of vibrations in a rotary flexible link and the resulting mode shapes. Any beam-like structure exhibits vibrations either due to external changing loads or due to reorientation via actuators. The experiment deals with the modeling and identification of modal frequencies and mode shapes of free vibrations for a rotary flexible link. This analysis is particularly useful to understand structural vibrations and modes and how to contain them in real-world applications like aircraft and spacecraft structures as well as robotic links/manipulators.

Equipment Required

Flexible link plant (Quanser FLEXGAGE module)
Power amplifier (Quanser VoltPaq-X2)
Data acquisition board (Quanser Q2-USB)
Rotary servo plant (Quanser SRV02)
MATLAB and SIMULINK

Modeling and Vibration Analysis

This experiment involves system identification and modeling of the flexible link. The objective is to find the stiffness of the flexible link and conduct a frequency sweep across a range to determine the structural frequencies and mode shapes of the link.

Rotary Flexible Link Model

This experiment is performed using the Quanser Rotary Flexible Link mounted on an SRV-02 servo motor. This system, shown in Fig. 1, consists of an electromechanical plant, where a flexible link is rotated using a servo motor. The base of the flexible link is mounted on the load gear of the servo motor system. The servo angle, $\theta$ , increases positively when it rotates counter-clockwise (CCW). The servo (and thus the link) turn in the CCW direction when the control voltage is positive, i.e., $V_m>0$ .

The main components of the setup are labeled in Figs. 2 and 3 and are listed in Table 1.

Table 1. Setup Components

No.

Component

SRV02 Plant (Servo motor)

FLEXGAGE Module

FLEXGAGE Link

Strain Gauge

Strain Gauge Circuit

Thumbscrews

Sensor Connector

OFFSET Potentiometer

GAIN Potentiometer

The FLEXGAGE module consists of the strain gauge, the strain gauge circuitry, and a sensor connector. The flexible link is attached to this module and the strain gauge is fixed at the root of the link. The module is mounted onto the servo motor, which is the actuator for this system. The strain gauge sensor produces an analog signal proportional to the deflection of the link tip.

The link can be schematically represented as shown in Fig. 4. The flexible link has a total length of $L_l$ , a mass of $m_l$ , and its moment of inertia about the pivoted end is $J_l$ . The deflection angle of the link is denoted as $\alpha$ and increases positively when rotated CCW.

The complete flexible link system can be represented by the diagram shown in Fig. 5. The control variable is the input servo motor voltage, $V_m$ which is proportional to the angular rate of the servo motor. This generates a torque $\tau$ , at the load gear of the servo that rotates the base of the link, which is given by\tau=\displaystyle{\frac{\eta_g K_g\eta_mk_t(V_m-K_gk_m\dot{\theta})}{R_m}}=C_1V_m-C_2\dot{\theta} \qquad \qquad \qquad\tag{1.1}

\tau=\displaystyle{\frac{\eta_g K_g\eta_mk_t(V_m-K_gk_m\dot{\theta})}{R_m}}=C_1V_m-C_2\dot{\theta} \qquad \qquad \qquad\tag{1.1}

where the various constants are SRV02 parameters which are mentioned in Table 2.

The viscous friction coefficient of the servo is denoted by $B_{\rm eq}$ . This is the friction that opposes the torque being applied at the servo load gear. $J_{\rm eq}$ represents the moment of inertia of the SRV02 when there is no load. The friction acting on the link is represented by the viscous damping coefficient $B_l$ . The flexible link is modeled as a linear spring with the stiffness $K_s$ and with moment of inertia $J_l$ .

Table 2. Setup Parameters

Mass of flexible link

Length of flexible link

Width or breadth of flexible link

Thickness of flexible link

Young’s modulus of flexible link

Area moment of inertia of link cross-section

High-gear viscous damping coefficient of SRV02

Equivalent high-gear moment of inertia of SRV02 (no load)

Motor armature resistance

Motor torque constant

Motor efficiency

Back-emf constant

High-gear total gearbox ratio

Gearbox efficiency

Viscous damping coefficient of flexible link

To be calculated

Moment of inertia of flexible link about pivoted end

To be calculated

Stiffness of flexible link

To be calculated

Mass per unit length of flexible link

To be calculated

Equations of Motion

The servo and the flexible link can be modeled as a lumped mass system separated by an equivalent spring and damper, which represent the stiffness and damping coefficient of the flexible link, respectively. The equations that describe the motion of the servo and the link with respect to the servo motor torque, i.e., the dynamics, can be obtained using free body diagram (FBD) analysis of the lumped mass moments of inertia ( $J_{\rm eq}$ and $J_l$ ).

The torque balance on $J_{\rm eq}$ yield Eq. (1.2) and the torque balance on $J_l$ yield Eq. (1.3).

J_{\rm eq}\ddot{\theta}+B_{\rm eq}\dot{\theta}-B_l\dot{\alpha}-K_s\alpha=\tau \qquad \qquad \qquad \tag{1.2}

J_l\ddot{\theta}+J_l\ddot{\alpha}+B_l\dot{\alpha}+K_s\alpha=0 \qquad \qquad \qquad \tag{1.3}

On rearranging Eq. (1.3) to obtain an expression for $B_l\dot{a}+K_s\alpha$ and substituting it in Eq. (1.2), the equations of motion (EOM) for the rotary flexible link system can be obtained as

(J_{\rm eq}+J_l)\ddot{\theta}+J_l\ddot{\alpha}+B_{\rm eq}\dot{\theta}=\tau \qquad \qquad \qquad \tag{1.4}

J_l\ddot{\alpha}+J_l\ddot{\theta}+B_l\dot{\alpha}+K_s\alpha=0 \qquad \qquad \qquad \tag{1.5}

Stiffness Determination

The stiffness of the flexible link can be determined from the free oscillation of the link using a second-order model. The free-oscillatory equation of motion of this second-order system is obtained by setting the $\ddot{\theta}$ term to zero in Eq. (1.5), i.e., by holding $\theta$ constant, which is shown in Fig. 6. The resulting equation will be

J_l\ddot{\alpha}+B_l\dot{\alpha}+K_s\alpha=0 \qquad \qquad \qquad \tag{1.6}

Assuming the initial conditions $\alpha(0)=\alpha_0$ and $\dot{\alpha}(0)=0$ , the Laplace transform of Eqn. 1.6 yields:

\displaystyle{ A(s)={\frac{\displaystyle{\frac{\alpha_0}{J_l}}}{\displaystyle{s^2+\frac{B_l}{J_l}s+\frac{K_s}{J_l}}}}} \qquad \qquad \qquad \tag{1.7}

The prototype characteristic polynomial for a second-order system is defined as

s^2+2\zeta\omega_ns+\omega_n^2

where $\zeta$ is the damping ratio and $\omega_n$ is the natural frequency. Equating this to the characteristic polynomial (denominator) in Eq. (1.7) yields

\displaystyle{\omega_n^2=\frac{K_s}{J_l}} \qquad \qquad \qquad \tag{1.8}

\displaystyle{2\zeta\omega_n = \frac{B_l}{J_l}} \qquad \qquad \qquad \tag{1.9}

where $J_l$ represents the moment of inertia of the link about the pivot. This can be calculated approximately by considering the link as a rod rotating about a pivot at one edge $\displaystyle{(J_l=\frac{m_lL_l^2}{3})}$ . Equations (1.8) and (1.9) can be used to determine the stiffness and damping of the flexible link once the natural frequency and damping ratio are known.

The damping ratio of this second-order system can be found from its response (underdamped system) using the subsidence or decrement ratio given by

\displaystyle{\delta=\frac{1}{n-1}\ln{\frac{O_1}{O_n}}} \qquad \qquad \qquad \tag{1.10}

where $O_1$ is the peak of the first oscillation and $O_n$ is the peak of the nth oscillation. Note that $O_1>O_n$ , as this is a decaying response (positive damping).

The damping ratio is defined as

\displaystyle{\zeta=\frac{1}{\sqrt{1+(\displaystyle{\frac{2\pi}{\delta})^2}}}} \qquad \qquad \qquad \tag{1.11}

The period of oscillation in a system response can be found using the equation

\displaystyle{T_{\rm osc}=\frac{t_n-t_1}{n-1}} \qquad \qquad \qquad \tag{1.12}

where $t_n$ is the time of the $nth$ oscillation, $t_1$ is the time of the first peak, and $n$ is the number of oscillations considered.

From this, the damped natural frequency (in rad/s) is

\displaystyle{\omega_d=\frac{2\pi}{T_{\rm osc}}} \qquad \qquad \qquad \tag{1.13}

and the undamped natural frequency is

\displaystyle{\omega_n=\frac{\omega_d}{\sqrt{1-\zeta^2}}} \qquad \qquad \qquad \tag{1.14}

Part 1

Experimental Procedure

Mount the flexible link onto the calibration bench.
Download and open FlexLink_FreeOsc_Q2_USB.slx. This is the block diagram for this part of the experiment. Change the simulation time to 5-10 seconds.
Open the scope alpha.
Turn on the power supply.
Check you have good data and save the link deflection angle data for the free oscillation using the filename FreeOsc_1 to your folder.
Repeat the steps 7 and 8 two more times for different perturbation locations along the link (for example, around the middle and near the base) or different perturbation angles.

Analysis

Plot the measured angular deflection of the link vs. time for each case. Select peaks that are smooth (focus on the region after initial transients) as shown in the figure below. Ensure that there are at least 4-5 peaks in between the two chosen peaks. Option: Usage offindpeak MATLAB function might be helpful, but not required. If you use this function, verify the peaks identified.
From the peaks selected, determine the time period of oscillation. Use this information to determine the damped frequency, the undamped frequency, the damping ratio, the stiffness, the viscous damping coefficient and the moment of inertia of the flexible link.
Average the stiffness and the viscous damping coefficient values of the link for the three sets of data to get a single value of the link stiffness and damping coefficient.

The flexible link can be considered as a thin continuous (uniform) cantilever beam anchored at one end and free at the other end. Using the Euler-Bernoulli beam theory, the equation of motion can be written as

\displaystyle{EI\frac{\mathrm{\partial}^4Y(x,t)}{\mathrm{\partial }x^4}+m\frac{\partial ^2Y(x,t)}{\partial t^2}=q(x,t)} \qquad \qquad \qquad \tag{1.15}

where $E$ is the modulus of rigidity of beam material (assumed constant), $I$ is the area moment of inertia of the beam cross-section (assumed constant), $Y(x,t)$ is the displacement in $y$ direction at a distance $x$ from the fixed end at time $t$ , $\omega$ is the circular natural frequency, $m$ is the mass per unit length ( $m=\rho A$ , $\rho$ is the material density, $A$ is the cross-section area), $x$ is the distance measured from the fixed end and $q$ is the external applied force per unit length.

Also, the angle of deflection $\alpha$ is related to the displacement as

\displaystyle{\alpha=\frac{\partial Y}{\partial x}} \qquad \qquad \qquad \tag{1.16}

The general solution to Eq. (1.15) can be obtained using separation of variables, as in Eq. (1.17) below.

Y(x,t)=v(x)s(t) \qquad \qquad \qquad \tag{1.17}

Substituting (1.17) in (1.15), setting $q(x,t)=0$ and rearranging gives

\displaystyle{\frac{\mathrm{d}^4 v(x)}{\mathrm{d}x^4}s(t)+v(x)\frac{m}{EI}\frac{\mathrm{d}^2s(t)}{\mathrm{d}t^2}=0} \qquad \qquad \qquad \tag{1.18}

Equation (1.18) can be rewritten as

\displaystyle{\frac{\displaystyle{\frac{\mathrm{d}^4 v(x)}{\mathrm{d}x^4}}}{v(x)}=-\frac{m}{EI}\frac{\displaystyle{\frac{\mathrm{d}^2s(t)}{\mathrm{d}t^2}}}{s(t)}} \qquad \qquad \qquad \tag{1.19}

Since the left side of Eq. (1.19) is only a function of $x$ and the right side of Eq. (1.19) is only a function of $t$ , they both must equal a constant. Let this constant be $\beta^4$ . Thus, Eq. (1.19) can be written as the following two equations

\displaystyle{\frac{\mathrm{d}^4v(x)}{\mathrm{d}x^4}-\beta^4v(x)=0} \qquad \qquad \qquad \tag{1.20}

\displaystyle{\frac{\mathrm{d}^2s(t)}{\mathrm{d}t^2}+\frac{m}{EI}\beta^4s(t)=0} \qquad \qquad \qquad \tag{1.21}

The solution of Eq. (1.20) gives the displacement $v$ (as a function of $x$ ), which will be of the form

v(x)=a\sin(\beta x)+b\cos(\beta x)+c\sinh(\beta x)+d\cosh(\beta x) \qquad \qquad \tag{1.22}

where $a,b,c$ and $d$ are unknown constants. The general solution of Eq. (1.21) is given by

s(t)=g\sin(\omega t)+h\cos(\omega t) \qquad \qquad \qquad \tag{1.23}

where $\omega=\displaystyle{\sqrt{\frac{EI}{m}}\beta^2}$ , $g$ and $h$ are unknown constants.

The constants in Eq. (1.22) are determined from four boundary conditions, while the constants in Eq. (1.23) are determined from two initial conditions.

For a clamped-free or cantilever beam, the geometric boundary conditions are

Y(x=0, t)=0\to v(x=0)=0

\displaystyle{\left.\frac{\partial Y}{\mathrm{\partial} x}\right|_{(x=0,t)}=0\to \left.\frac{\mathrm{d}v}{\mathrm{d}x}\right|_{x=0}=0}

and the natural boundary conditions are

\displaystyle{M(x=L,t)=0\to EI\left.\frac{\partial^2Y}{\partial x^2}\right|_{x=L,t}=0\to \left.\frac{\mathrm{d}^2v}{\mathrm{d}x^2}\right|_{x=L}=0}

\displaystyle{V(x=L,t)=0\to -EI\left.\frac{\partial^3Y}{\partial x^3}\right|_{(x=L,t)}=0\to \left.\frac{\mathrm{d}^3v}{\mathrm{d}x^3}\right|_{(x=L)}=0}

where $M$ represents the bending moment and $V$ represents the shear force.

On substitution of the geometric boundary conditions at $x=0$ in Eq. (1.22) and its derivative, the following relations can be obtained

v(x=0)=0\to d=-b

\displaystyle{\left.\frac{\mathrm{d}v}{\mathrm{d}x}\right|_{(x=0)}=0\to c=-a}

Hence, Eq. (1.22) becomes

v(x)=a[\sin(\beta x)-\sinh(\beta x)]+b[\cos(\beta x)-\cosh(\beta x)] \qquad \qquad \qquad \tag{1.24}

On further substitution of the natural boundary conditions at $x=L$ in the derivatives of Eq. (1.24) yields

\left.\frac{\mathrm{d}^2v}{\mathrm{d}x^2}\right|_{(x=L)}=0\to a[\sin(\beta L)+\sinh(\beta L)]+b[\cos(\beta L)+\cosh(\beta L)]=0

\displaystyle{\left.\frac{\mathrm{d}^3v}{\mathrm{d}x^3}\right|_{(x=L)}=0\to a[\cos(\beta L)+\cosh(\beta L)]-b[\sin(\beta L)-\sinh(\beta L)]=0}

\displaystyle{\begin{bmatrix}\sin(\beta L)+\sinh(\beta L) & \cos(\beta L)+\cosh(\beta L)\\\cos(\beta L)+\cosh(\beta L) & -\sin(\beta L)+\sinh(\beta L)\end{bmatrix} \begin{bmatrix}a\\b\end{bmatrix}=0}

For a non-trivial solution, the determinant of the above matrix must be 0, i.e.

\displaystyle{\mathrm{det}\begin{bmatrix}\sin(\beta L)+\sinh(\beta L) & \cos(\beta L)+\cosh(\beta L)\\\cos(\beta L)+\cosh(\beta L) & -\sin(\beta L)+\sinh(\beta L)\end{bmatrix}=0}

which gives the following characteristic equation

[\sin(\beta L)+\sinh(\beta L)][-\sin(\beta L)+\sinh(\beta L)]-[\cos(\beta L)+\cosh(\beta L)]^2=0

The above equation simplifies to

\cos(\beta L)\cosh(\beta L)=-1

There are infinite solutions to this characteristic equation, which are given by

\beta L=1.875, 4.694, 7.855,...\space...

\beta_{1,2,3,...}=\frac{1.875}{L},\frac{4.694}{L}, \frac{7.855}{L},...\space...

Thus, the mode shapes are

\displaystyle{v_i(x)=[\cosh(\beta _ix)-\cos(\beta_ix)]+\frac{\cos(\beta_iL)+\cosh(\beta_iL)}{\sin(\beta_iL)+\sinh(\beta_iL)}[\sin(\beta_ix)-\sinh(\beta_ix)] \qquad \qquad \qquad \tag{1.25}}

Since $\omega = \sqrt{\frac{EI}{m}}\beta^2$ , the modal frequencies $\omega_i$ are given by

\displaystyle{\omega_i=\sqrt{\frac{EI}{m}}\beta_i^2,\space\space i=1,2,3,...\space...} \qquad \qquad \qquad \tag{1.26}

\omega_1=\sqrt{\frac{EI}{m}}\left(\frac{1.875}{L}\right)^2

\omega_2=\sqrt{\frac{EI}{m}}\left(\frac{4.694}{L}\right)^2

\omega_3=\sqrt{\frac{EI}{m}}\left(\frac{7.855}{L}\right)^2

$m$ is mass per unit length

The mode shapes of a uniform cantilever beam are shown in Fig. 7.

In order to determine the modal frequencies and mode shapes, the corresponding frequencies can be excited by providing a sinusoidal input to the link via external means. When the frequency of the input coincides with either the fundamental frequency or higher frequency modes, the corresponding modes will be excited due to resonance and their mode shapes can be observed physically. Hence, the following experiment involves a frequency sweep across a range provided as input to the flexible link via the servo motor to identify the frequencies and observe the corresponding mode shapes.

Calculations to be done before Part 2

Using the system parameters provided in Table 2, determine the mass per unit length of the beam $(m)$ .
Using Eqn. 1.26 and the values of necessary system parameters, calculate the first and second modal frequencies.

Part 2

Experimental Procedure

Mount the flexible link onto the rotary servo base.
Download and open FlexLink_ExciteMode.mdl
Open the scope alpha.
Turn on the power supply.
Ensure that the manual switch is connected to the Chirp signal input. This signal will provide a sinusoidal signal of fixed amplitude, with frequency increasing at a linear rate with time.
Open the Chirp signal command block and make sure that the Initial frequency is 0.1 Hz, the target time is 0.25 s and the Frequency at target time is 0.2 Hz. This will allow the frequency sweep to take place at a reasonable rate and ensure that the relevant frequencies are covered within the span of time.
Save the data using the format Osc_ChirpSignal into your folder.

Analysis

Plot the measured angular deflection of the link vs. time. Using the time domain plot, perform frequency analysis using Fast Fourier Transform (FFT) [refer to the Appendix for the code] or a similar technique to identify the number of dominant frequencies and their magnitudes present in the signal.

Part 3

Experimental Procedure

Again, open the FlexLink_ExciteMode.mdl Simulink file used in part 2.
Open the scope alpha.
Turn on the power supply.
Connect the manual switch to the Sine wave signal input. This signal will provide a sinusoidal signal of fixed amplitude and fixed frequency.
Open the Sine wave signal command block and make sure that the Amplitude is 3 and Phase is 0 rad. Enter the first modal frequency determined from Q2 of the calculations in rad/sec.
Tune the frequency value by increasing or decreasing in steps of 1 rad/sec until the first mode shape is clearly visible.
Observe the corresponding mode in the link and note down the number of nodes and their locations with respect to length of the link.
Save the data using the file name Osc_Sineinput into your folder.
Repeat steps 6 to 10 for the second modal frequency also identified from Q2 of the calculations.

Analysis

Plot the mode shapes using Eq. (1.25) in MATLAB (plot v(x) vs x/L) and record the node locations. Plot both theoretical shape and experiment shape. For theoretical, use $\beta_1$ and $\beta_2$ from $\beta_{1,2}=\frac{1.875}{L},\frac{4.694}{L}$ and for experimental, use two mode frequencies from FFT and $\omega = \sqrt{\frac{EI}{m}}\beta^2$ to find $\beta_{1,2}$ .
Record the number of nodes and their locations for each mode. Determine the locations as a ratio of the link length and compare them with the values obtained from the mode shape plots.

Results for Report

(A) From Part 1

Plot of the measured angular deflection of the link vs. time with the chosen peaks marked for each dataset, i.e., three plots corresponding to three perturbation locations. Plot as three subplots in one figure.
Equations for time period of oscillation, damped frequency, undamped frequency, damping ratio, stiffness and viscous damping coefficient.
Values of the damped frequency $\omega_d$ , the undamped frequency $\omega_n$ , the damping ratio $\zeta$ , the stiffness $K_s$ and the viscous damping coefficient $B_l$ of the link for each dataset, along with the corresponding average $\omega_d$ , $\omega_n$ , $\zeta$ , $K_s$ and $B_l$ value. Use Table B.1 to document these values. $\zeta$
Calculation for mass per unit length of the link $(m)$ and moment of inertia of inertia of the link ( $J_l$ )
Calculations and values for the theoretical first and second modal frequencies (calculated using Eqn. 1.26).

(B) From Part 2

Plot of measured angular deflection of the link vs. time
Fast Fourier Transform (FFT) plot of link deflection with dominant frequencies and their magnitudes identified (either marked on plot or mentioned in writing).
Compare the first dominant frequency (from FFT) with the natural frequency $\omega_n$ of the flexible link and state your observations. Explain any similarities or differences observed and the reasons behind them.
Compare the first and second dominant frequencies obtained from the FFT with the corresponding calculated first and second modal frequencies (Result A.4) and explain your observation.

(C) From Part 3

1st and 2nd Mode shapes (theoretical and experimental) plotted on the same graph with node locations marked. (plot v(x) vs x/L)
Number of nodes and their locations for each mode (theoretical and experimental). Use Table B.2 to document the values.
Comparison of experimental and theoretical mode shapes/node locations given and state the reason(s) for why they differ.

Table B.1 System Identification Parameters

Parameter

Dataset 1

Dataset 2

Dataset 3

Average

null

Table B.2 Locations of Nodes

Mode

Node Number

Location (Theoretical)

Location (Experimental)

First

...

First

Second

...

Second

Questions for Report

Compare the damped and undamped frequencies of the link and report your observation. What does this signify with respect to the flexible link damping?

Appendix

The following code performs FFT analysis on the flexible link angle response

load data.mat %data that requires to use FFT
dT = time(2)-time(1);
Fs = 1/dT;
L = length(output_alpha)-1;
f = Fs*(0:(L/2))/L;
Ya = fft(output_alpha);
Pa2 = abs(Ya/L);
Pa1 = Pa2(1:L/2+1);
Pa1(2:end-1) = 2*Pa1(2:end-1);
% Frequency Response (FFT)
plot(f,Pa1)
xlim([0 150])
xlabel('Frequency (Hz)')
ylabel('|Amplitude {\alpha}|')
title('Single-Sided Amplitude Spectrum of \alpha')
% Angular Deflection vs. Time
figure
plot(time,output_alpha)
xlabel('Time (s)')
ylabel('Deflection Angle, \alpha (deg)')
title('Flexible Link Angle vs Time')

B & C. Control Design and Controller Implementation (Week 2)

Objective

The objective is to design a controller that allows rapid reorientation of the link with very little vibrations from the link flexibility. This analysis is a preliminary approach to containing structural vibrations in real-world applications like aircraft and spacecraft structures as well as robotic links/manipulators

B. Control Design

The methodology used is to design the control law assuming full state feedback to meet a given set of specifications. This will be followed by studying the loss of controller performance with servo feedback alone. The method we will use for full state feedback is based on the linear quadratic regulator (LQR) theory, which requires the given system to be controllable.

State Space Model

The equations that describe the motions of the servo and the link with respect to the servo motor voltage are obtained using the Euler-Lagrange method and can be written as

(J_{eq} + J_l)\ddot{\theta} + J_l\ddot{\alpha} + B_{eq}\dot{\theta} = \tau \quad \quad \tag{1}

J_{l}\ddot{\theta} + J_l\ddot{\alpha} + B_{l}\dot{\theta} + K_s\alpha= 0 \quad \quad \tag{2}

Assuming the viscous damping of the link is negligible, i.e. $B_l$ = 0 and substituting the torque $\tau$ using the expression,

the EOM can be rearranged as

The various system parameters can be found in Table 2 of Part A: System Identification.

Equations 3 and 4 can be represented in state-space form

In the output equation, only the position of the servo and link angles are being measured.

As such, the velocities of the servo and link angles can be computed in the digital controller if required, e.g., by taking the derivative and filtering the result though a high-pass filter. Hence, with all the states of the system as outputs, the C and D matrices are

Linear Quadratic Regulator (LQR) Design

Controllability

equals the number of states of the system,

LQR controller

Assuming the time-domain design specifications are

Then one can use the Bryson’s rule to choose the weighting matrices as

In our case,

The reference state is defined as

and the controller is

Controller Specifications:

Controller Design Procedure:

Note: The results from the following steps have to be included in the Analysis section of the report.

>> [K,S,E] = lqr(A,B,Q,R);
Using Figure 1 as a guide, develop a Simulink model of the system.
Make sure to include servo position command input, saturation blocks, and appropriate outputs. Use scopes and To Workspace blocks so that reference, output, and response signals can be tracked. (Make sure to maintain consistency between deg and rad in SIMULINK)
- Check if all the given specifications are met.
- Check your plots with the TA.
- Check if the specifications are met (no need to tune Q and R matrices.)

C. Controller Implementation and Evaluation

In this experiment, the servo position is controlled while minimizing the link deflection using the LQR based control found. Measurements will then be taken to ensure that the specifications are satisfied.

Warning: In closed loop experiments, the link may swing around rather quickly. Stay clear of it during closed loop experiments.

Experiment

Turn the amplifier power on.
Download the q_flex.zip and open the Simulink file q_flex for controller implementation. This is the block diagram for this part of the experiment. Make sure you have K gain variable in the MATLAB workspace. Note: The Smooth Signal Generator block generates a 0.33 Hz square wave (with amplitude of 1) that is passed through a Rate Limiter block to smooth the signal (Rising slew rate of 100 and falling slew rate of -100). The Amplitude (deg) gain block is used to change the desired servo position command. The state-feedback gain K is set in the LQR Control gain block and is read from the MATLAB workspace.
Align the flexible link along the 0 mark on the servo stand.
Double click on the scopes and open them.
Ensure that the manual switch is connected to the full state feedback.
A periodic square wave is given as input voltage to the servo motor, causing the link to move left and right alternately for 10 seconds.
Copy the response data to your group folder. The response data is saved in variables data_theta, data_alpha, and data_vm. DO NOT DELETE or SAVE THE SIMULINK MODEL.
Connect the manual switch to the partial state feedback and repeat steps 4-9 again. Make sure your files name is different from full state feedback.
Turn the power off.

Results & Questions for Report

Control Design

Results

Include the results from the controller design procedure Steps 1-6.

Questions

Estimate the bandwidth for your full state and partial state feedback controllers.
Estimate the stability margins (gain and phase margins) of your full state and partial state feedback controllers.

Controller Implementation

Results

Plots of servo angle, link deflection angle and input voltage vs time for full state and partial state feedback. Compare the results and explain the differences.
Run the Simulink that you built in control design with the same command input used in the experiment. Compare the simulated response with the experimental response for both full state and partial state feedback and explain any differences between them. Note 1: Generate two sets of figures with one set containing the plots of simulation and experiment for full state feedback and the second set containing the plots of simulation and experiment for full state feedback. For each feedback case, the simulated response and corresponding experimental response must be on the same graph. Note 2: Average the first few seconds of the experimental data to obtain initial conditions for your state-space block in Simulink.

Questions

Explain why the link deflection angle goes opposite to servo rotation initially.
For a controllable system, can partial state feedback guarantee stability like full state feedback? Briefly explain.

B & C. Control Design and Controller Implementation (Week 2)

Objective

B. Control Design

State Space Model

The equations that describe the motions of the servo and the link with respect to the servo motor voltage are obtained using the Euler-Lagrange method and can be written as

(J_{eq} + J_l)\ddot{\theta} + J_l\ddot{\alpha} + B_{eq}\dot{\theta} = \tau \quad \quad \tag{1}

J_{l}\ddot{\theta} + J_l\ddot{\alpha} + B_{l}\dot{\theta} + K_s\alpha= 0 \quad \quad \tag{2}

Assuming the viscous damping of the link is negligible, i.e. $B_l$ = 0 and substituting the torque $\tau$ using the expression,

\tau = \displaystyle\frac{\eta_gK_g\eta_mk_t(V_m - K_gk_m\dot\theta)}{R_m} = C_1V_m - C_2\dot\theta

the EOM can be rearranged as

\ddot{\theta} = -(\frac{B_{eq} + C_2}{J_eq})\dot{\theta} + \frac{K_s}{J_{eq}}\alpha + \frac{C_1}{J_{eq}}V_m \quad \quad \tag{3}

{\ddot{\alpha}} = (\frac{B_{eq}+C_2}{J_{eq}}){\dot{\theta}} - K_s(\frac{J_{eq}+J_l}{J_{eq}J_l}){\alpha }- \frac{C_1}{J_{eq}}V_m \quad \quad \tag{4}

The various system parameters can be found in Table 2 of Part A: System Identification.

Equations 3 and 4 can be represented in state-space form

\dot{x} = Ax + Bu \\ y = Cx + Du

with $x=[\theta \;\; \alpha \;\; \dot{\theta} \;\; \dot{\alpha}]^T$ and $u = V_m$ , where

A = \begin{bmatrix} 0 &0 & 1 &0 \\ 0 &0 & 0 &1 \\ 0 &\frac{K_s}{J_{eq}} & -(\frac{B_{eq}+C_2}{J_{eq}}) &0 \\ 0 &-K_s(\frac{J_{eq}+J_l}{J_{eq}J_l}) & (\frac{B_{eq}+C_2}{J_{eq}}) & 0 \end{bmatrix} \quad \quad \tag{5}

B = \begin{bmatrix} 0 \\ 0\\ \frac{C_1}{J_{eq}} \\ -\frac{C_1}{J_{eq}} \end{bmatrix} \quad \quad \tag{6}

In the output equation, only the position of the servo and link angles are being measured.

C = \begin{bmatrix} 1 & 0& 0& 0\\ 0 & 1& 0& 0\\ 0 & 0& 1& 0\\ 0 & 0& 0& 1\ \end{bmatrix} \qquad \qquad \tag{7}

D = \begin{bmatrix} 0\\ 0 \\0 \\ 0 \end{bmatrix} \qquad \qquad \tag{8}

Linear Quadratic Regulator (LQR) Design

Controllability

If the control input $'u'$ of a system can take each state variable, $x_i$ where $i = 1 . . . n$ , from an initial state to a final state in finite time, then the system is controllable, otherwise it is uncontrollable. The controllability of a system can be determined using the ‘Rank Test’. According to the Rank Test, the system is controllable if the rank of its controllability matrix

T=[B \;\; AB \;\;A^2B… A^{n−1}B]

equals the number of states of the system,

rank(T)=n

LQR controller

If (A,B) are controllable, then the Linear Quadratic Regular optimization method can be used to find full state feedback controller gains. Given the plant model in state-space form where $x$ represents the state vector, this technique finds a control input ‘ $u$ ’ that minimizes the cost function

J= \int_0^\infty(x^TQx+ u^TRu)dt \qquad \qquad \tag{9}

where Q and R are the user-defined weighting matrices $Q = Q^T$ and $R = R^T$ . It is assumed that Q is positive semi-definite (all the eigenvalues of $Q$ are non-negative) and $R$ is positive definite (all the eigenvalues of $R$ are positive). The weighting matrices affect how LQR minimizes the function and are, essentially, tuning variables. Here, $x$ is the state vector $[\theta \;\; \alpha \;\; \dot{\theta} \;\; \dot{\alpha}]^T$ and $u$ is the servo motor voltage input $V_m$ .

Given the control law $u = −Kx$ , the state-space of the closed loop system becomes

\dot{x} = Ax+B(−Kx)=(A−BK)x \quad \quad \tag{10}

where $K$ is given by

K=[k_{\theta} \;\;k_{\alpha}\;\;k_{\dot{\theta}}\;\;k_{\dot{\alpha}}] \qquad \qquad \tag{11}

is the controller gain vector. Note that the input voltage $Vm$ is scalar, and therefore, $R$ is also scalar. In practice, the weighting matrices $Q$ and $R$ are chosen to be diagonal. A simple way for choosing $Q$ and $R$ is using the Bryson's rule. This method works best when control design specifications regarding the system response and control input are given in terms of limits on excursions from equilibrium.

Assuming the time-domain design specifications are

|x_i| \leq x_{i_{\max}} \; i = 1, 2, ..., n \; \; \; \\ |u_j| \leq u_{j_{\max}} \; j = 1, 2, ..., m

Then one can use the Bryson’s rule to choose the weighting matrices as

Q=\displaystyle {diag\left (\frac{1}{x_{i_{\max}}^2}\right) \; \; i = 1, 2, ..., n} \quad \quad (12) \\ R = diag\left (\frac{1}{u_{j_{\max}}^2}\right) \; \; j = 1, 2, ..., m \quad \quad (13)

Further tuning of $Q$ and $R$ may be required for improving the controller performance. Note that increasing $R$ means increasing the penalty on the control input (using less control effort). Similarly, increasing $Q$ means increasing the penalty on the states (states reaching equilibrium position quickly and with less overshoot). Matrix $Q$ sets the weight on the states and determines how $u$ will minimize $J$ (and hence how it generates gain $K$ ).

In our case,

Q = \begin{bmatrix} \frac{1}{\theta_{max}^2} & 0& 0& 0\\ 0 & \frac{1}{\alpha_{max}^2}& 0& 0\\ 0 & 0& \frac{1}{\dot{\theta}_{max}^2}& 0\\ 0 & 0& 0& \frac{1}{\dot{\alpha}_{max}^2}\ \end{bmatrix} \qquad \qquad \tag{14}

R = \frac{1}{V^2_{m, max}} \quad \quad (15)

Substituting $Q$ and $R$ in Eqn. 9 gives

J= ∫_0^\infty (q_{θ}\theta^2+q_{α}\alpha^2+ q_{\dot{θ}}{\dot{\theta}}^2+q_{\dot{α}}{\dot{\alpha}}^2+RV_m^2)dt \quad \quad \tag{16}

where $q_\theta= \displaystyle\frac{1}{\theta_{max}^2}$ , $q_\alpha= \displaystyle\frac{1}{\alpha_{max}^2}$ , $q_{\dot{\theta}}= \displaystyle\frac{1}{\dot{\theta}_{max}^2}$ , $q_{\dot{\alpha}}= \displaystyle\frac{1}{\dot{\alpha}^2}$ , $R=\displaystyle\frac{1}{V_{m,max}^2}$

On observing Eq. 12, it can be seen that increasing $q_i$ causes control input $V_m$ to work harder to minimize state $x_i$ , which will predominately increase $k_i$ . If $q_{\theta}$ is increased, then $k_{\theta}$ will increase to compensate for the larger weight placed on state $\theta$ . Depending on the model, changing a single $q_i$ can affect multiple $k_i$ gains because each state is not independent. If $R$ is increased, then control input $V$ has to work less to minimize $J$ . In that case, LQR would generate a lower overall value of $K$ . A block diagram of the system with fill state feedback controller is shown in Figure 1.

The feedback control loop in Figure 1 is designed to stabilize the servo to a desired position, $θ_d$ , while minimizing the deflection of the flexible link, where $x$ is the state column vector and $K$ is the gain row vector as shown earlier. Note that $C$ is $4×4$ identity matrix and $D$ is $4×1$ column of zeros in Figure 1.

The reference state is defined as

x_d=[θ_d \;\; 0\;\; 0\;\;0]^T

and the controller is

u=K(x_d−x)

Controller Specifications:

The following specifications in the time domain are to be met when the rotary arm is tracking a $30°$ step change in its angular position.

Specification 1: Servo angle 2% settling time: $t_s≤0.5s$

Specification 2: Servo angle percentage overshoot: $PO≤5 \%$

Specification 3: Maximum link angle deflection: $|α|_{max}≤10°$

Specification 4: Maximum control effort (voltage): $|V_{m}|_{max}≤10 V$

Controller Design Procedure:

Note: The results from the following steps have to be included in the Analysis section of the report.

Obtain the $A$ and $B$ matrices of the system using the parameters from Part A: System identification. Check your A and B matrices with the TA.
Using Bryson’s rule (Eq. 12-15) calculate $Q$ and $R$ matrices, based on the following suggested maximum values of state and control excursions:
$0.05 rad ≤|θ_{max}|≤0.1 rad \;\;\;\;\;\; |α_{max}|\simeq0.1 rad$ $∣\dot{θ}_{max}∣\simeq1rad/s \;\;\;\; \; \; 0.4rad/s≤|\dot{α}_{max}|≤0.6rad/s$ $0.8 V ≤ |V_{m,max}|≤1.2 V$ Check if $Q \ge 0$ and $R > 0$ (by finding eigenvalues of the matrix).
Using the lqr command in MATLAB, calculate the optimal feedback gain vector $K$ . The lqr command solves the Algebraic Riccati equation which is required to determine the gain vector $K$ is
>> [K,S,E] = lqr(A,B,Q,R);
where $K$ is the gain matrix, $S$ is the solution to the Riccati equation and E are closed loop system poles. Check your gain matrix K with the TA.
Using Figure 1 as a guide, develop a Simulink model of the system.
Note: SIMULINK references the current MATLAB workspace when setting variables. Use the variable names for the system matrices $A$ and $B$ and the control gain matrix $K$ from your MATLAB script file for the LQR design.
Make sure to include servo position command input, saturation blocks, and appropriate outputs. Use scopes and To Workspace blocks so that reference, output, and response signals can be tracked. (Make sure to maintain consistency between deg and rad in SIMULINK)
Using the controller gains found in step 3, simulate the closed loop system for a step servo position command of $30\degree$ over a 10s duration.
- Record (and plot for lab report) the servo angle ( $\theta$ ) link angle ( $\alpha$ ) and voltage input $V_m$ .
- Check if all the given specifications are met.
- If not, keep adjusting the $Q$ and $R$ weighting matrices parameters, till all specifications are met. Record the final Q and R matrices, the servo angle ( $\theta$ ), link angle ( $\alpha$ ) and voltage input $V_m$ plots.
Study in simulation the degradation, if any, in closed loop system performance (in meeting the given specifications) using partial feedback. For partial feedback consider feedback of servo angle ( $\theta$ ) and servo angle rate ( $\dot{\theta}$ ).
- This is easily done by retaining the gain value of the chosen feedback loops and zeroing out the remaining gains in the gain matrix $K$ used in your Simulink model.
- Check your plots with the TA.
- Check if the specifications are met (no need to tune Q and R matrices.)
- Record (and plot for lab report) the servo angle ( $\theta$ ) link angle ( $\alpha$ ) and voltage input $V_m$

C. Controller Implementation and Evaluation

Warning: In closed loop experiments, the link may swing around rather quickly. Stay clear of it during closed loop experiments.

Experiment

Turn the amplifier power on.
Download the q_flex.zip and open the Simulink file q_flex for controller implementation. This is the block diagram for this part of the experiment. Make sure you have K gain variable in the MATLAB workspace. Note: The Smooth Signal Generator block generates a 0.33 Hz square wave (with amplitude of 1) that is passed through a Rate Limiter block to smooth the signal (Rising slew rate of 100 and falling slew rate of -100). The Amplitude (deg) gain block is used to change the desired servo position command. The state-feedback gain K is set in the LQR Control gain block and is read from the MATLAB workspace.
To build the model, click down arrow on Monitor & Tune under Hardware tab and then Build for monitoring . This generates the controller code.
Align the flexible link along the 0 mark on the servo stand.
Double click on the scopes and open them.
Ensure that the manual switch is connected to the full state feedback.
Click Connect button under Monitor & Tune and then click Start
A periodic square wave is given as input voltage to the servo motor, causing the link to move left and right alternately for 10 seconds.
Copy the response data to your group folder. The response data is saved in variables data_theta, data_alpha, and data_vm. DO NOT DELETE or SAVE THE SIMULINK MODEL.
Connect the manual switch to the partial state feedback and repeat steps 4-9 again. Make sure your files name is different from full state feedback.
Turn the power off.

Results & Questions for Report

Control Design

Results

Include the results from the controller design procedure Steps 1-6.

Questions

Estimate the bandwidth for your full state and partial state feedback controllers.
Estimate the stability margins (gain and phase margins) of your full state and partial state feedback controllers.

Controller Implementation

Results

Plots of servo angle, link deflection angle and input voltage vs time for full state and partial state feedback. Compare the results and explain the differences.
Run the Simulink that you built in control design with the same command input used in the experiment. Compare the simulated response with the experimental response for both full state and partial state feedback and explain any differences between them. Note 1: Generate two sets of figures with one set containing the plots of simulation and experiment for full state feedback and the second set containing the plots of simulation and experiment for full state feedback. For each feedback case, the simulated response and corresponding experimental response must be on the same graph. Note 2: Average the first few seconds of the experimental data to obtain initial conditions for your state-space block in Simulink.

Questions

Explain why the link deflection angle goes opposite to servo rotation initially.
For a controllable system, can partial state feedback guarantee stability like full state feedback? Briefly explain.

A. System Identification (Week 1)

Objective

Equipment Required

Flexible link plant (Quanser FLEXGAGE module)
Power amplifier (Quanser VoltPaq-X2)
Data acquisition board (Quanser Q2-USB)
Rotary servo plant (Quanser SRV02)
MATLAB and SIMULINK

Modeling and Vibration Analysis

Rotary Flexible Link Model

The main components of the setup are labeled in Figs. 2 and 3 and are listed in Table 1.

Table 1. Setup Components

No.

Component

SRV02 Plant (Servo motor)

FLEXGAGE Module

FLEXGAGE Link

Strain Gauge

Strain Gauge Circuit

Thumbscrews

Sensor Connector

OFFSET Potentiometer

GAIN Potentiometer

\tau=\displaystyle{\frac{\eta_g K_g\eta_mk_t(V_m-K_gk_m\dot{\theta})}{R_m}}=C_1V_m-C_2\dot{\theta} \qquad \qquad \qquad\tag{1.1}

where the various constants are SRV02 parameters which are mentioned in Table 2.

Table 2. Setup Parameters

Mass of flexible link

Length of flexible link

Width or breadth of flexible link

Thickness of flexible link

Young’s modulus of flexible link

Area moment of inertia of link cross-section

High-gear viscous damping coefficient of SRV02

Equivalent high-gear moment of inertia of SRV02 (no load)

Motor armature resistance

Motor torque constant

Motor efficiency

Back-emf constant

High-gear total gearbox ratio

Gearbox efficiency

Viscous damping coefficient of flexible link

To be calculated

Moment of inertia of flexible link about pivoted end

To be calculated

Stiffness of flexible link

To be calculated

Mass per unit length of flexible link

To be calculated

Equations of Motion

The torque balance on $J_{\rm eq}$ yield Eq. (1.2) and the torque balance on $J_l$ yield Eq. (1.3).

J_{\rm eq}\ddot{\theta}+B_{\rm eq}\dot{\theta}-B_l\dot{\alpha}-K_s\alpha=\tau \qquad \qquad \qquad \tag{1.2}

J_l\ddot{\theta}+J_l\ddot{\alpha}+B_l\dot{\alpha}+K_s\alpha=0 \qquad \qquad \qquad \tag{1.3}

On rearranging Eq. (1.3) to obtain an expression for $B_l\dot{a}+K_s\alpha$ and substituting it in Eq. (1.2), the equations of motion (EOM) for the rotary flexible link system can be obtained as

(J_{\rm eq}+J_l)\ddot{\theta}+J_l\ddot{\alpha}+B_{\rm eq}\dot{\theta}=\tau \qquad \qquad \qquad \tag{1.4}

J_l\ddot{\alpha}+J_l\ddot{\theta}+B_l\dot{\alpha}+K_s\alpha=0 \qquad \qquad \qquad \tag{1.5}

Stiffness Determination

J_l\ddot{\alpha}+B_l\dot{\alpha}+K_s\alpha=0 \qquad \qquad \qquad \tag{1.6}

Assuming the initial conditions $\alpha(0)=\alpha_0$ and $\dot{\alpha}(0)=0$ , the Laplace transform of Eqn. 1.6 yields:

\displaystyle{ A(s)={\frac{\displaystyle{\frac{\alpha_0}{J_l}}}{\displaystyle{s^2+\frac{B_l}{J_l}s+\frac{K_s}{J_l}}}}} \qquad \qquad \qquad \tag{1.7}

The prototype characteristic polynomial for a second-order system is defined as

s^2+2\zeta\omega_ns+\omega_n^2

where $\zeta$ is the damping ratio and $\omega_n$ is the natural frequency. Equating this to the characteristic polynomial (denominator) in Eq. (1.7) yields

\displaystyle{\omega_n^2=\frac{K_s}{J_l}} \qquad \qquad \qquad \tag{1.8}

\displaystyle{2\zeta\omega_n = \frac{B_l}{J_l}} \qquad \qquad \qquad \tag{1.9}

The damping ratio of this second-order system can be found from its response (underdamped system) using the subsidence or decrement ratio given by

\displaystyle{\delta=\frac{1}{n-1}\ln{\frac{O_1}{O_n}}} \qquad \qquad \qquad \tag{1.10}

where $O_1$ is the peak of the first oscillation and $O_n$ is the peak of the nth oscillation. Note that $O_1>O_n$ , as this is a decaying response (positive damping).

The damping ratio is defined as

\displaystyle{\zeta=\frac{1}{\sqrt{1+(\displaystyle{\frac{2\pi}{\delta})^2}}}} \qquad \qquad \qquad \tag{1.11}

The period of oscillation in a system response can be found using the equation

\displaystyle{T_{\rm osc}=\frac{t_n-t_1}{n-1}} \qquad \qquad \qquad \tag{1.12}

where $t_n$ is the time of the $nth$ oscillation, $t_1$ is the time of the first peak, and $n$ is the number of oscillations considered.

From this, the damped natural frequency (in rad/s) is

\displaystyle{\omega_d=\frac{2\pi}{T_{\rm osc}}} \qquad \qquad \qquad \tag{1.13}

and the undamped natural frequency is

\displaystyle{\omega_n=\frac{\omega_d}{\sqrt{1-\zeta^2}}} \qquad \qquad \qquad \tag{1.14}

Part 1

Experimental Procedure

Mount the flexible link onto the calibration bench.
Download and open FlexLink_FreeOsc_Q2_USB.slx. This is the block diagram for this part of the experiment. Change the simulation time to 5-10 seconds.
To build the model, click the down arrow on Monitor & Tune under the Hardware tab and then click Build for monitoring . This generates the controller code.
Open the scope alpha.
Turn on the power supply.
Press Connect button under Monitor & Tune and hold on to the base to prevent any rotation at the root.
Press Start and immediately perturb the flexible link (for example, at the tip). Keep holding the base until the data is collected for the complete run.
Check you have good data and save the link deflection angle data for the free oscillation using the filename FreeOsc_1 to your folder.
Repeat the steps 7 and 8 two more times for different perturbation locations along the link (for example, around the middle and near the base) or different perturbation angles.

Analysis

Plot the measured angular deflection of the link vs. time for each case. Select peaks that are smooth (focus on the region after initial transients) as shown in the figure below. Ensure that there are at least 4-5 peaks in between the two chosen peaks. Option: Usage offindpeak MATLAB function might be helpful, but not required. If you use this function, verify the peaks identified.
From the peaks selected, determine the time period of oscillation. Use this information to determine the damped frequency, the undamped frequency, the damping ratio, the stiffness, the viscous damping coefficient and the moment of inertia of the flexible link.
Average the stiffness and the viscous damping coefficient values of the link for the three sets of data to get a single value of the link stiffness and damping coefficient.

\displaystyle{EI\frac{\mathrm{\partial}^4Y(x,t)}{\mathrm{\partial }x^4}+m\frac{\partial ^2Y(x,t)}{\partial t^2}=q(x,t)} \qquad \qquad \qquad \tag{1.15}

Also, the angle of deflection $\alpha$ is related to the displacement as

\displaystyle{\alpha=\frac{\partial Y}{\partial x}} \qquad \qquad \qquad \tag{1.16}

The general solution to Eq. (1.15) can be obtained using separation of variables, as in Eq. (1.17) below.

Y(x,t)=v(x)s(t) \qquad \qquad \qquad \tag{1.17}

Substituting (1.17) in (1.15), setting $q(x,t)=0$ and rearranging gives

\displaystyle{\frac{\mathrm{d}^4 v(x)}{\mathrm{d}x^4}s(t)+v(x)\frac{m}{EI}\frac{\mathrm{d}^2s(t)}{\mathrm{d}t^2}=0} \qquad \qquad \qquad \tag{1.18}

Equation (1.18) can be rewritten as

\displaystyle{\frac{\displaystyle{\frac{\mathrm{d}^4 v(x)}{\mathrm{d}x^4}}}{v(x)}=-\frac{m}{EI}\frac{\displaystyle{\frac{\mathrm{d}^2s(t)}{\mathrm{d}t^2}}}{s(t)}} \qquad \qquad \qquad \tag{1.19}

\displaystyle{\frac{\mathrm{d}^4v(x)}{\mathrm{d}x^4}-\beta^4v(x)=0} \qquad \qquad \qquad \tag{1.20}

\displaystyle{\frac{\mathrm{d}^2s(t)}{\mathrm{d}t^2}+\frac{m}{EI}\beta^4s(t)=0} \qquad \qquad \qquad \tag{1.21}

The solution of Eq. (1.20) gives the displacement $v$ (as a function of $x$ ), which will be of the form

v(x)=a\sin(\beta x)+b\cos(\beta x)+c\sinh(\beta x)+d\cosh(\beta x) \qquad \qquad \tag{1.22}

where $a,b,c$ and $d$ are unknown constants. The general solution of Eq. (1.21) is given by

s(t)=g\sin(\omega t)+h\cos(\omega t) \qquad \qquad \qquad \tag{1.23}

where $\omega=\displaystyle{\sqrt{\frac{EI}{m}}\beta^2}$ , $g$ and $h$ are unknown constants.

The constants in Eq. (1.22) are determined from four boundary conditions, while the constants in Eq. (1.23) are determined from two initial conditions.

For a clamped-free or cantilever beam, the geometric boundary conditions are

Y(x=0, t)=0\to v(x=0)=0

\displaystyle{\left.\frac{\partial Y}{\mathrm{\partial} x}\right|_{(x=0,t)}=0\to \left.\frac{\mathrm{d}v}{\mathrm{d}x}\right|_{x=0}=0}

and the natural boundary conditions are

\displaystyle{M(x=L,t)=0\to EI\left.\frac{\partial^2Y}{\partial x^2}\right|_{x=L,t}=0\to \left.\frac{\mathrm{d}^2v}{\mathrm{d}x^2}\right|_{x=L}=0}

\displaystyle{V(x=L,t)=0\to -EI\left.\frac{\partial^3Y}{\partial x^3}\right|_{(x=L,t)}=0\to \left.\frac{\mathrm{d}^3v}{\mathrm{d}x^3}\right|_{(x=L)}=0}

where $M$ represents the bending moment and $V$ represents the shear force.

On substitution of the geometric boundary conditions at $x=0$ in Eq. (1.22) and its derivative, the following relations can be obtained

v(x=0)=0\to d=-b

\displaystyle{\left.\frac{\mathrm{d}v}{\mathrm{d}x}\right|_{(x=0)}=0\to c=-a}

Hence, Eq. (1.22) becomes

v(x)=a[\sin(\beta x)-\sinh(\beta x)]+b[\cos(\beta x)-\cosh(\beta x)] \qquad \qquad \qquad \tag{1.24}

On further substitution of the natural boundary conditions at $x=L$ in the derivatives of Eq. (1.24) yields

\left.\frac{\mathrm{d}^2v}{\mathrm{d}x^2}\right|_{(x=L)}=0\to a[\sin(\beta L)+\sinh(\beta L)]+b[\cos(\beta L)+\cosh(\beta L)]=0

\displaystyle{\left.\frac{\mathrm{d}^3v}{\mathrm{d}x^3}\right|_{(x=L)}=0\to a[\cos(\beta L)+\cosh(\beta L)]-b[\sin(\beta L)-\sinh(\beta L)]=0}

\displaystyle{\begin{bmatrix}\sin(\beta L)+\sinh(\beta L) & \cos(\beta L)+\cosh(\beta L)\\\cos(\beta L)+\cosh(\beta L) & -\sin(\beta L)+\sinh(\beta L)\end{bmatrix} \begin{bmatrix}a\\b\end{bmatrix}=0}

For a non-trivial solution, the determinant of the above matrix must be 0, i.e.

\displaystyle{\mathrm{det}\begin{bmatrix}\sin(\beta L)+\sinh(\beta L) & \cos(\beta L)+\cosh(\beta L)\\\cos(\beta L)+\cosh(\beta L) & -\sin(\beta L)+\sinh(\beta L)\end{bmatrix}=0}

which gives the following characteristic equation

[\sin(\beta L)+\sinh(\beta L)][-\sin(\beta L)+\sinh(\beta L)]-[\cos(\beta L)+\cosh(\beta L)]^2=0

The above equation simplifies to

\cos(\beta L)\cosh(\beta L)=-1

There are infinite solutions to this characteristic equation, which are given by

\beta L=1.875, 4.694, 7.855,...\space...

\beta_{1,2,3,...}=\frac{1.875}{L},\frac{4.694}{L}, \frac{7.855}{L},...\space...

Thus, the mode shapes are

\displaystyle{v_i(x)=[\cosh(\beta _ix)-\cos(\beta_ix)]+\frac{\cos(\beta_iL)+\cosh(\beta_iL)}{\sin(\beta_iL)+\sinh(\beta_iL)}[\sin(\beta_ix)-\sinh(\beta_ix)] \qquad \qquad \qquad \tag{1.25}}

Since $\omega = \sqrt{\frac{EI}{m}}\beta^2$ , the modal frequencies $\omega_i$ are given by

\displaystyle{\omega_i=\sqrt{\frac{EI}{m}}\beta_i^2,\space\space i=1,2,3,...\space...} \qquad \qquad \qquad \tag{1.26}

\omega_1=\sqrt{\frac{EI}{m}}\left(\frac{1.875}{L}\right)^2

\omega_2=\sqrt{\frac{EI}{m}}\left(\frac{4.694}{L}\right)^2

\omega_3=\sqrt{\frac{EI}{m}}\left(\frac{7.855}{L}\right)^2

$m$ is mass per unit length

The mode shapes of a uniform cantilever beam are shown in Fig. 7.

Calculations to be done before Part 2

Using the system parameters provided in Table 2, determine the mass per unit length of the beam $(m)$ .
Using Eqn. 1.26 and the values of necessary system parameters, calculate the first and second modal frequencies.

Part 2

Experimental Procedure

Mount the flexible link onto the rotary servo base.
Download and open FlexLink_ExciteMode.mdl
To build the model, click the down arrow on Monitor & Tune under the Hardware tab and then click Build for monitoring . This generates the controller code.
Open the scope alpha.
Turn on the power supply.
Ensure that the manual switch is connected to the Chirp signal input. This signal will provide a sinusoidal signal of fixed amplitude, with frequency increasing at a linear rate with time.
Open the Chirp signal command block and make sure that the Initial frequency is 0.1 Hz, the target time is 0.25 s and the Frequency at target time is 0.2 Hz. This will allow the frequency sweep to take place at a reasonable rate and ensure that the relevant frequencies are covered within the span of time.
Press the Connect button under Monitor & Tune and click on Start . Run the servo motor with the chirp input voltage for 60 seconds.
Save the data using the format Osc_ChirpSignal into your folder.

Analysis

Plot the measured angular deflection of the link vs. time. Using the time domain plot, perform frequency analysis using Fast Fourier Transform (FFT) [refer to the Appendix for the code] or a similar technique to identify the number of dominant frequencies and their magnitudes present in the signal.

Part 3

Experimental Procedure

Again, open the FlexLink_ExciteMode.mdl Simulink file used in part 2.
To build the model, click the down arrow on Monitor & Tune under the Hardware tab and then click Build for monitoring . This generates the controller code.
Open the scope alpha.
Turn on the power supply.
Connect the manual switch to the Sine wave signal input. This signal will provide a sinusoidal signal of fixed amplitude and fixed frequency.
Open the Sine wave signal command block and make sure that the Amplitude is 3 and Phase is 0 rad. Enter the first modal frequency determined from Q2 of the calculations in rad/sec.
Press Connect button under Monitor & Tune and Press Start . Run the servo motor with the sine wave for at least 20 seconds.
Tune the frequency value by increasing or decreasing in steps of 1 rad/sec until the first mode shape is clearly visible.
Observe the corresponding mode in the link and note down the number of nodes and their locations with respect to length of the link.
Save the data using the file name Osc_Sineinput into your folder.
Repeat steps 6 to 10 for the second modal frequency also identified from Q2 of the calculations.

Analysis

Plot the mode shapes using Eq. (1.25) in MATLAB (plot v(x) vs x/L) and record the node locations. Plot both theoretical shape and experiment shape. For theoretical, use $\beta_1$ and $\beta_2$ from $\beta_{1,2}=\frac{1.875}{L},\frac{4.694}{L}$ and for experimental, use two mode frequencies from FFT and $\omega = \sqrt{\frac{EI}{m}}\beta^2$ to find $\beta_{1,2}$ .
Record the number of nodes and their locations for each mode. Determine the locations as a ratio of the link length and compare them with the values obtained from the mode shape plots.

Results for Report

(A) From Part 1

Plot of the measured angular deflection of the link vs. time with the chosen peaks marked for each dataset, i.e., three plots corresponding to three perturbation locations. Plot as three subplots in one figure.
Equations for time period of oscillation, damped frequency, undamped frequency, damping ratio, stiffness and viscous damping coefficient.
Values of the damped frequency $\omega_d$ , the undamped frequency $\omega_n$ , the damping ratio $\zeta$ , the stiffness $K_s$ and the viscous damping coefficient $B_l$ of the link for each dataset, along with the corresponding average $\omega_d$ , $\omega_n$ , $\zeta$ , $K_s$ and $B_l$ value. Use Table B.1 to document these values. $\zeta$
Calculation for mass per unit length of the link $(m)$ and moment of inertia of inertia of the link ( $J_l$ )
Calculations and values for the theoretical first and second modal frequencies (calculated using Eqn. 1.26).

(B) From Part 2

Plot of measured angular deflection of the link vs. time
Fast Fourier Transform (FFT) plot of link deflection with dominant frequencies and their magnitudes identified (either marked on plot or mentioned in writing).
Compare the first dominant frequency (from FFT) with the natural frequency $\omega_n$ of the flexible link and state your observations. Explain any similarities or differences observed and the reasons behind them.
Compare the first and second dominant frequencies obtained from the FFT with the corresponding calculated first and second modal frequencies (Result A.4) and explain your observation.

(C) From Part 3

1st and 2nd Mode shapes (theoretical and experimental) plotted on the same graph with node locations marked. (plot v(x) vs x/L)
Number of nodes and their locations for each mode (theoretical and experimental). Use Table B.2 to document the values.
Comparison of experimental and theoretical mode shapes/node locations given and state the reason(s) for why they differ.

Table B.1 System Identification Parameters

Parameter

Dataset 1

Dataset 2

Dataset 3

Average

null

Table B.2 Locations of Nodes

Mode

Node Number

Location (Theoretical)

Location (Experimental)

First

...

First

Second

...

Second

Questions for Report

Compare the damped and undamped frequencies of the link and report your observation. What does this signify with respect to the flexible link damping?

Appendix

The following code performs FFT analysis on the flexible link angle response

load data.mat %data that requires to use FFT
dT = time(2)-time(1);
Fs = 1/dT;
L = length(output_alpha)-1;
f = Fs*(0:(L/2))/L;
Ya = fft(output_alpha);
Pa2 = abs(Ya/L);
Pa1 = Pa2(1:L/2+1);
Pa1(2:end-1) = 2*Pa1(2:end-1);
% Frequency Response (FFT)
plot(f,Pa1)
xlim([0 150])
xlabel('Frequency (Hz)')
ylabel('|Amplitude {\alpha}|')
title('Single-Sided Amplitude Spectrum of \alpha')
% Angular Deflection vs. Time
figure
plot(time,output_alpha)
xlabel('Time (s)')
ylabel('Deflection Angle, \alpha (deg)')
title('Flexible Link Angle vs Time')

Lab 3: Rotary Flexible Link

A. System Identification (Week 1)

Objective

Equipment Required

Modeling and Vibration Analysis

Rotary Flexible Link Model

Equations of Motion

Stiffness Determination

Part 1

Experimental Procedure

Analysis

Modal Frequencies and Mode Shapes

Calculations to be done before Part 2

Part 2

Experimental Procedure

Analysis

Part 3

Experimental Procedure

Analysis

Results for Report

(A) From Part 1

(B) From Part 2

(C) From Part 3

Table B.1 System Identification Parameters

Table B.2 Locations of Nodes

Questions for Report

Appendix

B & C. Control Design and Controller Implementation (Week 2)

Objective

B. Control Design

State Space Model

Linear Quadratic Regulator (LQR) Design

Controllability

LQR controller

Controller Specifications:

Controller Design Procedure:

C. Controller Implementation and Evaluation

Experiment

Results & Questions for Report

Control Design

Controller Implementation

Lab 3: Rotary Flexible Link

B & C. Control Design and Controller Implementation (Week 2)

Objective

B. Control Design

State Space Model

Linear Quadratic Regulator (LQR) Design

Controllability

LQR controller

Controller Specifications:

Controller Design Procedure:

C. Controller Implementation and Evaluation

Experiment

Results & Questions for Report

Control Design

Controller Implementation

A. System Identification (Week 1)

Objective

Equipment Required

Modeling and Vibration Analysis

Rotary Flexible Link Model

Equations of Motion

Stiffness Determination

Part 1

Experimental Procedure

Analysis

Modal Frequencies and Mode Shapes

Calculations to be done before Part 2

Part 2

Experimental Procedure

Analysis

Part 3

Experimental Procedure

Analysis

Results for Report

(A) From Part 1

(B) From Part 2

(C) From Part 3

Table B.1 System Identification Parameters

Table B.2 Locations of Nodes