EXAMPLE 009

Cessna 172, ISA1976 integrated with Flat Earth (Euler angles).

Evolution of the aircraft after a pitch perturbation (delta doublet applied to the elevator) at t=2, a roll perturbation (delta doublet applied to the ailerons) at t=5 and a yaw perturbation (delta doublet applied to the rudder) at t=10.

Initially trimmed to a stationary, horizontal, symmetric, wings level flight.


In [1]:
# -*- coding: utf-8 -*-

Import python libreries needed.


In [2]:
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

In [3]:
import plotly.offline as offline
import plotly.graph_objs as go
import plotly.tools as tls

# notebook mode - inline
offline.init_notebook_mode()


Import PyFME classes


In [4]:
from pyfme.aircrafts import Cessna172
from pyfme.environment.environment import Environment
from pyfme.environment.atmosphere import ISA1976
from pyfme.environment.gravity import VerticalConstant
from pyfme.environment.wind import NoWind
from pyfme.models.systems import EulerFlatEarth
from pyfme.simulator import BatchSimulation
from pyfme.utils.trimmer import steady_state_flight_trimmer
from pyfme.utils.input_generator import doublet

Initialize variables


In [5]:
aircraft = Cessna172()
atmosphere = ISA1976()
gravity = VerticalConstant()
wind = NoWind()
environment = Environment(atmosphere, gravity, wind)

Initial conditions


In [6]:
TAS = 45  # m/s
h0 = 2000  # m
psi0 = 1  # rad
x0, y0 = 0, 0  # m
turn_rate = 0.0  # rad/s
gamma0 = 0.00  # rad

Define system


In [7]:
system = EulerFlatEarth(lat=0, lon=0, h=h0, psi=psi0, x_earth=x0, y_earth=y0)

not_trimmed_controls = {'delta_elevator': 0.0,
                        'delta_aileron': 0.01 * np.sign(turn_rate),
                        'delta_rudder': 0.01 * np.sign(turn_rate),
                        'delta_t': 0.5}

controls2trim = ['delta_elevator', 'delta_aileron', 'delta_rudder', 'delta_t']

In [8]:
trimmed_ac, trimmed_sys, trimmed_env, results = steady_state_flight_trimmer(
    aircraft, system, environment, TAS=TAS, controls_0=not_trimmed_controls,
    controls2trim=controls2trim, gamma=gamma0, turn_rate=turn_rate, verbose=1)


`gtol` termination condition is satisfied.
Function evaluations: 6, initial cost: 2.7864e+01, final cost 3.9052e-19, first-order optimality 1.78e-09.

Steady state flight trimmer results


In [9]:
print('delta_elevator = ',"%8.4f" % np.rad2deg(results['delta_elevator']), 'deg')
print('delta_aileron = ', "%8.4f" % np.rad2deg(results['delta_aileron']), 'deg')
print('delta_rudder = ', "%8.4f" % np.rad2deg(results['delta_rudder']), 'deg')
print('delta_t = ', "%8.4f" % results['delta_t'], '%')
print()
print('alpha = ', "%8.4f" % np.rad2deg(results['alpha']), 'deg')
print('beta = ', "%8.4f" % np.rad2deg(results['beta']), 'deg')
print()
print('u = ', "%8.4f" % results['u'], 'm/s')
print('v = ', "%8.4f" % results['v'], 'm/s')
print('w = ', "%8.4f" % results['w'], 'm/s')
print()
print('psi = ', "%8.4f" % np.rad2deg(psi0), 'deg')
print('theta = ', "%8.4f" % np.rad2deg(results['theta']), 'deg')
print('phi = ', "%8.4f" % np.rad2deg(results['phi']), 'deg')
print()
print('p =', "%8.4f" % results['p'], 'rad/s')
print('q =', "%8.4f" % results['q'], 'rad/s')
print('r =', "%8.4f" % results['r'], 'rad/s')


delta_elevator =   -3.3200 deg
delta_aileron =    0.0000 deg
delta_rudder =   -0.0000 deg
delta_t =    0.5913 %

alpha =    4.9304 deg
beta =    0.0000 deg

u =   44.8335 m/s
v =    0.0000 m/s
w =    3.8675 m/s

psi =   57.2958 deg
theta =    4.9304 deg
phi =    0.0000 deg

p =  -0.0000 rad/s
q =   0.0000 rad/s
r =   0.0000 rad/s

Initialise simulation


In [10]:
my_simulation = BatchSimulation(trimmed_ac, trimmed_sys, trimmed_env)

In [11]:
tfin = 15  # seconds
N = tfin * 100 + 1
time = np.linspace(0, tfin, N)
initial_controls = trimmed_ac.controls

In [12]:
controls = {}
for control_name, control_value in initial_controls.items():
    controls[control_name] = np.ones_like(time) * control_value
  • Elevator, aileron and rudder doublets
  • Elevator max travel: +28º/-26º
  • Aileron max travel: +20º/-15º
  • Rudder max travel: +16º/-16º

In [13]:
amplitude_elev = np.deg2rad(26)
controls['delta_elevator'] = initial_controls['delta_elevator'] + \
                             doublet(t_init=2,
                                     T=2,
                                     A=amplitude_elev,
                                     time=time,
                                     offset=np.deg2rad(0.0))

amplitude_aile = np.deg2rad(15)
controls['delta_aileron'] = doublet(t_init=5,
                                    T=2,
                                    A=amplitude_aile,
                                    time=time,
                                    offset=np.deg2rad(0.0))

amplitude_rud = np.deg2rad(32)
controls['delta_rudder'] = doublet(t_init=10,
                                   T=2,
                                   A=amplitude_rud,
                                   time=time,
                                   offset=np.deg2rad(0.0))

In [14]:
my_simulation.set_controls(time, controls)

In [15]:
par_list = ['x_earth', 'y_earth', 'height',
            'psi', 'theta', 'phi',
            'u', 'v', 'w',
            'v_north', 'v_east', 'v_down',
            'p', 'q', 'r',
            'alpha', 'beta', 'TAS',
            'F_xb', 'F_yb', 'F_zb',
            'M_xb', 'M_yb', 'M_zb']

Run Simulation


In [16]:
my_simulation.set_par_dict(par_list)
my_simulation.run_simulation()

Plot results


In [17]:
plt.style.use('ggplot')

In [18]:
# for ii in range(len(par_list) // 3):
#     three_params = par_list[3 * ii:3 * ii + 3]
#     fig, ax = plt.subplots(3, 1, sharex=True)
#     for jj, par in enumerate(three_params):
#         ax[jj].plot(time, my_simulation.par_dict[par])
#         ax[jj].set_ylabel(par)
#         ax[jj].set_xlabel('time (s)')

# fig = plt.figure()
# ax = Axes3D(fig)
# ax.plot(my_simulation.par_dict['x_earth'],
#         my_simulation.par_dict['y_earth'],
#         my_simulation.par_dict['height'])

# ax.plot(my_simulation.par_dict['x_earth'],
#         my_simulation.par_dict['y_earth'],
#         my_simulation.par_dict['height'] * 0)
# ax.set_xlabel('x_earth')
# ax.set_ylabel('y_earth')
# ax.set_zlabel('z_earth')

# plt.show()

In [19]:
camera = dict(
    up=dict(x=2, y=1, z=1),
    center=dict(x=0, y=0, z=0),
    eye=dict(x=2, y=1, z=0.7)
)
layout = {
    'autosize': False,
    'width': 900,
    'height': 500,
    'xaxis': {
        'range': [-2, 4.5],
        'zeroline': False,
    },
    'yaxis': {
        'range': [-2, 4.5]
    },
    'width': 800,
    'height': 800,
    'scene':{
        'camera': camera
    }
}
trace_with_height = go.Scatter3d(
    x=my_simulation.par_dict['x_earth'],
    y=my_simulation.par_dict['y_earth'],
    z=my_simulation.par_dict['height'],
    mode='lines')
trace_without_height = go.Scatter3d(
    x=my_simulation.par_dict['x_earth'],
    y=my_simulation.par_dict['y_earth'],
    z=my_simulation.par_dict['height'] * 0,
    mode='lines')
data = [trace_with_height, trace_without_height]
fig = {
    'data': data,
    'layout': layout,
}
offline.iplot(fig)



In [20]:
for ii in range(len(par_list) // 3):
    three_params = par_list[3 * ii:3 * ii + 3]
    trace1 = go.Scatter(
        x=time,
        y=my_simulation.par_dict[three_params[0]]
    )
    trace2 = go.Scatter(
        x=time,
        y=my_simulation.par_dict[three_params[1]]
    )
    trace3 = go.Scatter(
        x=time,
        y=my_simulation.par_dict[three_params[2]]
    )

    fig = tls.make_subplots(rows=3, cols=1, subplot_titles=('Plot 1', 'Plot 2', 'Plot 3'))
    # add subplots
    fig.append_trace(trace1, 1, 1)
    fig.append_trace(trace2, 2, 1)
    fig.append_trace(trace3, 3, 1)

    # Edit the layout
    # All of the axes properties here: https://plot.ly/python/reference/#XAxis
    fig['layout']['xaxis3'].update(title='time (s)', showgrid=True)

    # All of the axes properties here: https://plot.ly/python/reference/#YAxis
    fig['layout']['yaxis1'].update(title=three_params[0], showgrid=True)
    fig['layout']['yaxis2'].update(title=three_params[1], showgrid=True)
    fig['layout']['yaxis3'].update(title=three_params[2], showgrid=True)
    
    fig['layout'].update(title='Output subplots')

    offline.iplot(fig)


This is the format of your plot grid:
[ (1,1) x1,y1 ]
[ (2,1) x2,y2 ]
[ (3,1) x3,y3 ]

This is the format of your plot grid:
[ (1,1) x1,y1 ]
[ (2,1) x2,y2 ]
[ (3,1) x3,y3 ]

This is the format of your plot grid:
[ (1,1) x1,y1 ]
[ (2,1) x2,y2 ]
[ (3,1) x3,y3 ]

This is the format of your plot grid:
[ (1,1) x1,y1 ]
[ (2,1) x2,y2 ]
[ (3,1) x3,y3 ]

This is the format of your plot grid:
[ (1,1) x1,y1 ]
[ (2,1) x2,y2 ]
[ (3,1) x3,y3 ]

This is the format of your plot grid:
[ (1,1) x1,y1 ]
[ (2,1) x2,y2 ]
[ (3,1) x3,y3 ]

This is the format of your plot grid:
[ (1,1) x1,y1 ]
[ (2,1) x2,y2 ]
[ (3,1) x3,y3 ]

This is the format of your plot grid:
[ (1,1) x1,y1 ]
[ (2,1) x2,y2 ]
[ (3,1) x3,y3 ]