PID Control for Nonlinear Processes

Abstract

This chapter presents a proportional-integral-derivative (PID) Takagi-Sugeno fuzzy system controller that can be trained by the particle swarm optimization-cuckoo search (PSOCS) technique to control nonlinear multi-input multi-output (MIMO) systems. Instead of the standard methods that are widely used in the literature, the PSOCS is used to adjust all of the PID parameters by the minimization of a given objective function. A nonlinear MIMO system has been selected to be controlled by this controller. The simulation results show the notable control accuracy and generalization ability of this MIMO controller. Finally, a comparative study with a PSO algorithm and CS algorithm shows the superiority of the PSOCS over these two optimization methods in terms of guaranteeing the desired performance.

1. Introduction

The proportional-integral-derivative (PID) control technique is used in this chapter to achieve the control objective of making the output signals follow the desired trajectory or arrive at the required places accurately and fast. However, there are a number of characteristics in big system models that prevent the direct application of PID control approaches. The largest (and unknown) model order, hazy connections between subsystems, wide parameter fluctuation, and complex organizational structure stand out among these characteristics. There has been a lot of research done on the stability problem with fuzzy control systems. However, there are frequently ambiguities in many real-world systems, which are a source of instability. Thus, the robust fuzzy stabilization problem for uncertain nonlinear systems has received considerable interest [1, 2]. We can use fuzzy logic theory to break down the effort of modeling and control design into a collection of simpler local tasks by using qualitative, linguistic information about a complex nonlinear system. Additionally, it offers a method for combining these local activities to produce the overall model and control design. On the other hand, improvements in the theory of linear systems have led to the development of numerous potent design tools. As a result, the analysis and controller synthesis of the nonlinear system may be done using the productive linear system theory based on the linear Takagi Sugeo (TS) fuzzy model [3]. Although the local design structure is used to build the fuzzy controller, the global design conditions should be used to calculate the feedback gains. To ensure global stability and control performance, the global design conditions are required.

PID control with feedback signals accessible at the site of each controlled device is typically the most advantageous for the nonlinear problem due to the ease of the practical implementation of the controllers [4, 5]. The difficulty of precisely and automatically tweaking the gains of PID controller poses a significant obstacle to the efficient application of this method. Because, it is a computationally expensive combinational optimization issue, and because it may be difficult, time-consuming, and process-specific to extract the right set of static benefits for every subsystem. In order to overcome these drawbacks, several optimization algorithms are used, for example, genetic algorithm (GA), equilibrium optimizer algorithm (EOA) [6], particle swarm optimization (PSO), arithmetic optimizer (AO) [7], and cuckoo search (CS).

For the best tuning of the PID gains in each control region and to accelerate the convergence of the algorithms in this work, the particle swarm optimization-cuckoo search (PSOCS) technique, a mix of the PSO and the cs algorithm, is utilized in this chapter. PSOCS, a unique population-based metaheuristic, has become an effective tool for engineering optimization because it makes use of the swarm intelligence produced by the collaboration and competition among the particles in a swarm. Additionally, it has proven to be effective in resolving issues with nonlinearity, non-differentiability, and high dimensionality [8].

2. Identifying a MIMO system

2.1 Representation of MISO systems

Consider a MIMO system with $n$ inputs and $n$ outputs. After decomposition, each MISO system will be described by:

$$y_i\left(k+1\right)=f_i\left(x\left(k\right)\right),i=1,2,\dots ,n$$E1

with the regression vector given by:

$$x_i\left(k\right)=\left[y_i{\left(k\right)}_0^n,u_1{\left(k-d_{1i}\right)}_0^n,\dots ,u_n{\left(k-d_{\mathit{ni}}\right)}_0^n\right]$$E2

here $k$ denotes the discrete time sample, $n$ is an integer relative to the order of the system, and $d_{\mathit{ij}}$ is the pure delay. MISO systems are independently identified, for simplicity of notation the index i will be ignored. The output of the system is written as:

$$y\left(k+1\right)=\mathit{Ay}\left(k\right)+\mathit{Bu}\left(k\right)+\alpha$$E3

where $\alpha$ is a constant. A fuzzy TS model makes an attempt to approximate the unknown function $f\left(.\right)$. By using the Fuzzy C-Means (FCM) algorithm, a set of fuzzy regions, in this case characterized by Gaussian membership functions, are created in the data space, and the ensuing portions describe how the system behaves in these areas. The MISO system’s governing body will now be

$$R_l:\mathit{\text{if}}x\left(k\right)\mathit{\text{is}}{\mathrm{\Omega }}_l\mathit{\text{then}}{y}^l\left(k+1\right)=A_{l}y\left(k\right)+B_{l}u\left(k\right)+{\alpha }_l,l=1,\dots ,r$$E4

where ${\mathrm{\Omega }}_l$ is the antecedent fuzzy set of the $l^{\mathit{th}}$ rule, $A_l=\left[A_{l1},A_{l2},\dots ,A_{\mathit{ln}}\right]$ and $B_l=\left[B_{l1},B_{l2},\dots ,B_{\mathit{ln}}\right]$ are the vectors of the two polynomials $A_l$ and $B_l$, and $r$ is the number of rules.

Rule (4) can then be written as:

$$\begin{array}{c}{R}_l:\mathit{\text{if}}{x}_1\left(k\right)\mathit{\text{is}}{\mathrm{\Omega }}_{l1}\mathit{\text{and}}{x}_p\left(k\right)\mathit{\text{is}}{\mathrm{\Omega }}_{\mathit{lp}}\mathit{\text{then}}\\ y^l\left(k+1\right)=A_{l}y\left(k\right)+B_{l}u\left(k\right)+{\alpha }_l,l=1,\dots ,r\end{array}$$E5

with $p=2n^2+1$

The output of the TS model is then evaluated by:

$$y\left(k+1\right)=\frac{{\sum }_{i=1}^r{\mu }_i\left(x\left(k\right)\right)y^i\left(k+1\right)}{{\sum }_{i=1}^r{\mu }_i\left(x\left(k\right)\right)}$$E6

By asking:

$${\mathrm{\Phi }}_j\left(x,c_i,{\sigma }_i\right)=\frac{{\mu }_j\left(x\left(k\right)\right)}{\sum _{i=1}^r{\mu }_i\left(x\left(k\right)\right)}$$E7

with ${\mathrm{\Phi }}_j\left(x,c_i,{\sigma }_i\right)$ is the validation function of the Gaussian function having as parameters the centers $c_i$ and the variances ${\sigma }_i$.

$${\mu }_i\left(x\left(k\right)\right)=\mathit{exp}\left(\frac{-{\left(x_i-c_{i1}\right)}^2}{2{\sigma }_{i1}^2}\right)\dots \mathit{exp}\left(\frac{-{\left(x_i-c_{\mathit{ip}}\right)}^2}{2{\sigma }_{\mathit{ip}}^2}\right)$$E8

The formula (6) becomes

$$y\left(k+1\right)=\sum _{i=1}^r{y}^i\left(k+1\right){\mathrm{\Phi }}_j\left(x,c_i,{\sigma }_i\right)$$E9

A Gaussian function’s center and variance ($\mathit{\text{sigmai}}$) in this situation, as well as the linear parameters of the consequents, are determined in the first stage of the identification of MISO systems, which is often done offline. The second stage, which is online, applies the recursive least squares technique to update the local models’ parameters.

2.2 Offline fuzzy model identification

The data set, denoted $Z$, is constructed by concatenating the regression matrix $X$ and the regressing vector $Y$:

$$X=\left[\begin{array}{c}\dots \\ x\left(k\right)\\ \dots \\ x\left(N-1\right)\end{array}\right],Y=\left[\begin{array}{c}\dots \\ y\left(k\right)\\ \dots \\ y\left(N-1\right)\end{array}\right],Z^T=\left[XY\right]$$E10

where $N$ is the number of observations.

The data set $Z$ will be divided into $\mathit{Nc}$ subfuzzy sets by using the fuzzy classification. Several algorithms, including the C-means, Gatha-Geva, and Gustafson Kessel (GK) algorithms, execute this process. We shall employ these algorithms in the following [9].

Data group membership values will be described by a fuzzy partition matrix $U={\left[{\mu }_{\mathit{ik}}\right]}_{N_c\times N}$ with ${\mu }_{\mathit{ik}}\in \left[01\right]$ represents the degree to which observation $x_k$ belongs to group $i$. Each group is characterized by a center $c_i$, where $C=\left[c_1,\dots ,C_{N_c}\right]$ is the vector of the centers. And a covariance matrix $F=\left[F_1,\dots ,F_{\mathit{Nc}}\right]$ describes the variance of the data in the group $F_i$.

The Gaussian-type membership functions chosen in the context of this chapter are given by:

$${\mathrm{\Omega }}_{\mathit{ij}}\left(x_j\left(k\right)\right)=\mathit{exp}\left(\frac{-1}{2}\frac{\left(x_j\right)-c_{\mathit{ij}}^2}{{\sigma }_{\mathit{ij}}^2}\right)$$E11

The parameters of the consequents ${\theta }_i=\left[A_i,B_i,C_i\right]$ are estimated separately by the recursive least squares algorithm by minimizing the following objective function:

$${\mathit{min}}_{{\theta }_i}\frac{1}{N}{\left[Y-{\mathit{\xi \theta }}_i\right]}^T{Q}_i\left[Y-{\mathit{\xi \theta }}_i\right]$$E12

with $\xi =\left[X1\right]$ is the regression matrix augmented by adding a unit column vector and $Q_i$ is a matrix containing the values of the validity functions ${\mathrm{\Phi }}_i$ of the $i^{\mathit{th}}$ local linear model of each data group:

$$Q=\left[\begin{array}{cccc}{\mathrm{\Phi }}_{-}\left(x\left(1\right),c_i,{\sigma }_i\right)& 0& \dots & 0\\ 0& {\mathrm{\Phi }}_{-}\left(x\left(2\right),c_i,{\sigma }_i\right)& \dots & 0\\ ⋮& ⋮& ⋮& ⋮\\ 0& 0& \dots & {\mathrm{\Phi }}_{-}\left(x\left(N\right),c_i,{\sigma }_i\right)\end{array}\right],$$E13

The least squares estimate of the consequent parameters, (${\theta }_i=\left[a_{i1}\dots a_{\mathit{im}},b_{i1}\dots b_{\mathit{im}},c_i\right]$) is given by:

$${\theta }_i={\left[{\xi }^T{Q}_i\xi \right]}^{-1}{\xi }^T{Q}_{i}Y$$E14

2.3 Online fuzzy model identification

The parameters of any real system’s model change over time in accordance with the conditions of the experiment and thus necessitate their adaptation. The recursive least squares with normalization and projection adaption algorithm is as follows [10]:

$${\varphi }_i\left(k+1\right)=\left[{\mu }_{i}y\left(k-1\right),\dots {\mu }_{i}y\left(k-m\right),{\mu }_{i}u\left(k-1\right),\dots {\mu }_{i}u\left(k-m\right),1\right]$$

$${\theta }_i=\left[a_{i1}\dots a_{\mathit{im}},b_{i1}\dots b_{\mathit{im}},c_i\right]$$

$$\varphi \left(k-1\right)=\left[\varphi {\left(k-1\right)}_1^T{\varphi }_2^T\left(k-1\right)\dots {\varphi }_r^T\left(k-1\right)\right]$$

$$\theta ={\left[{\theta }_1^T{\theta }_2^T\dots {\theta }_r^T\right]}^T$$

Obtaining the estimated parameter vector $\widehat{\theta }$ from the system parameter vector $\theta$ is given by the following least squares algorithm:

$$\widehat{\theta }\left(k\right)=\widehat{\theta }\left(k-1\right)+\frac{p\left(k-1\right)\overline{\varphi }\left(k-1\right){\overline{e}}_1\left(k-1\right)}{\lambda \left(k\right)+{\overline{\varphi }}^T\left(k-1\right)p\left(k-1\right)\overline{\varphi }\left(k-1\right)}$$E15

With:

• ${\overline{e}}_1\left(k\right)=\overline{y}\left(k\right)-{\overline{\varphi }}^T\left(k-1\right)\theta \left(k-1\right)$

• $p\left(k\right)=p\left(k-1\right)-\frac{p\left(k-1\right)\overline{\varphi }\left(k\right)\overline{\varphi }\widehat{T}\left(k\right)p\left(k-1\right)}{\lambda \left(k\right)+{\overline{\varphi }}^T\left(k-1\right)p\left(k-1\right)\overline{\varphi }\left(k-1\right)}$

• $\lambda \left(k\right)={\lambda }_0\lambda \left(k-1\right)+\left(1-{\lambda }_0\right)$, with ${\lambda }_0=0.95,p\left(0\right)=\mathit{\epsilon I},0<\epsilon <\infty$

• $\overline{y}\left(k\right)=\frac{y\left(k\right)}{\eta \left(k-1\right)}$, $\overline{\varphi }\left(k\right)=\frac{\varphi \left(k-1\right)}{\eta \left(k-1\right)}$, $\eta \left(k\right)=\mathit{max}1,\parallel \varphi \left(k\right)\parallel$, $\eta \left(k\right)$ is a normalization signal.

3. PID controller

3.1 Continous PID controller

The structure of the control law of a proportional-integral-derivative (PID) controller results from the sum of three actions:

• Proportional action (P) is as follows:

$$u_p\left(t\right)=K_p\epsilon \left(t\right)$$E16

where $\epsilon \left(t\right)=y_c\left(t\right)-y\left(t\right)$ represents the error signal with $y_c\left(t\right)$ the set point and $y\left(t\right)$ the output of the system to be controlled.

• Integral action (I):

$$u_I\left(t\right)=\frac{K_p}{T_i}{\int }_0^t\epsilon \left(\tau \right)$$E17

where $T_i$ denotes the integration constant.

• Derivative action (D) is as follows:

$$u_D\left(t\right)=K_p{T}_d\frac{\mathit{d\epsilon }\left(t\right)}{\mathit{dt}}$$E18

where $T_d$ denotes the derivative constant.

Hence, the command signal is as follows:

$$u\left(t\right)=u_p\left(t\right)+u_i\left(t\right)+u_d\left(t\right)$$E19

The transfer function of a theoretical PID regulator will be given by:

$$R\left(p\right)=K_p\left(1+\frac{1}{T_{i}p}+t_{d}p\right)$$E20

In practice, it is not possible to achieve a perfect derived action. We often use a filtered version:

$$T_{d}p\to \frac{T_{d}p}{1+\frac{T_d}{N_{f}p}}\mathit{\text{with}}{N}_f\ge 5$$E21

He then just wrote the transfer function of a filtered PID as follows:

$$R\left(p\right)=K_p\left(1+\frac{1}{T_{i}p}+\frac{T_{d}p}{1+\frac{T_d}{N_{f}p}}p\right)$$E22

3.2 Discrete version of the PID controller

To obtain the discrete version, we replace:

• the integral by a discrete sum:

$$u_I\left(t\right)=\frac{K_p}{T_i}{\int }_0^t\epsilon \left(\tau \right)\mathit{d\tau }\to u_I\left(k\right)=\frac{K_p{T}_e}{T_i}\sum _{i=0}^k\epsilon \left(k\right)$$E23

where $\mathit{Te}$ represents the sampling period and k is an integer.

• The derivative by a difference:

$$u_D\left(t\right)=K_p{T}_d\frac{\mathit{d\epsilon }\left(t\right)}{\mathit{dt}}\to u_D\left(k\right)=K_p{T}_d\frac{\epsilon \left(k\right)-\epsilon \left(k-1\right)}{T_e}$$E24

The control law corresponding to the discrete PID regulator also results from the sum of three terms:

$$u_{\mathit{PID}}\left(k\right)=u_p\left(k\right)+u_i\left(k\right)+u_d\left(k\right)=K_p\left(\epsilon \left(k\right)+\frac{T_e}{T_i}\sum _{i=0}^k\epsilon \left(k\right)+\frac{T_d}{T_e}\left(\epsilon \left(k\right)-\epsilon \left(k-1\right)\right)\right)$$E25

This method of discretization is known by the method of upper rectangles. The correspondence between the Laplace transform and the $z$ transform of the integral and derivative terms is given by:

• For the integral term:

$$\frac{1}{T_{i}p}\to \frac{T_e}{T_i}\frac{1}{1-z^{-1}}=\frac{T_e}{T_i}\frac{z}{z-1}$$

• For the derived term:

$$T_{d}p\to \frac{T_d}{T_e}\left(1-z^{-1}\right)=\frac{T_d}{T_e}\frac{z-1}{z}$$

This results in the transfer function of a discrete PID controller:

$$R\left(z^{-1}\right)=K_p\left(1+\frac{T_e}{T_i}\frac{1}{1-z^{-1}}+\frac{N_f{T}_d}{\left(N_f{T}_e+T_d\right)}\frac{1-z^{-1}}{1-\frac{T_d}{N_f{T}_e+T_d}{z}^{-1}}\right)$$

3.3 Different forms of a numerical PID

3.3.1 Parallel form

The transfer function of the digital PID regulator with derivative filtering in parallel form is given by:

$$R_{\mathit{Par}}\left(z^{-1}\right)=K_p+K_i\frac{1}{1-z^{-1}}+K_d\frac{1-z^{-1}}{1-\iota z^{-1}}$$E26

with $K_i=K_p\frac{T_e}{T_i}$, $K_d=K_p\frac{N_f{T}_d}{N_f{T}_e+T_d}$ and $\iota =\frac{T_d}{N_f{T}_e+T_d}$

Note 1: if there is no filter on the derivative, that is to say that $N_f\to \infty$, then we will have

$$\iota \to 0\mathit{\text{and}}{K}_d\to K_p\frac{T_d}{T_e}$$

3.3.2 Mixed form

The transfer function of the digital PID regulator with derivative filtering in mixed form is given by:

$$R_{\mathit{Mix}}\left(z^{-1}\right)=k_p\left(1+k_i\frac{1}{1-z^{-1}}+k_d\frac{1-z^{-1}}{1-\iota z^{-1}}\right)$$E27

with $k_p=K_p$, $k_i=\frac{T_e}{T_i}$ and $d_d=\frac{N_f{T}_d}{N_f{T}_e+T_d}$

3.3.3 Series form

The transfer function of the digital PID regulator with filtering of the derivative in series form is written in the following form:

$$R_{\mathit{Ser}}\left(z^{-1}\right)=\frac{r_0+r_1{z}^{-1}+r_2{z}^{-2}}{\left(1-z^{-1}\right)\left(1+s_1{z}^{-1}\right)}$$E28

with

$$\begin{array}{c}{r}_0=K_p\left(1+\frac{T_e}{T_i}-N_f{s}_1\right)\\ r_1=K_p\left(s_1\left(1+\frac{T_e}{T_i}+2N_f\right)-1\right)\\ r_2=-K_p{s}_1\left(1+N_f\right)\\ s_1=\frac{-T_d}{N_f{T}_e+\mathit{td}}\end{array}$$

Note 2: the continuous equivalent does not always exist!

The existing condition requires that:

$$-1\le s_1\le 1\mathit{\text{thatistosay}}\frac{T_d}{N_f}>0$$

Note 3: In this chapter, we will use the transfer function of the digital PID regulator with filtering of the derivative in series form.

3.3.4 RST structure

The standard form of RST controller is presented as follows (Figure 1):

Figure 1.

Standard form of RST.

In this chapter, we consider only two branches $R$ and $S$, as shown in the diagram below, that is to say:

$$T\left(z^{-1}\right)=R\left(z^{-1}\right)$$

The transfer function in a closed loop is given by:

$$H_{\mathit{Bf}}\left(z^{-1}\right)=\frac{B\left(z^{-1}\right)}{A\left(z^{-1}\right)S\left(z^{-1}\right)+B\left(z^{-1}\right)R\left(z^{-1}\right)}$$E29

again:

$$H_{\mathit{Bf}}\left(z^{-1}\right)=\frac{B\left(z^{-1}\right)}{P{\text{'}}^{\prime }\left(z^{-1}\right)}$$E30

with:

$$\begin{array}{c}e\left(k\right)=r\left(k\right)-y\left(k\right)\\ B\left(z^{-1}\right)=b_1{z}^{-1}+b_2{z}^{-2}\\ A\left(z^{-1}\right)=1+a_1{z}^{-1}+a_2{z}^{-2}\\ R\left(z^{-1}\right)=r_0+r_1{z}^{-1}+r_2{z}^{-2}\\ S\left(z^{-1}\right)=\left(1-z^{-1}\right)\left(1+s_1{z}^{-1}\right)\end{array}$$

From the after scheme (Figure 2), the equation of regulator is written by:

Figure 2.

RST structure.

$$S\left(z^{-1}\right)u\left(z^{-1}\right)=R\left(z^{-1}\right)e\left(z^{-1}\right)$$E31

Hence, the control law of the system can be written in this form:

$$u\left(k\right)=\left(1-s_1\right)u\left(k-1\right)+s_{1}u\left(k-2\right)+r_{0}e\left(k\right)+r_{1}e\left(k-1\right)+r_{2}e\left(k-2\right)$$E32

4. Method of determining the parameters of PID

How to proceed to the synthesis of a digital PID controller under form RST is as follows :

• Performance specifications;

• A method of calculation of the parameters $r_0$, $r_1$, $r_2$, and $s_1$;

• Obtaining parameters of PID controller.

4.1 Pole placement

With the help of this technique, it is feasible to swiftly acquire results that are satisfactory for the temporal criteria related to the step response of the corrected system. It involves using a PID controller to create a system whose corrected $\mathit{Hbf}$ possesses the following features:

• A damping coefficient $m$ for the dominant mode;

• A maximum steady-state error (for the step response);

• An overshoot $D_1$ and a peak time $T_p$.

For the choice of the sampling period $T_e$, we consider that the studied systems have a behavior close to that of a second order so:

$$0,25<w_0{T}_e<1,5$$

with $w_0$ represents the undamped pulsation. For the choice of pulsation and damping coefficient, it is important to satisfy the following conditions :

$$0,25<w_0{T}_e<1,5$$

$$0,7<m<1$$

Because, the two settings $w_0$ and $m$ directly affect the behavior of the system, rise time $t_m$, settling time $t_s$, and the maximal overshoot $D_1$. The sampling period is smaller than $m$ is small.

Given the desired characteristic polynomial system:

$$P\left(z^{-1}\right)=1+p_1{z}^{-1}+p_2{z}^{-2}$$E33

where $p_1$ and $p_2$ depend on the imposed specifications (rise time, overshoot, and damping), where:

$$\begin{array}{c}{p}_1=-2e^{-{\mathit{mw}}_0{T}_e}\mathit{cos}\left(w_0{T}_e\mathit{\text{sqrt}}\left(1-m^2\right)\right)\\ p_2=e^{-{\mathit{mw}}_0{T}_e}\end{array}$$

where equation Diophantine is written as:

$$P\left(z^{-1}\right)=A\left(z^{-1}\right)B\left(z^{-1}\right)+B\left(z^{-1}\right)R\left(z^{-1}\right)$$E34

such as $A\left(z^{-1}\right)$, $B\left(z^{-1}\right)$, and $P\left(z^{-1}\right)$ are the known polynomials, and $S\left(z^{-1}\right)$ and $R\left(z^{-1}\right)$ are the wanted polynomials. We recall that:

$$\frac{B\left(z^{-1}\right)}{A\left(z^{-1}\right)}=\frac{b-1z^{-1}+b_2{z}^{-2})}{1+a_1{z}^{-1}+a_2{z}^{-2}}$$E35

whence (34) becomes

$$\begin{array}{c}P\left(z^{-1}\right)=\left(1+a_1{z}^{-1}+a_2{z}^{-2}\right)\left(1-z^{-1}\right)\left(1+s_1{z}^{-1}\right)\\ +(b-1z^{-1}+b_2{z}^{-2}))\left(r_0+r_1{z}^{-1}+r_2{z}^{-2}\right)\end{array}$$E36

By identification, we can write

$$\left\{\begin{array}{l}{a}_1-1+s_1+r_0{b}_1=p_1\\ \left(a_1-1\right)s_1-a_1+a_2+r_0{b}_2+r_1{b}_1=p_2\\ \left(a_2-a_1\right)s_1-a_2+r_1{b}_2+r_2{b}_1=0\\ -a_2{s}_1+r_2{b}_2=0\end{array}\right.$$E37

where $p_1$ and $p_2$ are the coefficients of the characteristic polynomial $P\left(z^{-1}\right)$:

$$\left\{\begin{array}{l}+s_1+r_0{b}_1=p_1-a_1+1\\ \left(a_1-1\right)s_1+r_0{b}_2+r_1{b}_1=p_2-a_2+a_1\\ \left(a_2-a_1\right)s_1+r_1{b}_2+r_2{b}_1=a_2\\ -a_2{s}_1+r_2{b}_2=0\end{array}\right.$$E38

We can write equations in a matrix form, we then have $\mathit{MQ}=P$ with:

$$M=\left[\begin{array}{cccc}1& b_1& 0& 0\\ a_1-1& b_2& b_1& 0\\ a_2-a_1& 0& b_2& b_1\\ -a_2& 0& 0& b_2\end{array}\right]Q=\left[\begin{array}{c}{s}_1\\ r_0\\ r_1\\ r_2\end{array}\right]P=\left[\begin{array}{c}{p}_1-a_1+1\\ p_2-a_2+a_1\\ a_2\\ 0\end{array}\right],$$E39

The matrix $M$ is called Sylvester matrix. To arrive at the values of the parameters RST controller, first the inversibility of the matrix $M$ must be ensured. If this last holds, then:

$$Q=M^{-1}P$$E40

The structure of the corrector is the standard structure discretized by approximating above rectangles, and to find the values of PID parameters we will take the following expression:

$$\begin{array}{c}{K}_p=\frac{\left(r_0{s}_1-r_1-\left(2+s_1\right)r_2\right)}{{\left(1+s_1\right)}^2}\\ T_i=T_e\frac{k_p\left(1+s_1\right)}{r_0+r_1+r_2}\\ T_d=T_e\frac{r_0{s}_1^2-r_1{s}_1+r_2}{K_p{\left(1+s_1\right)}^3}\\ \frac{T_d}{N_f}=\frac{-s_1{T}_e}{1+s_1}\end{array}$$E41

4.2 PID controller based on PSOCS

4.2.1 PSO

This approach was introduced in 1995 [11]. The swarm $X_p,p=1,2,\dots ,N_p$ which has $N_p$ particles is considered in the standard PSO. In D-dimensional space, the $p^{\mathit{th}}$ particle $X_p=x_{p1},\dots ,x_{\mathit{pD}}$ is a potential solution to the researched problem. The velocity of the $p^{\mathit{th}}$ particle is $V_p=v_{p1},\dots ,v_{\mathit{pD}}$. For the $p^{\mathit{th}}$ particle, the best position in the $p^{\mathit{th}}$ step is expressed as ${\mathit{\text{pbest}}}_p\left(t\right)={\mathit{\text{pbest}}}_{p1}\left(t\right),\dots ,{\mathit{\text{pbest}}}_{\mathit{pD}}\left(t\right)$. Meanwhile, the best position of the entire swarm is defined as $\mathit{\text{Gbest}}\left(t\right)={\mathit{\text{Gbest}}}_1\left(t\right),\dots ,{\mathit{\text{Gbest}}}_D\left(t\right)$. Thus, at the ${\left(t+1\right)}^{\mathit{th}}$ step, the new position of the $p^{\mathit{th}}$ particle is calculated as follows:

$$x_{\mathit{pj}}\left(t+1\right)=x_{\mathit{pj}}\left(t\right)+v_{\mathit{pj}}\left(t+1\right)$$E42

$$v_{\mathit{pj}}\left(t+1\right)=\omega v_{\mathit{pj}}\left(t\right)+c_1\times {\mathit{rand}}_1\times \left[{\mathit{\text{Pbest}}}_{\mathit{pj}}\left(t\right)-x_{\mathit{pj}}\left(t\right)\right]+c_2\times {\mathit{rand}}_2\times \left[{\mathit{\text{Gbest}}}_p\left(t\right)-x_{\mathit{pj}}\left(t\right)\right]$$E43

where $x_{\mathit{pj}}\left(t\right)$ and $v_{\mathit{pj}}\left(t\right)$ represent the position and velocity of the $p^{\mathit{th}}$ particle with respect to the $j^{\mathit{th}}$ dimension, $\omega$ is the inertia weight factor, $c_1$ and $c_2$ are two positive constants called acceleration coefficients, and ${\mathit{rand}}_1$ and ${\mathit{rand}}_2$ are two uniformly distributed random values in the range $\left[0,1\right]$.

4.2.2 Cuckoo search and Lévy flights

4.2.2.1 Key step of Cuckoo search

A type of metaheuristic algorithm inspired by nature, CS, is motivated by the aggressive reproduction tactics of particular cuckoo species. The last rule involves adding some additional random solutions to the algorithm [12]. Three idealised rules are delineated. The friction $p_a$ of the $n_p$ host nests to create new nests can be used to estimate it.

The fundamental processes of the CS can be ascertained by studying the cuckoo breeding behaviour, which is described in [13] and was condensed in [14]. The objective function $f\left(x\right)$, $x=x_1,\dots ,x_D$ in the D-dimensional space is how the optimization problem is defined. We normalised the CS variables with PSO in order to elucidate on the PSOCS algorithm’s constructional aspects. The given search space contains $N_p$ host nests $x_p,p=1,\dots ,N_p$. Each nest $x_p$ represents a potential answer to the optimization problem that needs to be solved. The representations are the same as the particle $x_p$ in PSO. Finding the new population of nests $x_p\left(t+1\right)$ is one of the CS’s crucial phases. Furthermore, Lévy flying is used to obtain the new nests:

$$x_{\mathit{pj}}\left(t+1\right)=x_{\mathit{pj}}\left(t\right)+\alpha \oplus \mathit{\text{Levy}}\left(\lambda \right)$$E44

where $\alpha$ is the step size, $\oplus$ represents the entry-wise multiplication operation, and $\lambda$ is a Lévy flight parameter.

4.2.2.2 Lévy flights

A type of random walk with a step length that has the Lévy distribution is called an Lévy flight. It is added to the CS algorithm to obtain an intermittent scale-free search pattern akin to Lévy’s flight. It is shown in [15] that Lévy flights can maximize the effectiveness of resource searches in the ambiguous setting. A definition of the Lévy distribution is as follows:

$$L\left(s,\gamma ,\nu \right)=\left\{\begin{array}{l}\sqrt{\frac{\gamma }{2\mathrm{\Pi }}}\mathit{exp}\left[-\frac{\gamma }{2\left(s-\nu \right)}\frac{1}{{\left(s-\nu \right)}^{3/2}}\right],0<\nu <s<1\\ 0\mathit{\text{Otherwise}}\end{array}\right..$$E45

where $\nu >0$ is a minimum step and $\gamma$ is a scale parameter. In Mantegna’s algorithm, the step length $s$ can be calculated by:

$$s=\frac{u}{{\left|v\right|}^{1/\beta }}$$E46

where $u$ and $V$ are drawn from normal distribution.That is,

$$\left(\begin{array}{c}u\sim N\left(0,{\sigma }_u^2\right),\\ v\sim N\left(0,{\sigma }_v^2\right),\\ {\sigma }_u={\frac{\tau \left(1+\beta \right)\left(\mathit{sin}\left(\mathit{\pi \beta }/2\right)\right)}{\tau \left[\left(1+\beta \right)/2\right]\beta 2^{\left(1+\beta \right)/2}}}^{1/\beta }\\ {\sigma }_v=1\end{array}\right..$$E47

where $\tau \left(z\right)$ is the Gamma function $\tau \left(z\right)={\int }_0^{\infty }{t}^{z-1}{e}^{-t}\mathit{dt}$. In the case when $z=n_p$ is an integer, $\tau \left(n_p\right)=\left(n_p-1\right)!$

4.2.3 Particle Swarm Optimization-Cuckoo search algorithm with Lévy flights

Although the entry-wise product $\oplus$ employed in PSO is comparable, the random walk via Lévy flight is more effective at navigating the search space since its step length is significantly greater over time. The Lévy flight can exhibit self-similarity and fractal behavior in flight patterns because a power-law distribution is frequently associated with some scale-free properties [16]. Naturally, the Lévy flight is thought to replace the random searching approach of the conventional PSO algorithm in order to improve PSO’s performance. Thus, PSOCS is the name given to the modified PSO algorithm.

The searching ability of PSO is influenced by random variables ${\mathit{rand}}_1$ and ${\mathit{rand}}_2$, when we set the parameters $\omega$, $c_1$, and $c_2$ as fixed values [17]. We introduced the Lévy flight for the change of random step length. Thus the formed PSO-CS algorithm is detailed as follows:

• Step 1. The priori values of parameters are initialized, such as population size of swarm $N_p$, minimum and maximum weights (${\omega }_{\mathit{min}}$, ${\omega }_{\mathit{max}}$), and acceleration coefficients ($c_1$, $c_2$).

• Step 2. The lower and upper bounds for each particle and the particles’ velocities are specified in different neighborhood.

• Step 3. The first generation of particles is randomly initialized within the specified space, $X_1=x_{11},\dots ,x_{1D}$.

while$t<\mathit{\text{Maxgeneratio}}$ or other stop criterion.

• Step 4. The fitness function of each particle $f\left(X_p\right)=f_p$ is evaluated. The best position $\mathit{\text{Pbest}}\left(p\right)\left(t\right)$ and the best position of the whole swarm $\mathit{\text{Gbest}}\left(t\right)$ are found (similar to PSO):

  $${\mathit{\text{Pbest}}}_p\left(t+1\right)=\left(\begin{array}{c}{x}_p\left(t+1\right)\mathit{\text{if}}f\left(x_p\left(t+1\right)\right)<f\left(x_p\left(t\right)\right)\\ {\mathit{\text{Pbest}}}_{p}2\mathit{\text{cmOtherwise}}\end{array}\right.$$E48
  $$\mathit{\text{Gbest}}\left(t\right)={\mathit{\text{minPbest}}}_1\left(t\right),{\mathit{\text{Pbest}}}_2\left(t\right),\dots ,{\mathit{\text{Pbest}}}_{N_p}\left(t\right)$$

• Step 5. A friction $p_a$ of the worst performing particle is chosen in terms of the fitness function. The selected particles should be abandoned, and then the replacement of randomly generated ones is undertaken within the specified search space (similar to CS).

• Step 6. The inertia weight $\omega$ is updated as:

$$\omega ={\omega }_{\mathit{max}}-\frac{{\omega }_{\mathit{max}}-{\omega }_{\mathit{min}}}{T}\times t$$E49

The velocity $v_p$ and position $X_p$ of each particle are updated according to Eqs. (42) and (43). Parameters ${\mathit{rand}}_1$ and ${\mathit{rand}}_2$vary with Lévy flight pattern following (45) to (47),which is different from the former PSO (hybrid of PSO and CS):

$$\begin{array}{c}{v}_{\mathit{pj}}\left(t+1\right)=\omega v_{\mathit{pj}}\left(t\right)+\left(c_1\oplus \mathit{\text{Levy}}\left(\lambda \right)\right)\times \left[{\mathit{\text{Pbest}}}_{\mathit{pj}}\left(t\right)-x_{\mathit{pj}}\left(t\right)\right]\\ +\left(c_2\oplus \mathit{\text{Levy}}\left(\lambda \right)\right)\times \left[{\mathit{\text{Gbest}}}_j\left(t\right)-x_{\mathit{pj}}\left(t\right)\right]\end{array}$$E50

$$x_{\mathit{pj}}\left(t+1\right)=x_{\mathit{pj}}\left(t\right)+v_{\mathit{pj}}\left(t+1\right)$$E51

• Step 7. The iteration step increases ($t=t+1$). As with other metaheuristic algorithms, the terminating generation or the fitness function’s maximum value serves as the termination criterion. If the termination criterion is not met, go to Step 4. Else return $X_{\mathit{\text{best}}}$ as the final solution to the optimization problem. End while

• Step 8. List the optimization results. Meanwhile plot each value of the best fitness function in the optimization processing.

The problem’s dimension determines the size of the swarm $N_p$. Similar to the conventional PSO and CS algorithms, the hybrid algorithm’s parameters depend on the task at hand [17]. They are frequently established following study and improved through repeated calculations. A Monte Carlo method, which is employed in the research of classical particle swarm optimization [18], can be incorporated to acquire a set for fixed parameters for the proposed method’s performance analysis and parameter selection.

Besides, $\tau \left(z\right)={\int }_0^{\infty }{t}^{z-1}{e}^{-t}\mathit{dt}$ is realized by Gamma function $\tau =\gamma \left(z\right)$ using MATLAB. The inertia weight changes from ${\omega }_{\mathit{min}}$ to ${\omega }_{\mathit{max}}$.

The friction $p_a$ controls the launching of a random long-distance exploration strategy, reflecting the probability whether the nest will be abandoned or be updated [19]. The higher the value of this parameter, the closer the search process is to the random search.

5. PID-PSOCS algorithm

The PID-PSOCS algorithm is an algorithm for determining the parameters of the PID regulator which satisfies the stability conditions and is capable of guaranteeing good performance (overshoot, rise time, stabilization time, and static error). This is due to the minimization of an objective function. This function depends on all system performance indexes (static error quality, allowed overshoot, rise time, and settling time). This objective function to be minimized used in the PID-PSOCS algorithm is described by:

$$\begin{array}{c}J=\left(1-e^{\left(-\beta \right)}\right)\left(D+E_{\mathit{ss}}\right)+e^{-\beta }\left(t_s-t_m\right)\\ \mathit{minimize}J{\left(G\right)}_{G=r_0,r_1,r_2,s_1}\\ G_{\mathit{min}}<G<G_{\mathit{max}}\end{array}$$E52

Within the initial search space ($R_{\mathit{min}},R_{\mathit{max}}$), the controller parameters ($r_0,r_1,r_2,s_1$) will be optimized using the optimization algorithms. The graphic below shows the regulator block diagram for a nonlinear system using a PID controller and the PSOCS optimization algorithm (Figure 3):

Figure 3.

Block diagram of a PID-PSOCS algorithm.

The PID-PSOCS algorithm is summarized by the following steps

Algorithm 1 : Proposed AFMPC-PSOCS

• Phase 1. Construction of the fuzzy model using the GK algorithm.

• Phase 2. Optimization of the PID regulator parameters

 1. Step 1. Set parameters $\mathit{Np}$, $c_1$, $c_2$. Initialize the position and velocity of each particle.

  Repeat for h1=1,2,...

 2. Step 2. For $t=1\to t_{\mathit{max}}$ do

  For each particle do:

 3. Step 3. Calculate the control law $u_i\left(k\right)$.

 4. Step 4. Determine the parameters $r_0$, $r_1$, $r_2$ and $s_1$

 5. Step 5. Calculate the value of the fitness function from (52).

 6. Step 6. Update the value of $\mathit{\text{Pbest}}$ and $\mathit{\text{Gbest}}$.

 7. Step 7. Update the position and velocity of each particle from Eqs. (??).

  End for

  End for

 8. Until, the stability conditions are satisfied, then stop. Otherwise, $h1=h1+1$ and return to step 2.

 9. Step 8. Find the vector of optimum parameters $G_b\mathit{est}=\left[r_0{r}_1{r}_2{s}_1\right]$

6. Application of the proposed PID method

This section presents an example of simulation on a Continuous Stirred Tank Reactor (CSTR) nonlinear multivariate system to illustrate the efficiency of the proposed algorithm. We take as an example the system proposed by [20]:

$$\left\{\begin{array}{l}{C}_a\left(k+1\right)=C_a\left(k\right)+T_e\left(\frac{1}{\upsilon }q\left(k\right)\right)\left(C_{\mathit{Af}}-C_a\left(k\right)-k_0{C}_a\left(k\right)e^{\frac{-E}{\mathit{RT}\left(k\right)}}\right)\\ T\left(k+1\right)=\begin{array}{l}T\left(k\right)+T_e\left(\frac{1}{\upsilon }q\left(k\right)\right)\left(T_f-T\left(k\right)\right)+k_1{C}_a\left(k\right)e^{\frac{-E}{\mathit{RT}\left(k\right)}}\\ +k_2{q}_c\left(k\right)\left(1-e^{-\left(\frac{k_3}{q_c\left(k\right)}\right)}\left(T_{\mathit{cf}}-T\left(k\right)\right)\right)\end{array}\end{array}\right.$$E53

with $k_1=-\frac{\mathrm{\Delta }{\mathit{Hk}}_0}{\rho C_p}$, $k_2=\frac{{\rho }_c{C}_{\mathit{pc}}}{\rho C_p\upsilon }$ and $k_3=\frac{h_a}{{\rho }_c{C}_{\mathit{pc}}}$ The output variables are the measured product concentration Ca and the reactor temperature $T$. The input variables are the system flow rate q and the coolant flow rate $q_c$. The nominal system parameters are listed in Table 1. We applied the PID-PSOCS algorithm to optimize the PID controller parameters for the $C_a$ and $T$ outputs.

Measured product concentration	Ca	0,1 $\mathit{mol}{L}^{-1}$
Reactor temperature	T	438,51 $K$
Coolant flow	qc	103,41 ${\mathit{min}}^{-1}$
Process flow	q	100 L ${\mathit{min}}^{-1}$
Focus feeding	CAf	1 $\mathit{mol}$$L^{-1}$
temperature supply	Tf	350 $K$
Coolant temperature	Tcf	350 $K$
Volume of CSTRs	v	100 $L$
heat transfer	q	$7\times {10}^5$cal ${\mathit{min}}^{-1}K$
Constant rate of reaction	CAf	02 $\times {10}^{10}$${\mathit{min}}^{-1}$
Activation energy	Tf	${10}^4$$K$
Heat of reaction	Tcf	$-2\times {10}^5$$\mathit{cal}{\mathit{mol}}^{-1}$
Liquid densities	${\rho }_c$,$\rho$	1 $\times {10}^{3}g{L}^{-1}$
Specific heats	ha	1 $\mathit{cal}$$g^{-1}$$K^{-1}$

Table 1.

CSTR system parameters.

The PID-PSOCS algorithm’s performance was compared to that of the PID-PSO algorithm and the PID-CS algorithm. In fact, we collected 1000 points to use in the GK algorithm’s identification of the model, which uses input variables like $\left[\mathit{Ca}\left(k-1\right)\mathit{Ca}\left(k-2\right)q_c\left(k-1\right)q_c\left(k-2\right)q\left(k-1\right)\right]$ and $\left[T\left(k-1\right)T\left(k-2\right)q\left(k-1\right)q\left(k-2\right)q_c\left(k-1\right)\right]$ and four fuzzy rules.

Under the same circumstances, a comparison investigation was conducted. The three algorithms PID-PSOCS, PID-CS, and PID-PSO have the following adjustment parameters: $\mathit{\text{tmax}}=150$, $\mathit{Np}=10$, $w_{\mathit{max}}=0.9$, $w_{\mathit{min}}=0.4$, $c_1=1.5$, and $c_2=2.5$. The three algorithms’ optimal parameters for the two outputs ($C_a$ and $T$) are provided in Tables 2 and 3, respectively.

Algorithm	kp	Ti	Td
PID-PSO	1.345	0.024	0.0017
PID-CS	0.856	0.013	0.0013
PID-PSOCS	2.1405	0.081	0.0011

Table 2.

Different values of PID parameters for $C_a$ output.

Algorithm	kp	Ti	Td
PID-PSO	1.235	0.023	0.0035
PID-CS	0.746	0.012	0.0045
PID-PSOCS	2.1345	0.072	0.0029

Table 3.

Different values of PID parameters for $T$ output.

Tables 4 and 5 provide a full breakdown of the performance indices for the two algorithms in terms of rise time $\mathit{tm}$, stabilization time $\mathit{ts}$, and overshoot D($\%$). The PID-PSOCS algorithm outperformed the PIDPSO and PID-CS algorithms, according to the Tables 4 and 5.

Algorithm	PID-PSO	PID-CS	PID-PSOCS
tm (s)	8	3	3
ts (s)	9,1	6	5
D%	30%	14%	00%

Table 4.

Performance of different methods for $C_a$ output.

Algorithm	PID-PSO	PID-CS	PID-PSOCS
tm (s)	8	6	7
ts(s)	10	8	9
D%	27.5%	10,9%	00%

Table 5.

Performance of different methods for $T$ output.

Figures 4 and 5 shows the evolution of the two system outputs $\mathit{Ca}$ and $T$ at their reference setpoints. The PID-PSOCS algorithm superimposed the other two algorithms in terms of the other performance indices namely $\mathit{tm}$, $\mathit{ts}$, and $D\left(\%\right)$.

Figure 4.

Concentration evolution.

Figure 5.

Temperature evolution.

7. Conclusions

The proposed solution addresses the PID control of MIMO nonlinear systems issue. Each MISO nonlinear system is given a fuzzy PID control. An approximation of the alleged unknown nonlinearity, coupling between inputs, unmodeled dynamics, and disturbances is provided by the TS fuzzy model. The latter makes it possible to convert a nonlinear problem to a linear one. This approximation allows the decoupling to be released and the nonlinearity problem to be solved via the control synthesis. A numerical example shows the robustness and adaptability of the proposed method to a large class of MISO nonlinear systems.