Exam 2
Part A
(a.)
We first created the “trim” file to determine the required control (elevator and thurst setting) to achieve specific equilibrium flight by minimizing the left hand side of equations (2.5-34) which are equations for the time derivatives of the velocity, angle of attack, pitch rate, and pitch angle. The “transp” function is used in the “trim” file to determine the rate of changes for these specified controls. Once trimmed flight is achieved, the non-dimensional characteristics of this flight are used to determine the dimensional derivatives needed to set up the state space model for the trimmed flight. Once the dimensional derivatives are computed, the state space angle of attack, pitch rate, velocity, and pitch can be set up. Several characteristics of both the short period and phugoid of this state space model can be seen in Table 1.
Table 1: A table of system characteristics for both the short term period and phugoid Characteristic | Short Period | Phugoid | Eigenvalues | -0.6252 +/- 1.3188i | -0.0021 +/- 0.0780i | Time Period (s) | 1.5995 | 486.9 | Frequency (s^-1) | 1.3188 | 0.078 | Settling Time (s) | 6.398 | 1947.6 | Period of Oscillation (s) | 4.7645 | 80.5 | Natural Frequency (s^-1) | 1.4594 | 0.0781 | Damping Ratio | 0.4284 | 0.0263 |

Part B
We were tasked with plotting the Argand diagrams associated with the individual phugoid and short-period modes. An argand diagram is a scaled representation of the eigenvectors, represented in a complex plane. To do this, we took our original eigenvectors, and found the maximum value across each column. We excluded the eigenvector representation for velocity since it was significantly larger than the others, and we wanted to better visualize the other eigenvectors. The result of this can be seen below

Figure [ 1 ]. Argand Diagrams showing eigenvectors for the short-period and phugoid modes.
Some very useful conclusions can be made using these graphs. The first, is that based on the scaled output, the short-period response is most influenced by the Pitch Rate, while the phugoid response was most influenced by the pitch-angle, theta. Other conclusions can then be derived by analyzing the magnitude and direction of the other plotted vectors. A larger real component will represent a more direct response, where a larger imaginary component will represent an oscillatory motion.

Part C

The brief analysis conducted in part B was simple and could give us some information about the response, but we now wanted to plot the actual responses to conduct a more full and detailed analysis. Using our A, B, C, and D matrices as found in Part A, along with matlab’s “ss( )” function, we were able to generate the system characteristics. We then used Matlab’s “initial( )” function, to deliver the response of the system for the desired 10% initial perturbation. Below is the overall time response of the system generated by this function.

Figure [ 2 ]. Response of the system.
This shows the generic oscillatory motion across all modes and across all states, however it does not show any information regarding responses in individual phugoid and short-period modes. The function does also deliver the response in individual states; however we wanted to break those responses down even further into their modal components. To remedy this, we first broke the eigenvectors down into their real and imaginary components. We then took the real components of the eigenvectors and used them to help determine the individual short-period and phugoid responses, using modal participation. xτ=TeΛτT-1x0 (1)
Equation 1 above allows us to generate the time response for each individual state, but also allows us to isolate the individual contributions of the phugoid and short-period modes, based on the desired x0, the real component of each mode’s eigenvectors. In theory, this should be a straightforward procedure, where we would be able to plot each state’s modal response against the component state responses generated from the original “initial( )” function above. However due to an unknown error, our graphs for the individual modes were scaled incorrectly. Ultimately, this would generate the following:

Figure [ 3 ]. Example for response of individual mode contribution, for pitch-angle 'q'.
To view the responses for the other states, please view the attached Matlab script and its output, located in the appendix. As stated above, the expected results should have been such that the addition of the phugoid and short-period modal responses should sum to equal the overall total response.
Part D
The procedure to accomplish this section was very similar to the one described above, in section (c.), except that there is no 10% perturbation as an input into the system. Instead, both the short period and phugoid modes are excited separately by using their corresponding eigenvectors as initial conditions. The responses of the system from exciting the short period and phugoid modes can be seen in Figures 4 and 5. Note that the angle of attack and pitch angle have nearly identical responses, resulting in level flight.

Figure [ 4 ]: Time response of the system through exciting the short period mode

Figure [ 5 ]: Time response of the system through exciting the phugoid mode

Part E
This section asks to define the transfer functions, so that we could input a desired throttle setting and output velocity, or input elevator angle and output angle-of-attack. To do this, we first redefined our state space model. We needed to change our C matrix, since we only want a desired output of either velocity or angle of attack, and our B matrix to isolate the desired input. We then used functions like “zpkdata( )” and “ss2tf( )” to get a visual representation of the these transfer functions. This would output the poles and zeros, as well as the general transfer function as we wanted defined. Using vTδT as an example:
Transfer function: 8.188 s^3 + 10.18 s^2 + 16.31 s - 0.5668 ------------------------------------------------- s^4 + 1.254 s^3 + 2.141 s^2 + 0.01637 s + 0.01298

Z =

-0.6383 + 1.2758i
-0.6383 - 1.2758i
0.0340

P =

-0.6252 + 1.3188i
-0.6252 - 1.3188i
-0.0021 + 0.0780i
-0.0021 - 0.0780i

By analyzing the poles and zeros, we can say that the first two zeros are fairly close to the pole location. Through this assumption, we can perform a pole-zero cancellation to help simplify the system. To do this, the function “mineral( )” was used, and it gave the following results.

Zapprox =

0.0093

Papprox =

-0.0020 + 0.0782i -0.0020 - 0.0782i

To ensure that this pole-zero cancellation was appropriate, we then modeled response of both systems. This would cancel the effect of the short-period modes.

Figure [ 1 ]. A plot showing the side-by-side comparison of the approximated and the full system responses for vTδT.
Although the steady-state output changes by performing this simplification, it can be said that the overall response of the system is pretty accurate, and therefore the pole-zero cancellation was appropriate. However, for αδe, this was not the case.

Figure [ 2 ]. A side-by-side comparison of αδe, showing that pole-zero approximation is impossible.

Part F
For this step, the disturbances are no longer a change in the initial x vector. They are now a disturbance to the initial “u” vector. Technically, to determine the change in the system’s response due to δT and δe, we should just be able to multiply the transfer functions and by the disturbances. However, we were unsure of how to do this using Matlab, so we instead noticed that instead of scaling the “u” vector in the state space, we could scale the “B” matrix instead, since x=Ax+Bu. By doing this, we were able to use the state space equations to generate the desired time responses and plots. Graphs for this part can be visualized in the Appendix.
Part G
We attempted to do this step inside of part (c.) first. VT was broken up into x and z components over every time interval through the multiplication of cosine of angle of attack and sine of angle of attack terms. We did not have enough time to complete this step, but we anticipated finding the flight path of the vehicle through finding the area under the velocity curve. Since the velocity components were known, they could be multiplied by dt for every time interval and then summed to achieve the flight path data in the inertial frame. The flight path data in the relative flight frame could be extracted from the inertial flight data through the subtraction of the equilibrium flight path (i.e. VTe*dt). The flight paths for part (c.) can be seen in Figure 3, but since the code could not be finished, the paths are not correct.

Figure 3: Flight paths for part (c.) in the inertial and relative frames, which are not correct yet