Tuning of controller parameters using Pythagorean fuzzy similarity measure for stable and time delayed unstable plants

Grunwald–Letnikov, Caputo, and Riemann–Liouville fractional definitions

The Grunwald–Letnikov derivative of function f(t) with fractional order p, here m is an integer value satisfying the condition m < p < m + 1 and a and t are the limit values, is given as Eq. (1) (Podlubny, 1999b): (1)${\displaystyle a^{D_t^p}f\left(t\right)=\sum _{k=0}^m\frac{f^{\left(k\right)}\left(a\right){\left(t-a\right)}^{-p+k}}{\Gamma \left(-p+k+1\right)}+\frac{1}{\Gamma \left(-p+m+1\right)}{\int }_a^t{\left(t-\tau \right)}^{m-p}{f}^{\left(m+1\right)}\left(\tau \right)d\tau .}$

Here, the Gamma (Γ) function is defined by Eq. (2): (2)${\displaystyle \Gamma \left(x\right)={\int }_0^{\infty }{t}^{x-1}{e}^{-t}dt,\mathfrak{R}\left(Z\right)>0.}$

The Caputo derivative of function f(t) with fractional order p is given as Eq. (3) (Shah & Agashe, 2016): (3)${\displaystyle a^{D_t^p}f\left(t\right)=\frac{1}{\Gamma \left(m-p\right)}{\int }_a^t\frac{f^{\left(m\right)}\left(\tau \right)}{{\left(t-\tau \right)}^{p-m+1}}d\tau .}$

Here, m − 1 < p < m, m ∈ N, a,  and t are the limit values of integration.

The Riemann–Liouville derivative of function f(t) with fractional order p, here m is an integer value satisfying the condition m − 1 < p < m, m ∈ N, is given as Eq. (4) (Shah & Agashe, 2016): (4)${\displaystyle a^{D_t^{p}f\left(t\right)}=\frac{1}{\Gamma \left(m-p\right)}\frac{d^m}{d{t}^m}{\int }_a^t\frac{f\left(\tau \right)}{{\left(t-\tau \right)}^{p-m+1}}d\tau .}$

Integer approximation methods for fractional order operators

Fractional calculus has significant advantages in modeling and control system performance, but it brings higher computational complexity than integer order calculations because the fractional derivative operator is a non-local operator that includes all the past values of the function. This means that more memory elements are needed to store all the past values. To overcome this computational complexity, integer approximation methods have been developed to approximate the responses of fractional order elements within limited operating ranges (Deniz et al., 2020). Continuous fraction expansion (CFE) (Khovanskii, 1963), Outstalup (Oustaloup et al., 2000), and Matsuda (Matsuda & Fujii, 1993) are examples of well-known integer-order models in the literature.

In the continued fraction expansion method (CFE), s^α is determined by expanding the expression (1 + x)^α. The expansion of the expression (1 + x)^α satisfying the condition 0 < α < 1 is determined by Eq. (5) (Zawadzki & Włodarczyk, 2017): (5)${\displaystyle {\left(1+x\right)}^{\alpha }=\frac{1}{1-\frac{\alpha x}{1+\frac{\left(1+\alpha \right)x}{2+\frac{\left(1-\alpha \right)x}{3+\frac{\left(2+\alpha \right)x}{2+\frac{\left(2-a\right)x}{5+\dots }}}}}}}$

The expansion of (1 + x)^α can alternatively be expressed by Eq. (6): (6)${\displaystyle {\left(1+x\right)}^{\alpha }=\frac{1}{1-}\frac{\alpha x}{1+}\frac{\left(1+\alpha \right)x}{2+}\frac{\left(1-\alpha \right)x}{3+}\frac{\left(2+\alpha \right)x}{2+}\frac{\left(2-\alpha \right)x}{5+\dots }.}$

Here, the 1st-degree continued fraction expansion of s^α is expressed as Eq. (7) by writing s − 1 instead of x: (7)${\displaystyle S^{\alpha }\cong \frac{\left(1+\alpha \right)s+\left(1-\alpha \right)}{\left(1-\alpha \right)s+\left(1+\alpha \right)}.}$

The Oustaloup method gives an integer-order approximation of the fractional-order operator s^α using integer-order transfer functions in a certain frequency range [ωl, ωh] in specified lower and upper limits (Marushchak et al., 2022). Oustaloup approximation of the fractional-order operator s^α given as Eq. (8): (8)${\displaystyle s^{\alpha }={\left(\frac{{\omega }_u}{{\omega }_h}\right)}^{\alpha }\prod _{k=-N}^{k=N}\frac{1+s/{\omega }_k^{\prime }}{1+s/{\omega }_k}.}$

Here, N is the approximation order, ${\omega }_u=\sqrt{{\omega }_l{\omega }_h}$ and ${\omega }_k^{\prime }$, ω_k are the zeros and poles of the equivalent integer order transfer functions, respectively. ${\omega }_k^{\prime }$ and ω_k are given by Eqs. (9) and (10): (9)${\displaystyle {\omega }_k^{{}^{\prime }}={\omega }_l{\left(\frac{{\omega }_h}{{\omega }_l}\right)}^{\left(k+N+0.5-0.5\cdot \alpha \right)/\left(2N+1\right)}}$ (10)${\displaystyle {\omega }_k={\omega }_l{\left(\frac{{\omega }_h}{{\omega }_l}\right)}^{\left(k+N+0.5+0.5\cdot \alpha \right)/\left(2N+1\right)}.}$

The coefficient of the transfer function of the 2N + 1 degree approximation is expressed by Eq. (11): (11)${\displaystyle k_{\Pi }={\left(\frac{{\omega }_u}{{\omega }_h}\right)}^{\alpha }.}$

The integer order Oustaloup approximation of s^α can be expressed by Eq. (12): (12)${\displaystyle s^{\alpha }\cong k_{\Pi }\prod _{k=-N}^N\frac{s+{\omega }_k^{{}^{\prime }}}{s+{\omega }_k}=k_{\Pi }\frac{B_n{s}^n+B_{n-1}{s}^{n-1}+\dots +B_{1}s+B_0}{A_n{s}^n+A_{n-1}{s}^{n-1}+\dots +A_{1}s+A_0}.}$

Based on the CFE method, the Matsuda method approximates the fractional-order operator s^α to integer order by calculating the gain at logarithmically spaced frequencies. Logarithmically spaced frequencies, here n approximation order and ω_l and ω_h lower and upper-frequency values, are given as Eq. (13) (Koseoglu et al., 2021): (13)${\displaystyle {\omega }_k=\left\{\begin{array}{c}{\omega }_l{\begin{array}{c}\hfill \left(\frac{{\omega }_h}{{\omega }_l}\right)\hfill \\ \hfill \hfill \end{array}}^{\frac{k-1}{n-1}},n>1.\hfill \\ {\omega }_h,n=1\hfill \end{array}\right.}$

The gains corresponding to these frequencies are expressed by Eqs. (14) and Eq. (15) (Krishna, 2011): (14)${\displaystyle d_o\left({\omega }_k\right)=\left|G\left(j{\omega }_k\right)\right|}$ (15)${\displaystyle d_n\left({\omega }_k\right)=\frac{{\omega }_k-{\omega }_{n-1}}{d_{n-1}\left({\omega }_k\right)-d_{n-1}\left({\omega }_{n-1}\right)}.}$

Then, gain values are placed in a (N + 1)(N + 1) upper triangular matrix. The obtained upper triangular matrix is expressed as Eq. (16) (Krishna, 2011): (16)${\displaystyle D=\left[\begin{array}{ccccc}\hfill d_0\left({\omega }_0\right)\hfill & \hfill d_0\left({\omega }_1\right)\hfill & \hfill d_0\left({\omega }_2\right)\hfill & \hfill \cdots \hfill & \hfill d_0\left({\omega }_n\right)\hfill \\ \hfill \hfill & \hfill d_1\left({\omega }_1\right)\hfill & \hfill d_1\left({\omega }_2\right)\hfill & \hfill \cdots \hfill & \hfill d_1\left({\omega }_n\right)\hfill \\ \hfill \hfill & \hfill \hfill & \hfill d_2\left({\omega }_2\right)\hfill & \hfill \cdots \hfill & \hfill d_2\left({\omega }_n\right)\hfill \\ \hfill \cdots \hfill & \hfill \cdots \hfill & \hfill \cdots \hfill & \hfill \cdots \hfill & \hfill \cdots \hfill \\ \hfill \cdots \hfill & \hfill \cdots \hfill & \hfill \cdots \hfill & \hfill \cdots \hfill & \hfill d_n\left({\omega }_n\right)\hfill \end{array}\right].}$

The approximate integer order model’s gain coefficients can be expressed as Eq. (17) using the diagonal elements D_kk of matrix D (Krishna, 2011): (17)${\displaystyle {\alpha }_k=D_{kk}=\left\{\begin{array}{c}\left|G\left(j\omega \right)\right|ifk=0\hfill \\ \frac{{\omega }_k-{\omega }_{k-1}}{d_{k-1}\left({\omega }_k\right)-d_{k-1}\left({\omega }_{k-1}\right)}ifk=1,2,\dots ,n\hfill \end{array}\right.}$

Then, using the CFE method, the integer degree approximation of the fractional order operator s^α is expressed by Eq. (18). (18)${\displaystyle s^{\alpha }={\alpha }_0+\frac{s-{\omega }_0}{{\alpha }_1+\frac{s-{\omega }_1}{{\alpha }_2+\frac{s-{\omega }_2}{{\alpha }_3+\dots }}}={\alpha }_0+\frac{s-{\omega }_0}{{\alpha }_1+}\frac{s-{\omega }_1}{{\alpha }_2+}\cdots }$

Fractional PID and PI-PD controllers

The fractional PID controller has five parameter values with λ (integrator) and µ (derivative) parameters in addition to the K_p, K_i, K_d parameters that exist in conventional PID. These additional fractional controller parameters allow the controller to operate over a wide range. Also, the FOPID controller is more robust than conventional PID controllers when the parameters of the system being controlled change (Agarwal et al., 2019). The block diagram of the FOPID controller is shown in Fig. 1.

The expressions of the FOPID controller output in the time domain and s-plane are given by Eqs. (19) and (20), respectively: (19)${\displaystyle u\left(t\right)=K_{p}e\left(t\right)+K_i\frac{d^{-\lambda }}{d{t}^{-\lambda }}e\left(t\right)+K_d\frac{d^{\mu }}{d{t}^{\mu }}e\left(t\right)}$ (20)${\displaystyle u\left(s\right)=K_{p}e\left(s\right)+K_i{s}^{-\lambda }e\left(s\right)+K_d{s}^{\mu }e\left(s\right),\left(\lambda ,\mu >0\right)}$

In the PI-PD controller, which is obtained by combining a forward PI controller and an internal PD controller with feedback, the PD controller part transforms an open-loop unstable system into an open-loop stable system (Tan, 2009), whereas the PI controller part has a significant impact on enhancing the stability process (Kaya, 2003). The block diagram of the system including the PI-PD controller is shown in Fig. 2.

Four parameter values are K_p, K_i, K_d, and K_f for the PI–PD controller. The transfer functions of PI-PD controller components are given as Eqs. (21) and (22): (21)${\displaystyle C_{PI}\left(s\right)=K_p+\frac{K_i}{s}=\frac{K_{p}s+K_i}{s}}$ (22)${\displaystyle C_{PD}\left(s\right)=K_f+K_{d}s}$

Pythagorean fuzzy similarity measure

Zadeh (1965) considering the uncertainty in decision-making, proposed the concept of fuzzy sets, and great success has been achieved in various areas involving uncertainty. The fuzzy set A, defined in the universal set X = x₁, x₂, x₃, …, can be expressed by Eq. (23): (23)${\displaystyle A=〈x,{\mu }_A\left(x\right)〉|x\in X}$

Here, μ_A(x) is a membership function in the range [0, 1] and is a measure of how compatible an element x is with the fuzzy set A. The most common types of membership functions used in fuzzy applications are triangle, trapezoidal, and Gaussian (Can & Ozguven, 2017). Researchers have developed different fuzzy set extensions to deal with the increasing fuzziness and uncertainty of real-world data. For this purpose, intuitionistic fuzzy sets (IFS) were introduced by Atanassov (1999) as an extension of the concept of fuzzy sets. The intuitionistic fuzzy set A defined in X is expressed as Eq. (24): (24)${\displaystyle A=〈x,{\mu }_A\left(x\right),{\nu }_A\left(x\right)〉|x\in X}$ μ_A(x) and ν_A(x) are membership functions showing membership and non-membership degrees, respectively, and take values between [0, 1]. Also, the condition $0\le {\mu }_A\left(x\right)+{\nu }_A\left(x\right)\le 1$ must be satisfied and ${\pi }_A\left(x\right)$ (degree of hesitation) is expressed by Eq. (25): (25)${\displaystyle {\pi }_A\left(x\right)=1-{\mu }_A\left(x\right)-{\nu }_A\left(x\right).}$

There may be cases where the sum of ${\mu }_A\left(x\right)+{\nu }_A\left(x\right)$ in IFS is greater than 1. Pythagorean fuzzy set (PFS) is proposed by Yager (2014) to deal with uncertainty satisfying the conditions ${\mu }_A\left(x\right)+{\nu }_A\left(x\right)\le 1$ or ${\mu }_A\left(x\right)+{\nu }_A\left(x\right)\ge 1$. A Pythagorean fuzzy set defined in X is given as Eq. (26). It must also satisfy the $0\le {\mu }_A^2\left(x\right)+{\nu }_A^2\left(x\right)\le 1$ condition. ${\pi }_A\left(x\right)$ (degree of hesitation) is expressed by Eq. (27) In this way, in decision-making problems; PFS significantly increased the working range compared to IFS (Bryniarska, 2020): (26)${\displaystyle A=〈x,{\mu }_A\left(x\right),{\nu }_A\left(x\right)〉|x\in X}$ (27)${\displaystyle {\pi }_A\left(x\right)=\sqrt{1-{\mu }_A^2\left(x\right)-{\nu }_A^2\left(x\right)}.}$

The similarity measure is used to determine the degree of similarity between two or more sets. There are distance-based, probability-based, fuzzy set theory-based, and graph theory-based approaches to similarity measurement (Can & Ozguven, 2017). The similarity measure is usually expressed as a numeric value. It is higher when the data samples are more similar.

X = x₁, x₂, x₃, … is a universal set, A and B are two fuzzy sets defined in X, the Pythagorean distance measure between fuzzy sets A and B can be expressed by Eq. (28) (Peng, 2019): (28)${\displaystyle d\left(A,B\right)=\sqrt[P]{\frac{1}{2n{\left(t+1\right)}^P}\sum _{i=1}^n\begin{array}{c}\hfill \text{(}{\left|\left(t+1-a\right)\left({\mu }_A^2\left(x_i\right)-{\mu }_B^2\left(x_i\right)\right)-a\left(v_A^2\left(x_i\right)-v_B^2\left(x_i\right)\right)\right|}^P+\hfill \\ \hfill {\left|\left(t+1-b\right)\left(v_A^2\left(x_i\right)-v_B^2\left(x_i\right)\right)-b\left({\mu }_A^2\left(x_i\right)-{\mu }_B^2\left(x_i\right)\right)\right|}^P\text{)}\hfill \end{array}}}$ p is the L_p standard and t, a,  and b are parameters that satisfy the condition a + b ≤ t + 1 , 0 < a, b ≤ t + 1, t > 0, and show the uncertainty level. The similarity measure between A and B is expressed by Eq. (29): (29)${\displaystyle S\left(A,B\right)=1-d\left(A,B\right).}$

Optimization methods

Ant colony optimization for continuous domains

Ant colony optimization is a population-based metaheuristic optimization algorithm inspired by the foraging behavior of ants (Dorigo, 1992). When an ant finds a food source, it evaluates its quality and quantity and carries some of it back to its nest. On its way back to its nest, it leaves a substance called chemical pheromone trails on the road. The amount of pheromone left is related to the quantity or quality of the food found. These trails guide other ants to the same food source. Communication between ants through pheromone trails helps them find the shortest path (Ojha, Abraham & Snášel, 2015). Initially developed for discrete optimization problems, the ant colony algorithm was extended to solve continuous optimization problems by Socha & Dorigo (2008). Continuous-time ant colony optimization (ACO_R) is a population-based metaheuristic algorithm that iteratively generates solutions.

In the ACO_R algorithm, the pheromone information is stored in a table called a solution archive. The solution archive is shown in Fig. 3. The solution archive contains k solutions, each of n dimensions, initially initialized with random values. A solution s_j in the archive represents a solution vector. The ith variable of the jth solution is denoted by $s_j^i.f\left(s_j\right)$ represents the value of the cost function. These solutions are added sequentially according to their cost functions as f(s₁) ≤ f(s₂) ≤ … ≤ f(s_k).

In Fig. 3, ω represents the weight vector. The weight vector ω is a Gaussian function. The best solution obtains the maximum weight. The weight value ω of the solution s_j is calculated by Eq. (32). q is a tuning parameter that affects the convergence rate of the algorithm. If it is small, the area around the best-ranked solutions is searched and an early convergence may occur (Afshar & Madadgar, 2008): (32)${\displaystyle {\omega }_j=\frac{1}{qk\sqrt{2\pi }}{e}^{\frac{-{\left(rank\left(j\right)-1\right)}^2}{2q^2{k}^2}}.}$

To select and update a solution in the solution table, the probability of selecting each row in the solution table is calculated given the weight vector ω (Eq. (33)): (33)${\displaystyle p_j=\frac{{\omega }_j}{\sum _{r=1}^k{\omega }_r}.}$

After the selection is made based on the solution probability value, the parameters of the solution are updated based on the σ value given in Eq. (34): (34)${\displaystyle {\sigma }_j^i=\xi \sum _{r=1}^k\frac{\left|s_r^i-s_j^i\right|}{k-1}.}$ σ is the average distance between the selected solution and the other solutions. ξ is the parameter of the algorithm known as the evaporation rate, and as ξ increases, the convergence of the algorithm decreases. In the final stage of ACO_R, m newly generated solutions are added to the initial k solutions. The cost function values of k + m solutions are ranked in ascending order, and in the next step, m worst solutions are removed. Thus, the solution archive is updated with the values of k best solutions. The algorithm is executed until the termination criterion is met and the optimum result is obtained. The flowchart of the ABC and ACO_R algorithm is given in Fig. 4.

The proposed cost function for tuning of controller parameters

When determining the controller parameters based on the optimization, it usually uses cost functions based on the control error variable e. However, e’s minimization alone is not enough for its effective performance. In addition to minimizing e, the system’s transient and steady-state properties (undershoot, rise time, peak time, overshoot, settling time, and steady-state error value) must also be kept at a minimum value. Here, each of the undershoot, rise time, peak time, overshoot, settling time, and steady-state error value values are criteria that must be brought to a certain value. Considering all the values that the controller parameters will take into account, these mentioned parameters will also take a large number of values. In this respect, determining the best system response according to the values of all the transient and steady-state properties of the system is a multi-criteria decision-making problem (Can & Ozguven, 2017). There are many different methods for MCDM, but fuzzy similarity measure-based MCDM is a preferred method for the tuning of controller parameters in recent years. In this method, the ideal solution must be known. The optimum solution is reached by minimizing the distance between the instant solution and the ideal solution.

In the study, the Pythagorean fuzzy similarity measure was used for the information on the distance between the ideal solution and the instantaneous solution, as it allows evaluation the positive and negative distances together. This distance information was defined as a cost function and its minimization was performed with optimization algorithms. GWO, PSO, ABC, and ACO_R optimization algorithms were tried, but the desired results could not be obtained with GWO and PSO. On the other hand, the results obtained with ABC and ACO _R could be obtained at a satisfactory level. For this reason, ABC and ACO_R were preferred in the study. For the ABC algorithm, the number of food sources is 30, and the trial limit parameter is 90. For the ACOR algorithm, the archive size is 30, the sample size is 40, the q parameter value is 0.5, and Pythagorean fuzzy similarity measure parameters are p = 1, a = 1, b = 2, t = 3 (Peng, 2019). In both optimization algorithms, the number of iterations is 100 and the number of populations is 30.

The most commonly used cost functions in Table 1 were selected to be used in these optimization algorithms, and the cost function obtained with the proposed method was compared with these cost functions. In the equations in Table 1, e (t) represents the difference (error) sign between the reference value applied to the controlled system and the output of the system, and t is the time variable.

The first step in the methodology of the study is to determine the ideal solution for MCDM. The ideal solution is determined by the opinion of an expert. The study, it is aimed to determine the controller parameters that will provide the most suitable transient and steady-state properties for the system output with MCDM. For this reason, the ideal solution is the transient and steady-state properties of the system determined according to expert opinion. Since the Pythagorean fuzzy similarity is used in the study, the ideal solution and the instantaneous solutions that occur during the optimization should be in the Pythagorean fuzzy space. This is done by fuzzification of the real values. Fuzzification is done using Pythagorean fuzzy membership functions for each transient and steady-state feature. The Pythagorean fuzzy similarity between the ideal solution (or fuzzy set) and the fuzzy instantaneous solution (or instantaneous fuzzy set) is calculated and the instantaneous distance is found. With this distance information, a cost function is defined in the study. Minimization of this cost function was performed with ABC and ACO_R optimization algorithms. The parameters found at the end of the optimization are the controller parameters sought. The methodology of the study is described below.

First, based on expert opinion, Pythagorean membership functions are created individually for parameters such as rise time, settling time, overshoot, undershoot, peak time, and steady-state error, etc. to satisfy the condition $0\le {\mu }^2\left(x\right)+{\nu }^2\left(x\right)\le 1$. After, K_p, K_i, K_d, λ, μ, K_f values are initially generated at random. These generated values are provided to the controller selected based on the controller type. Then, the step input is applied to the system, and the transient and steady-state response characteristics (rise time, settling time, overshoot, undershoot, peak time, and steady-state error) are obtained. Each of the obtained system response parameters is passed to membership functions. The real (or instantaneous) Pythagorean fuzzy set B, which consists of member and non-member degrees, is the fuzzy set A, which consists of targeted ideal values. ${\displaystyle B_1=〈x,{\mu }_{B-RiseTime}\left(x\right),{\nu }_{B-RiseTime}\left(x\right)〉|x\in X}$ ${\displaystyle B_2=〈x,{\mu }_{B-SettlingTime}\left(x\right),{\nu }_{B-SettlingTime}\left(x\right)〉|x\in X}$ ${\displaystyle B_3=〈x,{\mu }_{B-Peak}\left(x\right),{\nu }_{B-Peak}\left(x\right)〉|x\in X}$

The resulting real Pythagorean fuzzy set vector B is obtained as B = [ B₁, B₂, B₃ …]. The ideal set vector A is obtained as A = [ A₁, A₂, A₃], here A₁, A₂, and A₃ represent the ideal sets created for each parameter. Here, A₁ = 〈x, 1, 0〉|x ∈ X, A₂ = 〈x, 1, 0〉|x ∈ X, A₃ = 〈x, 1, 0〉|x ∈ X. The Pythagorean distance measure between sets A and B is determined by Eq. (28). Then, using the distance measure obtained, the Pythagorean fuzzy similarity measure between these two sets is determined by Eq. (29). In the proposed method, the value of the Pythagorean similarity measure is used as the cost function that optimization algorithms attempt to minimize and cost function is given as Eq. (35): (35)${\displaystyle cost=1-S\left(A,B\right)}$

The optimal K_p, K_i, K_d, λ, μ, K_f values for the controller are determined by running optimization algorithms until the cost value reaches the target value or up to a determined number of iterations. The method of the study is summarized below in the form of processing steps.

Step 1: Randomly determine the values of K_p, K_i, K_d, λ, μ,  and K_f within the specified lower and upper bounds as defined by the optimization algorithm.

Step 2: Provide these values to the controller.

Step 3: Apply a step input to the system.

Step 4: Obtain the response parameters from the system: rise time, settling time, overshoot, undershoot, peak time, steady-state error, etc.

Step 5: Each parameter in the response is evaluated using the corresponding membership function from the B vector, and fuzzy sets are obtained for each parameter.

Step 6: Calculate the Pythagorean distance between sets A and B using Eq. (28).

Step 7: Using this distance, calculate the Pythagorean similarity measure between sets A and B by Eq. (29).

Step 8: Utilize this similarity measure as a cost function which optimization algorithms aim to minimize.

Step 9: Iteratively tune the values of K_p, K_i, K_d, λ, μ starting from Step 2 until the value of the cost function falls below a pre-determined acceptable target.

Z-type and S-type membership functions generated for an example rise time are shown in Table 2 and Fig. 5. µis the Z-type membership function denoting the membership degrees and ν is the S-type membership function denoting the non-membership degrees.

For example, if the rise time value is obtained to be 1.8 s after applying the step input to the system, the membership degree of this value to the µmembership function is calculated as 0.08 and the membership degree to the ν membership function is calculated as 0.18 according to the ranges determined in Table 2. In this case, the real fuzzy Pythagorean set B of the rise time is obtained as B = (0.08, 0.18). The membership degree of any value corresponding to the red-orange area in Fig. 5 constitutes the ideal set A. These membership degrees also correspond to the blue area. In other words, a rise time value of 1 s or less is the ideal target. In these ranges, the degree of membership of the rise time to the function µis 1, the degree of membership to the function ν is 0, and the ideal set is shown as A = (1, 0). The similarity measure between sets A and B is calculated by Eq. (29) and is found as 0.3668. In another case, if the rise time value is 1.1 s, the real fuzzy set B is obtained as (0.98.0), and the similarity measure of set B to set A is found as 0.9752. These values are a measure of how similar set B is to set A. The closer the value of the similarity measure is to 1, the higher the similarity between the two sets. In this example, for set B to approach the ideal Pythagorean fuzzy set A, the incoming rise time value must be less than or close to 1 s. In this case, the real fuzzy set B will be close to A = (1.0). If desired, this process is repeated for other system response characteristics (settling time, overshoot, undershoot, peak, etc.), and a total similarity measure is obtained with the help of Eq. (29). This similarity measure value is created a cost function with Eq. (35).

Example 1

In this example, a second-order non-minimum phase system is used (Bingul & Karahan, 2018). The transfer function of the controlled system is expressed by Eq. (36). This system is difficult to control due to the undershoot and time delay caused by zero in the right half s plane. (36)${\displaystyle G\left(s\right)=\frac{1-5s}{\left(1+10s\right)\left(1+20s\right)}.}$

For the control of this system, the 4th-integer Matsuda approximation method is used in the FOPID controller. In the Matsuda approximation, the lower and upper frequencies are chosen as [10⁻², 10²]. The lower and upper boundary values used in the ABC and ACO_R algorithms to be used in setting the FOPID controller parameters are shown in Table 3.

Membership functions were created separately for steady-state error, settling time, peak, and undershoot. The range values of the membership functions created for the system are shown in Table 4. The graphical representations corresponding to these range values are given in Fig. 6. It is desirable that the settling time is less than or close to 9 s, the steady-state error is 0.020 or less, the peak value is between 1 and 1.016, and the undershoot value is between 150 and 200.

After optimization, the FOPID controller parameters were obtained, as shown in Table 5, using ITAE, ITSE, IAE, ISE, and the proposed method. The step response characteristics corresponding to the controller parameters are given in Table 6, and the graphical representation of the step response is given in Fig. 7.

The study for Example 1, it is aimed to reduce the undershoot using the proposed method. In addition, the overshoot, settling time, and rise time are desired to be kept as small as possible. Table 6 shows that the IAE cost function stands out among the classical cost functions used in both optimization algorithms. It can be said that the proposed method significantly reduces the undershoot compared to the IAE cost function (156.5210% for ABC and 154.0753% for ACO _R). In terms of rise time and settling time, the proposed method is very close to the IAE cost function. Overall, it can be said that the proposed method provides a better control performance.

Example 2

In this example, a first-order unstable system with time delay is used (Zheng, Huang & Zhang, 2019). The transfer function of the controlled system is expressed by Eq. (37). Here is a time delay of 2 s in the system. (37)${\displaystyle G\left(s\right)=\frac{4}{4s-1}{e}^{-2s}.}$

PI-PD controller is used to control this system. The lower and upper boundary values used in ABC and ACO_R algorithms are shown in Table 7.

Membership functions were created separately for rise time, settling time, peak, and steady-state error. The range values of the membership functions created for the system are given in Table 8 and the graphical representations corresponding to these range values are given in Fig. 8. According to Table 8, it is desired that the rise time value is between 2.7 and 3.2 s, the settling time value is between 10 and 12.3 s, the peak value is between 1 and 1.004, and the steady-state error is below 0.020.

As a result of the optimizations, the PI-PD controller parameters were obtained as in Table 9. The step response characteristics corresponding to these controller parameters are given in Table 10, and Fig. 9 shows the step response graphically.

In Example 2, aimed to reduce the overshoot and keep the settling time and rise time as small as possible with the proposed method. Among the classical cost functions, the lowest overshoot was obtained by using the ITAE cost function in the ABC algorithm (0.5446%). When the proposed method is used in the ABC and ACO _R algorithms, much lower overshoots of 0.0350% and 0.0117% occur, respectively. The duration and amplitude of the oscillation are less when the proposed method is used compared to other cost functions. In addition, the target values for rise time and settling time are achieved with this method.

Example 3

In this example, a fractional-order unstable system with time delay is used (Zheng, Huang & Zhang, 2019). The transfer function of the controlled system is expressed by Eq. (38). A delay of 0.2 s exists in the system. The Matsuda 4th-integer approximation method was used for fractional-order expressions. In the Matsuda approximation, the lower and upper frequencies are chosen as [10⁻², 10²]. (38)${\displaystyle G\left(s\right)=\frac{1}{s^{1.2}-1}{e}^{-0.2s}.}$

PI-PD controller is used to control this system. The lower and upper boundary values used in ABC and ACO_R algorithms are given in Table 11.

Membership functions were created separately for rise time, settling time, peak value, and steady-state error. The range values of the membership functions created for the system are given in Table 12 and the graphical representations corresponding to these range values are given in Fig. 10. It is desired that the rise time value is between 0.2 and 0.3, the settling time value is between 1.2 and 1.4 s, the peak value is between 1.010 and 1.020, and the steady-state error is 0.020 or less.

As a result of the optimization, the PI-PD controller parameters were obtained as in Table 13. The step response characteristics corresponding to these controller parameters are shown in Table 14, also with their graphical representation in Fig. 11.

The proposed method is aimed to reduce the overshoot and settling time and to ensure the rapid rise of the system by keeping the rise time at a low value. In the study for Example 3, a high overshoot occurs in classical cost functions. Using the proposed method significantly reduced the overshoot (1.6146% for ABC and 1.9951% for ACO_R). When Table 14 is examined, the closest values to the proposed method are obtained with the IAE cost function. However, when this cost function is used, the overshoot and oscillation amplitude remain high. The proposed method is more successful than other cost functions in terms of both settling time and overshoot.