Simulation of a Mathematical Model of an Aircraft Using Parallel Techniques (MPI and GPU)

Abstract

This chapter is focused on creation, accuracy, and simulation of two-parameter control of a mathematical model of motion of aircraft in a flying simulator. We are discussing many of the important advances in applied aircraft modeling. Modeling on various computer architectures (central, distributed, parallel) has an impact on a structure of a mathematical model of aircraft. An important part is the way of description of a numerical method and its accuracy, use of distributed memory system, and shared memory system are presented in the chapter. Motivation of this research is implementation of the general-purpose message passing interface and graphics processing units as inexpensive arithmetic-processing units bring a relevant amount of computing power to desktop personal computers. The chapter is focused on exploitation of parallel techniques of simulation features, computation time of parallel methods of implementation, and improved simulation of a continuous mathematical model of aircraft motion in a flying simulator. The use and application of modeling methods and parallel simulation techniques determine the structure of the mathematical model used in the flying simulator. The effectiveness of our solution is confirmed by providing simulation results obtained by two-parameter control of the mathematical model of aircraft motion.

1. Introduction

Mathematical models are mostly used in natural sciences (physics, chemistry, earth science) and engineering disciplines (computer science, biological science, genetic engineering), as well as in social sciences. Substantial stages of implementation (creation, design, simulation, and visualization) of a continuous mathematical model of aircraft motion in a flying simulator were described in many papers. The one-parameter control of a distributed mathematical model of motion of aircraft with all phases of implementation in a flying simulator was described in the published paper [1, 2].

Motivation of this research is using a mathematical model of an object, which in practice and the real world often can be determined analytical otherwise. For this reason, it is necessary to design mathematical models characterized for this problem statement. The more accurate the mathematical model of the object we create, the more accurate the results of its simulation.

The mathematical model of motion of aircraft can be simulated using central computer architecture, which can be based on single-processor systems. The simulation of structured mathematical model of an aircraft can be implemented using central computer architecture, a distributed memory system (DMS) using a message interface (MPI) architecture; they are identified as node computers or shared memory system (SMS) and graphics processing unit (GPU) architectures [3]. The parallel computer architecture may be based on multiprocessors; each processor has a multicore architecture. Multiprocessor systems are computationally more efficient than such systems compared with the central computer architecture, see [4].

As known from specialized literature, creation of mathematical models of motion of aircraft consists of the following phases: definition of a physical base for creation of a mathematical model, selection of notation of a mathematical model, implementation of a mathematical model and others. The Laplace transformation, computation of parameters of aircraft for a selected flight phase, determination and computation of coefficients, numerical integration, programming, and simulation on a computer, etc. are used in our problem statement and proposed solution.

According to current knowledge, the effort made to create mathematical models of motion is aimed at improvement and extension of a field of view of simulation of mathematical models on the different computer architectures. The solution of large and complex mathematical models of aircraft problems leads to considerable demand for parallel computing and strategies [5]. We can confirm the most frequent use of the following phases: creation, design, implementation, and simulation of mathematical models. The parallel computing tools guarantee an impressive reduction of computing time [6].

A block structure mathematical model of motion of aircraft can be integrated numerically using parallel computer architecture. This architecture can be motivation for using a central processing unit (CPU)-based implementation with multiprocessor systems and a graphics processing unit (GPU); see [7]. Likewise, a computer system fitted with a GPU accelerator and a software application designed and created for this architecture is also suitable for simulation of a block mathematical model of aircraft motion. From our point of view, the following three principles are important for proposed solution:

• representation by solution architecture, i.e., the migration of codes of a mathematical model of motion of aircraft toward parallel computers.

• low-cost distributed platforms, combined suitable software tools that open very demanding perspectives.

• use GPU to speed up calculations.

Accordingly, in the chapter, we explain how to implement a compute unified device architecture (CUDA) algorithm of computation of a block structure mathematical model of motion of aircraft so as to achieve peak speed. The efficiency is demonstrated by discussing two-parameter control of aircraft, throttle control stick (fuel supply), and aircraft elevator control stick (angle of attack).

2. Description of aircraft mathematical model

The unstable characteristics such as aero elasticity impact, fuel density, changing geometry of aircraft, and some other parameters support complexity of our design mathematical models in simulation. The mathematical models of motion interact with intervening pilot’s control of aircraft and real equipment respond to the pilot’s interventions, data on equipment the flying aircraft are observed by the pilot. We can use a continuous simulation method to solve differential equations in mathematical models of motion of aircraft created this way. Aerospace engineers often use Newton’s laws of motion in design and creation of a mathematical model of motion of aircraft described by differential operators.

According to Newton’s third law, there is always the same opposite reaction to every action, or the interaction of two bodies on each other always leads to opposite parts. The basic system of equations has the form [8]:

$${\dot{x}}_1+f_i\left(x_1,\dots ,x_n;u_1,\dots ,u_m;{\xi }_1,\dots ,{\xi }_{\gamma }\right)=0,\left(i=0,1,\dots ,p\right)$$E1

where x₁, …, x_n are object coordinates, u₁, …, u_m are elements of control, and ζ₁, …, ζ_γ are failure functions. The problems of simulation and synchronization of mathematical models on computers is know, where the accuracy of a linear motion model is used [9]. The accuracy between the transient response results of nonlinear mathematical models compared with linear mathematical models from in terms of human accuracy can be neglected [10].

2.1 Equations numerical integration of mathematical model

The proposed solution based on the presented analysis of the state model space does not contain any discontinuity in either f_i(x, u, t) or in any of higher derivatives, x_i(t) itself is a continuous function of time. Currently used methods of numerical integration are used to find numerical approximations to solutions of ordinary differential equations (ODEs). The name, also known as numerical integration, represents a broad family of algorithms for calculating the numerical value of a particular integral. The methods of numerical integration can generally be described as combining evaluations of the integrand to get an approximation to the integral. The integrand is evaluated on a finite set of points called integration points, and the weighted sum of these values is used to approximate the integral [11]. This method of calculation is often used for polynomial functions, and these are suitably used on computers. Our proposed solution deals with creation of such functions, i.e., integration of mathematical models.

Two most basic algorithms, Euler (forward and backward) and Taylor series, are motivation for this research in our solution for numerical integration. A numerical stability domain is introduced as a pillar to characterize an integration algorithm and a general procedure to find the numerical stability domain of any integration scheme, see [11]. The mathematical model in state space has the following form:

$$\dot{\mathbf{x}}\left(t\right)=\mathbf{f}\left(\mathbf{x}\left(t\right),\mathbf{u}\left(t\right),t\right),$$E2

where x is a state vector, u is an input vector, and t represents time, with a set of initial conditions:

$${\mathbf{x}}_{\left(t=t_0\right)}={\mathbf{x}}_{\mathbf{0}}.$$E3

Let x_i(t) represents the ith state trajectory as a function of simulated time t. As long as the state-space model does not contain any discontinuity in either f_i(x, u, t) or any of higher derivatives, x_i(t) is itself a continuous function of time. Such function can be approximated by any desired precision by a Taylor series expansion of any given point along its trajectory. Let t^* denote a point in time, about which we wish to approximate the trajectory using Taylor series, and let t^* + h be the point in time, at which we wish to evaluate the approximation [11]. The trajectory at that point can then be given as follows:

$$x_i\left(t^{\ast }+h\right)=x_i\left(t^{\ast }\right)+\frac{{\mathit{dx}}_i\left(t^{\ast }\right)}{\mathit{dt}}\cdot h+\frac{d^2{x}_i\left(t^{\ast }\right)}{{\mathit{dt}}^2}\cdot \frac{h^2}{2!}+\cdots .$$E4

Different integration algorithms differ in how they approximate the higher-order derivatives and in the number of Taylor series members that they take into account in the approximation [10].

If the term n + 1 of the Taylor series is considered, the approximation accuracy of the second derivative d²x_i(t^*)/d²t = df_i(t^*)/dt should be of order n − 2, since this factor is multiplied by h². The accuracy of the third state derivatives should be of the order n − 3, since this factor is multiplied by h³, etc. In this way, the approximation is correct up to hⁿ. N is therefore called the approximation order of the integration method, see [11].

Many engineering simulation applications require a global relative accuracy of approximately 0.002. If the local integration error is of size e_l, then the per-unit-step integration error assumes the value of e_p.u.s = e_l/h. The global integration error is proportional to the per-unit-step integration error, as long as the integration error does not accumulate excessively across multiple steps.

In accordance with the previously made observation, this corresponds to an algorithm with an approximation order of h⁴ for the local integration error. We should require, for example, a local accuracy of 0.0001 [11]. In a digital computer, a real number can only be represented with a finite precision, 32-bit numbers, and double precision for 64-bit numbers. This type of error is called round-off error. It occurs in one of the two main general formats that have become common and called floating point. The most common problems due to a rounding error occur when multiple steps are rounded in each step.

3. An appropriate structured mathematical model of motion in state space

These practical requirements are determined to using linear models in a process of analyzing general processes [12]. For writing of mathematical models of motion of aircraft in a simulator, we can use state-space description. We define a linear controlled unobserved dynamic system based on the current state, see [13]:

$$\begin{array}{l}\dot{\mathbf{x}}\left(t\right)=\mathbf{Ax}\left(t\right)+\mathbf{Bu}\left(t\right),\\ \mathbf{y}\left(t\right)=\mathbf{Cx}\left(t\right),\\ \mathbf{x}\left(t\right)={\mathbf{x}}_0\end{array}$$E5

where A, B, C, x, u, and y have dimensions (n x n), (r x n), (l x n), (n x l), (m x l), and (r x l) are defined matrices, see the description of items [2]. The task easier that we will focus on the control object, the first and second equation from the system (5) can have a general shape:

$$\begin{array}{l}\dot{\mathbf{x}}=\mathbf{Ax}+\mathbf{Bu},\\ \mathbf{y}=\mathbf{Cx}.\end{array}$$E6

Eq. (6) in matrix form represents a complex aircraft dynamic system of a mathematical model of aircraft motion comprising in a process of simulation of 11 state variable sensors of information, 18 state variables-coordinates of performing elements in the system. These are two parts and the remaining 38 state variables are divided into variables that represent unmeasured noise and sensor disturbances [9]. The differential equations of mathematical model of aircraft motion have two-parameter control the form [14]:

$$\begin{array}{l}\mathrm{\Delta }\dot{V}+a_x^V\mathrm{\Delta }V+a_x^{\alpha }\mathrm{\Delta }\alpha +a_x^{\theta }\mathrm{\Delta }\theta =a_x^{\mathit{\delta M}}\mathrm{\Delta }{\delta }_M\\ \mathrm{\Delta }\theta +a_y^V\mathrm{\Delta }V+a_y^{\alpha }\mathrm{\Delta }\alpha +a_y^{\theta }\mathrm{\Delta }\theta =a_y^{\mathit{\delta M}}\mathrm{\Delta }{\delta }_V\\ \mathrm{\Delta }{\dot{\omega }}_z+a_{\mathit{mz}}^V\mathrm{\Delta }V+a_{\mathit{mz}}^{\alpha }\mathrm{\Delta }\alpha +a_{\mathit{mz}}^{\dot{\alpha }}\mathrm{\Delta }\theta +a_{\mathit{mz}}^{\mathit{\omega z}}\mathrm{\Delta }{\omega }_z=a_{\mathit{mz}}^{\mathit{\delta V}}\mathrm{\Delta }{\delta }_V\\ \mathrm{\Delta }\upsilon ={\omega }_z,\mathrm{\Delta }\upsilon =\mathrm{\Delta }\theta +\mathrm{\Delta }\alpha .\end{array}$$E7

Four parts represent a state matrix of mathematical models of motion of aircraft expressed by the state vector n = 4. Coefficients c_i and e_j represent aerodynamic parameters, their computation, in accordance with (Eq. (7)), we get [14]:

$$\begin{array}{l}\mathrm{\Delta }\dot{V}+a_x^V\mathrm{\Delta }V+a_x^{\alpha }\mathrm{\Delta }\alpha +a_x^{\theta }\mathrm{\Delta }\theta +a_x^H\mathrm{\Delta }H=a_x^{\mathit{\delta M}}\mathrm{\Delta }{\delta }_M\\ \mathrm{\Delta }\theta +a_y^V\mathrm{\Delta }V+a_y^{\alpha }\mathrm{\Delta }\alpha +a_y^{\theta }\mathrm{\Delta }\theta +a_y^H\mathrm{\Delta }H=a_y^{\mathit{\delta M}}\mathrm{\Delta }{\delta }_V\\ \mathrm{\Delta }\dot{H}-sin\left({\theta }_0\right)\mathrm{\Delta }V-cos\left({\theta }_0\right)\mathrm{\Delta }V=0;\mathrm{\Delta }\dot{\upsilon }={\omega }_z,\mathrm{\Delta }\upsilon =\mathrm{\Delta }\theta +\mathrm{\Delta }\alpha .\end{array}$$E8

Items ΔV, Δα, Δθ, ΔH, Δδ_M, and Δδ_V are defined in [2]. The following items: a^V_x ΔV, a^α_x Δα, a^θ_x Δθ, a^H_x ΔH for flight speed ΔV, a^V_y ΔV, a^α_y Δα, a^θ_y Δθ, a^H_y ΔH for thumb angle, and other parameters of aircraft motion. Transfer functions have the form [14]:

$$\begin{array}{l}\mathrm{\Delta }V\left(s\right)=-G_{V/\mathit{\delta M}}\left(s\right)\mathrm{\Delta }{\delta }_M\left(s\right)-G_{V/\mathit{\delta V}}\left(s\right)\mathrm{\Delta }{\delta }_V\left(s\right)\\ \mathrm{\Delta }\alpha \left(s\right)=-G_{\alpha /\mathit{\delta M}}\left(s\right)\mathrm{\Delta }{\delta }_M\left(s\right)-G_{\alpha /\mathit{\delta V}}\left(s\right)\mathrm{\Delta }{\delta }_V\left(s\right)\\ G_{V/\mathit{\delta M}}\left(s\right)=a_x^{\mathit{\delta M}}\left({\mathrm{\Delta }}_{11/}\mathrm{\Delta }\right);G_{V/\mathit{\delta V}}\left(s\right)=a_y^{\mathit{\delta M}}\left({\mathrm{\Delta }}_{21/}\mathrm{\Delta }\right)+a_{\mathit{mz}}^{\mathit{\delta V}}\left({\mathrm{\Delta }}_{31/}\mathrm{\Delta }\right);\\ G_{\alpha /\mathit{\delta M}}\left(s\right)=a_x^{\mathit{\delta M}}\left({\mathrm{\Delta }}_{12/}\mathrm{\Delta }\right);G_{\alpha /\mathit{\delta V}}\left(s\right)=a_y^{\mathit{\delta V}}\left({\mathrm{\Delta }}_{22/}\mathrm{\Delta }\right)+a_{\mathit{mz}}^{\mathit{\delta V}}\left({\mathrm{\Delta }}_{32/}\mathrm{\Delta }\right).\end{array}$$E9

where Δ_ij is the subdeterminant of the ith row and the jth column. Our problem and the proposed solution are described in the following. The first member (Eq. (9)) of the system is employed to describe the change of the aircraft speed ΔV; also the equation describing displacement of the speed aircraft depending on the displacement of throttle engine lever G_V/δM (s) and the displacement of elevator angle G_V/δV (s) is employed. The items of second (Eq. (9)) describe the change of aircraft attack of angle of the system G_α/δM (s) and G_α/δV (s). The third row of the system (Eq. (9)) defines how to compute these changes of two parameters for a change of aircraft speed, and fourth row of the system (Eq. (9)) defines how to compute these changes of two parameters for a change of an aircraft angle of attack. For the meaning of coefficients a^δM_x, a^δM_y, a^δV_y, a^δV_mz and ∆_ij(s), see the paragraph after (Eq. (20)). In our problem solution, the flight of the aircraft is steady and without any random interferences (wind, storm, or other outer interferences) [12].

3.1 Structure-defined mathematical model of an aircraft motion in a simulator

At present, matrix notation of a mathematical model in the state space, which we mentioned in the previous section, is very often used. We should express the first equation from (Eq. (6)) next [15]:

$$\left(\begin{array}{c}{\dot{x}}_1\\ {\dot{x}}_2\\ {\dot{x}}_3\\ {\dot{x}}_4\end{array}\right)=\left(\begin{array}{cccc}{a}_{11}& a_{12}& a_{13}& a_{14}\\ a_{21}& a_{22}& a_{23}& a_{24}\\ a_{31}& a_{32}& a_{33}& a_{34}\\ a_{41}& a_{42}& a_{43}& a_{44}\end{array}\right)\left(\begin{array}{c}{x}_1\\ x_2\\ x_3\\ x_4\end{array}\right).$$E10

The described solution is relatively simple due to description of the analyzed system using equations (Eq. (10)), see [5] for more. Let us divide the mathematical model of aircraft motion into four subsystems. The first will be a subsystem of state variable sensors, the second and third will be subsystems of power elements (two systems), the fourth will measure noise and system faults. Thus, the subsystems have the order n₁ = 1, n₂ = 1, n₃ = 1, n₄ = 1. The state space is divided into four parts:

$$\mathbf{x}={\left(x_1,x_2,x_3,x_4\right)}^T=\left(x_1{x}_1,x_1{x}_2,x_1{x}_3,\dots ,x_4{x}_4\right),$$E11

where x_ix_j represent the items of a state vector. If i represents order relevancy n, i.e., the number of the subsystem, then j – stands for sequential number of the item in the given subsystem. The research architecture of the discussed system matrix A in the state space can be formed after multiplying the following shape, see [15]:

$$\begin{array}{c}{\dot{x}}_{11}=a_{11}{x}_1{x}_1+a_{11}{x}_1{x}_2+a_{11}{x}_1{x}_3+a_{11}{x}_1{x}_4,\\ {\dot{x}}_{12}=a_{12}{x}_1{x}_1+a_{12}{x}_1{x}_2+a_{12}{x}_1{x}_3+a_{12}{x}_1{x}_4,\\ ⋮\\ {\dot{x}}_{44}=a_{44}{x}_1{x}_1+a_{44}{x}_1{x}_2+a_{44}{x}_1{x}_3+a_{44}{x}_1{x}_4,\end{array}$$E12

Then we divide 16 subsystems into four parts called isolated subsystems. The above decomposition process is performed by matrix and vector operations. The sixteen blocks A₁₁, A₁₂, …, A₄₃, A₄₄ of matrix A are marked in the order, where:

$$\begin{array}{c}{\dot{x}}_{11}=A_{11}{x}_{11}+A_{11}{x}_{12}+A_{11}{x}_{13}+A_{11}{x}_{14},\\ {\dot{x}}_{12}=A_{12}{x}_{11}+A_{12}{x}_{12}+A_{12}{x}_{13}+A_{12}{x}_{14},\\ ⋮\\ {\dot{x}}_{44}=A_{44}{x}_{11}+A_{44}{x}_{12}+A_{44}{x}_{13}+A_{44}{x}_{14}.\end{array}$$E13

The mathematical description of isolated subsystems has then the following form:

$${\dot{x}}_1=A_{11}^{\text{'}}{x}_1,{\dot{x}}_2=A_{22}^{\text{'}}{x}_2,{\dot{x}}_3=A_{33}^{\text{'}}{x}_3,{\dot{x}}_4=A_{44}^{\text{'}}{x}_4,$$E14

where:

$$A_{11}^{\text{'}}=\left(\begin{array}{c}{A}_{11}\\ A_{12}\\ A_{13}\\ A_{14}\end{array}\right),A_{22}^{\text{'}}=\left(\begin{array}{c}{A}_{21}\\ A_{22}\\ A_{23}\\ A_{24}\end{array}\right),A_{33}^{\text{'}}=\left(\begin{array}{c}{A}_{31}\\ A_{32}\\ A_{33}\\ A_{34}\end{array}\right),A_{44}^{\text{'}}=\left(\begin{array}{c}{A}_{41}\\ A_{42}\\ A_{43}\\ A_{44}\end{array}\right).$$E15

The mutual relations between the first and second isolated subsystems are described by l₁₂(x) meaning that the equation of the first and second isolated subsystem is:

$$l_{12}\left(x\right)=A_{11}^{\text{'}}\left(\begin{array}{c}{x}_1\\ x_2\\ 0\\ 0\end{array}\right)={\left(A_{11},A_{12},A_{13},A_{14}\right)}^T\left(\begin{array}{c}{x}_1\\ x_2\\ 0\\ 0\end{array}\right).$$E16

The state of solution method is known and well described in the reference [5]. The capability of being decomposed can be calculated by means of incidental matrices that can be utilized in cases when the mathematical model of the system is defined. The units in random matrix are positioned as its elements are permuted or transformed so that they can be placed diagonally.

4. Structured mathematical model of speed and angle of attack aircraft in a flight simulator

Current computer technologies allow formulation and solutions to a new intricate problem, which depends on construction of structured mathematical models of aircraft and methods of their solution. Achieving an accurate solution to this problem requires simultaneous solution of the whole complex of physical and geometric problems, which is based on considerable computational resources [16]. Our proposed solution is supported by the above approach, where there are, complex requirements for properties and parameters affecting the design process of a structured mathematical model of aircraft motion on computers.

To create a system of differential equations, it is necessary to know the aerodynamic coefficients, mathematical model of aircraft systems, and other parameters of the aircraft [1, 2]. The contribution in Laplace transform creates a structured mathematical model of aircraft motion in a flight simulator, see [5, 15].

Action of an aircraft elevator control stick angle and throttle control stick displacement on aircraft motion—two control parameters (Figure 1) are discussed below. Figure 1 is a block diagram of already mentioned two-parameter control of a structured mathematical model of motion of aircraft on a flight simulator, which will be the subject of research.

Figure 1.

Block diagram of a structured mathematical model of control of aircraft motion – speed and angle of attack.

The process of simulation experiments requires finding solutions for individual tasks defined at the beginning of the previous subchapter. The simulation program expresses the time factor in the form of an explicit simulation time and changes the state of model. Two opposite methods of working with simulation time are implemented in the programming technique. In the constant time step method, it is necessary to select a constant real-time interval, which is considered as a unit of simulation time. The structured mathematical models of all simulators need to process the simulation run—the simulating step during this time unit (Figure 2).

Figure 2.

Parallel simulation algorithm of simulators’ models.

4.1 Structure-defined mathematical model of the aircraft speed in a flight simulator

The input values are represented by a change of an aircraft elevator control stick angle and a throttle control stick of engine. The values are adjusted as required and forwarded to the input of the structured mathematical model of aircraft motion-speed. The aircraft speed displacement equation ΔV(s) defines a change in fuel supply and the change in the aircraft elevator angle [14]:

$$\mathrm{\Delta }V\left(s\right)=-G_{V/\mathit{\delta M}}\left(s\right)\mathrm{\Delta }{\delta }_M\left(s\right)-G_{V/\mathit{\delta V}}\left(s\right)\mathrm{\Delta }{\delta }_V\left(s\right),$$E17

where G_V/δM (s) defines a mathematical model of speed—transfer function for fuel supply, Δδ_M (s) is an input function for fuel supply in Laplace transform, G_V/δV (s) is a mathematical model of speed—transfer function for an aircraft lift, and Δδ_V (s) means the input function of an aircraft pitch angle in the Laplace transform. Derivation of equations of a structured mathematical model of speed increment is a suitable description of the problem. Information about mathematical solution of these equations is known [13]. The expression $A=s^4+1.134s^3+62.798s^2+28.659s+4.093,$ is the same in all equations. The items for l₁₂(x) defined by the Eq. (16) are as follows:

$$A_{11}=5\frac{s^3+1.12s^2+62.782s+25.32}{A},x_1=\mathrm{\Delta }{\delta }_M\left(s\right).$$E18

$$A_{12}=\frac{-0.11\cdot \left(9.81s+620.973\right)-0.42\cdot \left(-9.81s-10.006\right)}{A},x_2=\mathrm{\Delta }{\delta }_V\left(s\right).$$E19

In Eq. (17), G_V/δM (s) is the mathematical model of aircraft speed (a transfer function) depending on the fuel supply, δ_M(s) is the fuel supply, G_V/δV (s) is a mathematical model of aircraft speed depending on the elevator angle, δ_V(s) is the aircraft elevator angle. The mathematical model of aircraft speed in the longitudinal direction, the fuel supply, and angle of attack can be determined, see [14]:

$$G_{V/\mathit{\delta M}}\left(s\right)=-a_x^{\mathit{\delta M}}\frac{{\mathrm{\Delta }}_{11}}{\mathrm{\Delta }},G_{V/\mathit{\delta V}}\left(s\right)=-a_y^{\mathit{\delta V}}\frac{{\mathrm{\Delta }}_{21}}{\mathrm{\Delta }}-a_{\mathit{mz}}^{\mathit{\delta V}}\frac{{\mathrm{\Delta }}_{31}}{\mathrm{\Delta }},$$E20

where a^δM_x is a speed coefficient with respect to the fuel supply, a^δV_y is a pitch coefficient with respect to the elevator angle, a^δV_mz is the coefficient of the speed angle with respect to the elevator angle, ∆(s)—a determinant of the transfer function, ∆₁₁(s), ∆₂₁(s), ∆₃₁(s)—algebraic additions to the determinant ∆(s), see paragraph (Eq. (9)). The fuel supply transfer function (Eq. (17)) can be calculated using the following equation:

$$\begin{array}{ccc}\mathrm{\Delta }{V}_{V/\mathit{\delta M}}\left(s\right)& =& G_{V/\mathit{\delta M}}\left(s\right)\ast \mathrm{\Delta }{\delta }_M\left(s\right)\\ & =& 5\frac{s^3+1.12s^2+62.782s+25.32}{s^4+1.1338s^3+62.7975s^2+28.6585s+4.09291}\mathrm{\Delta }{\delta }_M\left(s\right).\end{array}$$E21

The transfer function of the elevator (Eq. (17)) is given by:

$$\begin{array}{rcl}\mathrm{\Delta }{V}_{V/\mathit{\delta V}}\left(s\right)& =& G_{V/\mathit{\delta V}}\left(s\right)\ast \mathrm{\Delta }{\delta }_V\left(s\right)\\ & =& \frac{-0.11\cdot \left(9.81s+620.973\right)-0.42\cdot \left(-9.81s-10.0062\right)}{s^4+1.1338s^3+62.7975s^2+28.6585s+4.09291}\mathrm{\Delta }{\delta }_V\left(s\right).\end{array}$$E22

4.2 Structure-defined mathematical model of the angle of aircraft attack in a flight simulator

From description angle of aircraft attack Δα(s), we can derive that the equation of angle of aircraft attack displacement defines a change aircraft throttle control (of fuel supply) and a change in the angle of aircraft elevator [10]:

$$\mathrm{\Delta }\alpha \left(s\right)=-G_{\alpha /\mathit{\delta M}}\left(s\right)\mathrm{\Delta }{\delta }_M\left(s\right)-G_{\alpha /\mathit{\delta V}}\left(s\right)\mathrm{\Delta }{\delta }_V\left(s\right),$$E23

where stability determined by zeroes of a characteristic equation is used as a numerator in the structured mathematical model of aircraft motion; see [2]. Next, we define permutation and transformation with regard to (Eq. (16)) and have coefficients for l₂₁(x):

$$A_{21}=5\frac{0.002s^2-0.252s-0.1}{A},x_1=\mathrm{\Delta }{\delta }_M\left(s\right).$$E24

$$\begin{array}{l}{A}_{22}=\frac{-0.11\cdot \left(-s^3+0.886s^2+0.0124s-2.453\right)-0.42\cdot \left(-s^2-0.414s-0.025\right)}{A},\\ x_2=\mathrm{\Delta }{\delta }_V\left(s\right).\end{array}$$E25

The equations Eq. (18) and Eq. (24) will use a step change of throttle control stick (fuel supply) in Laplace transformation ∆δ_M (s) = 1/s. Eqs. (19) and (25) will use a step change of aircraft elevator control (angle of attack) in Laplace transform ∆δ_V (s) = 1/s. The problem statement of a mathematical model of motion in a flying simulator is conditioned by identification of its stability. The roots of the characteristic equation, the denominator equation (Eq. (18)), (Eq. (19)), (Eq. (24)) and equation (Eq. (25)), see [1].

In Eq. (17), G_α/δV (s) is the calculated mathematical model angle of aircraft attack (a transfer function) depending on the fuel supply, δ_M(s) is the fuel supply, G_α/δV (s) is the calculated mathematical model angle of aircraft attack depending on the angle of the elevator, δ_V(s) is the angle of the elevator. A mathematical model of an angle of aircraft attack in the longitudinal direction, fuel supply, and the angle of attack can be determined; see [14]:

$$G_{\alpha /\mathit{\delta M}}\left(s\right)=-a_x^{\mathit{\delta M}}\frac{{\mathrm{\Delta }}_{12}}{\mathrm{\Delta }},G_{\alpha /\mathit{\delta V}}\left(s\right)=-a_y^{\mathit{\delta V}}\frac{{\mathrm{\Delta }}_{22}}{\mathrm{\Delta }}-a_{\mathit{mz}}^{\mathit{\delta V}}\frac{{\mathrm{\Delta }}_{32}}{\mathrm{\Delta }}.$$E26

The transfer function of fuel supply (Eq. (23)) can be calculated using the following equation:

$$\begin{array}{c}\mathrm{\Delta }{\alpha }_{\alpha /\mathit{\delta M}}\left(s\right)=G_{\alpha /\mathit{\delta M}}\left(s\right)\ast \mathrm{\Delta }{\delta }_M\left(s\right)\\ =5\frac{0.002s^2-0.2518s-0.1}{s^4+1.1338s^3+62.7975s^2+28.6585s+4.09291}\mathrm{\Delta }{\delta }_M\left(s\right).\end{array}$$E27

$$\begin{array}{r}\mathrm{\Delta }{\alpha }_{\alpha /\mathit{\delta V}}\left(s\right)=G_{\alpha /\mathit{\delta V}}\left(s\right)\ast \mathrm{\Delta }{\delta }_V\left(s\right)\\ =\frac{-0.11\cdot \left(-s^3+0.886s^2+0.0124s-2.453\right)-0.42\cdot \left(-s^2-0.414s-0.025\right)}{s^4+1.1338s^3+62.7975s^2+28.6585s+4.09291}\mathrm{\Delta }{\delta }_V\left(s\right).\end{array}$$E28

5. Visualization of results of parallel simulation of mathematical models

Initial or limiting restricting conditions in the given flight phase affect the form of equations of the system depending on which phase of aircraft motion they are calculated [17]. Visualization of results has an impact on quality simulation and simulation tries to obtain the information about the properties of a real system by means of an experiment, the so-called simulation model [18]. “Computer simulation of flying simulator is employed as enlargement or replacement of a structured mathematical model of aircraft motion for which an analytical solution is difficult of even impossible. [19]”.

The sequential program of the mathematical model program is characterized by the equations of simulation of aircraft motion in a single computer time in equidistant moments. The disadvantage of this method is the limitation of the power of the processor, which calculates mathematical models of aircraft motion [20]. For presentation of more accurate simulation results, we need a higher-quality simulation system and a visualization generator providing artificial surrounding of required quality; this surrounding is a three-dimensional scene, see Figure 3.

Figure 3.

Principle of a pilot’s activity and its visualization in a projection system.

The first term according to Eq. (17) or Eq. (23) represents the transfer function (structured mathematical model of aircraft motion) of speed displacement shift depending on the fuel supply of the aircraft and transfer function (structured mathematical model of aircraft motion) of the angle of attack depending on the aircraft fuel supply. In the polynomial expression of the transfer function, we derive the following form for the transfer function of increase in speed calculated from the change in fuel supply in meters per second and calculated displacement angle of attack shift depending on fuel supply in radians:

$$-G_{V/\mathit{\delta M}}\left(s\right)\mathrm{\Delta }{\delta }_M\left(s\right)=-5\frac{s^3+1.12s^2+62.782s+25.32}{A}\frac{M}{s}.$$E29

$$-G_{\alpha /\mathit{\delta M}}\left(s\right)\mathrm{\Delta }{\delta }_M\left(s\right)=-5\frac{0.002s^2-0.252s-0.1}{A}\frac{M}{s}.$$E30

If displacement of the speed is considered in Eqs. (29) and (30) and this is conditioned by the step of change in fuel supply (unit step).

The second member according to Eq. (17) or Eq. (23) represents the transfer function of the speed increment depending on the angle of the aircraft elevator and represents the transfer function of the angle of attack depending on the angle of the aircraft elevator. In polynomial expression of the transfer function, we derive the following form for the transfer function of the increase in speed calculated from the change of aircraft elevator in meters per second and of calculated displacement of angle of attack from the aircraft elevator in radians:

$$-G_{V/\mathit{\delta V}}\left(s\right)\mathrm{\Delta }{\delta }_V\left(s\right)=-\frac{-0.11\cdot \left(9.81s+620.973\right)-0.42\cdot \left(-9.81s-10.006\right)}{A}\frac{V}{s}.$$E31

$$\begin{array}{c}-G_{\alpha /\mathit{\delta V}}\left(s\right)\mathrm{\Delta }{\delta }_V\left(s\right)=\\ -\frac{-0.11\cdot \left(-s^3+0.886s^2+0.0124s-2.453\right)-0.42\cdot \left(-s^2-0.414s-0.025\right)}{A}\frac{V}{s}.\end{array}$$E32

The displacement of the elevation is considered in the Eq. (31) or Eq. (32) respectively, and this is conditioned by the step of elevator angle (unit step). The design and compilation of processes of a structured mathematical model of aircraft motion movement in the fuel supply to aircraft engines and angle of the aircraft elevator must be accurate. The simulation in our solution takes 30 s, and the intermediate data is sent in periodical time to the processor’s core or node processing recorded simulated data and creating a graphical form of calculated results after the end of simulation.

5.1 Parallel simulation of a structure-defined mathematical model of an aircraft using MPI

The sequential program of a mathematical model program is characterized by the calculation of equations in a single computer time. The code operations are performed sequentially in that order. The disadvantage of this method is the energy consumption of the processor that counts the models [20]. The state of the art of parallel system based on the standard MPI can be introduced as the first one.

5.1.1 Principle of parallel message passing interface

The simulation parallel program computers in distributed computer systems are referred to as node computers [21]. They usually consist of a primary message input queue, one or more equivalent processors, and necessary equipment for communicate over the interconnectors. The threads operate in either serial or parallel modes.

The MPI facilitates this approach by providing many ways to call open-source industrial MPI implementations such as MPICH and LAM-MPI. The send/receive commands are implemented to change messages in the source code application and are added to run on nodes. We use two basic functions to send and receive messages [22]:

• MPI_Send (parameters)

• MPI_Recv (parameters).

5.1.2 The program application created using MPI

There are n nodes, which consist of the processor P_i and the local memory M_j. The nodes communicate with each other using lines and an interconnection network. When executing a given program, the program is divided into concurrent processes, each of which is executed in a separate processor. This simultaneous execution of the same task on multiple processors is used to obtain results faster.

Distributed architecture was realized as connection of five nodes: one is a central computer, and the others are computing nodes (Figure 4). One node (N₁) is designed as a central computer, and the others are computational nodes, and each of them computes only one mathematical model. The MPI provides alternative methods for communication and movement of data among multiprocessors. There is no global memory, it is necessary to move data from one local memory to another by means of message passing [23].

Figure 4.

Message passing interface – architecture, Ni – node-computer, Pi – processor, Mi – local memory.

Where the MPI_Send() function on a side of the sender is responsible for sending messages. The corresponding MPI_Recv() function is inserted into a target process to receive messages. The simulation takes 30 seconds, and the intermediate data are sent in periodical time to the node that presents the received data in a graphical form.

5.1.3 The distributed modeling

The implementation MPICH2 is a portable, high-performance implementation of the entire MPI-2 standard and consists of a library of routines that can be called from the program. The TOOLKIT is an integrated set of tools that supports measurement, analysis, assignment, and presentation of application performance for sequential and parallel programs [24].

As follows from the expression, the structured mathematical model of aircraft defined by Eq. (18)A_{11 *}x₁ is simulated by the computer N₂ on the second node, (Eq. (19)) A_{12 *}x₂ is simulated by the computer N₃ on the third node, structured mathematical model of aircraft defined by Eq. (24)A_{21 *}x₁ is simulated by the computer N₄ on the fourth node, and Eq. (25)A_{22 *}x₂ is simulated by the computer N₅ on the fifth node.

From the graphical output of the central node, we get Figure 5 as a result of the simulation. The simulation of a structured mathematical model of aircraft motion using a flight simulator cluster technology is done according to Eqs. (18) and (19) or Eqs. (24) and (25).

Figure 5.

Simulation results of (Eq. (18)) – upper left, (Eq. (19)) – upper right, (Eq. (24)) – bottom left, (Eq. (25)) – bottom right.

5.1.4 Simulation results from MPI

The simulation results are in Figure 5, the upper left figure shows a graphical presentation of the aircraft speed increase depending on the fuel supply, and it is equal to 31.0192 [m/s], it is identical to the figure in Figure 6. The figure on the top right shows the aircraft speed increment depending on from the elevator, and it is equal to −15.6142 [m/s]. The figure at the bottom left shows a graphical representation of the increment angle of the aircraft attack depending on the fuel supply and is equal to −0.1237 [rad], it is identical to the presentation in Figure 7. The figure at the bottom right shows the increment angle of the pitch depending on the rudder and is equal to −0.0714 [rad]. The simulation time is set depending on the value of the integration error, which is less than 0.002 with these distributed simulation methods.

Figure 6.

A block diagram of simulation of a mathematical model of aircraft motion on the GPU computation.

Figure 7.

GeForce 560 GPU block scheme of Streaming Multiprocessors (SM), Load (LD), Store (ST), Special Function Unit (SFU).

Modeling of the parallel aspect of the decomposed flight simulator subsystems in the form of mathematical notation was performed in accordance with Eqs. (18) and (19) or Eqs. (24) and (25). The strengths of this model are the combination of both the advantages: efficiency (memory savings) and ease of programming of a shared-memory method and scalability of a distributed-memory method.

5.2 Parallel simulation of a structure-defined mathematical model of an aircraft using GPU

Computer simulation a structured mathematical model of aircraft motion is used for modeling of aircraft characteristics in such cases [25]. The numerical integration must be available in real time, in single computer time in equidistant moments. The above, two-parameter control of a block structured mathematical model of motion is performed by controlling the throttle lever of the engine and throttle rudder, which cause a change in the movement of the aircraft and the angle of attack of the aircraft. Parallel codes running on GPU hardware can yield results equivalent to the performance of dozens of traditional CPUs.

Input of block mathematical model represents control values of block structured mathematical model of aircraft motion—the fuel supply or the elevator angle. Output of block mathematical model represents simulated values of block structured mathematical model of aircraft motion—the speed or the attack angle. A given system can simulate multiple computers or processors. When several processors of a simulator system are connected (Block₁, Block₂, …, Block_n), they communicate with each other via SMS; see Figure 8.

Figure 8.

Outline of a revised host code ComputeFlightGPU ()

Control values of block structured mathematical models of aircraft motion—fuel supply and an aircraft elevator angle represent an input in a block diagram and an output in a block diagram, next figure.

5.2.1 Approach to CPU block mathematical model of an aircraft

The GPU blocks (Block₁, Block₂, …, Block₄) count and each of them numerically counts only one block of the mathematical model of aircraft motion. According to our solution, the algorithm of the block structured mathematical model of aircraft motion is the most suitable for implementation in a GPU-enabled manner [2, 26].

In proposed solution, CPU system communicates with GPU by means of an SMS. The CPU is designed with fewer processor cores that have higher clock speeds, allowing them to complete series of tasks very quickly. The CPU is the brain of operation that is responsible for providing instructions to the rest of system. On the other hand, a GPU has much larger number of cores and is designed for a different purpose; it has parallel processing architecture that allows it to perform multiple calculations across streams of data simultaneously. The GPU has a massive parallel architecture consisting of thousands of smaller, more efficient cores [26].

5.2.2 Achievements of implementation of structure-defined mathematical models of an aircraft using CUDA

In our solution, the user can monitor the simulation process in real time using the GPU (Figure 9). Parallel simulation is one of the areas of simulation optimization and global optimization, which are implemented by the method of estimating the differences between continuous and discrete versions. The division of the structured mathematical model of aircraft motion into blocks is the motivation of this research with the definition of the relationships between the created blocks of the mathematical model.

Figure 9.

The CUDA kernel of BlockSimulationModelGPU() function.

The design of a block structure of mathematical model of aircraft motion is realized by various methods of simulation software systems, and one of them is CUDA computing. The application programming interface (API) improves computing performance with a graphics card for numerical calculation in simulations, and it is standardized layer that allows applications to take advantage of software or hardware services and features. The CUDA includes C/C++ software development tools, function libraries, and a hardware abstraction mechanism that hides the GPU hardware from developers [22].

5.2.3 Basic software solution of computation on the GPU

The ComputeFlightGPU (…) function is essentially an outsourcing agent that sends input data to the device, activates the execution of the simulation on the device, and collects the results from the device [28]. To perform a simulation of a block structure mathematical model of aircraft motion on a GPU core, the programmer musts allocate the required memory by calling the cudaMalloc (…) function on the device. To transfer the relevant data from the host memory to the allocated device memory, we call cudaMemcpy (…) function, see Figure 10—Part 1.

Figure 10.

Screenshot of simulation results based on (Eq. (18)), horizontal axes – time, vertical axes – speed increment.

After completing the simulation, the programmer must transfer results from the device memory back to the host memory by calling function cudaMemcpy (…) and free the device memory by calling function cudaFree (…), see Figure 10 – Part 3. The CUDA kernel provides API simulation functions running as programmer-defined functions [27]. This core function determines the function of the block structure mathematical model of aircraft motion to be performed by all threads during the parallel simulation phase, see Figure 10—Part 2.

Parallel programming in C/C++ supports a CPU computer and the GPU card supports programming in C [28]. The inner loop in the CUDA implementation disappears, because the ComputeFlightGPU calculation is performed in parallel by CUDA threads, see previous figure. Our presented solution was simulated on a personal computer that consists of a CPU Intel Quad Core Q9650 processor with 4 cores, 3.00 GHz each, cache 12 MB L2, bus speed 1333 MHz (FSB), 4GB RAM DDR3, and GPU accelerator NVIDIA GeForce GTX 560 Ti.

The first thing to notice is the __global__ keyword, see Figure 11. This simply means that this function can be called either from host computer or the CUDA device. Each thread runs the same block mathematical model of aircraft motion, so the only way to differentiate yourself from other threads is to use their threadIdx and your blockIdx [29].

Figure 11.

Screenshot of simulation results based on (Eq. (24)), horizontal axes – time, vertical axes – angle of attack increment.

The index in the thread array is calculated and determined by the block and thread ID. Each thread uses its own threadIdx.x and threadIdx.y to identify the elements of a block mathematical model of aircraft motion defined by the system of equations (Eq. (17)) and (Eq. (23)). The expression blockIdx.y ˟ 2 ˟ 2 equals the number of threads in all grid rows above the current thread position. The expression blockIdx.x ˟ 2 equals the number of all columns in the current grid row. Finally, threadIdx.x equals the number of threads above the current block.

The function BlockSimulationModelGPU() calculates a unique ID in the Idx register variable, which is then used as an identifier in the array and in the aircraft motion simulation calculation model.

5.2.4 Simulation results from CUDA

As follows from the expression, the block of the structured mathematical model of the aircraft is defined by Eq. (18)A_{11 *}x₁ and is simulated by the GPU block Block₁, the block structure of the mathematical model of the aircraft is defined by Eq. (24)A_{21 *}x₁ and is simulated by the GPU block Block₂. From the obtained graphical output, Figures 6 and 7, which are the results of the numerical integration of the block structure of the mathematical models of aircraft motion.

Graphically, the course of Figure 6 shows the increase in aircraft speed as a function of fuel delivery not exceeding 31.0174 [m/s]. Figure 7 shows the angle of increment of the aircraft attack as a function of the change in the position of the rudder, at which it does not exceed −0.1247 [rad].

The CUDA memory model consists of different memory spaces, which differ significantly in latency times [30]. Results from comparing of the simulation of block structured mathematical model of aircraft motion on the CPU and on the GPU are in Table 1.

Number of block mathematical models	CPU; Runtime [ms]	GPU; Runtime [ms]	Speed-Up
1	1,60	0,07	22,72
2	3,10	0,14	22,01
5	7,08	0,33	21,74
10	13,82	0,65	21,22
20	25,37	1,24	20,44
50	47,39	2,37	20,02

Table 1.

The performance obtained from simulation different members of the mathematical models of aircraft motion.

In the number of block mathematical models column, the number means the number of block structure of mathematical model of aircraft motion simulated in parallel. The times achieved from the measurement of the simulation of mathematical aircraft models on the CPU are in CPU runtime column and on the GPU in the GPU runtime column. The last column is the acceleration of the calculation on the GPU compared on the CPU.

The table shows the performance acceleration obtained by simulating a different number of block structured mathematical models of motion of aircraft. The columns show the number of block mathematical models, the CPU run time, GPU run time, and acceleration. Decreasing the simulation speed of a block structure of mathematical models of aircraft motion on the GPU with an increasing number of simulated mathematical models is caused by an increasing the time required to read or retrieve data from or to the host computer.

This subchapter presents GPU calculation and CUDA programming with an explanation of how to implement an efficient GPU application. The subchapter aims to offer data obtained by simulating the block structure of a mathematical model of aircraft motion using the GPU. Even a non-optimized parallel implementation of the block structured mathematical model of aircraft motion on the GPU can lead to a significant reduction in computational time compared with the implementation to the CPU.

6. Conclusions

The results from simulation help to overcome the physical and architectural limitations of computational power that can be achieved with a single-processor system. It is then very convenient for the application to use a more powerful processor, a faster cache, and faster access to the RAM with higher transfer capacity. The results obtained in the simulation of structured mathematical models of aircraft motion and the use of this method appear to be effective and pragmatic according to the requirements of the integration algorithm given in the chapter.

In general, in accordance with the proposed solution, it is also possible to run the program faster due to factors such as the availability of better bandwidth for communication within node. Greater bandwidth provides additional communication with the interface to transfer message between nodes and the higher performance of the graphics processing unit.

The described design of a structured mathematical model of aircraft motion to perform simulation on the GPU, this simulation was chosen as an example of an algorithm confirming the proposed solution. The presented problem of simulation methods of block structure of mathematical models of aircraft motion was solved by predetermined known procedures. The advantages of the methods are a combination of both the features: the efficiency of computer system programming and the scalability of a GPU computational method.

Based on knowledge gained in design of the block structure of mathematical models of aircraft motion, we can say that the graphics card can be used in multiple and parallel simulation, which confirms high speed of computations. The advantage of faster simulation can sometimes be obtained at the earliest defined quality. These facts and accuracies are well suited for use in real-time simulation and a virtual reality imaging system. The results of the simulation of such designed structured mathematical model of aircraft motion show very fast calculation and confirm the effectiveness of such a structure of the mathematical model of aircraft for the proposed system of the simulator, which was the motivation of our research.

Acknowledgments

The part of solving the chapter was supported by the Slovak Grant Agency for Science, by the VEGA 4/0330/09 grant.