From 79d893d3b62ac2f487d3d1e1ed09efa06d5d0093 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?=C3=81lex=20S=C3=A1ez?= Date: Sat, 30 Jan 2016 19:33:04 +0100 Subject: [PATCH 01/17] Steady-state trimmer in progress. Created main function and constrain functions. Minimizer is still to be implemented. --- src/pyfme/utils/trimmer.py | 123 +++++++++++++++++++++++++++++++++++++ 1 file changed, 123 insertions(+) create mode 100644 src/pyfme/utils/trimmer.py diff --git a/src/pyfme/utils/trimmer.py b/src/pyfme/utils/trimmer.py new file mode 100644 index 0000000..8666535 --- /dev/null +++ b/src/pyfme/utils/trimmer.py @@ -0,0 +1,123 @@ +# -*- coding: utf-8 -*- +""" +Python Flight Mechanics Engine (PyFME). +Copyright (c) AeroPython Development Team. +Distributed under the terms of the MIT License. + +Trimmer +------- +This module solves the problem of calculating the values of the state and +control vectors that satisfy the state equations of the aircraft at the initial +condition (or any other given condition. This cannot be done analytically +because of the very complex functional dependence of the erodynamic data. +Instead, it must be done with a numerical algorithm which iteratively adjusts +the independent variables until some solution criterion is met. +""" + +from math import sqrt, sin, cos, tan, atan + + +def steady_state_flight_trim(aircraft, h, TAS, gamma=0, turn_rate=0): + """Finds a combination of values of the state and control variables that + correspond to a steady-state flight condition. Steady-state aircraft flight + can be defined as a condition in which all of the motion variables are + constant or zero. That is, the linear and angular velocity components are + constant (or zero), and all acceleration components are zero. + + Parameters + ---------- + aircraft : aircraft class + Aircraft class with methods get_forces, get_moments + h : float + Geopotential altitude for ISA (m). + TAS : float + True Air Speed (m/s). + gamma : float, optional + Flight path angle (rad). + turn_rate : float, optional + Turn rate, d(psi)/dt (rad/s). + + Returns + ------- + lin_vel : + + ang_vel : + + theta : + + phi : + + alpha : + + beta : + + control_vector : + + + Notes + ----- + See section 3.4 in [1] for the algorithm description. + See section 2.5 in [1] for the definition of steady-state flight condition. + + References + ---------- + .. [1] Stevens, BL and Lewis, FL, "Aircraft Control and Simulation", + Wiley-lnterscience. + """ + + pass + + return + + +def turn_coord_cons1(turn_rate, alpha, beta, TAS, gamma=0): + """Calculates phi for coordinated turn + """ + + g0 = 9.8 + G = turn_rate * TAS / g0 + + if gamma == 0: + phi = G * cos(beta) / (cos(alpha) - G * sin(alpha) * sin(beta)) + phi = atan(phi) + + else: + a = 1 - G * tan(alpha) * sin(beta) + b = sin(gamma) / cos(beta) + c = 1 + G**2 * cos(beta)**2 + + sq = sqrt(c * (1 - b**2) + G**2 * sin(beta)**2) + var1 = b**2 * (1 + c * tan(alpha)**2) + + temp = (a - b)**2 + b * tan(alpha) * sq / a**2 - b**2 * var1 + + phi = atan(G * cos(beta) / cos(alpha) * temp) + + return phi + + +def turn_coord_cons2(turn_rate, alpha, TAS): + """Calculates phi for coordinated given that gamma is equal to cero and + beta is small (beta << 1). + """ + + g0 = 9.8 + G = turn_rate * TAS / g0 + + phi = G / cos(alpha) + + return phi + + +def rate_of_climb_cons(gamma, alpha, beta, phi): + """ + """ + a = cos(alpha) * cos(beta) + b = sin(phi) * sin(beta) + cos(phi) * sin(alpha) * cos(beta) + + sq = sqrt(a**2 - sin(gamma)**2 + b**2) + + theta = (a * b + sin(gamma) * sq) / (a**2 - sin(gamma)**2) + theta = atan(theta) + + return theta From 7a59afba45845f72760c91511fd2c01fb8b74802 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?=C3=81lex=20S=C3=A1ez?= Date: Sun, 31 Jan 2016 00:15:22 +0100 Subject: [PATCH 02/17] Added funciton to be optimized (still to be tested) The evaluation of forces and moments suggests a change in the aircraft model. --- src/pyfme/utils/trimmer.py | 47 +++++++++++++++++++++++++++++++++++--- 1 file changed, 44 insertions(+), 3 deletions(-) diff --git a/src/pyfme/utils/trimmer.py b/src/pyfme/utils/trimmer.py index 8666535..8b120aa 100644 --- a/src/pyfme/utils/trimmer.py +++ b/src/pyfme/utils/trimmer.py @@ -15,6 +15,9 @@ """ from math import sqrt, sin, cos, tan, atan +import numpy as np + +from pyfme.utils.coordinates import wind2body def steady_state_flight_trim(aircraft, h, TAS, gamma=0, turn_rate=0): @@ -71,7 +74,7 @@ def steady_state_flight_trim(aircraft, h, TAS, gamma=0, turn_rate=0): def turn_coord_cons1(turn_rate, alpha, beta, TAS, gamma=0): - """Calculates phi for coordinated turn + """Calculates phi for coordinated turn. """ g0 = 9.8 @@ -97,8 +100,8 @@ def turn_coord_cons1(turn_rate, alpha, beta, TAS, gamma=0): def turn_coord_cons2(turn_rate, alpha, TAS): - """Calculates phi for coordinated given that gamma is equal to cero and - beta is small (beta << 1). + """Calculates phi for coordinated turn given that gamma is equal to cero + and beta is small (beta << 1). """ g0 = 9.8 @@ -121,3 +124,41 @@ def rate_of_climb_cons(gamma, alpha, beta, phi): theta = atan(theta) return theta + + +def func(trimmed_params, h, TAS, gamma, turn_rate, aircraft, dynamic_eqs): + """Function to optimize + """ + + # Fixme: if turn_rate != 0, p, q, r must also be calculated + + alpha = trimmed_params[0] + beta = trimmed_params[1] + + delta_e = trimmed_params[2] + delta_ail = trimmed_params[3] + delta_r = trimmed_params[4] + delta_t = trimmed_params[5] + + phi = turn_coord_cons1(turn_rate, alpha, beta, TAS, gamma) + theta = rate_of_climb_cons(gamma, alpha, beta, phi) + + # w = turn_rate * k_h + # k_h = sin(theta) i_b + sin(phi) * cos(theta) j_b + cos(theta) * sin(phi) + # w = p * i_b + q * j_b + r * k_b + p = - turn_rate * sin(theta) + q = turn_rate * sin(phi) * cos(theta) + r = turn_rate * cos(theta) * sin(phi) + + ang_vel = np.array([p, q, r]) + + lin_vel = wind2body((TAS, 0, 0), ) + + # FIXME: This implies to change the aircraft model... + forces, moments = aircraft.get_forces_and_moments() # Parametets...? + mass = aircraft.mass + inertia = aircraft.inertia + + vel = np.concatenate(lin_vel, ang_vel) + + return dynamic_eqs(0, vel, mass, inertia, forces, moments) From bd1deacce3def75b0c25942df380acf86fa6f542 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?=C3=81lex=20S=C3=A1ez?= Date: Sun, 31 Jan 2016 15:02:00 +0100 Subject: [PATCH 03/17] Finished function to be optimized during the trimming processs. * A new function was introduced in cessna_310 to get total forces & moments. * Minor corrections in cessna_310. * Only trimming main function is still to be implemented. TESTS ARE PENDING. --- src/pyfme/aircrafts/cessna_310.py | 57 +++++++++++++++++++++++-------- src/pyfme/utils/trimmer.py | 19 +++++++---- 2 files changed, 55 insertions(+), 21 deletions(-) diff --git a/src/pyfme/aircrafts/cessna_310.py b/src/pyfme/aircrafts/cessna_310.py index ae7cfb0..0392696 100644 --- a/src/pyfme/aircrafts/cessna_310.py +++ b/src/pyfme/aircrafts/cessna_310.py @@ -13,9 +13,10 @@ import numpy as np from pyfme.utils.anemometry import calculate_dynamic_pressure +from pyfme.utils.coordinates import hor2body -def geometric_Data(): +def geometric_data(): """ Provides the value of some geometric data. @@ -42,7 +43,7 @@ def geometric_Data(): return Sw, c, span -def mass_and_Inertial_Data(): +def mass_and_inertial_data(): """ Provides the value of some mass and inertial data. @@ -78,7 +79,7 @@ def mass_and_Inertial_Data(): return mass, inertia -def long_Aero_Coefficients(): +def long_aero_coefficients(): """Assigns the value of the coefficients of stability in cruise conditions and order them in a matrix. @@ -141,7 +142,7 @@ def long_Aero_Coefficients(): return long_coef_matrix -def lat_Aero_Coefficients(): +def lat_aero_coefficients(): """Assigns the value of the coefficients of stability in cruise conditions and order them in a matrix. @@ -243,8 +244,8 @@ def get_aerodynamic_forces(TAS, rho, alpha, beta, delta_e, ih, delta_ail, chapter 3 and 4 """ - long_coef_matrix = long_Aero_Coefficients() - lat_coef_matrix = lat_Aero_Coefficients() + long_coef_matrix = long_aero_coefficients() + lat_coef_matrix = lat_aero_coefficients() CL_0, CL_a, CL_de, CL_dih = long_coef_matrix[0, :] CD_0, CD_a, CD_de, CD_dih = long_coef_matrix[1, :] @@ -254,7 +255,7 @@ def get_aerodynamic_forces(TAS, rho, alpha, beta, delta_e, ih, delta_ail, CD_full = CD_0 + CD_a * alpha + CD_de * delta_e + CD_dih * ih CY_full = CY_b * beta + CY_da * delta_ail + CY_dr * delta_r - Sw = geometric_Data()[0] + Sw = geometric_data()[0] aerodynamic_forces = calculate_dynamic_pressure(rho, TAS) *\ Sw * np.array([-CD_full, CY_full, -CL_full]) # N @@ -301,8 +302,8 @@ def get_aerodynamic_moments(TAS, rho, alpha, beta, delta_e, ih, delta_ail, chapter 3 and 4 """ - long_coef_matrix = long_Aero_Coefficients() - lat_coef_matrix = lat_Aero_Coefficients() + long_coef_matrix = long_aero_coefficients() + lat_coef_matrix = lat_aero_coefficients() Cm_0, Cm_a, Cm_de, Cm_dih = long_coef_matrix[2, :] Cl_b, Cl_da, Cl_dr = lat_coef_matrix[1, :] @@ -312,9 +313,9 @@ def get_aerodynamic_moments(TAS, rho, alpha, beta, delta_e, ih, delta_ail, Cl_full = Cl_b * beta + Cl_da * delta_ail + Cl_dr * delta_r Cn_full = Cn_b * beta + Cn_da * delta_ail + Cn_dr * delta_r - span = geometric_Data()[2] - c = geometric_Data()[1] - Sw = geometric_Data()[0] + span = geometric_data()[2] + c = geometric_data()[1] + Sw = geometric_data()[0] aerodynamic_moments = calculate_dynamic_pressure(rho, TAS) * Sw\ * np.array([Cl_full * span, Cm_full * c, Cn_full * span]) @@ -337,7 +338,7 @@ def get_engine_force(delta_t): ------- engine_force : float - Trust (N) + Thrust (N) References ---------- @@ -349,11 +350,37 @@ def get_engine_force(delta_t): q_cruise = 91.2 * 47.880172 # Pa - Sw = geometric_Data()[0] + Sw = geometric_data()[0] engine_force = Ct * Sw * q_cruise # N else: raise ValueError('delta_t must be between 0 and 1') - return engine_force \ No newline at end of file + return engine_force + + +def get_forces_and_moments(TAS, rho, alpha, beta, delta_e, ih, delta_ail, + delta_r, delta_t, attitude): + """ + """ + # As stability axes are coincident with wind axes at the moment of + # linearization (with alpha=0 and beta=0), stability axes are parallel to + # body axes. + aero_forces = get_aerodynamic_forces(TAS, rho, alpha, beta, delta_e, 0, + delta_ail, delta_r) + thrust = get_engine_force(delta_t) + # Assuming that all the engine thrust is prodecued in the x_body direction + engine_force = np.array([thrust, 0, 0]) + + g0 = 9.81 # m/s² + theta, phi, psi = attitude + gravity_force = hor2body((0, 0, g0), theta, phi, psi) + + forces = aero_forces + engine_force + gravity_force + + # It is assumed that the engine moments are zero. + moments = get_aerodynamic_moments(TAS, rho, alpha, beta, delta_e, ih, + delta_ail, delta_r) + + return forces, moments diff --git a/src/pyfme/utils/trimmer.py b/src/pyfme/utils/trimmer.py index 8b120aa..54db6dd 100644 --- a/src/pyfme/utils/trimmer.py +++ b/src/pyfme/utils/trimmer.py @@ -18,6 +18,7 @@ import numpy as np from pyfme.utils.coordinates import wind2body +from pyfme.environment.isa import atm def steady_state_flight_trim(aircraft, h, TAS, gamma=0, turn_rate=0): @@ -77,7 +78,7 @@ def turn_coord_cons1(turn_rate, alpha, beta, TAS, gamma=0): """Calculates phi for coordinated turn. """ - g0 = 9.8 + g0 = 9.81 G = turn_rate * TAS / g0 if gamma == 0: @@ -104,7 +105,7 @@ def turn_coord_cons2(turn_rate, alpha, TAS): and beta is small (beta << 1). """ - g0 = 9.8 + g0 = 9.81 G = turn_rate * TAS / g0 phi = G / cos(alpha) @@ -154,10 +155,16 @@ def func(trimmed_params, h, TAS, gamma, turn_rate, aircraft, dynamic_eqs): lin_vel = wind2body((TAS, 0, 0), ) - # FIXME: This implies to change the aircraft model... - forces, moments = aircraft.get_forces_and_moments() # Parametets...? - mass = aircraft.mass - inertia = aircraft.inertia + # FIXME: This implied some changes in the aircraft model. + # psi angle does not influence the attitude of the aircraft for gravity + # force projection. So it is set to 0. + attitude = np.array([theta, phi, 0]) + _, _, rho = atm(h) + forces, moments = aircraft.get_forces_and_moments(TAS, rho, alpha, beta, + delta_e, 0, delta_ail, + delta_r, delta_t, + attitude) + mass, inertia = aircraft.mass_and_inertial_data() vel = np.concatenate(lin_vel, ang_vel) From aaa3f3962c506471e933b77e63df4ee8e3a0492a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?=C3=81lex=20S=C3=A1ez?= Date: Sun, 31 Jan 2016 18:16:10 +0100 Subject: [PATCH 04/17] Changed atm doc in isa to include also sound velocity as an output --- src/pyfme/environment/isa.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/pyfme/environment/isa.py b/src/pyfme/environment/isa.py index b866bbc..1c990bf 100644 --- a/src/pyfme/environment/isa.py +++ b/src/pyfme/environment/isa.py @@ -51,7 +51,9 @@ def atm(h): p : float Pressure (Pa). rho : float - Density (kg/m³) + Density (kg/m³). + a : float + Sound speed for given altitude (m/s). Raises ------ From 0618f415f2c8660ea8d4318b10d125c5a691f3f2 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?=C3=81lex=20S=C3=A1ez?= Date: Sun, 31 Jan 2016 18:19:43 +0100 Subject: [PATCH 05/17] First attempt to trim the aircraft: DOES NOT SEEM TO CONVERGE TO THE RIGHT SOLUTION * Added example: trim_example_001.py. * Another minimization algorithm should be tried... this one in converginf to another solution. Maybe is a good idea to choose one that accepts constrains (ie. delta_t <= 1). --- examples/trim_example_001.py | 34 ++++++++++++++++ src/pyfme/utils/trimmer.py | 75 +++++++++++++++++++++++++++++------- 2 files changed, 95 insertions(+), 14 deletions(-) create mode 100644 examples/trim_example_001.py diff --git a/examples/trim_example_001.py b/examples/trim_example_001.py new file mode 100644 index 0000000..cf15561 --- /dev/null +++ b/examples/trim_example_001.py @@ -0,0 +1,34 @@ +# -*- coding: utf-8 -*- +""" +Python Flight Mechanics Engine (PyFME). +Copyright (c) AeroPython Development Team. +Distributed under the terms of the MIT License. + + + +""" +from pyfme.aircrafts import cessna_310 +from pyfme.utils.trimmer import steady_state_flight_trim +from pyfme.utils.coordinates import wind2body + +TAS = 120 # m/s +h = 1000 + +results = steady_state_flight_trim(cessna_310, h, TAS, gamma=0, turn_rate=0) +# lin_vel, ang_vel, theta, phi, alpha, beta, control_vector + +alpha = results[3] +beta = results[4] +lin_vel = wind2body((TAS, 0, 0), alpha, beta) + +print('') +print('Results:') +print("lin_vel = {}".format(lin_vel)) +print("ang_vel = {}".format(results[0])) +print("theta = {}".format(results[1])) +print("phi = {}".format(results[2])) +print("alpha = {}".format(results[3])) +print("beta = {}".format(results[4])) +print("control = {}".format(results[5])) +print('') + diff --git a/src/pyfme/utils/trimmer.py b/src/pyfme/utils/trimmer.py index 54db6dd..a740a87 100644 --- a/src/pyfme/utils/trimmer.py +++ b/src/pyfme/utils/trimmer.py @@ -13,15 +13,17 @@ Instead, it must be done with a numerical algorithm which iteratively adjusts the independent variables until some solution criterion is met. """ - -from math import sqrt, sin, cos, tan, atan +from warnings import warn import numpy as np +from math import sqrt, sin, cos, tan, atan +from scipy.optimize import fmin from pyfme.utils.coordinates import wind2body from pyfme.environment.isa import atm -def steady_state_flight_trim(aircraft, h, TAS, gamma=0, turn_rate=0): +def steady_state_flight_trim(aircraft, h, TAS, gamma=0, turn_rate=0, + dyn_eqs=None): """Finds a combination of values of the state and control variables that correspond to a steady-state flight condition. Steady-state aircraft flight can be defined as a condition in which all of the motion variables are @@ -69,9 +71,51 @@ def steady_state_flight_trim(aircraft, h, TAS, gamma=0, turn_rate=0): Wiley-lnterscience. """ - pass + if dyn_eqs is None: + from pyfme.models.euler_flat_earth import linear_and_angular_momentum_eqs + dynamic_eqs = linear_and_angular_momentum_eqs + + # TODO: try to look for a good inizialization + initial_guess = np.zeros([6]) + 0.01 + + args = (h, TAS, gamma, turn_rate, aircraft, dynamic_eqs) + + results = fmin(func, x0=initial_guess, args=args, full_output=True, + disp=False) + trimmed_params, fopt = results[0:2] + + if abs(fopt) > 1e-3: + warn("Trim process did not converge. fopt={}".format(fopt), + RuntimeWarning) + + alpha = trimmed_params[0] + beta = trimmed_params[1] + + delta_e = trimmed_params[2] + delta_ail = trimmed_params[3] + delta_r = trimmed_params[4] + delta_t = trimmed_params[5] + + # What happens if we have two engines and all that kind of things...? + control_vector = delta_e, delta_ail, delta_r, delta_t + + if abs(turn_rate) < 1e-8: + phi = 0 + else: + phi = turn_coord_cons1(turn_rate, alpha, beta, TAS, gamma) - return + theta = rate_of_climb_cons(gamma, alpha, beta, phi) + + # w = turn_rate * k_h + # k_h = sin(theta) i_b + sin(phi) * cos(theta) j_b + cos(theta) * sin(phi) + # w = p * i_b + q * j_b + r * k_b + p = - turn_rate * sin(theta) + q = turn_rate * sin(phi) * cos(theta) + r = turn_rate * cos(theta) * sin(phi) + + ang_vel = np.array([p, q, r]) + + return ang_vel, theta, phi, alpha, beta, control_vector def turn_coord_cons1(turn_rate, alpha, beta, TAS, gamma=0): @@ -81,7 +125,7 @@ def turn_coord_cons1(turn_rate, alpha, beta, TAS, gamma=0): g0 = 9.81 G = turn_rate * TAS / g0 - if gamma == 0: + if abs(gamma) < 1e-8: phi = G * cos(beta) / (cos(alpha) - G * sin(alpha) * sin(beta)) phi = atan(phi) @@ -130,9 +174,6 @@ def rate_of_climb_cons(gamma, alpha, beta, phi): def func(trimmed_params, h, TAS, gamma, turn_rate, aircraft, dynamic_eqs): """Function to optimize """ - - # Fixme: if turn_rate != 0, p, q, r must also be calculated - alpha = trimmed_params[0] beta = trimmed_params[1] @@ -141,7 +182,11 @@ def func(trimmed_params, h, TAS, gamma, turn_rate, aircraft, dynamic_eqs): delta_r = trimmed_params[4] delta_t = trimmed_params[5] - phi = turn_coord_cons1(turn_rate, alpha, beta, TAS, gamma) + if abs(turn_rate) < 1e-8: + phi = 0 + else: + phi = turn_coord_cons1(turn_rate, alpha, beta, TAS, gamma) + theta = rate_of_climb_cons(gamma, alpha, beta, phi) # w = turn_rate * k_h @@ -153,19 +198,21 @@ def func(trimmed_params, h, TAS, gamma, turn_rate, aircraft, dynamic_eqs): ang_vel = np.array([p, q, r]) - lin_vel = wind2body((TAS, 0, 0), ) + lin_vel = wind2body((TAS, 0, 0), alpha, beta) # FIXME: This implied some changes in the aircraft model. # psi angle does not influence the attitude of the aircraft for gravity # force projection. So it is set to 0. attitude = np.array([theta, phi, 0]) - _, _, rho = atm(h) + _, _, rho, _ = atm(h) forces, moments = aircraft.get_forces_and_moments(TAS, rho, alpha, beta, delta_e, 0, delta_ail, delta_r, delta_t, attitude) mass, inertia = aircraft.mass_and_inertial_data() - vel = np.concatenate(lin_vel, ang_vel) + vel = np.concatenate((lin_vel[:], ang_vel[:])) + + output = dynamic_eqs(0, vel, mass, inertia, forces, moments) - return dynamic_eqs(0, vel, mass, inertia, forces, moments) + return np.sum(np.dot(output, output)) From 0ba7045a9065f4bc110afde8d18cf56d950c922b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?=C3=81lex=20S=C3=A1ez?= Date: Sat, 6 Feb 2016 18:00:24 +0100 Subject: [PATCH 06/17] Using least_squares from scipy v0.17 provides a better trim --- examples/trim_example_001.py | 5 ++--- src/pyfme/utils/trimmer.py | 21 +++++++++++++-------- 2 files changed, 15 insertions(+), 11 deletions(-) diff --git a/examples/trim_example_001.py b/examples/trim_example_001.py index cf15561..3193cff 100644 --- a/examples/trim_example_001.py +++ b/examples/trim_example_001.py @@ -7,11 +7,11 @@ """ -from pyfme.aircrafts import cessna_310 +from pyfme.aircrafts import cessna_310 from pyfme.utils.trimmer import steady_state_flight_trim from pyfme.utils.coordinates import wind2body -TAS = 120 # m/s +TAS = 80 # m/s h = 1000 results = steady_state_flight_trim(cessna_310, h, TAS, gamma=0, turn_rate=0) @@ -31,4 +31,3 @@ print("beta = {}".format(results[4])) print("control = {}".format(results[5])) print('') - diff --git a/src/pyfme/utils/trimmer.py b/src/pyfme/utils/trimmer.py index a740a87..9f712f7 100644 --- a/src/pyfme/utils/trimmer.py +++ b/src/pyfme/utils/trimmer.py @@ -16,7 +16,7 @@ from warnings import warn import numpy as np from math import sqrt, sin, cos, tan, atan -from scipy.optimize import fmin +from scipy.optimize import least_squares from pyfme.utils.coordinates import wind2body from pyfme.environment.isa import atm @@ -80,13 +80,18 @@ def steady_state_flight_trim(aircraft, h, TAS, gamma=0, turn_rate=0, args = (h, TAS, gamma, turn_rate, aircraft, dynamic_eqs) - results = fmin(func, x0=initial_guess, args=args, full_output=True, - disp=False) - trimmed_params, fopt = results[0:2] + lower_bounds = (-np.inf, -np.inf, -np.inf, -np.inf, -np.inf, 0) + upper_bounds = (+np.inf, +np.inf, +np.inf, +np.inf, +np.inf, 1) - if abs(fopt) > 1e-3: - warn("Trim process did not converge. fopt={}".format(fopt), - RuntimeWarning) + results = least_squares(func, x0=initial_guess, args=args, verbose=2, + bounds=(lower_bounds, upper_bounds)) + + trimmed_params = results['x'] + fun = results['fun'] + cost = results['cost'] + + if cost > 1e-7 or any(abs(fun) > 1e-3): + warn("Trim process did not converge", RuntimeWarning) alpha = trimmed_params[0] beta = trimmed_params[1] @@ -215,4 +220,4 @@ def func(trimmed_params, h, TAS, gamma, turn_rate, aircraft, dynamic_eqs): output = dynamic_eqs(0, vel, mass, inertia, forces, moments) - return np.sum(np.dot(output, output)) + return output From fc580e369d6fc827e4b700c916a5d68650ffda3e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?=C3=81lex=20S=C3=A1ez?= Date: Fri, 12 Feb 2016 20:01:47 +0100 Subject: [PATCH 07/17] First 'working' version --- examples/example_003.py | 126 ++++++++++++++++++++++++++++++ src/pyfme/aircrafts/cessna_310.py | 12 +-- src/pyfme/utils/trimmer.py | 30 +++++-- 3 files changed, 151 insertions(+), 17 deletions(-) create mode 100644 examples/example_003.py diff --git a/examples/example_003.py b/examples/example_003.py new file mode 100644 index 0000000..3776d67 --- /dev/null +++ b/examples/example_003.py @@ -0,0 +1,126 @@ +# -*- coding: utf-8 -*- +""" +Python Flight Mechanics Engine (PyFME). +Copyright (c) AeroPython Development Team. +Distributed under the terms of the MIT License. + + + +""" + +# -*- coding: utf-8 -*- +# -*- coding: utf-8 -*- +""" +Dummy example + +@author: asaez +""" + +import numpy as np +import matplotlib.pyplot as plt + +from pyfme.aircrafts import cessna_310 +from pyfme.models.system import System +from pyfme.utils.trimmer import steady_state_flight_trim +from pyfme.environment.isa import atm + + +if __name__ == '__main__': + + # Case params + mass, inertia = cessna_310.mass_and_inertial_data() + + # initial conditions + TAS = 85 # m/s + h = 3000 # m + psi_0 = 3 # rad + x_0, y_0 = 0, 0 # m + + trim_results = steady_state_flight_trim(cessna_310, h, TAS, gamma=0.5, + turn_rate=0) + + lin_vel, ang_vel, theta, phi, alpha, beta, control_vector = trim_results + # Time + t0 = 0 + tf = 30 + dt = 1e-1 + time = np.arange(t0, tf, dt) + + # Results initialization + results = np.zeros([time.size, 12]) + + velocities = np.empty([time.size, 6]) + results[0, 0:3] = lin_vel + results[0, 3:6] = ang_vel + results[0, 6:9] = theta, phi, psi_0 + results[0, 9:12] = x_0, y_0, h + + # Linear Momentum and Angular Momentum eqs + equations = System(integrator='dopri5', model='euler_flat_earth', jac=False) + u, v, w = lin_vel + p, q, r = ang_vel + equations.set_initial_values(u, v, w, p, q, r, theta, phi, psi_0, x_0, y_0, + h) + + _, _, rho, _ = atm(h) + + delta_e, delta_ail, delta_r, delta_t = control_vector + attitude = theta, phi, psi_0 + + for ii, t in enumerate(time[1:]): + forces, moments = cessna_310.get_forces_and_moments(TAS, rho, alpha, beta, + delta_e, 0, delta_ail, delta_r, delta_t, + attitude) + + results[ii+1, :] = equations.propagate(mass, inertia, forces, moments, + dt) + + lin_vel = results[ii+1, 0:3] + ang_vel = results[ii+1, 3:6] + attitude = results[ii+1, 6:9] + position = results[ii+1, 9:12] + + _, _, rho, _ = atm(position[2]) + + TAS = np.sqrt(np.sum(lin_vel*lin_vel)) + + +velocities = results[:, 0:6] +attitude_angles = results[:, 6:9] +position = results[:, 9:12] + +plt.close('all') +plt.style.use('ggplot') + +# TODO: improve graphics: legend title.... +plt.figure('pos') +plt.plot(time, position[:, 0], label='x') +plt.plot(time, position[:, 1], label='y') +plt.plot(time, position[:, 2], label='z') +plt.xlabel('time (s)') +plt.ylabel('position (m)') +plt.legend() + +plt.figure('angles') +plt.plot(time, attitude_angles[:, 0], label='theta') +plt.plot(time, attitude_angles[:, 1], label='phi') +plt.plot(time, attitude_angles[:, 2], label='psi') +plt.xlabel('time (s)') +plt.ylabel('attitude (rad)') +plt.legend() + +plt.figure('velocities') +plt.plot(time, velocities[:, 0], label='u') +plt.plot(time, velocities[:, 1], label='v') +plt.plot(time, velocities[:, 2], label='w') +plt.xlabel('time (s)') +plt.ylabel('velocity (m/s)') +plt.legend() + +plt.figure('ang velocities') +plt.plot(time, velocities[:, 3], label='p') +plt.plot(time, velocities[:, 4], label='q') +plt.plot(time, velocities[:, 5], label='r') +plt.xlabel('time (s)') +plt.xlabel('angular velocity (rad/s)') +plt.legend() \ No newline at end of file diff --git a/src/pyfme/aircrafts/cessna_310.py b/src/pyfme/aircrafts/cessna_310.py index 0392696..3db46bb 100644 --- a/src/pyfme/aircrafts/cessna_310.py +++ b/src/pyfme/aircrafts/cessna_310.py @@ -50,17 +50,9 @@ def mass_and_inertial_data(): Returns ------ mass : float - mass (lb * 0.453592 = kg) - Ixx_b : float - Moment of Inertia x-axis ( slug * ft2 * 1.3558179 = Kg * m2) - Iyy_b : float - Moment of Inertia y-axis ( slug * ft2 * 1.3558179 = Kg * m2) - Izz_b : float - Moment of Inertia z-axis ( slug * ft2 * 1.3558179 = Kg * m2) - Ixz_b : float - Product of Inertia xz-plane ( slug * ft2 * 1.3558179 = Kg * m2) + mass (lb * 0.453592 = kg) inertia : array_like - I_xx_b,I_yy_b,I_zz_b] + Inertia tensor (Kg * m^2) References ---------- diff --git a/src/pyfme/utils/trimmer.py b/src/pyfme/utils/trimmer.py index 9f712f7..5feb2f7 100644 --- a/src/pyfme/utils/trimmer.py +++ b/src/pyfme/utils/trimmer.py @@ -23,7 +23,7 @@ def steady_state_flight_trim(aircraft, h, TAS, gamma=0, turn_rate=0, - dyn_eqs=None): + dyn_eqs=None, verbose=0): """Finds a combination of values of the state and control variables that correspond to a steady-state flight condition. Steady-state aircraft flight can be defined as a condition in which all of the motion variables are @@ -75,15 +75,28 @@ def steady_state_flight_trim(aircraft, h, TAS, gamma=0, turn_rate=0, from pyfme.models.euler_flat_earth import linear_and_angular_momentum_eqs dynamic_eqs = linear_and_angular_momentum_eqs - # TODO: try to look for a good inizialization - initial_guess = np.zeros([6]) + 0.01 + # TODO: try to look for a good inizialization method + alpha_0 = 0.05 * np.sign(gamma) + betha_0 = 0.001 * np.sign(turn_rate) + delta_e_0 = 0.05 + delta_ail_0 = 0.01 * np.sign(turn_rate) + delta_r_0 = 0.01 * np.sign(turn_rate) + delta_t_0 = 0.5 + + initial_guess = (alpha_0, + betha_0, + delta_e_0, + delta_ail_0, + delta_r_0, + delta_t_0) + args = (h, TAS, gamma, turn_rate, aircraft, dynamic_eqs) - lower_bounds = (-np.inf, -np.inf, -np.inf, -np.inf, -np.inf, 0) - upper_bounds = (+np.inf, +np.inf, +np.inf, +np.inf, +np.inf, 1) + lower_bounds = (-1, -0.5, -1, -1, -1, 0) + upper_bounds = (+1, +0.5, +1, +1, +1, 1) - results = least_squares(func, x0=initial_guess, args=args, verbose=2, + results = least_squares(func, x0=initial_guess, args=args, verbose=verbose, bounds=(lower_bounds, upper_bounds)) trimmed_params = results['x'] @@ -92,6 +105,8 @@ def steady_state_flight_trim(aircraft, h, TAS, gamma=0, turn_rate=0, if cost > 1e-7 or any(abs(fun) > 1e-3): warn("Trim process did not converge", RuntimeWarning) + if trimmed_params[5] > 0.99: + warn("Probably not enough power for demanded conditions") alpha = trimmed_params[0] beta = trimmed_params[1] @@ -119,8 +134,9 @@ def steady_state_flight_trim(aircraft, h, TAS, gamma=0, turn_rate=0, r = turn_rate * cos(theta) * sin(phi) ang_vel = np.array([p, q, r]) + lin_vel = wind2body((TAS, 0, 0), alpha, beta) - return ang_vel, theta, phi, alpha, beta, control_vector + return lin_vel, ang_vel, theta, phi, alpha, beta, control_vector def turn_coord_cons1(turn_rate, alpha, beta, TAS, gamma=0): From 8741e10043ba68f50566f10ca12a3e4c67a32f4e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?=C3=81lex=20S=C3=A1ez?= Date: Sun, 14 Feb 2016 17:02:15 +0100 Subject: [PATCH 08/17] Added example 003 to check if the trimmed condition is maintained --- examples/example_003.py | 171 +++++++++++++++++++++++----------------- 1 file changed, 98 insertions(+), 73 deletions(-) diff --git a/examples/example_003.py b/examples/example_003.py index 3776d67..2803c98 100644 --- a/examples/example_003.py +++ b/examples/example_003.py @@ -4,16 +4,10 @@ Copyright (c) AeroPython Development Team. Distributed under the terms of the MIT License. - - -""" - -# -*- coding: utf-8 -*- -# -*- coding: utf-8 -*- -""" -Dummy example - -@author: asaez +Example with trimmed aircraft. +The main purpose of this example is to check if the aircraft trimmed in a given +state maintains the trimmed flight condition. +Trimmed in stationary, horizontal, symmetric, wings level flight. """ import numpy as np @@ -22,55 +16,73 @@ from pyfme.aircrafts import cessna_310 from pyfme.models.system import System from pyfme.utils.trimmer import steady_state_flight_trim +from pyfme.utils.anemometry import calculate_alpha_beta_TAS from pyfme.environment.isa import atm if __name__ == '__main__': - # Case params + # Aircraft parameters. mass, inertia = cessna_310.mass_and_inertial_data() - # initial conditions - TAS = 85 # m/s - h = 3000 # m + # Initial conditions. + TAS_ = 312 * 0.3048 # m/s + h = 8000 * 0.3048 # m psi_0 = 3 # rad x_0, y_0 = 0, 0 # m - trim_results = steady_state_flight_trim(cessna_310, h, TAS, gamma=0.5, + # Trimming. + trim_results = steady_state_flight_trim(cessna_310, h, TAS_, gamma=0, turn_rate=0) - lin_vel, ang_vel, theta, phi, alpha, beta, control_vector = trim_results - # Time - t0 = 0 - tf = 30 - dt = 1e-1 + lin_vel, ang_vel, theta, phi, alpha_, beta_, control_vector = trim_results + + # Time. + t0 = 0 # s + tf = 30 # s + dt = 1e-2 # s + time = np.arange(t0, tf, dt) - # Results initialization + # Results initialization. results = np.zeros([time.size, 12]) - velocities = np.empty([time.size, 6]) results[0, 0:3] = lin_vel results[0, 3:6] = ang_vel results[0, 6:9] = theta, phi, psi_0 results[0, 9:12] = x_0, y_0, h - - # Linear Momentum and Angular Momentum eqs - equations = System(integrator='dopri5', model='euler_flat_earth', jac=False) + alpha = np.empty_like(time) + alpha[0] = alpha_ + beta = np.empty_like(time) + beta[0] = beta_ + TAS = np.empty_like(time) + TAS[0] = TAS_ + + # Linear Momentum and Angular Momentum eqs. + equations = System(integrator='dopri5', + model='euler_flat_earth', + jac=False) u, v, w = lin_vel p, q, r = ang_vel - equations.set_initial_values(u, v, w, p, q, r, theta, phi, psi_0, x_0, y_0, - h) + + equations.set_initial_values(u, v, w, + p, q, r, + theta, phi, psi_0, + x_0, y_0, h) _, _, rho, _ = atm(h) - delta_e, delta_ail, delta_r, delta_t = control_vector + # Define control vectors. + d_e, d_a, d_r, d_t = control_vector + attitude = theta, phi, psi_0 + # Rename function to make it shorter + forces_and_moments = cessna_310.get_forces_and_moments for ii, t in enumerate(time[1:]): - forces, moments = cessna_310.get_forces_and_moments(TAS, rho, alpha, beta, - delta_e, 0, delta_ail, delta_r, delta_t, - attitude) + + forces, moments = forces_and_moments(TAS[ii], rho, alpha[ii], beta[ii], + d_e, 0, d_a, d_r, d_t, attitude) results[ii+1, :] = equations.propagate(mass, inertia, forces, moments, dt) @@ -80,47 +92,60 @@ attitude = results[ii+1, 6:9] position = results[ii+1, 9:12] + alpha[ii+1], beta[ii+1], TAS[ii+1] = calculate_alpha_beta_TAS(*lin_vel) + _, _, rho, _ = atm(position[2]) - TAS = np.sqrt(np.sum(lin_vel*lin_vel)) - - -velocities = results[:, 0:6] -attitude_angles = results[:, 6:9] -position = results[:, 9:12] - -plt.close('all') -plt.style.use('ggplot') - -# TODO: improve graphics: legend title.... -plt.figure('pos') -plt.plot(time, position[:, 0], label='x') -plt.plot(time, position[:, 1], label='y') -plt.plot(time, position[:, 2], label='z') -plt.xlabel('time (s)') -plt.ylabel('position (m)') -plt.legend() - -plt.figure('angles') -plt.plot(time, attitude_angles[:, 0], label='theta') -plt.plot(time, attitude_angles[:, 1], label='phi') -plt.plot(time, attitude_angles[:, 2], label='psi') -plt.xlabel('time (s)') -plt.ylabel('attitude (rad)') -plt.legend() - -plt.figure('velocities') -plt.plot(time, velocities[:, 0], label='u') -plt.plot(time, velocities[:, 1], label='v') -plt.plot(time, velocities[:, 2], label='w') -plt.xlabel('time (s)') -plt.ylabel('velocity (m/s)') -plt.legend() - -plt.figure('ang velocities') -plt.plot(time, velocities[:, 3], label='p') -plt.plot(time, velocities[:, 4], label='q') -plt.plot(time, velocities[:, 5], label='r') -plt.xlabel('time (s)') -plt.xlabel('angular velocity (rad/s)') -plt.legend() \ No newline at end of file + velocities = results[:, 0:6] + attitude_angles = results[:, 6:9] + position = results[:, 9:12] + + # PLOTS + plt.close('all') + plt.style.use('ggplot') + + plt.figure('pos') + plt.plot(time, position[:, 0], label='x') + plt.plot(time, position[:, 1], label='y') + plt.plot(time, position[:, 2], label='z') + plt.xlabel('time (s)') + plt.ylabel('position (m)') + plt.legend() + + plt.figure('angles') + plt.plot(time, attitude_angles[:, 0], label='theta') + plt.plot(time, attitude_angles[:, 1], label='phi') + plt.plot(time, attitude_angles[:, 2], label='psi') + plt.xlabel('time (s)') + plt.ylabel('attitude (rad)') + plt.legend() + + plt.figure('velocities') + plt.plot(time, velocities[:, 0], label='u') + plt.plot(time, velocities[:, 1], label='v') + plt.plot(time, velocities[:, 2], label='w') + plt.plot(time, TAS, label='TAS') + plt.xlabel('time (s)') + plt.ylabel('velocity (m/s)') + plt.legend() + + plt.figure('ang velocities') + plt.plot(time, velocities[:, 3], label='p') + plt.plot(time, velocities[:, 4], label='q') + plt.plot(time, velocities[:, 5], label='r') + plt.xlabel('time (s)') + plt.ylabel('angular velocity (rad/s)') + plt.legend() + + plt.figure('aero angles') + plt.plot(time, alpha, label='alpha') + plt.plot(time, beta, label='beta') + plt.xlabel('time (s)') + plt.ylabel('angle (rad)') + plt.legend() + + plt.figure('2D Trajectory') + plt.plot(position[:, 0], position[:, 1]) + plt.xlabel('X (m)') + plt.ylabel('Y (m)') + plt.legend() From ec04e975b59e196f95f3b773fa719a12bbcdbf7b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?=C3=81lex=20S=C3=A1ez?= Date: Sun, 14 Feb 2016 17:07:27 +0100 Subject: [PATCH 09/17] Renamed example_003 to example_001 and removed example 1 --- examples/example_001.py | 163 ++++++++++++++++++++++++++++------------ examples/example_003.py | 151 ------------------------------------- 2 files changed, 116 insertions(+), 198 deletions(-) delete mode 100644 examples/example_003.py diff --git a/examples/example_001.py b/examples/example_001.py index 10fe420..2803c98 100644 --- a/examples/example_001.py +++ b/examples/example_001.py @@ -1,82 +1,151 @@ # -*- coding: utf-8 -*- """ -Dummy example - -@author: asaez +Python Flight Mechanics Engine (PyFME). +Copyright (c) AeroPython Development Team. +Distributed under the terms of the MIT License. + +Example with trimmed aircraft. +The main purpose of this example is to check if the aircraft trimmed in a given +state maintains the trimmed flight condition. +Trimmed in stationary, horizontal, symmetric, wings level flight. """ import numpy as np import matplotlib.pyplot as plt +from pyfme.aircrafts import cessna_310 from pyfme.models.system import System +from pyfme.utils.trimmer import steady_state_flight_trim +from pyfme.utils.anemometry import calculate_alpha_beta_TAS +from pyfme.environment.isa import atm +if __name__ == '__main__': -def forces(): - return np.array([0, 0, -9.8*500]) - - -def moments(): - return np.array([0, 0, 0]) - + # Aircraft parameters. + mass, inertia = cessna_310.mass_and_inertial_data() -if __name__ == '__main__': + # Initial conditions. + TAS_ = 312 * 0.3048 # m/s + h = 8000 * 0.3048 # m + psi_0 = 3 # rad + x_0, y_0 = 0, 0 # m - # Case params - mass = 500 - inertia = np.array([[1, 0, 0], - [0, 1, 0], - [0, 0, 1]]) + # Trimming. + trim_results = steady_state_flight_trim(cessna_310, h, TAS_, gamma=0, + turn_rate=0) - # initial conditions - u, v, w = 0, 0, 0 - p, q, r = 0, 0, 0 + lin_vel, ang_vel, theta, phi, alpha_, beta_, control_vector = trim_results - theta, phi, psi = 0, 0, 0 - x, y, z = 0, 0, 5000 + # Time. + t0 = 0 # s + tf = 30 # s + dt = 1e-2 # s - t0 = 0 - tf = 10 - dt = 1e-2 time = np.arange(t0, tf, dt) + # Results initialization. results = np.zeros([time.size, 12]) - velocities = np.empty([time.size, 6]) - results[0, 0:6] = u, v, w, p, q, r - results[0, 6:9] = theta, phi, psi - results[0, 9:12] = x, y, z + results[0, 0:3] = lin_vel + results[0, 3:6] = ang_vel + results[0, 6:9] = theta, phi, psi_0 + results[0, 9:12] = x_0, y_0, h + alpha = np.empty_like(time) + alpha[0] = alpha_ + beta = np.empty_like(time) + beta[0] = beta_ + TAS = np.empty_like(time) + TAS[0] = TAS_ + + # Linear Momentum and Angular Momentum eqs. + equations = System(integrator='dopri5', + model='euler_flat_earth', + jac=False) + u, v, w = lin_vel + p, q, r = ang_vel + + equations.set_initial_values(u, v, w, + p, q, r, + theta, phi, psi_0, + x_0, y_0, h) + + _, _, rho, _ = atm(h) + + # Define control vectors. + d_e, d_a, d_r, d_t = control_vector + + attitude = theta, phi, psi_0 + + # Rename function to make it shorter + forces_and_moments = cessna_310.get_forces_and_moments + for ii, t in enumerate(time[1:]): - # Linear Momentum and Angular Momentum eqs - equations = System(integrator='dopri5', model='euler_flat_earth', jac=False) - equations.set_initial_values(u, v, w, p, q, r, theta, phi, psi, x, y, z) + forces, moments = forces_and_moments(TAS[ii], rho, alpha[ii], beta[ii], + d_e, 0, d_a, d_r, d_t, attitude) + results[ii+1, :] = equations.propagate(mass, inertia, forces, moments, + dt) - for ii, t in enumerate(time[1:]): - results[ii+1, :] = equations.propagate(mass, inertia, forces(), moments(), - dt=dt) + lin_vel = results[ii+1, 0:3] + ang_vel = results[ii+1, 3:6] + attitude = results[ii+1, 6:9] + position = results[ii+1, 9:12] + + alpha[ii+1], beta[ii+1], TAS[ii+1] = calculate_alpha_beta_TAS(*lin_vel) + + _, _, rho, _ = atm(position[2]) velocities = results[:, 0:6] attitude_angles = results[:, 6:9] position = results[:, 9:12] - # TODO: improve graphics: legend title.... + # PLOTS + plt.close('all') + plt.style.use('ggplot') + plt.figure('pos') - plt.plot(time, position[:, 0]) - plt.plot(time, position[:, 1]) - plt.plot(time, position[:, 2]) + plt.plot(time, position[:, 0], label='x') + plt.plot(time, position[:, 1], label='y') + plt.plot(time, position[:, 2], label='z') + plt.xlabel('time (s)') + plt.ylabel('position (m)') + plt.legend() plt.figure('angles') - plt.plot(time, attitude_angles[:, 0]) - plt.plot(time, attitude_angles[:, 1]) - plt.plot(time, attitude_angles[:, 2]) + plt.plot(time, attitude_angles[:, 0], label='theta') + plt.plot(time, attitude_angles[:, 1], label='phi') + plt.plot(time, attitude_angles[:, 2], label='psi') + plt.xlabel('time (s)') + plt.ylabel('attitude (rad)') + plt.legend() plt.figure('velocities') - plt.plot(time, velocities[:, 0]) - plt.plot(time, velocities[:, 1]) - plt.plot(time, velocities[:, 2]) + plt.plot(time, velocities[:, 0], label='u') + plt.plot(time, velocities[:, 1], label='v') + plt.plot(time, velocities[:, 2], label='w') + plt.plot(time, TAS, label='TAS') + plt.xlabel('time (s)') + plt.ylabel('velocity (m/s)') + plt.legend() plt.figure('ang velocities') - plt.plot(time, velocities[:, 3]) - plt.plot(time, velocities[:, 4]) - plt.plot(time, velocities[:, 5]) + plt.plot(time, velocities[:, 3], label='p') + plt.plot(time, velocities[:, 4], label='q') + plt.plot(time, velocities[:, 5], label='r') + plt.xlabel('time (s)') + plt.ylabel('angular velocity (rad/s)') + plt.legend() + + plt.figure('aero angles') + plt.plot(time, alpha, label='alpha') + plt.plot(time, beta, label='beta') + plt.xlabel('time (s)') + plt.ylabel('angle (rad)') + plt.legend() + + plt.figure('2D Trajectory') + plt.plot(position[:, 0], position[:, 1]) + plt.xlabel('X (m)') + plt.ylabel('Y (m)') + plt.legend() diff --git a/examples/example_003.py b/examples/example_003.py deleted file mode 100644 index 2803c98..0000000 --- a/examples/example_003.py +++ /dev/null @@ -1,151 +0,0 @@ -# -*- coding: utf-8 -*- -""" -Python Flight Mechanics Engine (PyFME). -Copyright (c) AeroPython Development Team. -Distributed under the terms of the MIT License. - -Example with trimmed aircraft. -The main purpose of this example is to check if the aircraft trimmed in a given -state maintains the trimmed flight condition. -Trimmed in stationary, horizontal, symmetric, wings level flight. -""" - -import numpy as np -import matplotlib.pyplot as plt - -from pyfme.aircrafts import cessna_310 -from pyfme.models.system import System -from pyfme.utils.trimmer import steady_state_flight_trim -from pyfme.utils.anemometry import calculate_alpha_beta_TAS -from pyfme.environment.isa import atm - - -if __name__ == '__main__': - - # Aircraft parameters. - mass, inertia = cessna_310.mass_and_inertial_data() - - # Initial conditions. - TAS_ = 312 * 0.3048 # m/s - h = 8000 * 0.3048 # m - psi_0 = 3 # rad - x_0, y_0 = 0, 0 # m - - # Trimming. - trim_results = steady_state_flight_trim(cessna_310, h, TAS_, gamma=0, - turn_rate=0) - - lin_vel, ang_vel, theta, phi, alpha_, beta_, control_vector = trim_results - - # Time. - t0 = 0 # s - tf = 30 # s - dt = 1e-2 # s - - time = np.arange(t0, tf, dt) - - # Results initialization. - results = np.zeros([time.size, 12]) - - results[0, 0:3] = lin_vel - results[0, 3:6] = ang_vel - results[0, 6:9] = theta, phi, psi_0 - results[0, 9:12] = x_0, y_0, h - alpha = np.empty_like(time) - alpha[0] = alpha_ - beta = np.empty_like(time) - beta[0] = beta_ - TAS = np.empty_like(time) - TAS[0] = TAS_ - - # Linear Momentum and Angular Momentum eqs. - equations = System(integrator='dopri5', - model='euler_flat_earth', - jac=False) - u, v, w = lin_vel - p, q, r = ang_vel - - equations.set_initial_values(u, v, w, - p, q, r, - theta, phi, psi_0, - x_0, y_0, h) - - _, _, rho, _ = atm(h) - - # Define control vectors. - d_e, d_a, d_r, d_t = control_vector - - attitude = theta, phi, psi_0 - - # Rename function to make it shorter - forces_and_moments = cessna_310.get_forces_and_moments - for ii, t in enumerate(time[1:]): - - forces, moments = forces_and_moments(TAS[ii], rho, alpha[ii], beta[ii], - d_e, 0, d_a, d_r, d_t, attitude) - - results[ii+1, :] = equations.propagate(mass, inertia, forces, moments, - dt) - - lin_vel = results[ii+1, 0:3] - ang_vel = results[ii+1, 3:6] - attitude = results[ii+1, 6:9] - position = results[ii+1, 9:12] - - alpha[ii+1], beta[ii+1], TAS[ii+1] = calculate_alpha_beta_TAS(*lin_vel) - - _, _, rho, _ = atm(position[2]) - - velocities = results[:, 0:6] - attitude_angles = results[:, 6:9] - position = results[:, 9:12] - - # PLOTS - plt.close('all') - plt.style.use('ggplot') - - plt.figure('pos') - plt.plot(time, position[:, 0], label='x') - plt.plot(time, position[:, 1], label='y') - plt.plot(time, position[:, 2], label='z') - plt.xlabel('time (s)') - plt.ylabel('position (m)') - plt.legend() - - plt.figure('angles') - plt.plot(time, attitude_angles[:, 0], label='theta') - plt.plot(time, attitude_angles[:, 1], label='phi') - plt.plot(time, attitude_angles[:, 2], label='psi') - plt.xlabel('time (s)') - plt.ylabel('attitude (rad)') - plt.legend() - - plt.figure('velocities') - plt.plot(time, velocities[:, 0], label='u') - plt.plot(time, velocities[:, 1], label='v') - plt.plot(time, velocities[:, 2], label='w') - plt.plot(time, TAS, label='TAS') - plt.xlabel('time (s)') - plt.ylabel('velocity (m/s)') - plt.legend() - - plt.figure('ang velocities') - plt.plot(time, velocities[:, 3], label='p') - plt.plot(time, velocities[:, 4], label='q') - plt.plot(time, velocities[:, 5], label='r') - plt.xlabel('time (s)') - plt.ylabel('angular velocity (rad/s)') - plt.legend() - - plt.figure('aero angles') - plt.plot(time, alpha, label='alpha') - plt.plot(time, beta, label='beta') - plt.xlabel('time (s)') - plt.ylabel('angle (rad)') - plt.legend() - - plt.figure('2D Trajectory') - plt.plot(position[:, 0], position[:, 1]) - plt.xlabel('X (m)') - plt.ylabel('Y (m)') - plt.legend() From a284e661c1bc186c28c69519442cab53d79c6817 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?=C3=81lex=20S=C3=A1ez?= Date: Sun, 14 Feb 2016 17:25:03 +0100 Subject: [PATCH 10/17] Added more examples for different flight conditions and for longitudinal perturbations --- examples/example_002.py | 151 ++++++++++++++++++++++++++++++++++++ examples/example_003.py | 151 ++++++++++++++++++++++++++++++++++++ examples/example_004.py | 166 ++++++++++++++++++++++++++++++++++++++++ 3 files changed, 468 insertions(+) create mode 100644 examples/example_002.py create mode 100644 examples/example_003.py create mode 100644 examples/example_004.py diff --git a/examples/example_002.py b/examples/example_002.py new file mode 100644 index 0000000..40c8b9a --- /dev/null +++ b/examples/example_002.py @@ -0,0 +1,151 @@ +# -*- coding: utf-8 -*- +""" +Python Flight Mechanics Engine (PyFME). +Copyright (c) AeroPython Development Team. +Distributed under the terms of the MIT License. + +Example with trimmed aircraft. +The main purpose of this example is to check if the aircraft trimmed in a given +state maintains the trimmed flight condition. +Trimmed in stationary descent, symmetric, wings level flight. +""" + +import numpy as np +import matplotlib.pyplot as plt + +from pyfme.aircrafts import cessna_310 +from pyfme.models.system import System +from pyfme.utils.trimmer import steady_state_flight_trim +from pyfme.utils.anemometry import calculate_alpha_beta_TAS +from pyfme.environment.isa import atm + + +if __name__ == '__main__': + + # Aircraft parameters. + mass, inertia = cessna_310.mass_and_inertial_data() + + # Initial conditions. + TAS_ = 312 * 0.3048 # m/s + h = 8000 * 0.3048 # m + psi_0 = 3 # rad + x_0, y_0 = 0, 0 # m + + # Trimming. + trim_results = steady_state_flight_trim(cessna_310, h, TAS_, gamma=0.20, + turn_rate=0) + + lin_vel, ang_vel, theta, phi, alpha_, beta_, control_vector = trim_results + + # Time. + t0 = 0 # s + tf = 120 # s + dt = 1e-2 # s + + time = np.arange(t0, tf, dt) + + # Results initialization. + results = np.zeros([time.size, 12]) + + results[0, 0:3] = lin_vel + results[0, 3:6] = ang_vel + results[0, 6:9] = theta, phi, psi_0 + results[0, 9:12] = x_0, y_0, h + alpha = np.empty_like(time) + alpha[0] = alpha_ + beta = np.empty_like(time) + beta[0] = beta_ + TAS = np.empty_like(time) + TAS[0] = TAS_ + + # Linear Momentum and Angular Momentum eqs. + equations = System(integrator='dopri5', + model='euler_flat_earth', + jac=False) + u, v, w = lin_vel + p, q, r = ang_vel + + equations.set_initial_values(u, v, w, + p, q, r, + theta, phi, psi_0, + x_0, y_0, h) + + _, _, rho, _ = atm(h) + + # Define control vectors. + d_e, d_a, d_r, d_t = control_vector + + attitude = theta, phi, psi_0 + + # Rename function to make it shorter + forces_and_moments = cessna_310.get_forces_and_moments + for ii, t in enumerate(time[1:]): + + forces, moments = forces_and_moments(TAS[ii], rho, alpha[ii], beta[ii], + d_e, 0, d_a, d_r, d_t, attitude) + + results[ii+1, :] = equations.propagate(mass, inertia, forces, moments, + dt) + + lin_vel = results[ii+1, 0:3] + ang_vel = results[ii+1, 3:6] + attitude = results[ii+1, 6:9] + position = results[ii+1, 9:12] + + alpha[ii+1], beta[ii+1], TAS[ii+1] = calculate_alpha_beta_TAS(*lin_vel) + + _, _, rho, _ = atm(position[2]) + + velocities = results[:, 0:6] + attitude_angles = results[:, 6:9] + position = results[:, 9:12] + + # PLOTS + plt.close('all') + plt.style.use('ggplot') + + plt.figure('pos') + plt.plot(time, position[:, 0], label='x') + plt.plot(time, position[:, 1], label='y') + plt.plot(time, position[:, 2], label='z') + plt.xlabel('time (s)') + plt.ylabel('position (m)') + plt.legend() + + plt.figure('angles') + plt.plot(time, attitude_angles[:, 0], label='theta') + plt.plot(time, attitude_angles[:, 1], label='phi') + plt.plot(time, attitude_angles[:, 2], label='psi') + plt.xlabel('time (s)') + plt.ylabel('attitude (rad)') + plt.legend() + + plt.figure('velocities') + plt.plot(time, velocities[:, 0], label='u') + plt.plot(time, velocities[:, 1], label='v') + plt.plot(time, velocities[:, 2], label='w') + plt.plot(time, TAS, label='TAS') + plt.xlabel('time (s)') + plt.ylabel('velocity (m/s)') + plt.legend() + + plt.figure('ang velocities') + plt.plot(time, velocities[:, 3], label='p') + plt.plot(time, velocities[:, 4], label='q') + plt.plot(time, velocities[:, 5], label='r') + plt.xlabel('time (s)') + plt.ylabel('angular velocity (rad/s)') + plt.legend() + + plt.figure('aero angles') + plt.plot(time, alpha, label='alpha') + plt.plot(time, beta, label='beta') + plt.xlabel('time (s)') + plt.ylabel('angle (rad)') + plt.legend() + + plt.figure('2D Trajectory') + plt.plot(position[:, 0], position[:, 1]) + plt.xlabel('X (m)') + plt.ylabel('Y (m)') + plt.legend() diff --git a/examples/example_003.py b/examples/example_003.py new file mode 100644 index 0000000..7a5c646 --- /dev/null +++ b/examples/example_003.py @@ -0,0 +1,151 @@ +# -*- coding: utf-8 -*- +""" +Python Flight Mechanics Engine (PyFME). +Copyright (c) AeroPython Development Team. +Distributed under the terms of the MIT License. + +Example with trimmed aircraft. +The main purpose of this example is to check if the aircraft trimmed in a given +state maintains the trimmed flight condition. +Trimmed in stationary, horizontal, turn. +""" + +import numpy as np +import matplotlib.pyplot as plt + +from pyfme.aircrafts import cessna_310 +from pyfme.models.system import System +from pyfme.utils.trimmer import steady_state_flight_trim +from pyfme.utils.anemometry import calculate_alpha_beta_TAS +from pyfme.environment.isa import atm + + +if __name__ == '__main__': + + # Aircraft parameters. + mass, inertia = cessna_310.mass_and_inertial_data() + + # Initial conditions. + TAS_ = 312 * 0.3048 # m/s + h = 8000 * 0.3048 # m + psi_0 = 3 # rad + x_0, y_0 = 0, 0 # m + + # Trimming. + trim_results = steady_state_flight_trim(cessna_310, h, TAS_, gamma=0.02, + turn_rate=0) + + lin_vel, ang_vel, theta, phi, alpha_, beta_, control_vector = trim_results + + # Time. + t0 = 0 # s + tf = 120 # s + dt = 1e-2 # s + + time = np.arange(t0, tf, dt) + + # Results initialization. + results = np.zeros([time.size, 12]) + + results[0, 0:3] = lin_vel + results[0, 3:6] = ang_vel + results[0, 6:9] = theta, phi, psi_0 + results[0, 9:12] = x_0, y_0, h + alpha = np.empty_like(time) + alpha[0] = alpha_ + beta = np.empty_like(time) + beta[0] = beta_ + TAS = np.empty_like(time) + TAS[0] = TAS_ + + # Linear Momentum and Angular Momentum eqs. + equations = System(integrator='dopri5', + model='euler_flat_earth', + jac=False) + u, v, w = lin_vel + p, q, r = ang_vel + + equations.set_initial_values(u, v, w, + p, q, r, + theta, phi, psi_0, + x_0, y_0, h) + + _, _, rho, _ = atm(h) + + # Define control vectors. + d_e, d_a, d_r, d_t = control_vector + + attitude = theta, phi, psi_0 + + # Rename function to make it shorter + forces_and_moments = cessna_310.get_forces_and_moments + for ii, t in enumerate(time[1:]): + + forces, moments = forces_and_moments(TAS[ii], rho, alpha[ii], beta[ii], + d_e, 0, d_a, d_r, d_t, attitude) + + results[ii+1, :] = equations.propagate(mass, inertia, forces, moments, + dt) + + lin_vel = results[ii+1, 0:3] + ang_vel = results[ii+1, 3:6] + attitude = results[ii+1, 6:9] + position = results[ii+1, 9:12] + + alpha[ii+1], beta[ii+1], TAS[ii+1] = calculate_alpha_beta_TAS(*lin_vel) + + _, _, rho, _ = atm(position[2]) + + velocities = results[:, 0:6] + attitude_angles = results[:, 6:9] + position = results[:, 9:12] + + # PLOTS + plt.close('all') + plt.style.use('ggplot') + + plt.figure('pos') + plt.plot(time, position[:, 0], label='x') + plt.plot(time, position[:, 1], label='y') + plt.plot(time, position[:, 2], label='z') + plt.xlabel('time (s)') + plt.ylabel('position (m)') + plt.legend() + + plt.figure('angles') + plt.plot(time, attitude_angles[:, 0], label='theta') + plt.plot(time, attitude_angles[:, 1], label='phi') + plt.plot(time, attitude_angles[:, 2], label='psi') + plt.xlabel('time (s)') + plt.ylabel('attitude (rad)') + plt.legend() + + plt.figure('velocities') + plt.plot(time, velocities[:, 0], label='u') + plt.plot(time, velocities[:, 1], label='v') + plt.plot(time, velocities[:, 2], label='w') + plt.plot(time, TAS, label='TAS') + plt.xlabel('time (s)') + plt.ylabel('velocity (m/s)') + plt.legend() + + plt.figure('ang velocities') + plt.plot(time, velocities[:, 3], label='p') + plt.plot(time, velocities[:, 4], label='q') + plt.plot(time, velocities[:, 5], label='r') + plt.xlabel('time (s)') + plt.ylabel('angular velocity (rad/s)') + plt.legend() + + plt.figure('aero angles') + plt.plot(time, alpha, label='alpha') + plt.plot(time, beta, label='beta') + plt.xlabel('time (s)') + plt.ylabel('angle (rad)') + plt.legend() + + plt.figure('2D Trajectory') + plt.plot(position[:, 0], position[:, 1]) + plt.xlabel('X (m)') + plt.ylabel('Y (m)') + plt.legend() diff --git a/examples/example_004.py b/examples/example_004.py new file mode 100644 index 0000000..7146ed9 --- /dev/null +++ b/examples/example_004.py @@ -0,0 +1,166 @@ +# -*- coding: utf-8 -*- +""" +Python Flight Mechanics Engine (PyFME). +Copyright (c) AeroPython Development Team. +Distributed under the terms of the MIT License. + +Example with trimmed aircraft. +The main purpose of this example is to see the evolution of the aircraft after +a longitudinal perturbation (delta doublet). +Trimmed in stationary, horizontal, symmetric, wings level flight. +""" + +import numpy as np +import matplotlib.pyplot as plt + +from pyfme.aircrafts import cessna_310 +from pyfme.models.system import System +from pyfme.utils.trimmer import steady_state_flight_trim +from pyfme.utils.anemometry import calculate_alpha_beta_TAS +from pyfme.environment.isa import atm + + +if __name__ == '__main__': + + # Aircraft parameters. + mass, inertia = cessna_310.mass_and_inertial_data() + + # Initial conditions. + TAS_ = 312 * 0.3048 # m/s + h = 8000 * 0.3048 # m + psi_0 = 3 # rad + x_0, y_0 = 0, 0 # m + + # Trimming. + trim_results = steady_state_flight_trim(cessna_310, h, TAS_, gamma=0, + turn_rate=0) + + lin_vel, ang_vel, theta, phi, alpha_, beta_, control_vector = trim_results + + # Time. + t0 = 0 # s + tf = 30 # s + dt = 1e-2 # s + + time = np.arange(t0, tf, dt) + + # Results initialization. + results = np.zeros([time.size, 12]) + + results[0, 0:3] = lin_vel + results[0, 3:6] = ang_vel + results[0, 6:9] = theta, phi, psi_0 + results[0, 9:12] = x_0, y_0, h + alpha = np.empty_like(time) + alpha[0] = alpha_ + beta = np.empty_like(time) + beta[0] = beta_ + TAS = np.empty_like(time) + TAS[0] = TAS_ + + # Linear Momentum and Angular Momentum eqs. + equations = System(integrator='dopri5', + model='euler_flat_earth', + jac=False) + u, v, w = lin_vel + p, q, r = ang_vel + + equations.set_initial_values(u, v, w, + p, q, r, + theta, phi, psi_0, + x_0, y_0, h) + + _, _, rho, _ = atm(h) + + # Define control vectors. + delta_e, delta_ail, delta_r, delta_t = control_vector + d_e = np.ones_like(time) * delta_e + d_e[np.where(time<2)] = delta_e * 1.30 + d_e[np.where(time<1)] = delta_e * 0.70 + d_a = np.ones_like(time) * delta_ail + d_r = np.ones_like(time) * delta_r + d_t = np.ones_like(time) * delta_t + + attitude = theta, phi, psi_0 + + # Rename function to make it shorter + forces_and_moments = cessna_310.get_forces_and_moments + for ii, t in enumerate(time[1:]): + + forces, moments = forces_and_moments(TAS[ii], rho, alpha[ii], beta[ii], + d_e[ii], 0, d_a[ii], d_r[ii], d_t[ii], + attitude) + + results[ii+1, :] = equations.propagate(mass, inertia, forces, moments, + dt) + + lin_vel = results[ii+1, 0:3] + ang_vel = results[ii+1, 3:6] + attitude = results[ii+1, 6:9] + position = results[ii+1, 9:12] + + alpha[ii+1], beta[ii+1], TAS[ii+1] = calculate_alpha_beta_TAS(*lin_vel) + + _, _, rho, _ = atm(position[2]) + + + velocities = results[:, 0:6] + attitude_angles = results[:, 6:9] + position = results[:, 9:12] + + # PLOTS + plt.close('all') + plt.style.use('ggplot') + + plt.figure('pos') + plt.plot(time, position[:, 0], label='x') + plt.plot(time, position[:, 1], label='y') + plt.plot(time, position[:, 2], label='z') + plt.xlabel('time (s)') + plt.ylabel('position (m)') + plt.legend() + + plt.figure('angles') + plt.plot(time, attitude_angles[:, 0], label='theta') + plt.plot(time, attitude_angles[:, 1], label='phi') + plt.plot(time, attitude_angles[:, 2], label='psi') + plt.xlabel('time (s)') + plt.ylabel('attitude (rad)') + plt.legend() + + plt.figure('velocities') + plt.plot(time, velocities[:, 0], label='u') + plt.plot(time, velocities[:, 1], label='v') + plt.plot(time, velocities[:, 2], label='w') + plt.plot(time, TAS, label='TAS') + plt.xlabel('time (s)') + plt.ylabel('velocity (m/s)') + plt.legend() + + plt.figure('ang velocities') + plt.plot(time, velocities[:, 3], label='p') + plt.plot(time, velocities[:, 4], label='q') + plt.plot(time, velocities[:, 5], label='r') + plt.xlabel('time (s)') + plt.ylabel('angular velocity (rad/s)') + plt.legend() + + plt.figure('aero angles') + plt.plot(time, alpha, label='alpha') + plt.plot(time, beta, label='beta') + plt.xlabel('time (s)') + plt.ylabel('angle (rad)') + plt.legend() + + plt.figure('2D Trajectory') + plt.plot(position[:, 0], position[:, 1]) + plt.xlabel('X (m)') + plt.ylabel('Y (m)') + plt.legend() + + plt.figure() + plt.plot(time, alpha, label='alpha') + plt.plot(time, attitude_angles[:, 0], label='theta') + plt.plot(time, d_e, label='delta_e') + plt.plot(time, velocities[:, 4], label='q') + plt.legend() From 53edaf70a44a96cf8df902000d11dc3f75e56e75 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?=C3=81lex=20S=C3=A1ez?= Date: Sun, 14 Feb 2016 17:49:42 +0100 Subject: [PATCH 11/17] Improved documentation --- src/pyfme/aircrafts/cessna_310.py | 97 +++++++++++++++++++------------ src/pyfme/utils/trimmer.py | 29 ++++----- 2 files changed, 74 insertions(+), 52 deletions(-) diff --git a/src/pyfme/aircrafts/cessna_310.py b/src/pyfme/aircrafts/cessna_310.py index 3db46bb..c9d3044 100644 --- a/src/pyfme/aircrafts/cessna_310.py +++ b/src/pyfme/aircrafts/cessna_310.py @@ -22,13 +22,12 @@ def geometric_data(): Returns ---- - Sw : float - Wing surface (ft2 * 0.3048 * 0.3048 = m2) + Wing surface (m2). c : foat - Mean aerodynamic Chord (ft * 0.3048 = m) + Mean aerodynamic Chord (m). span : float - Wing span (ft * 0.3048 = m) + Wing span (m). References ---------- @@ -98,7 +97,6 @@ def long_aero_coefficients(): Returns ------- - Long_coef_matrix : array_like [ [CL_0, CL_a, CL_de, CL_dih], @@ -206,29 +204,27 @@ def get_aerodynamic_forces(TAS, rho, alpha, beta, delta_e, ih, delta_ail, ---------- rho : float - density (kg/(m3)) + density (kg/(m3). TAS : float - velocity (m/s) - + velocity (m/s). alpha : float - attack angle (rad). + attack angle (rad). beta : float - sideslip angle (rad). + sideslip angle (rad). delta_e : float - elevator deflection (rad). + elevator deflection (rad). ih : float - stabilator deflection (rad). + stabilator deflection (rad). delta_ail : float - aileron deflection (rad). + aileron deflection (rad). delta_r : float - rudder deflection (rad). + rudder deflection (rad). Returns ------- - forces : array_like - 3 dimensional vector with (F_x_s, F_y_s, F_z_s) forces - in stability axes. + 3 dimensional vector with (F_x_s, F_y_s, F_z_s) forces in stability + axes. References ---------- @@ -262,31 +258,28 @@ def get_aerodynamic_moments(TAS, rho, alpha, beta, delta_e, ih, delta_ail, Parameters ---------- - rho : float - density (kg/m3) + density (kg/m3). TAS : float - velocity (m/s) - + velocity (m/s). alpha : float - attack angle (rad). + attack angle (rad). beta : float - sideslip angle (rad). + sideslip angle (rad). delta_e : float - elevator deflection (rad). + elevator deflection (rad). ih : float - stabilator deflection (rad). + stabilator deflection (rad). delta_ail : float - aileron deflection (rad). + aileron deflection (rad). delta_r : float - rudder deflection (rad). + rudder deflection (rad). - returns + Returns ------- - moments : array_like - 3 dimensional vector with (Mxs, Mys, Mzs) forces - in stability axes. + 3 dimensional vector with (Mxs, Mys, Mzs) forces + in stability axes. References ---------- @@ -316,21 +309,17 @@ def get_aerodynamic_moments(TAS, rho, alpha, beta, delta_e, ih, delta_ail, def get_engine_force(delta_t): - """ Calculates forces Parameters ---------- - delta_t : float - trust_lever (0 = 0 Newton, 1 = CT) - + trust_lever (0 = 0 Newton, 1 = CT). returns ------- - engine_force : float - Thrust (N) + Thrust (N). References ---------- @@ -354,7 +343,39 @@ def get_engine_force(delta_t): def get_forces_and_moments(TAS, rho, alpha, beta, delta_e, ih, delta_ail, delta_r, delta_t, attitude): - """ + """Return the total forces and moments including aerodynamics, thrust and + gravity. + + Parameters + ---------- + TAS : float + velocity (m/s). + rho : float + density (kg/m3). + alpha : float + attack angle (rad). + beta : float + sideslip angle (rad). + delta_e : float + elevator deflection (rad). + ih : float + stabilator deflection (rad). + delta_ail : float + aileron deflection (rad). + delta_r : float + rudder deflection (rad). + delta_t : float + Thrust level (between 0-1). + attitude : array_like + Attitude angles: (theta, phi, psi). + + Returns + ------- + forces : array_like + 3 dimensional vector with (F_x_s, F_y_s, F_z_s) forces in body axes. + (N) + moments : array_like + 3 dimensional vector with (Mxs, Mys, Mzs) forces in body axes. (N·m) """ # As stability axes are coincident with wind axes at the moment of # linearization (with alpha=0 and beta=0), stability axes are parallel to diff --git a/src/pyfme/utils/trimmer.py b/src/pyfme/utils/trimmer.py index 5feb2f7..3a1568c 100644 --- a/src/pyfme/utils/trimmer.py +++ b/src/pyfme/utils/trimmer.py @@ -45,20 +45,20 @@ def steady_state_flight_trim(aircraft, h, TAS, gamma=0, turn_rate=0, Returns ------- - lin_vel : - - ang_vel : - - theta : - - phi : - - alpha : - - beta : - - control_vector : - + lin_vel : float array + [u, v, w] air linear velocity body-axes (m/s). + ang_vel : float array + [p, q, r] air angular velocity body-axes (m/s). + theta : float + Pitch angle (rad). + phi : float + Bank angle (rad). + alpha : float + Angle of attack (rad). + beta : float + Sideslip angle (rad). + control_vector : array_like + [delta_e, delta_ail, delta_r, delta_t]. Notes ----- @@ -93,6 +93,7 @@ def steady_state_flight_trim(aircraft, h, TAS, gamma=0, turn_rate=0, args = (h, TAS, gamma, turn_rate, aircraft, dynamic_eqs) + # TODO: pass max deflection of the controls inside aircraft. lower_bounds = (-1, -0.5, -1, -1, -1, 0) upper_bounds = (+1, +0.5, +1, +1, +1, 1) From 69321554763f1c20ea051f17ca679f85fc20433a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?=C3=81lex=20S=C3=A1ez?= Date: Sun, 14 Feb 2016 19:52:50 +0100 Subject: [PATCH 12/17] Added tests for trimming constrians --- examples/example_003.py | 2 +- src/pyfme/utils/trimmer.py | 7 ++++--- 2 files changed, 5 insertions(+), 4 deletions(-) diff --git a/examples/example_003.py b/examples/example_003.py index 7a5c646..e4d22a9 100644 --- a/examples/example_003.py +++ b/examples/example_003.py @@ -33,7 +33,7 @@ # Trimming. trim_results = steady_state_flight_trim(cessna_310, h, TAS_, gamma=0.02, - turn_rate=0) + turn_rate=0.02) lin_vel, ang_vel, theta, phi, alpha_, beta_, control_vector = trim_results diff --git a/src/pyfme/utils/trimmer.py b/src/pyfme/utils/trimmer.py index 3a1568c..ce21540 100644 --- a/src/pyfme/utils/trimmer.py +++ b/src/pyfme/utils/trimmer.py @@ -157,11 +157,11 @@ def turn_coord_cons1(turn_rate, alpha, beta, TAS, gamma=0): c = 1 + G**2 * cos(beta)**2 sq = sqrt(c * (1 - b**2) + G**2 * sin(beta)**2) - var1 = b**2 * (1 + c * tan(alpha)**2) - temp = (a - b)**2 + b * tan(alpha) * sq / a**2 - b**2 * var1 + num = (a-b**2) + b * tan(alpha) * sq + den = a**2 - b**2 * (1 + c * tan(alpha)**2) - phi = atan(G * cos(beta) / cos(alpha) * temp) + phi = atan(G * cos(beta) / cos(alpha) * num / den) return phi @@ -175,6 +175,7 @@ def turn_coord_cons2(turn_rate, alpha, TAS): G = turn_rate * TAS / g0 phi = G / cos(alpha) + phi = atan(phi) return phi From 6482749e820807bd4ecff401ecef4df74e35d4e5 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?=C3=81lex=20S=C3=A1ez?= Date: Sun, 14 Feb 2016 19:57:20 +0100 Subject: [PATCH 13/17] Updated environment to take into account scipy dependency for least_squares --- environment.yml | 1 + 1 file changed, 1 insertion(+) diff --git a/environment.yml b/environment.yml index d00d77b..5c88359 100644 --- a/environment.yml +++ b/environment.yml @@ -1,3 +1,4 @@ name: pyfme dependencies: - numpy >=1.7 + - scipy >=0.17.0 From 9cf8870acdec59b874653a1a89c943715cbed3a1 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?=C3=81lex=20S=C3=A1ez?= Date: Sun, 14 Feb 2016 20:04:00 +0100 Subject: [PATCH 14/17] Actually added tests '=) --- src/pyfme/utils/tests/test_trimmer.py | 84 +++++++++++++++++++++++++++ 1 file changed, 84 insertions(+) create mode 100644 src/pyfme/utils/tests/test_trimmer.py diff --git a/src/pyfme/utils/tests/test_trimmer.py b/src/pyfme/utils/tests/test_trimmer.py new file mode 100644 index 0000000..d5583f6 --- /dev/null +++ b/src/pyfme/utils/tests/test_trimmer.py @@ -0,0 +1,84 @@ +# -*- coding: utf-8 -*- +""" +Python Flight Mechanics Engine (PyFME). +Copyright (c) AeroPython Development Team. +Distributed under the terms of the MIT License. + +Test functions for trimmer. +""" + +from numpy.testing import assert_almost_equal + +from pyfme.utils.trimmer import (rate_of_climb_cons, turn_coord_cons1, + turn_coord_cons2) + + +def test_rate_of_climb_cons(): + + alpha = 0.05 + beta = 0.01 + phi = 0.3 + gamma = 0.2 + + expected_theta = 0.25072488304787743 + + theta = rate_of_climb_cons(gamma, alpha, beta, phi) + + assert_almost_equal(theta, expected_theta) + + +def test_turn_coord_cons2(): + + turn_rate = 0.05 + alpha = 0.05 + TAS = 100 + + expected_phi = 0.4718708835828995 + + phi = turn_coord_cons2(turn_rate, alpha, TAS) + + assert_almost_equal(phi, expected_phi) + + +def test_turn_coord_cons1_against_2(): + + turn_rate = 0.05 + alpha = 0.05 + TAS = 100 + beta = 0 + gamma = 0 + + expected_phi = turn_coord_cons2(turn_rate, alpha, TAS) + phi = turn_coord_cons1(turn_rate, alpha, beta, TAS, gamma) + + assert_almost_equal(phi, expected_phi) + + +def test_turn_coord_cons1_small_gamma(): + + turn_rate = 0.05 + alpha = 0.05 + TAS = 100 + beta = 0.01 + gamma = 0 + + expected_phi = 0.4719539221164456 + + phi = turn_coord_cons1(turn_rate, alpha, beta, TAS, gamma) + + assert_almost_equal(phi, expected_phi) + + +def test_turn_coord_cons1_big_gamma(): + + turn_rate = 0.05 + alpha = 0.05 + TAS = 100 + beta = 0.01 + gamma = 0.2 + + expected_phi = 0.4766120129721381 + + phi = turn_coord_cons1(turn_rate, alpha, beta, TAS, gamma) + + assert_almost_equal(phi, expected_phi) From dab4fe07185efc6706aa0fecfebc59c93e11cdbf Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?=C3=81lex=20S=C3=A1ez?= Date: Sun, 14 Feb 2016 20:07:29 +0100 Subject: [PATCH 15/17] Modified travis configuration script to install scipy --- .travis.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.travis.yml b/.travis.yml index 8b82ea1..9ee456b 100644 --- a/.travis.yml +++ b/.travis.yml @@ -22,7 +22,7 @@ before_install: - source activate test-environment install: - - conda install numpy pytest + - conda install numpy pytest scipy - pip install -e . script: From 072a884df7c62eeb8288da988f5c9aa46a57df94 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?=C3=81lex=20S=C3=A1ez?= Date: Wed, 23 Mar 2016 10:45:27 +0100 Subject: [PATCH 16/17] Changed print style for results --- examples/trim_example_001.py | 26 ++++++++++++++++---------- 1 file changed, 16 insertions(+), 10 deletions(-) diff --git a/examples/trim_example_001.py b/examples/trim_example_001.py index 3193cff..5e2ccbd 100644 --- a/examples/trim_example_001.py +++ b/examples/trim_example_001.py @@ -21,13 +21,19 @@ beta = results[4] lin_vel = wind2body((TAS, 0, 0), alpha, beta) -print('') -print('Results:') -print("lin_vel = {}".format(lin_vel)) -print("ang_vel = {}".format(results[0])) -print("theta = {}".format(results[1])) -print("phi = {}".format(results[2])) -print("alpha = {}".format(results[3])) -print("beta = {}".format(results[4])) -print("control = {}".format(results[5])) -print('') +print(""" +Results +------- +Linear velocity: {lin_vel} (m/s) +Angular velocity: {ang_vel} (rad/s) +Theta, Phi: {angles} (rad) +Alpha, Beta: {wind_angles} (rad) +Control: {control} +""".format( + lin_vel=lin_vel, + ang_vel=results[1], + angles=results[2:4], + wind_angles=results[4:6], + control=results[6] + ) + ) From 2f0eba8a7e2e8cc8ed3ee0af9cb61c4ad8f0e84a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?=C3=81lex=20S=C3=A1ez?= Date: Wed, 23 Mar 2016 10:48:19 +0100 Subject: [PATCH 17/17] Changed name of constrain functions for turns --- src/pyfme/utils/tests/test_trimmer.py | 22 +++++----- src/pyfme/utils/trimmer.py | 58 +++++++++++++-------------- 2 files changed, 41 insertions(+), 39 deletions(-) diff --git a/src/pyfme/utils/tests/test_trimmer.py b/src/pyfme/utils/tests/test_trimmer.py index d5583f6..507f1c4 100644 --- a/src/pyfme/utils/tests/test_trimmer.py +++ b/src/pyfme/utils/tests/test_trimmer.py @@ -9,8 +9,8 @@ from numpy.testing import assert_almost_equal -from pyfme.utils.trimmer import (rate_of_climb_cons, turn_coord_cons1, - turn_coord_cons2) +from pyfme.utils.trimmer import (rate_of_climb_cons, turn_coord_cons, + turn_coord_cons_horizontal_and_small_beta) def test_rate_of_climb_cons(): @@ -27,7 +27,7 @@ def test_rate_of_climb_cons(): assert_almost_equal(theta, expected_theta) -def test_turn_coord_cons2(): +def test_turn_coord_cons_horizontal_and_small_beta(): turn_rate = 0.05 alpha = 0.05 @@ -35,7 +35,7 @@ def test_turn_coord_cons2(): expected_phi = 0.4718708835828995 - phi = turn_coord_cons2(turn_rate, alpha, TAS) + phi = turn_coord_cons_horizontal_and_small_beta(turn_rate, alpha, TAS) assert_almost_equal(phi, expected_phi) @@ -48,13 +48,15 @@ def test_turn_coord_cons1_against_2(): beta = 0 gamma = 0 - expected_phi = turn_coord_cons2(turn_rate, alpha, TAS) - phi = turn_coord_cons1(turn_rate, alpha, beta, TAS, gamma) + expected_phi = turn_coord_cons_horizontal_and_small_beta(turn_rate, + alpha, + TAS) + phi = turn_coord_cons(turn_rate, alpha, beta, TAS, gamma) assert_almost_equal(phi, expected_phi) -def test_turn_coord_cons1_small_gamma(): +def test_turn_coord_cons_small_gamma(): turn_rate = 0.05 alpha = 0.05 @@ -64,12 +66,12 @@ def test_turn_coord_cons1_small_gamma(): expected_phi = 0.4719539221164456 - phi = turn_coord_cons1(turn_rate, alpha, beta, TAS, gamma) + phi = turn_coord_cons(turn_rate, alpha, beta, TAS, gamma) assert_almost_equal(phi, expected_phi) -def test_turn_coord_cons1_big_gamma(): +def test_turn_coord_cons_big_gamma(): turn_rate = 0.05 alpha = 0.05 @@ -79,6 +81,6 @@ def test_turn_coord_cons1_big_gamma(): expected_phi = 0.4766120129721381 - phi = turn_coord_cons1(turn_rate, alpha, beta, TAS, gamma) + phi = turn_coord_cons(turn_rate, alpha, beta, TAS, gamma) assert_almost_equal(phi, expected_phi) diff --git a/src/pyfme/utils/trimmer.py b/src/pyfme/utils/trimmer.py index ce21540..c225d97 100644 --- a/src/pyfme/utils/trimmer.py +++ b/src/pyfme/utils/trimmer.py @@ -32,33 +32,33 @@ def steady_state_flight_trim(aircraft, h, TAS, gamma=0, turn_rate=0, Parameters ---------- - aircraft : aircraft class - Aircraft class with methods get_forces, get_moments - h : float - Geopotential altitude for ISA (m). - TAS : float - True Air Speed (m/s). - gamma : float, optional - Flight path angle (rad). - turn_rate : float, optional - Turn rate, d(psi)/dt (rad/s). + aircraft : aircraft class + Aircraft class with methods get_forces, get_moments + h : float + Geopotential altitude for ISA (m). + TAS : float + True Air Speed (m/s). + gamma : float, optional + Flight path angle (rad). + turn_rate : float, optional + Turn rate, d(psi)/dt (rad/s). Returns ------- - lin_vel : float array - [u, v, w] air linear velocity body-axes (m/s). - ang_vel : float array - [p, q, r] air angular velocity body-axes (m/s). - theta : float - Pitch angle (rad). - phi : float - Bank angle (rad). - alpha : float - Angle of attack (rad). - beta : float - Sideslip angle (rad). - control_vector : array_like - [delta_e, delta_ail, delta_r, delta_t]. + lin_vel : float array + [u, v, w] air linear velocity body-axes (m/s). + ang_vel : float array + [p, q, r] air angular velocity body-axes (m/s). + theta : float + Pitch angle (rad). + phi : float + Bank angle (rad). + alpha : float + Angle of attack (rad). + beta : float + Sideslip angle (rad). + control_vector : array_like + [delta_e, delta_ail, delta_r, delta_t]. Notes ----- @@ -123,7 +123,7 @@ def steady_state_flight_trim(aircraft, h, TAS, gamma=0, turn_rate=0, if abs(turn_rate) < 1e-8: phi = 0 else: - phi = turn_coord_cons1(turn_rate, alpha, beta, TAS, gamma) + phi = turn_coord_cons(turn_rate, alpha, beta, TAS, gamma) theta = rate_of_climb_cons(gamma, alpha, beta, phi) @@ -140,7 +140,7 @@ def steady_state_flight_trim(aircraft, h, TAS, gamma=0, turn_rate=0, return lin_vel, ang_vel, theta, phi, alpha, beta, control_vector -def turn_coord_cons1(turn_rate, alpha, beta, TAS, gamma=0): +def turn_coord_cons(turn_rate, alpha, beta, TAS, gamma=0): """Calculates phi for coordinated turn. """ @@ -166,7 +166,7 @@ def turn_coord_cons1(turn_rate, alpha, beta, TAS, gamma=0): return phi -def turn_coord_cons2(turn_rate, alpha, TAS): +def turn_coord_cons_horizontal_and_small_beta(turn_rate, alpha, TAS): """Calculates phi for coordinated turn given that gamma is equal to cero and beta is small (beta << 1). """ @@ -181,7 +181,7 @@ def turn_coord_cons2(turn_rate, alpha, TAS): def rate_of_climb_cons(gamma, alpha, beta, phi): - """ + """Calculates theta for the given ROC, wind angles, and roll angle. """ a = cos(alpha) * cos(beta) b = sin(phi) * sin(beta) + cos(phi) * sin(alpha) * cos(beta) @@ -208,7 +208,7 @@ def func(trimmed_params, h, TAS, gamma, turn_rate, aircraft, dynamic_eqs): if abs(turn_rate) < 1e-8: phi = 0 else: - phi = turn_coord_cons1(turn_rate, alpha, beta, TAS, gamma) + phi = turn_coord_cons(turn_rate, alpha, beta, TAS, gamma) theta = rate_of_climb_cons(gamma, alpha, beta, phi)