Design optimization of multi-objective proportional–integral–derivative controllers for enhanced handling quality of a twin-engine,propeller-driven airplane

Abstract

Herein, the design optimization of multi-objective controllers for the lateral–directional motion using proportional–integral–derivative controllers for a twin-engine, propeller-driven airplane is presented. The design optimization has been accomplished using the genetic algorithm and the main goal was to enhance the handling quality of the aircraft. The proportional–integral–derivative controllers have been designed such that not only the stability of the lateral–directional motion was satisfied but also the optimum result in longitudinal trim condition was achieved through genetic algorithm. Using genetic algorithm optimization, the handling quality was improved and placed in level 1 from level 2 for the proposed aircraft. A comprehensive sensitivity analysis to different velocities, altitudes and centre of mass positions is presented. Also, the performance of the genetic algorithm has been compared to the case where the particle swarm optimization tool is implemented. In this work, the aerodynamic coefficients as well as the stability and control derivatives were predicted using analytical and semi-empirical methods validated for this type of aircraft.

Introduction

Propeller effects and power contribution to the trim condition and stability of a propeller-driven airplane are significant. Since the beginning of the airplane introduction in 1909, this has always been of interest where a propeller provided significant force in the presence of side wind.^1–3 Initial analytical behaviour of this phenomenon was investigated by Harris4 in 1918 where he showed the resulting pitching moment due to the side wind effect. Later, Glauert^5,6 expanded this investigation to the other stability derivatives. In the United States, studies on propeller behaviour began in the 1920s and enhanced over World War II where new and more powerful engines were implemented. Since then, most of the methods were based on the linear classic momentum theory.⁷ Later in 1972, the National Aeronautics and Space Administration (NASA) of the United States did great research to investigate the aerodynamic coefficients and stability and control derivatives of twin-engine propeller-driven airplanes. The resulting investigations were published as longitudinal and lateral–directional aerodynamic characteristics of twin-engine propeller-driven airplanes.^8,9 Recently, several studies have been done to enhance semi-empirical methods to investigate aerodynamic characteristics of propeller-driven airplanes and their effects on the handling quality.^10–15

Generally, in single-engine light propeller-driven airplanes, in order to decrease the asymmetrical propeller thrust, one can use a small side angle (i.e. typically to the right-hand side for a clockwise rotating propeller from a pilot field of view). However, in the case of twin-engine propeller-driven airplanes with propellers turning in the same direction (i.e. clockwise from a pilot field of view), using a side angle results in an asymmetry in the design and should be avoided. Consequently, in order to compensate the propeller effects in these types of aircraft, usually a very small permanent deflection (∼1°–2°) is considered for the fixed part of the vertical stabilizer.^16,17

In the previous work, an enhanced semi-empirical method is developed to calculate longitudinal and lateral–directional stability and control derivatives of light airplanes.10 In this study, using the previous tool to estimate the aerodynamic characteristics of the proposed airplane which is a twin-engine propeller-driven aircraft, in the first place, the aircraft has been trimmed for the longitudinal motion and then the lateral–directional stability has been guaranteed. Due to the significant effect of propellers and corresponding power contribution on aerodynamic characteristics, the overall handling quality of the aircraft decreased and placed in level 2 for the Dutch-Roll mode. As mentioned, in twin-engine aircraft, in order to compensate the propeller effects, usually a very small permanent deflection is considered for the vertical tail. However, the proposed permanent deflection cannot guarantee the handling quality level 1 in all cruise flight conditions as it is designed for the conventional cruise flight condition. Consequently, pilots require to use controller surfaces to compensate the propeller effect in other flight conditions and this will increase their workload.^16,17 On the other side, passenger will also not fill very comfortable throughout the flight if handling quality level 1 is not guaranteed. Accordingly, the main goal of the proportional–integral–derivative (PID) controllers designed and optimized in this study has been to improve the overall handling quality of the proposed aircraft. In the following, the modelling process and the mathematical method are presented in detail and the corresponding results are investigated and compared.

Mathematical model

The airplane equations of motion are modelled in accordance with small perturbations equations as follows¹⁸

$$\text{m}\left(\dot{\mathrm{v}}+{\text{U}}_1\text{r}\right)=\text{mg}\mathrm{\varnothing }{\text{ cosθ}}_1+{\overline{q}}_1\text{S}\left\{{\text{C}}_{{\text{y}}_{\text{β}}}\text{β}+{\text{C}}_{{\text{y}}_{\text{p}}}\frac{\text{pb}}{2{\text{U}}_1}+{\text{ C}}_{{\text{y}}_{\text{r}}}\frac{\text{rb}}{2{\text{U}}_1}+{\text{C}}_{{\text{y}}_{{\text{δ}}_{\text{a}}}}{\text{δ}}_{\text{a}}+{\text{C}}_{{\text{y}}_{{\text{δ}}_{\text{r}}}}{\text{δ}}_{\text{r}}\right\}$$

(1)

$${\text{I}}_{\text{xx}}\dot{\mathrm{p}}-{\text{I}}_{\text{xz}}\dot{\mathrm{r}}={\overline{q}}_1\text{Sb}\left\{{\text{C}}_{{\text{l}}_{\text{β}}}\text{β}+{\text{C}}_{{\text{l}}_{\text{p}}}\frac{\text{pb}}{2{\text{U}}_1}+{\text{C}}_{{\text{l}}_{\text{r}}}\frac{\text{rb}}{2{\text{U}}_1}+{\text{C}}_{{\text{l}}_{{\text{δ}}_{\text{a}}}}{\text{δ}}_{\text{a}}+{\text{ C}}_{{\text{l}}_{{\text{δ}}_{\text{r}}}}{\text{δ}}_{\text{r}}\right\}$$

(2)

$${\text{ I}}_{\text{zz}}\dot{\mathrm{r}}-{\text{I}}_{\text{xz}}\dot{\mathrm{p}}={\overline{q}}_1\text{Sb}\left\{{\text{C}}_{{\text{n}}_{\text{β}}}\text{β}+{\text{C}}_{{\text{n}}_{\text{p}}}\frac{\text{pb}}{2{\text{U}}_1}+{\text{C}}_{{\text{n}}_{\text{r}}}\frac{\text{rb}}{2{\text{U}}_1}+{\text{C}}_{{\text{n}}_{{\text{δ}}_{\text{a}}}}{\text{δ}}_{\text{a}}+{\text{C}}_{{\text{n}}_{{\text{δ}}_{\text{r}}}}{\text{δ}}_{\text{r}}\right\}$$

(3)where $\text{p}=\dot{\varnothing },\text{ r}=\dot{\psi }$ , $\text{v}={\text{U}}_1\text{β}$ and ${\overline{q}}_1=0.5{\text{ρv}}^2$ . ${\text{I}}_{\text{xx}}$ and ${\text{ I}}_{\text{yy}}$ are airplane moments of inertia about X- and Y-axis, respectively, and ${\text{I}}_{\text{xz}}$ is the airplane product of inertia about Z-axis.

Considering the fact that here the main goal is to implement the PID controller to control the lateral–directional motion, the corresponding equations for ${\text{C}}_{{\text{y}}_{\text{β}}}$ , ${\text{C}}_{{\text{l}}_{\text{β}}}$ and ${\text{C}}_{{\text{n}}_{\text{β}}}$ with respect to the propeller effects and power on contribution are presented as follows¹⁸

$${\text{C}}_{{\text{y}}_{\text{β}}}={\text{C}}_{{\text{y}}_{{\text{β}}_{\text{propOff}}}}+{\text{C}}_{{\text{y}}_{{\text{β}}_{\text{propOn}}}}$$

(4)

$${\text{C}}_{{\text{l}}_{\text{β}}}={\text{C}}_{{\text{l}}_{{\text{β}}_{\text{propOff}}}}+{\text{C}}_{{\text{l}}_{{\text{β}}_{\text{propOn}}}}$$

(5)

$${\text{C}}_{{\text{n}}_{\text{β}}}={\text{C}}_{{\text{n}}_{{\text{β}}_{\text{propOff}}}}+{\text{C}}_{{\text{n}}_{{\text{β}}_{\text{propOn}}}}$$

(6)where ${\text{C}}_{{\text{y}}_{\text{β}}}$ is the variation of airplane side force coefficient with sideslip angle, ${\text{C}}_{{\text{l}}_{\text{β}}}$ is the variation of airplane rolling moment coefficient with angle of sideslip and ${\text{C}}_{{\text{n}}_{\text{β}}}$ is the variation of airplane yawing moment coefficient with sideslip angle.

On the other hand, the airplane moments of inertia are defined in the body-fixed axis system. Since the equations of motion are defined in the stability axis system, these moments are also needed to be defined in the same axis system. Accordingly, the transformation can be defined as follows¹⁸

$$\left\{\begin{array}{c}{\text{I}}_{{\text{xx}}_{\text{s}}}\\ {\text{I}}_{{\text{zz}}_{\text{s}}}\\ {\text{I}}_{{\text{xz}}_{\text{s}}}\end{array}\right\}=\left[\begin{array}{ccc}{\cos }^2{\text{α}}_1& {\sin }^2{\text{α}}_1& -{\sin }^2{\text{α}}_1\\ {\sin }^2{\text{α}}_1& {\cos }^2{\text{α}}_1& {\sin }^2{\text{α}}_1\\ \frac{1}{2}{\sin }^2{\text{α}}_1& -\frac{1}{2}{\sin }^2{\text{α}}_1& {\cos }^2{\text{α}}_1\end{array}\right]\left\{\begin{array}{c}{\text{I}}_{{\text{xx}}_{\text{B}}}\\ {\text{I}}_{{\text{zz}}_{\text{s}}}\\ {\text{I}}_{{\text{xz}}_{\text{B}}}\end{array}\right\}$$

(7)

With respect to the perturbed lateral–directional equations of motion provided in equations (1)–(3), the corresponding Laplace transforming equations for zero initial conditions can be defined as follows¹⁸

$$\begin{array}{c}\left({\text{sU}}_1-{\text{Y}}_{\text{β}}\right)\text{β}\left(\text{s}\right)-\left({\text{sY}}_{\text{p}}+\text{g}\text{cos}{\text{θ}}_1\right)\mathrm{\varphi }\left(\text{s}\right)\text{\hspace{0.17em}s}\left({\text{U}}_1-{\text{Y}}_{\text{r}}\right)\text{ψ}\left(\text{s}\right)={\text{Y}}_{\text{δ}}\text{δ}\left(\text{s}\right)\\ -{\text{L}}_{\text{β}}\text{β}\left(\text{s}\right)+\left({\text{s}}^2-{\text{L}}_{\text{p}}\text{s}\right)\mathrm{\varphi }\left(\text{s}\right)-\left({\text{s}}^2{\overline{A}}_1+{\text{sL}}_{\text{r}}\right)\text{ψ}\left(\text{s}\right)=\text{δ}\left(\text{s}\right)\\ -\left({\text{N}}_{\text{β}}+{\text{N}}_{{\text{T}}_{\text{β}}}\right)\text{β}\left(\text{s}\right)-\left({\text{s}}^2{\overline{B}}_1+{\text{N}}_{\text{p}}\text{s}\right)\mathrm{\varphi }\left(\text{s}\right)({\text{S}}^2-{\text{sN}}_{\text{r}})\text{ψ}\left(\text{s}\right)={\text{N}}_{\text{δ}}\text{δ}\left(\text{s}\right)\end{array}$$

(8)where $\text{β}(\text{s})$ , $\mathrm{\varphi }(\text{s})$ and $\text{ψ}(\text{s})$ are Laplace transformed of sideslip angle, bank angle and heading angle, respectively, and $\text{δ}(\text{s})$ is the Laplace transformed aileron or rudder input. In the equations above, ${\overline{A}}_1={\text{I}}_{\text{xz}}/{\text{I}}_{\text{xx}}$ and ${\overline{B}}_1={\text{I}}_{\text{xz}}/{\text{I}}_{\text{zz}}$ .

And writing equation (8) in the matrix and transfer function format will result in the following equations¹⁸

$$\left[\begin{array}{l}\begin{array}{cc}\left({\text{sU}}_1-{\text{Y}}_{\text{β}}\right)& \begin{array}{cc}-\left({\text{sY}}_{\text{p}}+\text{g}\text{cos}{\text{θ}}_1\right)& \text{s}\left({\text{U}}_1-{\text{Y}}_{\text{r}}\right)\end{array}\end{array}\\ \begin{array}{ccc}-{\text{L}}_{\text{β}}& \left({\text{s}}^2-{\text{L}}_{\text{p}}\text{s}\right)& -\left({\text{s}}^2{\overline{A}}_1+{\text{sL}}_{\text{r}}\right)\end{array}\\ \begin{array}{ccc}-{\text{N}}_{\text{β}}{\text{+N}}_{{\text{T}}_{\text{β}}}& -\left({\text{s}}^2{\overline{B}}_1+{\text{N}}_{\text{p}}\text{s}\right)& \left({\text{s}}^2-{\text{sN}}_{\text{r}}\right)\end{array}\end{array}\right]\left\{\begin{array}{c}\frac{\text{β}\left(\text{s}\right)}{\text{δ}\left(\text{s}\right)}\\ \frac{\mathrm{\varphi }\left(\text{s}\right)}{\text{δ}\left(\text{s}\right)}\\ \frac{\text{ψ}\left(\text{s}\right)}{\text{δ}\left(\text{s}\right)}\end{array}\right\}=\left\{\begin{array}{c}{\text{Y}}_{\text{δ}}\\ {\text{L}}_{\text{δ}}\\ {\text{N}}_{\text{δ}}\end{array}\right\}$$

(9)

Finally, the corresponding transfer functions of the aircraft can be defined as follows¹⁸

$$\frac{\text{β}\left(\text{s}\right)}{\text{δ}\left(\text{s}\right)}=\frac{\left[\begin{array}{ccc}{\text{Y}}_{\text{δ}}& -({\text{sY}}_{\text{p}}+\text{g}\text{cos}{\text{θ}}_1& \text{s}\left({\text{U}}_1-{\text{Y}}_{\text{r}}\right)\\ {\text{L}}_{\text{δ}}& \left({\text{s}}^2-{\text{L}}_{{\text{p}}_{\text{s}}}\right)& -\left({\text{s}}^2{\overline{A}}_1+{\text{sL}}_{\text{r}}\right)\\ {\text{N}}_{\text{δ}}& -\left({\text{s}}^2{\overline{B}}_1+{\text{N}}_{\text{p}}\text{s}\right)& \left({\text{s}}^2-{\text{sN}}_{\text{r}}\right)\end{array}\right]}{\left[\begin{array}{l}\begin{array}{cc}\left({\text{sU}}_1-{\text{Y}}_{\text{β}}\right)& \begin{array}{cc}-\left({\text{sY}}_{\text{p}}+\text{g}\text{cos}{\text{θ}}_1\right)& \text{s}\left({\text{U}}_1-{\text{Y}}_{\text{r}}\right)\end{array}\end{array}\\ \begin{array}{ccc}(-{\text{L}}_{\text{β}})& \left({\text{s}}^2-{\text{L}}_{\text{p}}\text{s}\right)& -\left({\text{s}}^2{\overline{A}}_1+{\text{sL}}_{\text{r}}\right)\end{array}\\ \begin{array}{ccc}(-{\text{N}}_{\text{β}}+\hspace{0.17em}\text{N}{\text{T}}_{\text{β}})& -\left({\text{s}}^2{\overline{B}}_1+{\text{N}}_{\text{p}}\text{s}\right)& \left({\text{s}}^2-{\text{sN}}_{\text{r}}\right)\end{array}\end{array}\right]}$$

(10)

$$\frac{\mathrm{\varphi }\left(\text{s}\right)}{\text{δ}\left(\text{s}\right)}=\frac{\left[\begin{array}{ccc}\left({\text{sU}}_1-{\text{Y}}_{\text{β}}\right)& {\text{Y}}_{\text{δ}}& \text{s}\left({\text{U}}_1-{\text{Y}}_{\text{r}}\right)\\ -{\text{L}}_{\text{β}}& {\text{L}}_{\text{δ}}& -\left({\text{s}}^2{\overline{A}}_1+{\text{sL}}_{\text{r}}\right)\\ -\left({\text{N}}_{\text{β}}+{\text{N}}_{{\text{T}}_{\text{β}}}\right)& {\text{N}}_{\text{δ}}& \left({\text{s}}^2-{\text{sN}}_{\text{r}}\right)\end{array}\right]}{\left[\begin{array}{l}\begin{array}{cc}\left({\text{sU}}_1-{\text{Y}}_{\text{β}}\right)& \begin{array}{cc}-\left({\text{sY}}_{\text{p}}+\text{g}\text{cos}{\text{θ}}_1\right)& \text{s}\left({\text{U}}_1-{\text{Y}}_{\text{r}}\right)\end{array}\end{array}\\ \begin{array}{ccc}(-{\text{L}}_{\text{β}})& \left({\text{s}}^2-{\text{L}}_{\text{p}}\text{s}\right)& -\left({\text{s}}^2{\overline{A}}_1+{\text{sL}}_{\text{r}}\right)\end{array}\\ \begin{array}{ccc}(-{\text{N}}_{\text{β}}{\text{+N}}_{{\text{T}}_{\text{β}}})& -\left({\text{s}}^2{\overline{B}}_1+{\text{N}}_{\text{p}}\text{s}\right)& \left({\text{s}}^2-{\text{sN}}_{\text{r}}\right)\end{array}\end{array}\right]}$$

(11)

$$\frac{\text{ψ}\left(\text{s}\right)}{\text{δ}\left(\text{s}\right)}=\frac{\left[\begin{array}{ccc}\left({\text{sU}}_1-{\text{Y}}_{\text{β}}\right)& -\left({\text{sY}}_{\text{p}}+\text{g}\text{cos}{\text{θ}}_1\right)& {\text{Y}}_{\text{δ}}\\ -{\text{L}}_{\text{β}}& \left({\text{s}}^2-{\text{L}}_{\text{p}}\text{s}\right)& {\text{L}}_{\text{δ}}\\ -\left({\text{N}}_{\text{β}}+{\text{N}}_{{\text{T}}_{\text{β}}}\right)& -\left({\text{s}}^2{\overline{B}}_1+{\text{N}}_{\text{p}}\text{s}\right)& {\text{N}}_{\text{δ}}\end{array}\right]}{\left[\begin{array}{l}\begin{array}{cc}\left({\text{sU}}_1-{\text{Y}}_{\text{β}}\right)& \begin{array}{cc}-\left({\text{sY}}_{\text{p}}+\text{g}\text{cos}{\text{θ}}_1\right)& \text{s}\left({\text{U}}_1-{\text{Y}}_{\text{r}}\right)\end{array}\end{array}\\ \begin{array}{ccc}(-{\text{L}}_{\text{β}})& \left({\text{s}}^2-{\text{L}}_{\text{p}}\text{s}\right)& -\left({\text{s}}^2{\overline{A}}_1+{\text{sL}}_{\text{r}}\right)\end{array}\\ \begin{array}{ccc}(-{\text{N}}_{\text{β}}{\text{+N}}_{{\text{T}}_{\text{β}}})& -\left({\text{s}}^2{\overline{B}}_1+{\text{N}}_{\text{p}}\text{s}\right)& \left({\text{s}}^2-{\text{sN}}_{\text{r}}\right)\end{array}\end{array}\right]}$$

(12)where $\text{β}(\text{s})/\text{δ}(\text{s})$ , $\mathrm{\varphi }(\text{s})/\text{δ}(\text{s})$ and $\text{ψ}(\text{s})/\text{δ}(\text{s})$ are the sideslip-to-aileron (or -rudder) transfer function, the bank-angle-to-aileron (or -rudder) and the heading-to-aileron (or -rudder), respectively. In the equations above, ${\overline{A}}_1={\text{I}}_{\text{xz}}/{\text{I}}_{\text{xx}}$ and ${\overline{B}}_1={\text{I}}_{\text{xz}}/{\text{I}}_{\text{zz}}$ .The open-loop transfer functions in response to an input can be shown as a block diagram in Figure 1 for the proposed airplane. Therefore, the resulting output can be expressed by

$$\mathrm{\varphi }\left(\text{s}\right)=\text{δ}\left(\text{s}\right)\left(\frac{\mathrm{\varphi }\left(\text{s}\right)}{\text{δ}\left(\text{s}\right)}\right)$$

(13)

Figure 1.

Schematic diagram of the open-loop transfer functions of the aircraft in the presence of a single input.

PID controller

The PID controller is a kind of controller usually used for practical purposes. Recently, over 85% of the controllers in dynamic systems are of this type. Each one of the P, I and D terms are representative of a controlling algorithm term used for a specific purpose. In some cases, one can use only one or two terms instead of the rest of the terms like PI, PD or even P controllers.^19–29

PID controllers can be schematically presented in Figure 2, where the control signal $\text{u}(\text{t})$ illustrates the sum of the three different components and each component is a function of the tracking error $\text{e}(\text{t})$ . ${\text{K}}_{\text{p}}$ is the proportional term, ${\text{K}}_{\text{i}}/\text{s}$ is the integral term and ${\text{K}}_{\text{d}}\text{s}$ is the derivative term. It should also be noted that each component of this system works independently of the others.¹⁰

Figure 2.

Schematic use of a PID controller in a closed-loop system.

In order to introduce each term, one can assume that the rest of the terms are zero. Therefore, assuming ${\text{K}}_{\text{i}}={\text{K}}_{\text{d}}=0$ ; $\text{u}(\text{t})$ is equal to ${\text{K}}_{\text{p}}\text{e }(\text{t})$ . Consequently, at each moment, the control is proportional to the corresponding error and is a function of the present value of the error.

Assuming ${\text{K}}_{\text{p}}={\text{K}}_{\text{d}}=0$ , u(t) is equal to ${\text{K}}_{\text{i}}{\displaystyle {\int }_0^{\text{t}}\text{e}}(\text{τ})\text{dτ}$ . The integral term can be presented by the proportional gain value as follows

$${\text{K}}_{\text{i}}=\frac{{\text{K}}_{\text{p}}}{{\text{τ}}_{\text{i}}}$$

(14)where ${\text{τ}}_{\text{i}}$ is the integral time.

Assuming ${\text{K}}_{\text{p}}={\text{K}}_{\text{i}}=0$ , u(t) is equal to ${\text{K}}_{\text{d}}({\text{d}}_{\text{e}}(\text{t})/\text{dt})$ . The differentiator term can be defined as a function of ${\text{τ}}_{\text{d}}$ , bandwidth limiter, as follows

$${\text{K}}_{\text{d}}={\text{τ}}_{\text{d}}{\text{K}}_{\text{p}}$$

(15)

The ‘non-interacting form’ PID controllers are defined based on the past, present and future values of the control value error in accordance with ${\text{K}}_{\text{i}}$ and ${\text{K}}_{\text{d}}$ as a function of ${\text{K}}_{\text{p}}$ as follows¹⁹

$$\text{u}\left(\text{t}\right)={\text{K}}_{\text{p}}\left(\text{e}\left(\text{t}\right)+\frac{1}{{\text{τ}}_{\text{i}}}{\displaystyle \underset{0}{\overset{\text{t}}{\int }}\text{e}}\left(\text{τ}\right)\text{dτ}+{\text{τ}}_{\text{d}}\frac{\text{de}\left(\text{t}\right)}{\text{dt}}\right)$$

(16)

Genetic algorithm

Overall, analytical optimizations in spite of being so powerful in solving analytical problems are facing several limitations to be used for real engineering problems. These methods are usually used for simple problems or in higher levels they have been used to provide the general solution of a real case. Conversely, numerical optimization methods have significant capabilities to solve complex problems. However, numerical optimizations are also suffering from providing analytical concepts. As numerical optimization methods have several advantages, different numerical optimizations are described in this section. Generally, in most scientific and technical problems, optimum product design depends on finding the best solution in available cases under specific conditions. Accordingly, a practical and suitable solution is required to find the optimum result.

Numerous evolutionary algorithms have been developed to solve global optimization problems. These algorithms are designed such a way that with precise parameters they can provide the general optimum working solution. However, unsuitable parameter selection can result in higher run time to find the optimum solution and consequently a reduction in the performance in terms of accuracy and time complexity.³⁰ Recently, several researches have been focused on extremal optimization (EO) which avoids elite collection while trying to change the bad members by mutation. Multi-objective population-based EO is a variation of EO that have been widely used by researchers in this regard.^31–33

More than a dozen stochastic algorithms have been developed over the past decades which in this study, the focus is on the genetic algorithm (GA) and more details will be provided in this regard. It is interesting to note that among current evolutionary algorithms, the best-known method is the GA. As this type of algorithm assesses numerous working solutions at the same time in the search space, it is more likely to find the global working solution of a certain problem. In addition, it has a simple scalar performance measure without derivative information enquiry which makes the GA as easy as to be used and implemented in most scenarios.³⁴

The GA is a computational method to obtain suitable working solutions while tracking future generations based on the structure of the biological organisms. This evolutionary search heuristic implements natural selection and nature-inspired operators to recognize optimum solutions.³⁵ The evolutionary algorithms idea was first suggested by Rechenberg, in the 1960s.^36,37 After that, GAs were developed by Holland38 in the 1970s and recently, they are very popular for solving optimization problems.

In the GA, the design environment is changing to the genetic environment. Therefore, GAs are working with the coded parameters. The main advantage of the coded parameters is that they are able to change the continuous variables to the discrete ones.^39–41 GAs are stochastic search algorithms that act on a population of possible solutions and leading the solutions towards the most optimum one.⁴² It should be noted that the GAs are not usually able to solve very complex problems where one is facing complex fitness functions or high scale of iterations. Therefore, in such cases, the GAs are facing a significant reduction in performance due to the complexity of time and lower accuracy.^43,44

In order to implement the GA, three concepts are important to be defined and implemented:

Cost functions;

Genetic environment;

GA functions.

The GA has been developed based on Darwin’s evolutionary theory which in the beginning starts with stochastic chromosomes called population. The possible working solutions are programmed as chromosomes and new working solutions can be achieved either by mutating members of the current population or by mating two answers together to produce a better working solution. The better answers are nominated to exist and mutate and the inferior ones are rejected. This procedure will be repeated until finding the best working solution.⁴²

Flying quality

As discussed earlier in the introduction, herein the aerodynamic characteristics of the airplane have been estimated through NAMAYEH.¹⁰ The investigated airplane is a six-place, low wing, twin-engine, propeller-driven general aviation airplane. General characteristics of the airplane are presented in Table 1. For this purpose, the cruise condition considering power on condition and propeller effects is investigated in this study. Properties of the investigated flight condition are presented in Table 2. Considering the cruise flight condition, the aircraft has been trimmed for the longitudinal motion. The longitudinal trim characteristics are provided for the proposed airplane in Table 3.

Table 1.

General characteristics of the investigated airplane.

Parameter	Value	Unit
b	11.8	$\text{m}$
c	1.61	$\text{m}$
S	17.96	${\text{m}}^2$
Aspect ratio	7.75	—
WTO	2500	kg

Table 2.

Properties of the investigated flight conditions.

Parameter	Cruise condition	Unit
M	0.28	—
m	2500	$\text{kg}$
${\overline{X}}_{\text{CG}}$	20	%
alt	3048	m
${\text{δ}}_{\text{flap}}$	0	°

Table 3.

Properties of the investigated trim condition.

Parameter	Value	Unit
${\text{α}}_{\text{trim}}$	0.9	°
${\text{δ}}_{{\text{e}}_{\text{trim}}}$	–1.0	°
Ptrim	398.1	Hp
${\text{C}}_{{\text{L}}_{\text{trim}}}$	0.36	—
${\text{C}}_{{\text{D}}_{\text{trim}}}$	0.038	kg

With respect to the trim condition provided in Table 3, the flying quality characteristics have been investigated and presented in Table. 4. As can be seen, the handling quality placed in level 2 for the Dutch-Roll mode. Accordingly, in the next step, the design and optimization of the PID controllers have been accomplished to increase the Dutch-Roll mode level while keeping all other flying quality characteristics in level 1.

Table 4.

Properties of the investigated longitudinal trim condition.

Flying quality characteristic	Level	$\text{ξ}$	${\text{ω}}_{\text{n}}$
Short period	1	0.76	5.26
Phugoid	1	0.15	0.13
Roll mode	1	1.00	2.79
Spiral mode	1	–1.00	0.00
Dutch-Roll mode	2	0.08	2.79

Controller design

Herein, controllers were designed for the cases where the aircraft is in longitudinal trim and calculations have been accomplished for the lateral–directional motions. Accordingly, two sets of control surfaces including the rudder and ailerons have been designed such that yaw and roll motions of the aircraft to be controlled, respectively. Simultaneously, disturbance from each surface has also been considered to increase the accuracy of the controller design. Besides, during the design process, the sideslip angle has also been optimized. Consequently, two PID controllers have been designed to control the lateral–directional motion in longitudinal trim which can be discussed as block diagrams provided in Figure 3.

Figure 3.

PID controllers design configuration as block diagrams.

Accordingly, the controller design procedure can be defined as follows in four phases:

Over phase 1, the controllers with respect to the desired phase margin and gain margin were designed such that the stability of the system was guaranteed. Consequently, the desired results using a multi-objective the GA were achieved to satisfy the GA function perquisites. In this phase, simple functions were used to investigate desired results with GA. Therefore, the corresponding cost function is as follows

$$\text{J}={\left(\frac{{\mathrm{\varphi }}_{\text{m}}}{{\text{c}}_2}\right)}^2+{\left(\frac{{\text{G}}_{\text{m}}}{{\text{c}}_2}\right)}^2+{\text{c}}_3$$

(17)

where ${\mathrm{\varphi }}_{\text{m}}$ is representative of the phase margin and ${\text{G}}_{\text{m}}$ is the gain margin. ${\text{c}}_1$ , ${\text{c}}_2$ and ${\text{c}}_3$ are constant values in equation (17).

2. Since results in phase 1 were based on the simple equations for GA functions, in phase 2 smarter equations were used. Therefore, this time the fuzzy equations were used as the cost function of the GA and the corresponding results showed lower overshoot. The corresponding cost function used in this phase is given as follows

$$\text{J}=\exp {\left({\mathrm{\varphi }}_{\text{m}}-{\text{c}}_1\right)}^2+\exp {\left({\text{G}}_{\text{m}}-{\text{c}}_2\right)}^2$$

(18)

3. In phase 3, as in Dutch-Roll mode, the proposed aircraft was in level 2 of handling quality, the controllers were designed such that the handling quality placed in level 1. Therefore, the related GA was defined to assure level 1 of handling quality for all lateral–directional flight modes. It should be noted that although the results achieved from phases 1 and 2 guaranteed stability of the system, they were not able to assure the performance of the aircraft. On the other hand, the non-coupled transfer function of the lateral–directional modes is just reliable in cases where the roots of the transfer functions stand below than 3.^45,46 Hence, the new cost function to guarantee stability and performance of the proposed aircraft is given as follows

$$\text{J}={\left(\frac{{\mathrm{\varphi }}_{\text{m}}}{{\text{c}}_2}\right)}^2+{\left(\frac{{\text{G}}_{\text{m}}}{{\text{c}}_2}\right)}^2+{\left({\text{l}}_d-1\right)}^2+{\left({\text{l}}_r-1\right)}^2+{\left({\text{l}}_s-1\right)}^2+{\text{c}}_3$$

(19)

where ${\text{l}}_d$ is the aircraft handling quality level for Dutch-Roll mode, ${\text{l}}_s$ is the aircraft handling quality in spiral mode and ${\text{l}}_r$ is the aircraft handling quality in roll mode. ${\text{c}}_1$ and ${\text{c}}_2$ are constant values.

4. Over phase 4, as the sideslip angle increased drastically in phase 3 comparing to the results from phases 1 and 2, the new cost functions for GA were defined such that the best results for sideslip to be achieved as well. Consequently, at the end of phase 4, all stability and performance terms compared to the first results placed in a better condition resulting in a higher handling quality (e.g. level 1) for all lateral–directional modes. The corresponding cost function used in this step can be expressed as follows

$$\text{J}=\exp {\left({\mathrm{\varphi }}_{\text{m}}-{\text{c}}_1\right)}^2+\exp {\left({\text{G}}_{\text{m}}-{\text{c}}_2\right)}^2+\exp {\left({\text{l}}_d-1\right)}^2+\exp {\left({\text{l}}_r-1\right)}^2+\exp {\left({\text{l}}_s-1\right)}^2$$

(20)

It should be noted that in all cases, the population size of the GA was considered equal to 100 and the average relative change in the best fitness function value over max stall generations was equal to 1e^–6.

Results and discussion

Herein, results for different controllers in different phases are provided and compared for the twin-engine propeller-driven aircraft. As can be seen in Figures 4–6, the first controller showed ∼50% and ∼30% overshoot for ϕ and ψ, respectively. Over the second phase, the corresponding overshoot of ϕ decreased and reached ∼20% but in the case of ψ remained at ∼30%. However, considering the handling quality of the lateral–directional modes, the results were not desirable. Hence, in phase 3, the main goal was to compensate for the handling quality issue and reach level 1 in all lateral–directional modes. However, the corresponding results showed a drastic increase of ∼40% in the overshoot of ψ. Therefore, in phase 4 with the new controller design, the overshoot for ψ decreased to less than 30% while the overshoot of ϕ remained at ∼20%.

Figure 4.

Changes in the results of each phase in $\text{ψ}$ considering the corresponding controller.

Figure 5.

Changes in the results of each phase in $\mathrm{\varphi }$ considering the corresponding controller.

Figure 6.

Changes in the results of each phase in sideslip angle considering the corresponding controller.

In Figure 7, root locus plots of the results for $\text{ψ}$ with and without the PID controller are presented for the closed-loop system to illustrate the changes in poles and zeros of both cases. As can be seen, positive real part of the pole moved to the negative side using PID controller to enhance the stability.

Figure 7.

Root locus plot for $\text{ψ}$ in the presence of the PID controller and without the PID controller.

Finally, with respect to Figure 8, root locus plots of the results for $\mathrm{\varphi }$ show the changes in poles and zeros from the case without the PID controller to the one in the presence of the PID controller.

Figure 8.

Root Locus plot for $\mathrm{\varphi }$ in the presence of the PID controller and without the PID controller.

Also, in order to compare the performance of the GA, the particle swarm optimization (PSO) tool has been used and the corresponding results for ψ, ϕ and β using the cost function provided in phase 4 are illustrated. PSO is a direct global optimization solver and is based on a model for the flocking behaviour of birds, where each bird is representative of a ‘particle’ in the solution space.^47,48 As can be seen in Figure 9, using the PSO tool, one can design the desired PID controllers capable of guaranteeing system stability, performance and handling quality. In comparison, the PSO tool was almost three times faster than the GA tool. However, results for $\mathrm{\varphi }$ showed higher fluctuations comparing to the results achieved by the GA in phase 4.Table 5 provides information about different terms of the tuned PID controllers designed at each step.

Figure 9.

Comparison of the results for $\text{ψ}$ , $\mathrm{\varphi }$ and β using genetic algorithm and particle swarm optimization technique.

Table 5.

Proportional, integral and differential terms corresponding to each controller at different phases.

Phase	PID controller 1			PID controller 2
	${\text{K}}_p$	${\text{K}}_i$	${\text{K}}_d$	${\text{K}}_p$	${\text{K}}_i$	${\text{K}}_d$
1	0.8952	0.0101	0.5948	–0.274	–0.526	0.063
2	0.7041	0.998	2.999	–0.999	–0.5001	–3.993
3	0.9959	0.9992	2.6883	–0.50	–0.5009	–3.00
4	0.9957	0.9961	2.6890	–0.50	–1.0096	–3.00
PSO	0.7457	0.9421	3.00	–0.50	–0.501	–3.00

PID: proportional–integral–derivative; PSO: particle swarm optimization.

Sensitivity analysis

Also, in order to check the performance of the designed PID controllers in other flight conditions, a sensitivity analysis has been done using the PID controller designed in phase 4. The sensitivity analysis presented here contains all cruise flight conditions available for the aircraft type used in this study. Accordingly, different velocities as well as the highest and lowest altitudes and most forward and most aft centre of gravity (CG) positions have been studied.

In the following, Figure 10 shows the variation of the results for $\text{ψ}$ , $\mathrm{\varphi }$ and $\mathrm{\beta }$ with different Mach numbers. As can be seen, increasing the Mach number will result in a faster response from the designed PID controller. This is while, the PID controllers delay in response by decreasing the Mach number.

Figure 10.

Variation of the results for $\text{ψ}$ , $\mathrm{\varphi }$ and $\mathrm{\beta }$ with different Mach numbers.

With respect to Figure 11, increasing or decreasing the altitude will not affect the controller response for $\text{ψ}$ and $\mathrm{\varphi }$ much. However, regarding the results for $\mathrm{\beta }$ , the graph shows that at lower altitudes the designed PID controllers’ response faster.

Figure 11.

Variation of the results for $\text{ψ}$ , $\mathrm{\varphi }$ and $\mathrm{\beta }$ with different altitudes.

Finally, variation of the CG position will not affect the PID controllers’ response. Here, the most forward CG position and most aft position have been shown in Figure 12.

Figure 12.

Variation of the results for $\text{ψ}$ , $\mathrm{\varphi }$ and $\mathrm{\beta }$ with different CG values (e.g. most forward and most aft positions).

Considering all flight scenarios presented in Figures 10–12, one can conclude that the designed PID controllers with the help of GA is capable of providing handling quality of level 1 at all cruise flight conditions for the proposed aircraft.

Conclusion

In this article, it is demonstrated that with the use of multi-objective GAs, one can optimize and regulate the PID controller particularly in cases where several controllers are used and perturbation of each one affects the other ones. Also, a comparison between the GA and PSO tools presented which showed a higher accuracy of GA against the higher speed of PSO. The advantages and disadvantages of each controller used in each phase as well as their time response are provided in Table. 6.

Table 6.

Advantages and disadvantages of different controllers designed in each phase.

Phase	Advantages	Disadvantages	Run time (min)
1	Guaranteed system stabilityLower run time compared to the other phasesLower sideslip angle comparing to phases 2 and 3	Performance failureDrastic overshoot ( $\mathrm{\varphi }:44.42\text{ \%}$ , $\text{ψ}:28.22\text{ \%}$ )	25
2	Guaranteed system stabilityLower overshoot for $\mathrm{\varphi }$ ( $\mathrm{\varphi }:17.84\text{\%}$ )Lower run time comparing to phases 3 and 4	Performance failureHigher sideslip angle compared to phase 1Higher overshoot for $\text{ψ}$ ( $\text{ψ}:31.79\text{ \%}$ )	30
3	Guaranteed system stabilityGuaranteed system performance and better handling qualityLower overshoot for $\mathrm{\varphi }$ ( $\mathrm{\varphi }:15.07\text{ \%}$ )	Higher run time compared to the other phasesHigher overshoot for $\text{ψ}$ ( $\text{ψ}:37.07\text{ \%}$ )Higher sideslip angle comparing to phase 1	40
4	Guaranteed system stabilityGuaranteed system performance and better handling qualityLower overshoot for $\mathrm{\varphi }$ and $\text{ψ}$ ( $\mathrm{\varphi }:17.16\text{ \%}$ , $\text{ψ}:30.66\text{ \%}$ )Lower sideslip angle compared to phases 2 and 3	Higher run time compared to the other phases	40
PSO	Guaranteed system stabilityGuaranteed system performance and handling qualityLower sideslip angle compared to phases 2 and 3Lower run time comparing to the rest	Higher overshoot for $\mathrm{\varphi }$ and $\text{ψ}$ ( $\mathrm{\varphi }:17.98\text{ \%}$ , $\text{ψ}:31.2\text{ \%}$ )Higher fluctuation for $\mathrm{\varphi }$	14

PSO: particle swarm optimization.

Also, the sensitivity analysis to different flight conditions showed performance of the proposed controllers for various scenarios. All in all, results showed that the proposed PID controllers using GA technique is capable of providing handling quality level 1 in all cruise flight conditions.

Future work

As mentioned, several evolutionary algorithms have been developed so far to solve global optimization problems, especially for multi-objective problems. Another evolutionary algorithm that can be suggested for future work is multi-objective population-based EO which is an improved multi-objective population-based EO algorithm.^31–33

Multi-objective GA is essential a stochastic algorithm. Another area that can be of interest for future investigation is to compare the results from GA algorithms to some non-parametric statistical tests. Several works have been performed in this area on the basis of EO to how to extend the multi-objective problems to more complex practical engineering systems.^49–52