User:Egm6341.s11.team4.KTH/HW7

=Problem 7.2: Optimal Solution of Aircraft Bunt Maneuver =

Refer to lecture slides [[media:nm1.s11.mtg40.djvu|40-1]], [[media:egm6341.s10.mtg40.djvu|40-3]] and [[media:egm6341.s10.mtg41.djvu|41-1]] for the problem statement.

Given
The parameters for this problem are given in both the problem statement in the lectures as well as Subchan and Zbikowski's paper Computational optimal control of the terminal bunt manoeuvre—Part 1: Minimum altitude case.

From Subchan and Zbikowski's paper we have the following:

Equations of Motion:


 * {| style="width:70%" border="0" align="center"

(2.1)
 * $$\displaystyle \dot \gamma = \frac{T - D}{mV}sin(\alpha) + \frac{L}{mV}cos(\alpha) - \frac{g cos(\gamma)}{V}$$
 * 
 * 
 * }
 * }


 * {| style="width:70%" border="0" align="center"

(2.2)
 * $$\displaystyle \dot V = \frac{T - D}{m}cos(\alpha) - \frac{L}{m}sin(\alpha) - g sin(\gamma)$$
 * 
 * 
 * }
 * }


 * {| style="width:70%" border="0" align="center"

(2.3)
 * $$\displaystyle \dot x = V cos(\gamma)$$
 * 
 * 
 * }
 * }


 * {| style="width:70%" border="0" align="center"

(2.4)
 * $$\displaystyle \dot y = V sin(\gamma)$$
 * 
 * 
 * }
 * }

And,


 * {| style="width:70%" border="0" align="center"

\begin{bmatrix} \gamma \\ v \\ x \\ y \end{bmatrix} $$
 * $$\mathbf{z} =
 * $$\mathbf{z} =
 * }
 * }

Where &gamma;, V, x and y are state variables defined such that &gamma; is the flight path angle, V is the speed, x is the horizontal position and y is the altitude.

Also from Subchan and Zbikowski we have:


 * {| style="width:70%" border="0" align="center"


 * $$\displaystyle D(h,V,\alpha) = \frac{1}{2}C_d \rho V^2 S_{ref}, \,\,\, C_d = A_1 \alpha^2 + A_2\alpha + A_3$$
 * }
 * }
 * }


 * {| style="width:70%" border="0" align="center"


 * $$\displaystyle L(h,V,\alpha) = \frac{1}{2}C_l \rho V^2 S_{ref}, \,\,\, C_l = B_1\alpha + B_2$$
 * }
 * }
 * }


 * {| style="width:70%" border="0" align="center"


 * $$\displaystyle \rho = C_1y^2 + C_2y + C_3$$
 * }
 * }
 * }

Where,


 * {| style="width:70%" border="0"

\begin{align} & m = 1005 \; kg \\ & g = 9.81 \; m/s^2 \\ & S_{ref} = 0.3376 \; m^2 \\ & A_{1} = -1.9431 \\ & A_{2} = -0.1499 \\ & A_{3} = 0.2359 \\ & B_{1} = 21.9 \\ & B_{2} = 0 \\ & C_{1} = 3.312 \times 10^{-9} \; kg/m^5 \\ & C_{2} = -1.142 \times 10^{-4} \; kg/m^4 \\ & C_{3} = 1.224 \; kg/m^3 \end{align} $$
 * $$\displaystyle
 * $$\displaystyle
 * }
 * }

The initial conditions are given as:


 * {| style="width:70%" border="0"

\begin{align} & \gamma_0 = 0 \, rad\\ & V_0 = 272 \, m/s \\ & x_0 = 0 \, m \\ & y_0 = 30 \, m \\ \end{align} $$
 * $$\displaystyle
 * $$\displaystyle
 * }
 * }

From the problem statement in the lecture notes we get the profiles for the two control inputs thrust, T(t) and angle of attack &alpha;(t):


 * {| style="width:70%" border="0"

T(t) = \begin{cases} 6000 \; N, & \mbox{for }0 \leq t < 27 \; \mbox{and} \; 33 \leq t \leq 40 \\ 1000 \; N, & \mbox{for }27 \leq t < 33 \end{cases} $$     (2.5)
 * $$\displaystyle
 * $$\displaystyle
 * 


 * }
 * }

For part 1:


 * {| style="width:70%" border="0"

\alpha(t) = \begin{cases} 0.03 \; rad, & \mbox{for }0 \leq t < 21 \\ \frac{0.13-0.09}{27-21}(t-21) + 0.09 \; rad, & \mbox{for }21 \leq t < 27 \\ \frac{-0.2+0.13}{33-27}(t-27) - 0.13 \; rad, & \mbox{for }27 \leq t < 33 \\ \frac{-0.13+0.2}{40-33}(t-33) - 0.2\; rad, & \mbox{for }33 \leq t \leq 40 \end{cases} $$     (2.6)
 * $$\displaystyle
 * $$\displaystyle
 * 
 * }
 * }

For part 2:


 * {| style="width:70%" border="0"

\alpha(t) = \begin{cases} 0.03 \; rad, & \mbox{for }0 \leq t < 21 \\ \frac{0.13-0.03}{27-21}(t-21) + 0.03 \; rad, & \mbox{for }21 \leq t < 27 \\ \frac{-0.2+0.13}{33-27}(t-27) - 0.13 \; rad, & \mbox{for }27 \leq t < 33 \\ - 0.2\; rad, & \mbox{for }33 \leq t \leq 40 \end{cases} $$     (2.7)
 * $$\displaystyle
 * $$\displaystyle
 * 
 * }
 * }

Objective
a.) Solve for z(t), for t contained on [0,40] such that successive numerical solutions are of the order 10-6, then plot the state variables as a function of time and compare them to the solution using the MatLab function "ode45".  Also compute


 * {| style="width:70%" border="0" align="center"


 * $$\displaystyle J_1 = \int_0^{t_f} y(t)dt \,\,\, s.t. \,\, y(t_f) = 0$$
 * }
 * }
 * }

Use given alpha for part 1.

b.) Solve for z(t), for t contained on [0,40] such that successive numerical solutions are of the order 10-6, then plot the state variables as a function of time and compare them to the solution using the MatLab function "ode45".


 * {| style="width:70%" border="0" align="center"


 * $$\displaystyle J_2 = \int_0^{t_f} y(t)dt \,\,\, s.t. \,\, y(t_f) = 0$$
 * }
 * }
 * }

Use given alpha for part 2.

Solution
Both parts of the problem will be solved using the same technique. The only difference will be the alpha function input into the Matlab code.

From equation (1) on Lecture Slide [[media:nm1.s11.mtg38.djvu|38-3]] we can derive:


 * {| style="width:70%" border="0" align="center"

$$     (2.8)
 * $$\displaystyle \mathbf{z}_{i+1} = \mathbf{z}_{i} + \frac{h/2}{3}\left[\mathbf{f}(\mathbf{z}_{i})+4\mathbf{f}(\mathbf{g}(\mathbf{z}_{i},\mathbf{z}_{i+1}))+\mathbf{f}(\mathbf{z}_{i+1})\right]
 * $$\displaystyle \mathbf{z}_{i+1} = \mathbf{z}_{i} + \frac{h/2}{3}\left[\mathbf{f}(\mathbf{z}_{i})+4\mathbf{f}(\mathbf{g}(\mathbf{z}_{i},\mathbf{z}_{i+1}))+\mathbf{f}(\mathbf{z}_{i+1})\right]
 * 
 * }
 * }

Where:


 * {| style="width:70%" border="0" align="center"

= \frac{1}{2}\left(\mathbf{z}_{i} + \mathbf{z}_{i+1}\right) + \frac{h}{8}\left(\mathbf{f}(\mathbf{z}_{i}) - \mathbf{f}(\mathbf{z}_{i+1})\right) $$     (2.9)
 * $$\displaystyle \mathbf{g}(\mathbf{z}_{i},\mathbf{z}_{i+1}) = \mathbf{z}_{i+1/2}
 * $$\displaystyle \mathbf{g}(\mathbf{z}_{i},\mathbf{z}_{i+1}) = \mathbf{z}_{i+1/2}
 * 
 * }
 * }

F(zi+1) can then be defined by subtracting the right side of eqn. (2.8) from the left side of eqn. (2.8).


 * {| style="width:70%" border="0" align="center"

$$     (2.10)
 * $$\displaystyle \mathbf{F}(\mathbf{z}_{i+1}) = \mathbf{z}_{i+1} - \mathbf{z}_{i} + \frac{h/2}{3}\left[\mathbf{f}(\mathbf{z}_{i})+4\mathbf{f}(\mathbf{g}(\mathbf{z}_{i},\mathbf{z}_{i+1}))+\mathbf{f}(\mathbf{z}_{i+1})\right] = \mathbf{0}
 * $$\displaystyle \mathbf{F}(\mathbf{z}_{i+1}) = \mathbf{z}_{i+1} - \mathbf{z}_{i} + \frac{h/2}{3}\left[\mathbf{f}(\mathbf{z}_{i})+4\mathbf{f}(\mathbf{g}(\mathbf{z}_{i},\mathbf{z}_{i+1}))+\mathbf{f}(\mathbf{z}_{i+1})\right] = \mathbf{0}
 * 
 * }
 * }

From Lecture Slide [[media:nm1.s11.mtg39.djvu|39-1]] we get the following:


 * {| style="width:70%" border="0" align="center"

$$     (2.11)
 * $$\displaystyle \mathbf{z}_{i+1}^{(k+1)} = \mathbf{z}_{i+1}^{(k)} - \left[\frac{\partial \mathbf{F}(\mathbf{z}_{i+1}^{(k)})}{\partial \mathbf{z}}\right]^{-1}\mathbf{F}(\mathbf{z}_{i+1}^{(k)})
 * $$\displaystyle \mathbf{z}_{i+1}^{(k+1)} = \mathbf{z}_{i+1}^{(k)} - \left[\frac{\partial \mathbf{F}(\mathbf{z}_{i+1}^{(k)})}{\partial \mathbf{z}}\right]^{-1}\mathbf{F}(\mathbf{z}_{i+1}^{(k)})
 * 
 * }
 * }

In order to solve the problem we need to calculate the Jacobian as follows:


 * {| style="width:70%" border="0" align="center"

$$     (2.12) The values in the Jacobian can be seen in the MatLab code below.  "Main Script" 
 * $$\displaystyle \frac{\partial \mathbf{F}}{\partial \mathbf{z}} = \left[\frac{\partial F_i}{\partial z_j}\right]
 * $$\displaystyle \frac{\partial \mathbf{F}}{\partial \mathbf{z}} = \left[\frac{\partial F_i}{\partial z_j}\right]
 * 
 * }
 * }

 Function "zdot" 

 Function "dfdz" 

Function "nonlinode1" and "nonlinode2" to be used with the ode45 function.  Function "nonlinode1"   Function "nonlinode2" 

a.) Part 1: Initial values for alpha










b.) Part 2: Modified values for alpha










c.) Compute $$\ J_1$$, $$\ J_2$$ and compare.
$$ \ J_1 = 16803.19 $$

$$ \ J_2 = 9024.61 $$

$$ \ J_1 > J_2 $$

So, using modified input makes better results.