Problem Set

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.

MATLAB & SIMULINK

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

75KB

Helpful MATLAB Stuff.pdf

pdf

577KB

SIMULINK Tutorial.pdf

pdf

Problem 1 (10 points)

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

$$\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 $\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,D] = eig(A)to return diagonal matrix D of eigenvalues and matrix V whose columns are the corresponding right eigenvectors, so that A*V = V*D.

(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 following block diagram of a feedback control system:

Given

$$G(s) = \displaystyle{\frac{1}{s(s+1)(s+2)}}$$

design a lead compensator

$$C(s) = K \frac{s+a}{s+b} \qquad \quad b > a >0$$

such that the phase margin is 60$°$, the gain margin is not less than 10 dB, and the static velocity error constant Kv is 0.75$\mathrm{s}^{-1}$. Your final answer should have the desired values of K, a and b. Plot the unit step response of the closed-loop system.

Note that the velocity error constant (refer to AE3531 notes) is defined as

$K_v = \lim_{s \to 0}~sG(s)C(s) = \displaystyle\frac{Ka}{2b}$

Hence, for selected values of a and b, the value of K is determined using the specified value of the velocity error constant.

Input the transfer functions G and C to MATLAB by typing: G=tf( [1] , [1 3 2 0] ) C=3*tf( [1 0.5] , [1 1] ) Alternatively, the following commands can be used. s=tf(‘s’) G=1/(s*(s+1)*(s+2)) C=3*(s+0.5)/(s+1)

Note that arbitrary initial values of $a = 0.5$ and $b = 1$ were chosen. The value of $K = 3$ was chosen so as to meet the specified value of the velocity error constant.

Launch the Control System Designer using:

controlSystemDesigner(G,C) [all MATLAB versions]

or

rltool(G,C) [only MATLAB R2020b or older].

Under the Control System tab, select "Tuning Methods", and subsequently choose the "Bode Editor" located under the "Graphical Tuning" section. Then, in the field "Select Response to Edit", choose "LoopTransfer_C" which corresponds to the open-loop transfer function and click "Plot". You should see in front of you, the root locus, and the Bode plot. If you double-click on ‘C’ in the controllers and fixed blocks, a small window will open with the controller transfer function.

Now try adjusting the pole and zero locations of the compensator C on the Bode plot so as to meet the specified values of phase and gain margins. Note that the Control System Designer uses the controller transfer function in the form $K_1 (\tau_1 s +1)/(\tau_2 s +1)$, which is displayed on the screen. When you change the locations of the pole and zero of the controller (using the cursor), you will observe that the gain and phase margins change in the Bode plot. Adjust the pole and zero locations to achieve the required design criteria. Once the desired criteria are achieved, you will need to convert the Control System Designer compensator C which is of the form $K_1 (\tau_1 s +1)/(\tau_2 s +1)$to the form $K (s+a)/(s+b)$to get $K, a,~\text{and}~b$ (you could also check these from the "Preview" tab when "C" is selected from the "Controllers And Fixed Blocks" tab). To plot the closed-loop system response to step input (if it does not automatically open when Control System Designer is launched), go to the "Data Browser" window, and under "Responses", click on "IOTransfer_r2y". Then, right-click on it and select "plot" -> "step". Alternatively, you can go to the "Analysis" section in the top menu bar, click on "New Plot" and select "New Step". Then, in the "Select Response to Plot" field, choose "IOTransfer_r2y" and click "Plot".

Print the Bode, root locus, and time response plots.

Problem 3 (10 points)

An engineering laboratory has presented a plan to operate an Earth-orbiting satellite that is to be controlled from a ground station. A block diagram of the proposed system is shown in the figure. It takes T seconds for the signal to reach the spacecraft from the ground station and an identical delay for the return signal. The proposed ground-based controller is a proportional-derivative (PD) controller. Assume no transmission delay (i.e., T=0), and design a controller to the following specifications: (i) percent overshoot of 8% to a unit step input, and (ii) time to peak roughly at 2 sec to a unit step input. Compute the closed loop system damping ratio and natural frequency and phase margin with the controller in the loop but assuming a zero transmission time delay. Estimate the amount of allowable time delay for a stable system from the phase margin calculation.

Similar to Problem 2, launch Control System Designer, with $G(s) = 1/ s^2$ , and $C(s) = K_d(s+K_p / K_d)$, using an arbitrary value for$K_p/K_d$, for example $K_p/K_d = 1.$

Note that the Delay blocks reduce to unity at this point. In Control System Designer, choose "IOTransfer_r2y". To plot the step response, right-click on "IOTransfer_r2y" and choose "plot" -> "step". A step response will be plotted. You can right-click on the graph, and choose “Characteristics” -> “Peak Response” to read the peak value (and thus the percent overshoot) and time to peak directly off the graph.

Now adjust the value of the gain $K_d$ and $K_p/K_d$(i.e., zero location of the controller) to achieve the required design criteria. Print the root locus and the time response plot. To find the closed-loop system poles, click on the purple squares on the root locus plot and at the bottom of the window, you will see pole location, damping ratio, and frequency. Find the phase margin of the system with the controller (still assuming no transmission delay) from the “Open Loop Bode” plot of the system. The maximum allowable time delay if there is just one time delay present can be calculated by the formula: Time delay = Phase margin in rad/Gain crossover frequency. The gain crossover frequency can be read off the Bode plot. With two time delay blocks, the maximum allowable time is reduced to half of its value with one time delay block.

Problem 4 (10 points)

Consider the differential equation given by:

$$\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=3$, $k=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

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

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

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

Note that the output of the first integrator is the velocity $(\dot{x})$, and the output of the second integrator is the position $(x)$. To start a new model in SIMULINK, simply choose “File” ->“New”->”Model” from the MATLAB menu. 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 5 (10 points)

The following differential equation describes the dynamics of a system:

$$\ddot{x} = u$$

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

$$u = -K_p (x - x_c) - K_d\dot{x}$$

Draw a SIMULINK block diagram for this system with $x_c$ as the input and $x$ as the output.

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

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

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

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

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 6 (10 points)

Develop the following block diagram in SIMULINK. The “Zero Order Hold” block accounts for the sampling time in digital control systems. 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. What is the effect of sampling time on system response?

Problem 7 (20 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$, i.e., $\delta_e = -2(\theta_c - \theta)$. Simulate the response over a time range of $0$ to $100 s$. Use a pitch angle pulse command $(\theta_c)$ input of magnitude of $0.1$ rad, pulse width of $10s$ and the pulse starting at $t = 5s$ and ending at $t = 15s$. Create plots of pitch angle command $(\theta_c)$, pitch attitude $(\theta)$ and altitude ($h$) versus time.

Note that climb rate ($\dot{h}$) is related to $\alpha$ and $\theta$ by the following linearized equation: $\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.

Problem 8 (20 points)

Your task is to design a position controller for a spacecraft. The spacecraft has a mass of 4000 kg, and you can idealize it as a point mass. There is a thruster on each axis, 3 perpendicular axes total. These thrusters allow both pushing and pulling forces (this way we don’t need 6). You can directly control how much thrust is applied to each thruster, positive or negative, in units of Newtons, where positive thrust will move the spacecraft in the positive axis direction. The thrust F = 1000u, where u is the thruster input (the input to your plant). Design a position controller using feedback from your IPS (interstellar position system), which outputs in units of meters. (Assume that when you start using your controller you will reset the system to output a value of 0m on each axis.) Use a PD controller.

i) First, do this by creating a state space model of your system. You should assume that the spacecraft’s IPS system is perfect, and simulate its output through the use of the state variables that you define. Remember that the only information the IPS provides is displacement. (Output from your plant should not include velocities.) You will need to show the A, B, C, and D matrices that you create. You will put this model into SIMULINK using the SS block.

First, obtain the equations of motion. With X, Y, Z as inertial position and U, V, W as inertial velocity components, then

$$\begin{align*} \dot{X} = U \\ \dot{Y} = V \\ \dot{Z} = W \\ m\dot{U} = F_X \\ m\dot{V} = F_Y \\ m\dot{W} = F_Z \end{align*}$$

​where $m$ is the spacecraft mass and $F_X,F_Y,F_Z$are the total external force components along the inertial X, Y, Z axes, respectively. If only forces acting on the spacecraft are the forces from the thrusters, then

$$\begin{align*} F_X = 1000\,u_X \\ F_Y = 1000\,u_Y \\ F_Z = 1000\,u_Z \end{align*}$$

​where $u_X, u_Y$and $u_Z$ are the control inputs to the thrusters. The above set of equations can be rearranged as

$$\begin{bmatrix} \dot{X} \\ \dot{Y} \\ \dot{Z} \\ \dot{U} \\ \dot{V} \\ \dot{W} \end{bmatrix} = \begin{bmatrix} 0&0&0&1&0&0 \\ 0&0&0&0&1&0 \\ 0&0&0&0&0&1 \\ 0&0&0&0&0&0 \\ 0&0&0&0&0&0 \\ 0&0&0&0&0&0 \end{bmatrix} \begin{bmatrix} X\\Y\\Z\\U\\V\\W \end{bmatrix} + \begin{bmatrix} 0&0&0\\0&0&0\\0&0&0\\ {1000}/{m} &0&0\\ 0 &{1000}/{m}&0\\ 0 &0&{1000}/{m} \end{bmatrix} \begin{bmatrix} u_X\\u_Y\\u_Z \end{bmatrix}$$

​By defining the state vector $x$ as

$$x = \begin{bmatrix} X&Y&Z&U&V&W \end{bmatrix} ^T$$

​the control vector $u$ as

$$u = \begin{bmatrix} u_X&u_Y&u_Z \end{bmatrix} ^T$$

and the output vector $y$ as

$$y = \begin{bmatrix} X&Y&Z \end{bmatrix} ^T$$

obtain the A, B, C and D matrices needed in the SS block of your SIMULINK model.

ii) Create your control system in SIMULINK. (You can only use variables that the IPS provides, or functions of those variables, as feedback to the control system. For now, assume that you can take the derivative of these variables, but note that in general taking the derivative of experimental digital variables is a BAD thing!) Then use it to determine gains that work. (You don’t have to optimize your gains, we will do that during the semester in lab, just pick something that gives a decent response. You may test your model using a step input on each axis if you wish.)

Use the mux block to combine $u_X,u_Y,u_Z$into the control vector u needed as input to the SS block. Likewise, use demux block to obtain position components X, Y, Z from the output vector y of the SS block.

iii) Next, simulate your control system for the following hypothetical situation. During a trip through interstellar space, an enemy ship intercepts your path and locks on with a tractor beam. Their ship has matched your velocity and is located at (-100, -200, -100), in meters, with respect to you. The tractor beam is pulling your ship towards theirs with 5000 N of force. Your response to this attack is to reposition your ship at coordinates (10, 5, 4) to get out of their firing path. (Assume that your movement won’t change the tractor beam orientation enough to matter. Remember that in interstellar space speed is relative, so it doesn’t matter how fast you were originally moving, and that your original position at the time of attack is the origin of your coordinate system.)

You need to turn in your simulation responses for each axis, include both the command input and the response on the same plot. Also, turn in your SIMULINK model, with labels, and all state space matrices, with the state, control, and output vectors defined, as well as any gains that you used. (Do not use the PID block in SIMULINK, ever. You should use the following blocks: state space, summation, constant source, step input, output to workspace, gain, mux, and demux.)

Note that for this case, the external forces considered in (i) need to be modified to include the force from the tractor beam. The magnitude of the force from the tractor beam is given as a constant equal to 5000 N with its direction aligned with the line-of-sight at the start of the simulation, i.e., the line joining the spacecraft position coordinates at t=0, i.e., the origin, and the position of the ship located at (-100, -200, -100). Then the external forces become $F_x = 1000\,u_X + \displaystyle{\frac{-100 \times 5000}{\sqrt{ (100)^2+ (200)^2+(100)^2 } }} \\[8pt] F_y = 1000\,u_Y + \frac{-200\times 5000}{\sqrt{(100)^2+ (200)^2+(100)^2 }} \\[8pt] F_z = 1000\,u_Z + \frac{-100\times 5000}{\sqrt{(100)^2+ (200)^2+(100)^2 } }$

​The force from the tractor beam can be included in the SIMULINK as a bias to your control inputs, which are proportional to forces, by amounts that represent the effect of the tractor beam, so that you do not have to modify your original model. Your SIMULINK diagram should basically have three copies of the same control system. Include the force bias term for each loop separately (Important Note: Do not use a vector representation to solve this problem). An example of the x-position loop is shown in the following figure:

Appendix