On the supporting nodes in the localized method of fundamental solutions for 2D potential problems with Dirichlet boundary condition

1.   Introduction

Nowadays, all hospitals around the world have confronted a huge surge of COVID-19 patients that require emergency health-care and services. This matter, however, has caused an enormous pressure on the available resources in many hospitals, especially with limited resources ^[1]. Meanwhile, worldwide knowledge and expertise have rapidly developed new and innovative ideas to reduce the spread of COVID-19 and treat infected patients. Acute Respiratory Distress Syndrome (ARDS) is one of the likely symptoms of coronavirus patients that can possibly cause death ^[2]. This intractable symptom can be overcome using one of the significant medical devices called an Artificial Ventilator (AV). In general, such a device is put forward to help individuals to regain healthy breathing when normal breathing is not possible. The principle of its operation is based on introducing oxygen into the lungs and withdrawing extra carbon dioxide from the body. This would enhance the gas exchange rate and hence provide comfortable breathing of sufferers. In view of technological advancement and continual modernization in introducing some technical procedures associated with control theory for medical devices, an AV system is one of these typical devices (see, e.g., ^{[3,4,5,6,7,8]}). Hence, for complementing this work, this paper proposes several optimal proportional-integral-derivative fractional-order controllers (or simply $PI^\rho D^\mu$-controllers), where $0 < \rho < 1$ and $0 < \mu < 1$ represent the orders of the integral and the differential parts of the controller, respectively. This is for the purpose of developing dynamic stabilization capabilities of the Volume-Controlled Artificial Ventilation (VCAV) system, see Figure 1.

It is known that the optimized $PI^\rho D^\mu$-controller has better performance than the traditional PID controller. As is known to all, it can offer more degrees of freedom by adding two parameters to the construction of the traditional one, as was reported in numerous references. For instance, it was demonstrated in ^[9] that one can obtain better simulation results with the use of the $PI^\rho D^\mu$-controller than with the integer-order PID controllers when we deal with some time-delay systems. In addition, the same result was presented in ^[10], when the authors designed a multivariable decoupling $PI^\rho D^\mu$-controller for variable air volume systems.

In this work, all proposed controllers have been optimally designed by carrying out two optimization techniques, namely, Particle Swarm Optimization (PSO) ^[11] and Bacteria Foraging Optimization (BFO) ^[12], for the purpose of fulfilling the high performance of the VCAV system. In particular, the key role of these two techniques is minimizing the objective function's value in light of specific constraints related to time or frequency domain, such as the system's rise time, overshoot and settling time. These constraints have a major task in quantifying the robustness of the controlled system. It is entirely normal that, once the two algorithms generate successfully the five optimal parameters ($\kappa_p$, $\kappa_i$, $\kappa_d$, $\rho$ and $\mu$) of the $PI^\rho D^\mu$-controller, the focus would be oriented towards the needs of approximating the so-called fractional-order Laplacian operator $s^\rho$ and/or $s^\mu$ (or simply $s^{\pm\gamma}$, where $\gamma = \{\rho, \ \mu\}$ and $0 < \gamma < 1$). Such approximation would be in the form of a finite integer-order rational transfer function due to the arbitrariness of $s^\gamma$ in its original form ^{[13,14,15,16,17,18]}. Actually, this new analytical function of $s^\gamma$, which permits one to design and analyze the system without a need to address some hard time-domain compositions, can be typically derived by applying multiple approximations, like Oustaloup's approximation ^[19,20], the Continued Fraction Expansion (CFE) ^[19,21] approximation and more recently the $1^{st}$- and $2^{nd}$-order El-Khazali approximations ^[22,23]. In this work, only these four approximations will be used in order to provide the operators $s^{\pm\gamma}$ with their corresponding transfer functions. Thus, we intend to tune eight $PI^\rho D^\mu$-BFO/PSO controllers by generating them through implementing the two aforesaid optimization algorithms. Numerical simulations of these eight controllers provide sufficient information to adopt the best performance of the controlled AV system, such as settling time, minimal overshoot and least rise time.

The rest of the paper is arranged as follows: An overview of the fractional-order linear time-invariant system is presented in the next section, followed by presenting some basic concepts of designing $PI^\rho D^\mu$-controllers in Section 2. Four approximations of finite-order rational transfer functions of the fractional-order integro-differential Laplacian operators are outlined in Section 3, while the final section includes the main results and all numerical outcomes of the proposed design methods.

2.   Fractional-order linear time-invariant systems

Fundamentally, the elementary principles of fractional calculus are typically utilized as a tool to transmit a class of control systems called the integer-order linear-time invariant (LTI) system into its fractional-order version, which is called the fractional-order LTI system (or simply FoLTI system) ^[24,25]. In light of some manifestations presented in ^[24,26], it has been confirmed that the FoLTI systems exceed the other integer-order counterparts by virtue of their flexibility in taking into account more additional parameters. In general, a FoLTI system can be expressed by the following fractional-order differential equation ^[24,27]:

where $y(t)$ and $u(t)$ are two variables over the time $t$ that indicate the control output and input of the system, respectively, and where $D^{\{\delta_i, \nu_k\}}$ represents the Caputo fractional differential operator of orders $\delta_i, i = 1, 2, 3, ..., n,$ and $\nu_k, k = 1, 2, 3, ..., m$, such that $n, m\in \mathbb{N}$.

Modeling system (2.1) can be accomplished by replacing $s^\pm\gamma$ by finite-order rational transfer functions, where $0 < \gamma < 1$. This would permit one to design and analyze the controlled system with finite-order controllers. The numerical method that generates such an approximate transfer function determines the bounds on the order of the controller. In a similar context, the frequency response of the LTI system can be represented by a transfer function which is the ratio of the Laplace transform of the system's output to the Laplace transform of its input for zero initial conditions ^[27], i.e.,

where $Y(s) = \mathscr{L}\{y(t)\}, \ U(s) = \mathscr{L}\{u(t)\}$ are the Laplace transforms of $y(t)$ and $u(t)$, respectively.

3.   Design of $PI^\rho D^\mu$ controllers

The principal structure of the $PI^\rho D^\mu$-controller was proposed by Podlubny et al. ^[28]. It was shown that the $PI^\rho D^\mu$-controller outperforms the traditional PID-controller. From this vantage point, many real-life technical applications have been improved by applying this controller. The structure of this controller is based upon appending two further parameters ($\rho$ and $\mu$) to the primary parameters $(\kappa_p, \kappa_i, \kappa_d)$ of the traditional PID-controller. Those two additional parameters would offer additional degrees of freedom in the controller algorithms. Nevertheless, the $PI^\rho D^\mu$-controller is definitely inferred from the following fractional-order integro-differential equation ^[24,28]:

where $e(t)$ is the error signal, $D^{\mu}$ is the Caputo operator of order $\mu$, and $I^{\rho}$ is the Riemann-Liouville operator of order $\rho$. The transfer function of the $PI^\rho D^\mu$-controller is given by

where $E(s) = \mathscr{L}\{e(t)\}$.

The subsequent mission concentrates on using the $PI^\rho D^\mu$-controller with one of the most significant industrial systems during this coronavirus time: the AV system. The key purpose of implementing such a controller is to further improve the process control of the VCAV system by enhancing its dynamic performance. This demands implementing a robust optimization algorithm to enhance the system's step response by means of optimally designing the five parameters of the $PI^\rho D^\mu$-controller. For this purpose, the BFO and PSO algorithms will be applied to determine the optimum values of these parameters through different approximations of $s^{\pm\gamma}$ that are encountered in (3.2), where $\gamma = \{\rho, \ \mu\}$, $0 < \gamma < 1$. It is necessary to establish the main objective function of the optimization algorithms. Notice that there are multiple standard objective functions that could be engaged for tuning the best parameters for the $PI^\rho D^\mu$-controller, like the Integral Absolute Error (IAE), Integral Time Square Error (ITSE), Integral Time-Absolute Error (ITAE) and Integral Square Error (ISE). For instance, minimizing any objective function value is the key goal of the chosen optimization algorithm with the aim of accomplishing the best values of the $PI^\rho D^\mu$-controller. The motivation of using the PSO and BFO algorithms comes back to their simple concepts, easy implementations, powerful abilities to control parameters and computational efficiencies when compared with other mathematical algorithms and other heuristic optimization techniques ^[29]. Nevertheless, to get a full overview about the BFO and PSO algorithms, one may consult ^{[27,30,31,32,33,34]}. The design procedure of the $PI^\rho D^\mu$-controller through the BFO and the PSO algorithms is described by the block diagram shown in Figure 2.

To obtain the best parameters of the $PI^\rho D^\mu$-controller, an objective function is established based on four terms to be minimized: the rise time, steady-state error, peak overshoot and the settling time of the controlled system. It takes the following form ^[32,35]:

where $e_{ss}$ is the steady state error, $T_s$ is the settling time, $T_r$ is the rise time, $M_p$ is the peak overshoot and $\beta$ is a scaling factor. It is relevant to note that although the scaling factor $\beta$ is typically chosen by a designer, it can definitely identify the roles of the four aforementioned items in the basic objective function ^[32,35]. In this work, this factor has been selected to be $0.5$ for the same reasons as reported in ^[32].

4.   Rational approximations of fractional-order Laplacian operators

This part briefly introduces four approximations formulated as finite-order rational transfer functions of the fractional-order integro-differential Laplacian operators, $s^{\pm\gamma}$, where $0 < \gamma < 1$. The need to employ such approximations promptly emerges after determining the best values of the $PI^\rho D^\mu$-controller by applying some optimization algorithms. Those four approximations permit one to design the target system without the need to treat some hard time-domain compositions ^[13,14].

4.1.   El-Khazali Approximations

El-Khazali proposed two practical approximations of $1^{st}$- and $2^{nd}$-orders, respectively. The $1^{st}$-order approximation is represented by a rational transfer function to replace $s^{\gamma}$, where $0 < \gamma < 1$, and given by ^[22]:

where $\omega_{cn}$ is a corner frequency, and

It is worth mentioning that the fractional-order integral Laplacian operator $s^{-\gamma}$ could be simply found from the reciprocal of (4.1) ^[22]. MATLAB Code C1, to calculate (4.1), is given in Appendix A for completeness.

In ^[23], a $2^{nd}$-order approximation was proposed for setting up a finite-order rational transfer function of the operators $s^{\gamma}$, where $0 < \gamma < 1$. The reciprocal of such approximation yields $s^{-\gamma}$. It provides an exact phase response at its center frequency and can be represented by the following modular structure ^[23]:

where $\omega_i, i = 1, 2, ..., n$, is the center frequency of each biquadratic module and $\omega_g = \sqrt[n]{\prod_{i = 1}^{n}\omega_i}$ is their geometric mean. If the leading center frequency $\omega_1$ of the beginning section is determined, then to detect the constant phase element, the next center frequencies of every section could be evaluated from the following recursive equation ^[19]:

where $\omega_x$ is the maximum real solution of the following equation:

where $\eta = \tan\frac{\gamma\pi}{4}$. Every biquadratic component in (4.3) can be determined by

where $i = 1, 2, 3, ...$, and

Note that (4.6) is the exclusive approximation that leads to $F_i(\frac{s}{\omega_i}) = (\frac{s}{\omega_i})$ as $\gamma\rightarrow1$. In other words,

Similarly, Code C2 in Appendix A simulates the approximation given by (4.6).

4.2.   Oustaloup's approximation

This approximation is the most widespread one among several approximations that are employed to generate approximate finite-order rational transfer functions for the operators $s^{\pm\gamma}$, where $0 < \gamma < 1$. The bandwidth could be defined to offer a proper adaptation to such operators by predetermining the frequency band. Typically, in order to obtain a finite-order approximation of $s^{\gamma}$ over the frequency range $(\omega_b, \omega_h)$, the following rational transfer function is implemented ^[36]:

where the gain, zeros and poles can be calculated from the following states:

It should be pointed out that the following equation allows one to compute the unity-gain geometric frequency $\omega_u$ ^[36]:

Having the form of the transfer function given in (4.9) in mind, one could observe that such a function will be always of odd order ($n = 2N+1$). Anyhow, Code C3, given in Appendix A, illustrates suitable MATLAB code for calculating a finite-order rational transfer function corresponding to $s^{\pm\gamma}$ using Oustaloup's approximation, where $0 < \gamma < 1$.

4.3.   The CFE approximation

This scheme is regarded as the principal mathematical method for offering approximate rational transfer functions for $s^{\pm\gamma}$, where $0 < \gamma < 1$. This approximation was set as follows ^[21]:

where $0 < \gamma < 1$, and $n\in\mathbb{N}$. For the sake of simplification, (4.12) can be rewritten in the following form ^[21]:

In order to find a proper finite-order approximation for $s^{\gamma}$ using (4.12), one may replace $z = (s^\gamma-1)$ in (4.12) or (4.13). This replacement step permits the $n^{th}$-order approximation of such operator to be found, about the center frequency $\omega_0 = 1\ rad/sec$, in the following form ^[26]:

where $\gamma_i\in(0, \ 1)$, which can be found in ^[37] for $i = 0, 1, ..., 5$. Two more observations should be made at the end of this part. The first one is associated with the operator $s^{-\gamma}$, which could be easily achieved by inverting upside down the formula given in (4.14). The second one is related to computing the finite-order rational transfer function corresponding to $s^{\pm\gamma}$ using the CFE approximation, which can be illustrated by tracking the MATLAB Code C4 given in Appendix A.

5.   Optimum design for VCAV System

It is common knowledge in the biomedical engineering field that the $1^{st}$-order lumped parameter model is the most streamlined model for the mechanics of the breathing system. In this model, the pressure in the upper airways is typically assumed to equal zero throughout spontaneous ventilation. The diaphragm, if it is active, may add subatmospheric pressure ^[38,39]. In general, there are commonly two actuators within many breathing systems of ventilators employed in ICU or anesthesia. These actuators are responsible for exhaling and inhaling processes. Since the AV system is just one of those breathing systems, its major components consist of these two actuators. In many circumstances, the exhaling actuator is actuated by a Positive End Expiratory Pressure (PEEP) valve, whereas the inhalator actuator, on the other hand, can be actuated by a piston drive mechanism. Also, the inhalator actuator system is linked with piston driven and lung mechanisms ^[38,39]. The overall system could be modeled by a $3^{rd}$-order differential equation by taking into account the linear friction and the electrical time constant of the motor (see ^[38,39]). Through such system, the two input variables $T_l$ and $V_a$, which respectively indicate the load torque and the applied motor voltage, are fed back to the subsystem of the lung mechanics. In particular, in view of possible estimations for all required parameters of the VCAV system, an open-loop transfer function of this system has been recently reported in ^[38] as follows:

where $Y(s) = \mathscr{L}\{y(t)\}$ and $V_i(s) = \mathscr{L}\{v_i(t)\}$ are respectively the Laplace transforms of the output variable, $y(t)$, which indicates the piston position, and the input variable of the system, $v_i(t)$, which indicates the voltage source of the motor.

The key target of this work is to enhance the dynamic behavior of the VCAV system by considering the transfer function given in (5.1). However, to show the efficiency of all feedback tuning schemes of $PI^\rho D^\mu$-controllers, the objective function given in (3.3) is minimized by running the BFO and PSO algorithms, which take into account the four approximations of $s^{\pm\gamma}$ in Section 3. The maximum value of iterations and the population size in both algorithms have been taken as 100 and 20, respectively. The optimum parameters of the $PI^\rho D^\mu$-controllers for $\beta = 0.5$ and $\beta = 0.9$, which are used to generate the corresponding approximation schemes, are given in Tables 1 and 2, respectively.

Remark 1. The scaling factor $\beta$ is chosen to define a trade-off between the overshoot and the steady-state error ($M_p+e_{ss}$) and the time difference between the settling and the rise times ($T_s-T_r$). For example, for $\beta > 0.7$ more emphasis will be given to ($M_p+e_{ss}$) than ($T_s-T_r$). The impact of $\beta$ is highlighted by several numerical results of the proposed system model for $\beta = 0.5$ and $\beta = 0.9$, reported below for completeness.

Clearly, Figures 3 and 4 show the performances of the closed-loop controlled system using both the PSO and BFO optimization algorithms by utilizing the four different approximations of the fractional-order operator repeated many times, coupled with assuming $\beta = 0.5$. In addition, we also plot in Figures 5 and 6 the performances of the closed-loop controlled system using both algorithms by assuming $\beta = 0.9$.

In both optimization techniques, the size of the controllers when using El-Khazali approximation is much less than that of the other two ones, as depicted in Tables 3 and 4. Moreover, there is a significant improvement in the step responses when using the PSO algorithm over that of the BFO algorithm. More precisely, one could notice that the minimum overshoot of the closed-loop VCAV system has been achieved when the $PI^\rho D^\mu$-controller was applying the PSO algorithm.

Tables 5 and 6, however, list the performance results of the unit-step responses of the controlled system using the four approximations for $\beta = 0.5$ and $\beta = 0.9$, respectively. It can be seen from Table 5 that the three approximations, i.e., Oustaloup, CFE, and the $2^{nd}$-order El-Khazali approximations, gave the same performance. However, the $2^{nd}$-order El-Khazali approximation yields a $4^{th}$-order $PI^\rho D^\mu$, while the other two yield $10^{th}$-order controllers. This is a significant reduction in the size of controller. A closer look at the rise time, settling time and the maximum overshoot justifies the use of the El-Khazali $2^{nd}$-order approximation. Similarly, when using the BFO optimization algorithm, the $2^{nd}$-order El-Khazali approximation was superior to the other three approximations. Therefore, both remarks satisfy the use of the biquadratic approximation described by Eqs (4.3) to (4.8).

Figures 7 and 8 shows the control signals of the controlled system using the PSO and the BFO optimization algorithms for $\beta = 0.5$, while Figures 9 and 10 show the control signals for $\beta = 0.9$. In addition, Table 7 shows the values of the objective function (3.3) using the two optimization algorithms for $\beta = 0.5$ and $\beta = 0.9$. One concludes that the PSO algorithm yields a better result than that of the BFO algorithm. Even though the initial magnitude of the control signal of the PSO algorithm is larger than that of the BFO algorithm, the value of the objective function $J$ described by (3.3) is smaller when using the PSO algorithm. This is also verified by the step responses of the controlled systems, as shown by Figures 5 and 6. It is clear that the $2^{nd}$-order El-Khazali approximation yields a smaller controller signal than the other approximations and provides a competitive steady-state error compared to its other counterparts in both algorithms and for both $\beta = 0.5$ and $\beta = 0.9$.

6.   Conclusions

A new $PI^\rho D^\mu$-controller is developed to control and improve the behavior of the Volume-Controlled Artificial Ventilation (VCAV) system. Two optimization algorithms, Bacteria Foraging Optimization (BFO) and Particle Swarm Optimization (PSO), are successfully used in conjunction with the use of four different approximations of the fractional-order integro-differential Laplacian operators, $s^{\pm\gamma}$, where $0 < \gamma < 1$. These approximations are the $1^{st}$- and the $2^{nd}$-order El-Khazali approximations, Oustaloup's approximation and the Continued Fractional Expansion (CFE) approximation. The PSO algorithm yields a better performance compared to that of the BFO algorithm, especially when using the $2^{nd}$-order El-Khazali approximations. It provided a significant improvement in terms of a smaller controller size, i.e., a $4^{th}$-order one, and in terms of the overall performance of the controlled system step response.

Acknowledgments

This work is supported by Ajman University grant 2020-COVID-19-08.

Conflict of interest

The authors declare no conflict of interest.

Appendix A

Code C1: MATLAB code for calculating the $1^{st}$-order El-Khazali approximation for $s^{\pm\gamma}$, where $0 < \gamma < 1$.

Code C2: MATLAB code for calculating the $2^{nd}$-order El-Khazali approximation for $s^{\pm\gamma}$, where $0 < \gamma < 1$.

Code C3: MATLAB code for calculating the Oustaloup approximation for $s^{\pm\gamma}$, where $0 < \gamma < 1$.

Code C4: MATLAB code for calculating the CFE approximation for $s^{\pm\gamma}$, where $0 < \gamma < 1$.