Lab 0: MATLAB & Simulink Problems

Important instructions for Lab 0 report: There is no need to prepare Lab 0 report using the specified lab report format. Just do it as homework, individually. Prepare your report using a word processor, convert it into a PDF file, and submit your work as a single PDF file on Canvas by the deadline. Begin each problem on a new sheet and answer the problems in order. Make sure to do the following: (i) Include all your work for each problem. If not, your work will not be graded and you will receive a zero. (ii) In each plot, include figure label, axes labels with units, legend, etc. (iii) Use a text font of at least 10 in all your plots and SIMULINK diagrams. (iv) Use different line types in case of multiple lines in the same plot and use a legend to clearly identify individual lines. (v) Include your MATLAB scripts and SIMULINK block diagrams in your report.

It is critical that you read following two pdf files to know how to use MATLAB and SIMULINK to perform well in this course.

Problem 1 (10 points)

The longitudinal dynamic model of a Boeing 747 aircraft flying at a speed 870 ft/s is obtained as

[u˙α˙q˙θ˙]=[0.0220.00200.0370.0840.392100.0871.560.53600010][uαqθ]+[00.0211.2090]δe\begin{bmatrix} \dot{u} \\ \dot{\alpha} \\ \dot{q} \\ \dot{\theta} \end{bmatrix} = \begin{bmatrix} -0.022 & 0.002 & 0 &-0.037\\ -0.084& -0.392 & 1 & 0\\ -0.087&-1.56&-0.536&0\\ 0&0&1&0 \end{bmatrix} \begin{bmatrix} u \\ \alpha \\ q\\ \theta \end{bmatrix} + \begin{bmatrix} 0\\-0.021\\-1.209\\0 \end{bmatrix} \delta_e

where u change in flight speed expressed as a fraction of the equilibrium speed (non- dimensional) α\alpha change in angle of attack from trim, rad q pitch rate, rad/sec θ\theta change in pitch attitude from trim, rad δe\delta_e change in elevator deflection from trim, rad

To Input this to MATLAB, type or copy and paste the following commands:

A = [-0.022 0.002 0 -0.037; -0.084 -0.392 1 0; -0.087 -1.56 -0.536 0 ; 0 0 1 0] 
B = [0; -0.021; -1.209; 0] 
C = [0 0 0 1] % Note: This selects pitch attitude as output
D = 0

Answer the following questions:

For Problem 1, a printout of your MATLAB script with the answers written in the script using comments is acceptable. Attach the plots that have been asked.

(i) Find the eigenvalues and eigenvectors of the system matrix. (1 point)

In MATLAB, use [V,S] = eig(A) to return diagonal matrix S of eigenvalues and matrix V whose columns are the corresponding right eigenvectors, so that A*V = V*S.

(ii) Obtain the transfer function between the pitch attitude and the elevator control input. (2 points)

In MATLAB, use [num,den]=ss2tf(A,B,C,D). num is the numerator of the transfer function, and den is the denominator. Then, to get a Transfer function, type G1=tf(num,den).

(iii) Obtain the poles of the transfer function of part (ii). Compare the poles of the transfer function with the eigenvalues of part (i). (2 points)

There are several ways to do this. One way is by typing zpk(G1). You get a transfer function in gain/zeroes/poles format. You can also type [z, p, k] = zpkdata(G1), which gives vectors containing zeroes (z), poles (p) and gains (k). Another is using the function pzmap(). Type ‘help pzmap’ to learn how it works.

(iv) Compute the damping ratios, natural frequencies, damped frequencies and time periods of the individual modes. (2 points)

There are several ways to get these values from eigenvalues (or alternatively from poles) of the system. Alternatively, you may use the ‘damp’ command in MATLAB. Use help to figure out how to use the damp command.

(v) Using the transfer function of part (ii), obtain frequency response plots (Bode plots). (1 point)

Type bode(G1). You get the Bode plot.

(vi) Plot the pitch attitude and pitch rate responses to a step elevator input of 0.02 rad. (2 points)

Remember that the C matrix of the system used so far implied that θ is the only output. To plot the response of both q and θ, C should be modified to: C = [0 0 1 0; 0 0 0 1]. Then, D becomes [0 ; 0]. Now, the new system can be defined by: G3 = ss(A,B,C,D) . Then, get the response to a step input of 0.02 rad using step(0.02*G3).

Problem 2 (10 points)

Consider the differential equation given by:

x¨(t)+cx˙(t)+kx(t)=r(t)\ddot{x}(t) + c\dot{x}(t) + kx(t) = r(t)

A SIMULINK model for this system is given in the following figure:

Using the figure as a guide, create your own Simulink model. For c=3c=3, k=10k=10, plot the position and the velocity over the time interval [0, 25] sec for the following cases. (i) r(t) is a unit step input, and zero initial conditions (1 point)

(ii) r(t) is zero, x(0)=1,x˙(0)=0x(0) = 1, \dot{x}(0) = 0 (1 point)

(iii) r(t) is zero, x(0)=0,x˙(0)=1x(0) = 0, \dot{x}(0) = 1 (1 point)

(iv) r(t) is a sinusoidal input of amplitude of 3 and frequency of 2 rad/s, i.e., r(t)=3sin(2t),r(t) = 3\sin(2t), and x(0)=0,x˙(0)=0x(0) = 0, \dot{x}(0) = 0. (2 points)

Note that the output of the first integrator is the velocity (x˙)(\dot{x}), and the output of the second integrator is the position (x)(x). To start a new model in SIMULINK, simply choose “File” ->“New”->”Model” from the MATLAB menu.

You can find the integrator blocks using SIMULINK Library browser under the ‘Continuous’ group and ‘Sum’ block and Gain block in the ‘Math Operations’ group. The step input is found in the “Sources” group.

To specify an initial condition, double-click on the integrator and input the initial condition in its text box. To output position and velocity, use the ’To Workspace’ block to export the data to MATLAB. The ‘To Workspace’ block can be found in the “Sinks” group. Remember to change the data type to array by double-clicking the block. The data can then be plotted against time (the tout variable) in MATLAB. Alternatively, you can also use ‘Scope’ from the “Sinks” group to create plots on the screen, and then use ‘print figure’ under ‘File’ in the ‘Scope’ menu. You can use the figure menu to include x-label, y-label, legend, etc.

Problem 3 (10 points)

The following differential equation describes the dynamics of a system:

x¨=u\ddot{x} = u

where ​xx is the system output, and uu is the control input. It is required to implement a proportional plus derivative controller (PD controller) with control law of the form

u=Kp(xxc)Kdx˙u = -K_p (x - x_c) - K_d\dot{x}

Draw a SIMULINK block diagram for this system with xcx_c as the input and xx as the output. (6 points)

Use Kp=10, and Kd=3.K_p=10,\text{ and } K_d =3.

Simulate the closed-loop system response to the following command inputs of xcx_c:

(a) pulse input of unit magnitude and pulse width of 2 seconds starting at t=1t = 1s. (2 points)

(b) step input of magnitude 2 starting at t=1t = 1 s. (2 points)

Print out your SIMULINK diagram, command input plots and response plots.

Substitute the given control law into the dynamic model of the system and use the resulting equation to create your SIMULINK model.

Problem 4 (10 points)

Develop the following block diagram in SIMULINK. The “Zero Order Hold” block accounts for the sampling time in digital control systems. (4 points)

Run the simulation for a sampling time of 0.2, 0.5, 0.8, 1.5 and 2 seconds, for a period of 30 sec. Plot these 5 graphs superimposed on the same axis. (5 points)

What is the effect of sampling time on system response? (1 point)

Problem 5 (10 points)

Refer to question 1 for the state space matrices A and B as well as the defined state vector. You will need to determine the C and D matrices as needed. Using SIMULINK, create a conventional proportional controller for a pitch angle (θ)(\theta) command system, using a gain of 2-2, i.e., δe=2(θcθ)\delta_e = -2(\theta_c - \theta). (5 points)

Simulate the response over a time range of 00 to 100s100 s. Use a pitch angle pulse command (θc)(\theta_c) input of magnitude of 0.10.1 rad, pulse width of 10s10s and the pulse starting at t=5st = 5s and ending at t=15st = 15s. Create plots of pitch angle command (θc)(\theta_c), pitch attitude (θ)(\theta) and altitude (hh) versus time. (5 points)

Note that climb rate (h˙\dot{h}) is related to α\alpha and θ\theta by the following linearized equation: h˙=V(θα)\dot{h} = V(\theta -\alpha) The flight speed V equals 870 ft/sec for this problem. With angle of attack (α\alpha) and pitch attitude (θ\theta) as outputs, obtain C and D matrices. Since there are two outputs in this problem, C is a 2x4 matrix and D is a 2x1 matrix. Use the Demux block to get α\alpha and pitch attitude θ\theta. To simulate a pulse input, combine step and negative step with a delay in starting times between them equaling the pulse width. Use the following SIMULINK diagram as a guide to creating your own SIMULINK model for this problem.

Show the SIMULINK model, C and D matrices, a plot of pitch attitude command and pitch attitude response (both on the same plot) versus time and a plot of climb rate response versus time. You should only need the following blocks in SIMULINK: step input, sum, gain, state space, demux and workspace blocks (alternatively, scopes).

Tip: If you use Scopes for your plots, you may use ‘Print to Figure’ under ‘File’ in the scope. Then, you may include x-label, y-label, legend, etc., using the figure options.

Appendix

The following material contains detailed information on how to use MATLAB for control that you might find useful throughout the semester.

Last updated

Was this helpful?