Preface

We are pleased to present this jointly edited book combining our respective expertise in academics and industry. Meeting the continuously evolving challenges in industrial process control requires efficient control strategies that can yield improved performance. Hence, this book covers a range of topics within the realm of proportional-integral-derivative (PID) controller design for various industrial processes.

Many researchers have reported performance improvement in PID-based industrial process control through the use of complex control strategies that require a large number of controllers and filter parameters. However, as simple and effective PID control schemes are more feasible in practice, such design approaches (for unstable processes, nonlinear systems and multi-input multi-output systems) are very much emphasized in Chapters 1-4 and 6 of this book.

An effective modelling technique is the basis of any model-based controller design approach, and in Chapter 2, a relay feedback-based modelling of a tank system is discussed.

Some PID controller designs based on linearized plant models have their own limitations when employed to control nonlinear systems in practical scenarios. Therefore, a PID controller design method for nonlinear processes is presented in Chapter 5. Advanced disturbance rejection control strategies are vital these days and are much more important than servo tracking in practice. Hence, Chapter 7 is dedicated to a linear quadratic regulator problem with emphasis on a ball and beam system.

We believe that the contents of this book will help engineers and researchers working in the domain of PID controller design for various applications. We thank all the authors of individual chapters for their contributions, and Sara Debeuc (Author Service Manager, IntechOpen), who was constantly in touch with us to ensure the smooth progress of the chapter review and the editing process. We also acknowledge our respective organizations (Billington Process Technology

and National Institute of Technology Patna) for providing us with the necessary facilities and infrastructure during this editorial task. Finally, it would have been impossible to successfully edit this book without the support and cooperation of our families.

### **Dr. Mohammad Shamsuzzoha**

Principal Process Engineer, Billington Process Technology, Sandvika, Norway

#### **Dr. G. Lloyds Raja**

**Chapter 1**

**1. Introduction**

Introductory Chapter: PID-Based

A PID controller is an instrument used in industrial control applications at the regulatory level to regulate process variables e.g., temperature, pressure, flow, etc. To meet the continuously evolving challenges in industrial process control, it is essential to formulate control strategies which can yield improved performance. Proportionalintegral-derivative (PID) controllers are still very much preferred in industries due to their simplicity and ability to yield reasonable closed-loop performance. A recent study has concluded that the preference for the PID, advanced and model predictive control in industries fall in the ratio 100:10:1 [1]. Another study states that about 90%

Majority of the control schemes use PID controllers in a unity feedback configuration [3]. However, unity feedback schemes are not suitable for plants having large time delays (LTD) and disturbances [4]. Hence, attempts have been made to design double-degree-of-freedom (DDOF) control schemes by adding

If an intermediate process output is available, cascaded control (CC) is more capable of giving better closed-loop performance compared to the DDOF control structures mentioned above. Based on the mode of operation, there are two varieties

Plant like boilers and reactors are often modeled as unstable processes having time delay [4]. In contrast, paper drum drier cans and boiler steam drums are of integrating type [4]. Having poles in the right half of the s-plane and origin makes unstable and integrating (UI) plants difficult to control. To control UI processes, modifications are required in single-loop, DDOF, CC and SP bases strategies [8, 13]. Hence, a lot of

of CC strategies: serial and parallel [13, 14]. In practice, time delays occur in transport and composition examination loops [4]. The control schemes reported in [15–17] fails to provide good servo response for processes with LTD. To compensate LTD, Smith predictor (SP) based schemes are reported in the literature [6]. However, SP based control strategies fails to yield satisfactory regulatory performance for processes having LTD in the presence of disturbances [18]. Hence, SP can be combined with cascade control to achieve both satisfactory servo and regulatory

research is still being carried out in the aforementioned domains.

Industrial Process Control

*Mohammad Shamsuzzoha and G. Lloyds Raja*

industrial controllers are of PID type [2] to meet the requirement.

**1.1 Literature review of PID control strategies**

additional controllers [4–12].

performance [18–20].

Assistant Professor, Department of Electrical Engineering, National Institute of Technology Patna, Patna, India

#### **Chapter 1**

## Introductory Chapter: PID-Based Industrial Process Control

*Mohammad Shamsuzzoha and G. Lloyds Raja*

#### **1. Introduction**

A PID controller is an instrument used in industrial control applications at the regulatory level to regulate process variables e.g., temperature, pressure, flow, etc. To meet the continuously evolving challenges in industrial process control, it is essential to formulate control strategies which can yield improved performance. Proportionalintegral-derivative (PID) controllers are still very much preferred in industries due to their simplicity and ability to yield reasonable closed-loop performance. A recent study has concluded that the preference for the PID, advanced and model predictive control in industries fall in the ratio 100:10:1 [1]. Another study states that about 90% industrial controllers are of PID type [2] to meet the requirement.

#### **1.1 Literature review of PID control strategies**

Majority of the control schemes use PID controllers in a unity feedback configuration [3]. However, unity feedback schemes are not suitable for plants having large time delays (LTD) and disturbances [4]. Hence, attempts have been made to design double-degree-of-freedom (DDOF) control schemes by adding additional controllers [4–12].

If an intermediate process output is available, cascaded control (CC) is more capable of giving better closed-loop performance compared to the DDOF control structures mentioned above. Based on the mode of operation, there are two varieties of CC strategies: serial and parallel [13, 14]. In practice, time delays occur in transport and composition examination loops [4]. The control schemes reported in [15–17] fails to provide good servo response for processes with LTD. To compensate LTD, Smith predictor (SP) based schemes are reported in the literature [6]. However, SP based control strategies fails to yield satisfactory regulatory performance for processes having LTD in the presence of disturbances [18]. Hence, SP can be combined with cascade control to achieve both satisfactory servo and regulatory performance [18–20].

Plant like boilers and reactors are often modeled as unstable processes having time delay [4]. In contrast, paper drum drier cans and boiler steam drums are of integrating type [4]. Having poles in the right half of the s-plane and origin makes unstable and integrating (UI) plants difficult to control. To control UI processes, modifications are required in single-loop, DDOF, CC and SP bases strategies [8, 13]. Hence, a lot of research is still being carried out in the aforementioned domains.

#### **1.2 Requirements for industrial process control**

It is essential that a control strategy must be capable of eliminating the load disturbances and tracking the reference input. Moreover, it must be robust towards uncertainties in process dynamics and noise that enters the system. Response of a control system to setpoint changes and disturbances are termed servo and regulatory responses, respectively. In process industries, changes in setpoint happen only when the production rate is altered. Mostly the production rate remains unaltered for years together. On the other hand, closed-loop performance is more frequently hindered by disturbances entering the system. Therefore, disturbance elimination is comparatively more vital than reference following [21]. The essential requirements for a PID control strategy are discussed below.

#### *1.2.1 Disturbance rejection*

The system output deviates from the desired value due to load disturbances which are of low frequency. Hence, rejecting such load disturbances is a primary task of a properly designed controlled system. The instantaneous error *e t*ð Þ is the deviation of setpoint (*r*1) from controlled output (*y*1) at time '*t*'. Using *e t*ð Þ, the performance of a closed-loop control system can be characterized by computing the following measures:

Integrated absolute error (IAE)

$$\text{IAE} = \int\_0^\infty |e(t)|dt\tag{1}$$

Integrated squared error (ISE)

$$\text{ISE} = \int\_0^\infty e(t)^2 dt \tag{2}$$

and Integrated time-weighted absolute error (ITAE)

$$\text{ITAE} = \int\_0^\infty t|e(t)|dt\tag{3}$$

Small values of (1) to (3) indicates better control performance.

#### *1.2.2 Setpoint tracking*

Whenever there is a change in the setpoint (reference input), it is expected that the system output should immediately follow the new reference value. The referencetracking capability of a closed-loop system is characterized by its rise-time (*t*r) and settling-time (*t*s). *t*<sup>r</sup> is the time consumed in system output raising from 10% to 90% of the expected value. Moreover, *ts* is the time consumed in system output to reach up to (and stay within) �2% or �5% of the final value. The system output is expected to have less overshoot, *t*r, *t*<sup>s</sup> and steady state error (error 'e' after reaching steady state) during a change in setpoint. In addition to the above, performance measures like IAE, ISE and ITAE are also used to characterize the servo performance.

#### *1.2.3 System robustness*

The plant model (*G*om) used to design controllers is an approximate version of the actual plant dynamics (*G*o). Therefore, it is important to ensure that the controller designed using *G*om to be robust enough to control *G*o. As per [22], the rule to achieve closed-loop robust stability is

$$\|l\_{\mathbf{m}}(\mathfrak{s})T\_{\mathbf{d}}(\mathfrak{s})\| < \mathbf{1} \forall \mathfrak{o} \in (-\infty,\infty) \tag{4}$$

Here, *T*dð Þ *s* ¼ *jω* denotes complementary sensitivity function. *l*mð Þ *s* ¼ *jω* is the multiplicative uncertainty as given below:

$$|l\_{\mathbf{m}}(\mathfrak{s}) = |\frac{\mathbf{G}\_{\mathbf{o}}(\mathfrak{s}) - \mathbf{G}\_{\mathbf{om}}(\mathfrak{s})}{\mathbf{G}\_{\mathbf{om}}(\mathfrak{s})}|\tag{5}$$

From the magnitude plots of *T*<sup>d</sup> and *l*m, the robust stability of a system is analyzed graphically. Furthermore, system robustness can be measured with maximum sensitivity (*M*s).*Ms* is defined as the inverse of the shortest distance from the Nyquist curve of the loop transfer function to the critical point '�1'. For an unity feedback system having a controller *G*<sup>c</sup> and process model *G*p, *M*<sup>s</sup> is obtained as follows:

$$M\_{\rm s} = \max\_{\omega} \frac{1}{\mathbf{1} + \mathbf{G}\_{\rm c}(j\omega)\mathbf{G}\_{\rm p}(j\omega)}\Big|\tag{6}$$

It is expected for *M*<sup>s</sup> to remain within 1.2 and 2 to ensure a good tradeoff between performance and robustness for stable plants with time delay [4].

#### *1.2.4 Control signal*

Softness of the control action *u*2ð Þ*t* is computed by its total variation (TV). Mathematically, TV is given as

$$\text{TV} = \sum\_{i=1}^{\infty} |u\_2(i+1) - u\_2(i)| \tag{7}$$

Moreover, the maximum magnitude of the control signal is given by *u*2 max ¼max{*u*2ð Þ*t* ∨}. TV and *u*2 max must remain as small as possible in practice.

#### **1.3 Motivation for this book**

The following observations are made from the contemporary works pertaining to PID-based industrial process control:


in practice. Hence, PID controller design for nonlinear processes have attracted good research attention in recent times [24].


Motivated by the above, the subsequent chapters of this book are presented to introduce the readers to some simple PID controller design strategies for unstable processes, nonlinear systems and MIMO systems. Also, considerable attention has been given to familiarize the reader with the concept of ADRC and relay-based auto tuning strategies.

#### **Author details**

Mohammad Shamsuzzoha<sup>1</sup> \* and G. Lloyds Raja<sup>2</sup>

1 Principal Process Engineer, Billington Process Technology, Sandvika, Norway

2 Department of Electrical Engineering, National Institute of Technology Patna, India

\*Address all correspondence to: smzoha@gmail.com

© 2022 The Author(s). Licensee IntechOpen. This chapter is distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

*Introductory Chapter: PID-Based Industrial Process Control DOI: http://dx.doi.org/10.5772/intechopen.109036*

#### **References**

[1] Kano M, Ogawa M. The state of the art in chemical process control in Japan: Good practice and questionnaire survey. Journal of Process Control. 2010;**20**(9): 969-982

[2] O'Dwyer A. Handbook of PI and PID Controller Tuning Rules. London: Imperial College Press; 2009

[3] Shamsuzzoha M, Lee M. IMCPID controller design for improved disturbance rejection of time-delayed processes. Industrial & Engineering Chemistry Research. 2007;**46**(7): 2077-2091

[4] Lloyds Raja G, Ali A. New PI-PD controller design strategy for industrial unstable and integrating processes with dead time and inverse response. Journal of Control, Automation and Electrical Systems. 2021;**32**(2):266-280

[5] Kumari S, Aryan P, Raja GL. Design and simulation of a novel FOIMC-PD/P double-loop control structure for CSTRs and bioreactors. International Journal of Chemical Reactor Engineering. 2021; **19**(12):1287-1303

[6] Kumar D, Aryan P, Raja GL. Design of a novel fractional-order internal model controller-based Smith predictor for integrating processes with large dead-time. Asia-Pacific Journal of Chemical Engineering. 2022;**17**(1):e2724

[7] Mukherjee D, Raja GL, Kundu P, Ghosh A. Improved fractional augmented control strategies for continuously stirred tank reactors. Asian Journal of Control. 2022:1-18. DOI: 10.1002/asjc.2887

[8] Kumar D, Raja GL. Unified fractional indirect IMC-based hybrid dual-loop

strategy for unstable and integrating type CSTRs. International Journal of Chemical Reactor Engineering. 2022. DOI: 10.1515/ijcre-2022-0120

[9] Kumar D, Aryan P, Raja GL. Decoupled double-loop FOIMC-PD control architecture for double integral with dead time processes. The Canadian Journal of Chemical Engineering. 2022;**100**(12):3691-3703. DOI: 10.1002/cjce.24355

[10] Chakraborty S, Singh J, Naskar AK, Ghosh S. A New Analytical Approach for Set-point Weighted 2DOF-PID Controller Design for Integrating Plus Time-Delay Processes: an Experimental Study. IETE Journal of Research. 2022:1-15

[11] Aryan P, Raja GL. A novel equilibrium optimized double-loop control scheme for unstable and integrating chemical processes involving dead time. International Journal of Chemical Reactor Engineering. 2022;**1**:20

[12] Kumari S, Aryan P, Kumar D, Raja GL. Hybrid dual-loop control method for dead-time second-order unstable inverse response plants with a case study on CSTR. International Journal of Chemical Reactor Engineering. 2022;**1**:11

[13] Raja GL, Ali A. Series cascade control: an outline survey. In: 2017 Indian Control Conference (ICC). IEEE; 2017. pp. 409-414

[14] Siddiqui MA, Anwar MN, Laskar SH. Cascade controllers design based on model matching in frequency domain for stable and integrating processes with time delay. COMPEL-The International Journal for Computation and Mathematics in Electrical and Electronic Engineering. 2022;**41**(5):1345-1375. DOI: 10.1108/COMPEL-06-2021-0185

[15] Ranganayakulu R, Seshagiri Rao A, Babu UB, G. Analytical design of fractional IMC filter–PID control strategy for performance enhancement of cascade control systems. International Journal of Systems Science. 2020;**51**(10): 1699-1713

[16] Siddiqui MA, Anwar MN, Laskar SH, Mahboob MR. A unified approach to design controller in cascade control structure for unstable, integrating and stable processes. ISA Transactions. 2021; **114**:331-346

[17] Raja GL, Ali A. Modified parallel cascade control strategy for stable, unstable and integrating processes. Isa Transactions. 2016;**65**: 394-406

[18] Raja GL, Ali A. Enhanced tuning of Smith predictor based series cascaded control structure for integrating processes. ISA Transactions. 2021;**114**: 191-205

[19] Raja GL, Ali A. Smith predictor based parallel cascade control strategy for unstable and integrating processes with large time delay. Journal of Process Control. 2017;**52**:57-65

[20] Vanavil B, Uma S, Rao AS. Smith predictor based parallel cascade control strategy for unstable processes with application to a continuous bioreactor. Chemical Product and Process Modeling. 2012;**7**(1):1-22

[21] Aryan P, Raja G. Equilibriumoptimized IMC-PD double-loop control strategy for industrial processes with dead time. In: Recent Advances in Mechanical Engineering. Singapore: Springer; 2023. pp. 37-50

[22] Morari M, Zafiriou E. Robust Process Control. Upper Saddle River: PTR Prentice Hall; 2002

[23] Siddiqui MA, Anwar MN, Laskar SH. Enhanced control of unstable cascade systems using direct synthesis approach. Chemical Engineering Science. 2021;**232**: 116322

[24] Chang XH, Jin X. Observer-based fuzzy feedback control for nonlinear systems subject to transmission signal quantization. Applied Mathematics and Computation. 2022;**414**:126657

[25] Wang C, Zhang C, Liu L, Ao L, He D. A finite-time adaptive fuzzy backstepping control for multiple-input multiple-output coupled nonlinear systems with tracking error constraints. International Journal of Adaptive Control and Signal Processing. 2022;**36** (11):2823-2837. DOI: 10.1002/acs.3485

[26] Ahmad S, Ali A. On active disturbance rejection control in presence of measurement noise. IEEE Transactions on Industrial Electronics. 2021;**69**(11):11600-11610

[27] Visioli A, Sánchez-Moreno J. A relayfeedback automatic tuning methodology of PIDA controllers for high-order processes. International Journal of Control. 2022. DOI: 10.1080/ 00207179.2022.2135019

#### **Chapter 2**

## Perspective Chapter: Modeling and Identification of Two Tank System Using Relay Feedback – An Experimental Approach

*Devarapalli Kishore and Vaska Lokesh*

### **Abstract**

An experimental approach for modeling and identification of two tank system was considered. For the system that we considered modeling has been carried out with theoretical values. The practically obtained values from series of experiments the used to correlate the simulated values, the process we have considered is the MIMO process which is difficult being the interaction exists, the interactions can be eliminated for better control action using the existing the methods. Simulations with MATLAB has been carried out and the simulated results are used in the real time experimental set using LAB VIEW. There is great degree of the Accuracy between the estimated values and simulated which can be practically demonstrated with the experimental setup. In this relay feedback approach is being considered for conducting the autocuing test to estimate the parameters.

**Keywords:** relay feedback, auto tuning, modeling, identification, MIMO

#### **1. Introduction**

Multivariable process systems are either time-lag dominated or dead-time dominated systems and can be modeled as First Order Plus Dead time (FOPDT) structures. The sequential identification and modified Ziegler-Nichols controller design method form the basic structure for the multivariable auto tuner. These methods consider desired response (output: for example, y1 in case of 2x2 system) and input (u1) for system identification and analysis. Also, they do not discuss about methods to reduce interactions nor do they give any exact analytical expressions for the relay responses that may help in analyzing the interaction behavior between input/output and may provide information regarding closed-loop parameters (PID using model based tuning rules) of the MIMO (Multi input and Multi output) system. Moreover, it is felt that exact model parameters and information on interactions can be better obtained/calculated from mathematical model of undesired relay responses for MIMO systems. As the offdiagonal closed loop transfer functions contain information on interactions, it is better to analyze the control system based on their time domain characteristics.

#### **1.1 Motivation**

Most of the chemical process are multivariable in nature.i.e. more than one input such as distillation column which are used for separation of the chemical components in the Industries. Identification and control of the process plays a crucial role in the Process industries to maintain the certain variable at set point. There are well defined Control strategies are there for identification and control of SISO [2014] (single input and single output process).

The control of MIMO (Multi Input and Multi out) is challenging due to presence of interaction exists among the manipulate variable and control variable. There are lot of methods available for identification and control of the process by reducing the interaction analysis.

Åström and Hägglund [1] were pioneered in the introducing the concept of autotuning to generate the sustained oscillations. Yu [2] has explored the relay feedback to investigate the shapes of various higher order systems. Tau and Liu [3] has presented the research gaps and further explored the further possible investigations. Shankar Prasad and Yaddanapudi Jaya [4] have proposed the new identifications for stable Process. Sujatha and Panda [5] has introduced the concept of estimating the parameters of MIMO systems. Ramesh and Panda [6] has explored the relay feedback and presented the report. Kalapana [7] has extended and developed the analytical method of identifying parameters. Chidambaram and Padma Sri [8] has discussed elaborately discussed the methods for unstable systems.

In this work an attempt is made by taking the two tank system (To control Liquid level) using the relay feedback (sequential tuning approach) and validation of the simulation results with the experimental results. The method used in this research work is sequential tuning approach which is widely acceptable for estimation and identification of the process parameters for proposed two tank (Liquid Level) system. Experimental transfer function is obtained from two tank system a simulation is carried out a comparison is made between experimental one and simulated one.

#### *1.1.1 Two tank interacting system*

The schematic diagram of coupled tank MIMO process is shown in **Figure 1**. The input flow for tank1 and tank2 are Fin1 and Fin2. The controlled variables are level h1 and h2 in the tank1 and tank2. The mass balance and Bernoulli's law yield.

$$A\_1 \frac{dh\_1}{dt} = k\_1 u\_1 - \beta\_1 a\_1 \sqrt{2gH\_1} - \beta\_x a\_{12} \sqrt{2g(H\_1 - H\_2)}\tag{1}$$

$$A\_2 \frac{dh\_2}{dt} = k\_2 u\_2 + \beta\_x a\_{12} \sqrt{2g(H\_1 - H\_2)} - \beta\_2 a\_2 \sqrt{2gH\_2} \tag{2}$$

A1, A2, cross sectional area of tank 1 and tank 2 cm2 ð Þ; a1,a1, cross sectional area of output pipe in tank 1 and tank 2 cm2 ð Þ; *<sup>a</sup>*12, cross sectional area of interaction pipe between tank 1 and tank 2 cm<sup>2</sup> ð Þ; h1 h2, water level of tank 1 and tank 2 (cm), Fin11 Fin22, inflow of tank 1 and tank 2 cm<sup>3</sup> ð Þ *<sup>=</sup>*<sup>s</sup> ; Fout11,Fout22, outflow of tank 1 and tank 2 cm<sup>3</sup> ð Þ *<sup>=</sup>*<sup>s</sup> ; u1,u2, input voltage to pump 1 and pump 2 (V); β1,β2, valve ratio at the outlet of tank 1 and tank 2;

*Perspective Chapter: Modeling and Identification of Two Tank System Using Relay Feedback… DOI: http://dx.doi.org/10.5772/intechopen.106568*

βx, valve ratio of jointed pipe between tank 1 and tank 2; k1,k2, gain of the pump 1 and pump 2 cm3 ð Þ *<sup>=</sup>*<sup>v</sup> � <sup>s</sup> ; g, gravity cm<sup>3</sup> ð Þ *<sup>=</sup>*<sup>s</sup> .

The parameter values of the coupled tank process are given in **Table 1**. The nominal operating conditions of the process are shown in **Table 2**.

The transfer function of coupled tank process is identified using system identification tool box. The levels in the tanks are initially maintained at the operating condition of 24.6cm and 14.4cm by giving the input voltage of 2.5 and 2 volts to the pump1 and pump2 respectively. Then the input to the pump1 is changed from 2.5 to 3 voltages by keeping pump2 input as constant and the level in tank1 and tank2 are recorded. The same procedure is repeated by changing the pump2 input from 2 to 2.5 volts by keeping the pump1 input as constant. The open loop response for the change in input1 and input2 are shown in **Figures 2** and **3**.

The experimentally identified transfer function model is

$$
\begin{bmatrix} h\_1 \\ h\_2 \end{bmatrix} = \begin{bmatrix} \frac{16.99e^{-12.89}}{214.7s+1} & \frac{6.69e^{-72.57s}}{204.93s+1} \\ \frac{9.23e^{-35.01s}}{256.44s+1} & \frac{11.38e^{-25.04s}}{169.15s+1} \end{bmatrix} \begin{bmatrix} u\_1 \\ u\_2 \end{bmatrix} \tag{3}
$$


#### **Table 1.**

*Parameters of coupled tank MIMO process.*


**Table 2.**

*Operating condition of coupled tank MIMO process.*

**Figure 2.** *Open loop response for input change in pump1.*

**Figure 3.** *Open loop responses for input change in pump2.*

#### **1.2 Relay test for the g12 interactive transfer function**

The relay feedback test is conducted on 2X2 MIMO (Two tanks) system. The undesired relay response of gp12 is obtained in the step2 of sequential auto tuning [9]. *Perspective Chapter: Modeling and Identification of Two Tank System Using Relay Feedback… DOI: http://dx.doi.org/10.5772/intechopen.106568*

The relay feedback test output thus obtained is shown in **Figure 4**. It is g12 the interactive transfer function relay response and it has the interaction of g11 transfer function. The time intervals are taken to derive the analytical expression for the above relay response.

#### *1.2.1 Relay test for the* **g21** *interactive transfer function*

The relay feedback test is conducted on 2X2 MIMO (Two tank) system. The undesired relay response of g21 is obtained in the step3 of sequential auto tuning. The relay feedback test output thus obtained is shown in **Figure 5**. It is g21 interactive transfer function relay response and it has the interaction of g22 transfer function. The time intervals are taken to derive the analytical expression for the above relay response.

Under decentralized PI control, with known pairing (y1–u1) and (y2–u2) 2x2 MIMO system is employed. First Relay is placed between y1 and u1, while loop 2 is on manual mode. Following the relay-feedback test, a controller can be designed from the ultimate gain and ultimate frequency. Next performed relay feedback test between y2 and u2, while loop 2 is on automatic mode. A controller is designed from the ultimate gain and ultimate frequency for the loop 2. Third performed relay feedback test between y1 and u1, while loop 2 is on automatic mode. Generally, new set of tuning constants are founded for the controller in loop1. Based on the concept of sequential auto tuning method each controller is designed in sequence. The controller's parameters are converged in 3–4 relay feedback test is shown in the below **Table 3**.

**Figure 4.** *Relay output response for g12 (step2).*

**Figure 5.** *Relay output response for g21.*


**Table 3.**

*Controller parameters computed from biased the relay feedback response of step 2 and step 3 by sequential autotuning.*

#### **1.3 Derivation of analytical expression for 2X2 MIMO system**

The ideal relay output consists of a series of step changes in manipulated variables. Hence, the stabilized output is a sum of infinite terms of step responses due to those step changes. The process output converges to the stationary oscillation in one period and the limit cycle oscillation for ideal relay test is characterized by the deriving analytical expressions. Analytical expressions are mathematical expressions for the stabilized relay feedback output responses and are useful for back calculation of exact process model parameters

### *1.3.1 Analytical expression for the* **g12** *interactive transfer function*

Analytical expression is derived for the undesired response obtained from the relay

*Perspective Chapter: Modeling and Identification of Two Tank System Using Relay Feedback… DOI: http://dx.doi.org/10.5772/intechopen.106568*

Feedback test. The biased relay response is assumed to be formed by n number of small

Step changes. Let <sup>μ</sup><sup>þ</sup> <sup>¼</sup> <sup>μ</sup><sup>0</sup> <sup>þ</sup> <sup>μ</sup> and <sup>μ</sup>‐ <sup>¼</sup> <sup>μ</sup>0–μ. The process input in the relay feedback test consists of a series of step changes with down amplitude, <sup>μ</sup>‐ and up amplitude μþ. At first interval (after synchronizing input with output by time shift), the response can be described.

as:

$$y\_1 = k\_{12} \mu\_0 \left[ 1 - e^{\frac{-t}{\tau\_{12}}} \right] - \frac{k\_{11} k\_{t1}}{\tau\_{i1}} \mu\_- \left[ t - (\tau\_{i1} + \tau\_{11}) e^{\frac{-t}{\tau\_{i1}}} \right] y\_1 (t - 1) \tag{4}$$

$$\begin{split} y\_{2} &= k\_{12} \mu\_{0} \left[ \left( 1 - e^{\frac{t + D\_{1}}{\sigma\_{1}}} \right) - 2 \left( 1 - e^{\frac{\omega}{\sigma\_{1}}} \right) \right] \\ &- \frac{k\_{11} k\_{c1}}{\tau\_{11}} \mu\_{-} \left[ (t + D\_{11}) - \left( \tau\_{11} + \tau\_{11} \left( 1 - e^{\frac{t + D\_{11}}{\sigma\_{11}}} \right) \right) \right] \nu\_{2} (t - 1) \left[ t - \left( \tau\_{11} + \tau\_{11} \left( 1 - e^{\frac{-t}{\sigma\_{11}}} \right) \right) \right] \end{split} \tag{5}$$

Where D21 and D11 are time delays of the individual transfer functions of the system. The above Eq. (5) can be simplified as follows:

$$\begin{aligned} \mu\_{2} &= k\_{12}\mu\_{0} \left\{ [\mathbf{1} - \mathbf{2}] - e^{\frac{-t}{\tau\_{21}}} \left[ e^{\frac{-D\_{t2}}{\tau\_{12}}} - 2 \right] \right\} - \frac{k\_{11}k\_{c1}}{\tau\_{i1}} \times \\ \mu\_{-} &\left\{ t[\mathbf{1} - \mathbf{2}] + \left( (\tau\_{i1} - \tau\_{11}) \left( \mathbf{1} - e^{\frac{-t}{\tau\_{i1}}} \right) \left[ e^{\frac{-D\_{t1}}{\tau\_{i1}}} \right] - 2 \right) \right\} \nu\_{2}(t - \mathbf{1}) \end{aligned} \tag{6}$$

Let p<sup>þ</sup> and p� be the positive and negative half cycle's periods. The third interval lags by an amount *<sup>D</sup>*<sup>21</sup>□*pu=*2 and *<sup>D</sup>*<sup>22</sup>□*pu=*2 from input can be given a

$$-\frac{k\_{11}k\_{c1}}{\tau\_{n1}}\mu\_{-}\left\{ [1-2+2]+\left( (r\_n+\tau\_{11})\left(1-e^{\frac{-t}{\tau\_{11}}}\right) \left[e^{\frac{-D\_{11}+\frac{D\_{11}}{2}}{\tau\_{11}}}-2e^{\frac{D\_{11}}{2}}+2\right] \right) \right\} \chi\_{3}(t-1)\tag{7}$$

The Eq. (7) can easily be simplified (8)

The Eq. (8) for *y*<sup>3</sup> slowly forms a series, as time tends to infinity the response becomes stabilized and it can be described as in given in Eq. (9)

The RHS of Eq. (9) has 4 parts:

$$y\_t = k\_{12}\mu\_0 \left\{ \left[ 1 - 2 + 2 - \cdots \right] \quad -e^{\frac{t}{T\_{12}}} \left[ e^{\frac{D\_{12} + \sum\_{k=1}^{n-1} \frac{t\_k}{2}}} - 2e^{\frac{\sum\_{k=1}^{n-1} \frac{t\_k}{2}}} + \cdots + 2e^{\frac{t\_1}{T\_{12}}} - 2 \right] \right\} \tag{8}$$

$$\frac{-k\_{11}k\_{c1}}{\tau\_{11}}\mu - x\left[ (\mathbf{1} - \mathbf{2} + \mathbf{2} \dots) + \left[ (\tau\_{11} - \tau\_{11}) \left( \mathbf{1} - \mathbf{e} \frac{-t}{\tau\_{11}} \right) \left[ \mathbf{e} \frac{D\_{11} + \sum\_{n=1}^{n-1} \frac{p\_{n}}{2} - 2 \frac{\sum\_{n=1}^{n} \tau\_{11}}{\tau} + \dots + 2 \mathbf{e} \frac{\rho}{\mathbf{n}} - 2 \mathbf{1} \right] \right] \right]$$

$$\left[ \mathbf{y}\_1(t - \mathbf{1}) + \mathbf{y}\_2(t - \mathbf{1}) + \mathbf{y}\_3(t - \mathbf{1}) + \dots + \mathbf{y}\_n(t - \mathbf{1}) \right] \tag{9}$$

The generalized analytical expression for g12 interactive transfer function of 2x2 MIMO process is given by: Similarly, the generalized analytical expression for G21 interactive transfer function of 2x2 MIMO process is given by:

$$y\_m = \left(k\_{21}\mu\_+ - 2k\_{21}\mu\_0 e^{\frac{r}{23}} \left(\frac{1 - e^{\frac{-p\_1}{p\_1}}}{1 + e^{\frac{-p\_1}{p\_2}}}\right)\right) - \left(\frac{k\_{22}k\_{e2}}{r\_{12}}(\mu\_{-2}t\_1 + \mu ur)\right) \tag{10}$$

$$-\frac{k\_{22}k\_{e2}}{r\_{12}}(r\_{12} + r\_{22})\right) \times \left[\mu\_+ - 2\mu e^{\frac{-t}{2x}}\right] \left(\frac{1 - e^{\frac{-p\_-}{p\_3}}}{1 + e^{\frac{-p\_-}{p\_3}}}\right)y\_n(t - 1)$$

The term ynð Þ t � 1 in the above Eqs. (9) and (10) is one step ahead prediction of ynð Þt .

#### **1.4 Boundary conditions to estimate the model parameters of 2X2 mimo system (two tanks)**

There are twelve parameters to be estimated in the 2x2 MIMO process. The four dead times D11,D12,D21,D22 can be obtained straight away from the initial relay response. The remaining eight parameters K11,K12,K22,K21,τ11,τ12τ<sup>22</sup> and τ<sup>21</sup> can be estimated by applying four different boundary conditions in Eqs. (9) and (10) and solving them.

First, we have to measure y1,t1,y2,t2y3,t3,y min and t min , as shown in **Figures 6** and **7**. The boundary conditions are as follows:

$$\mathbf{at}\,\mathbf{t}=\mathbf{t}\_1,\mathbf{y}=\mathbf{y}\_1\tag{11}$$

$$\mathbf{at}\,\mathbf{t}=\mathbf{t}\mathbf{z}, \mathbf{y}=\mathbf{y}\_2\tag{12}$$

$$\mathbf{at}\,\mathbf{t}=\mathbf{t}\_3,\mathbf{y}=\mathbf{y}\_3\tag{13}$$

$$\mathbf{at}\,\mathbf{t} = \mathbf{t}\_{\text{min}}\tag{14}$$

**Figure 6.** *Relay output response for g22.*

*Perspective Chapter: Modeling and Identification of Two Tank System Using Relay Feedback… DOI: http://dx.doi.org/10.5772/intechopen.106568*

**Figure 7.** *Boundary conditions from the biased* **g21** *relay response.*

#### **1.5 Procedure for parameter estimation using biased relay test**

Ideal relay feedback test is conducted on the 2x2 MIMO systems. The model Parameters are estimated as follows:


From the initial part of desirable relay response of 2x2 MIMO process.


#### *1.5.1 Parameter estimation of two tank system*

The dead times of the process is directly recorded from the responses. The information obtained during the stable oscillating condition for the process is given in **Table 4**.

Using this information, the process parameters are estimated by applying the boundary conditions given in Eqs. (10–13) in Eq. (10) and Eq. (11). The comparisons between the parameters and transfer function of true and estimated process are given **Tables 5** and **6** respectively.


#### **Table 4.**

*Parameters computed from biased relay feedback responses.*


#### **Table 5.**

*Comparison of process parameters of true process with estimated process.*


#### **Table 6.**

*Comparison of True process and estimated process.*

It is found that the estimated model parameters are very close to the true process Parameters.

#### **1.6 Interaction analysis using relative gain array (RGA)**

• From the estimated model parameters information on interaction is obtained by using RGA as follows: The steady (gain) model is expressed as in given in equation *Perspective Chapter: Modeling and Identification of Two Tank System Using Relay Feedback… DOI: http://dx.doi.org/10.5772/intechopen.106568*

$$k = |\mathbf{11}\begin{bmatrix} k & k \\ & \mathbf{12} \\ k & k \end{bmatrix} = |\begin{bmatrix} \mathbf{16.900} & \mathbf{6.5900} \\\\ \mathbf{9.200} & \mathbf{11.32} \end{bmatrix} \tag{15}$$

• The relative gain array for a 2 � 2 MIMO (Two tank) system can be expressed as

$$
\Lambda = \begin{bmatrix} 1.4639 & -0.4639 \\\\ -0.4639 & 1.4639 \end{bmatrix} \tag{16}
$$

• Pair the controlled and manipulated variables so that corresponding relative gains

Are positive and as close to one as possible. The input-output pairing of 2x2 MIMO systems is y1 and u1,y2 and u2.

#### **2. Real time implementation of two tank experimental setup using biased relay feedback test**

#### **2.1 Introduction**

Multivariable process systems are either time-lag dominated or dead-time dominated systems and can be modeled as First Order Plus Dead time (FOPDT) structures. The sequential identification and modified Ziegler-Nichols controller design method form the basic structure for the multivariable auto tuner. This method considers desired response (output: for example, y1 in case of 2x2 system) and input (u1) for system identification and analysis and they do not discuss about methods to reduce interactions nor do they give any exact analytical expressions for the relay responses that may help in analyzing the interaction behavior between input/output and may provide information regarding closed-loop parameters (PID using model based tuning rules) of the MIMO system. Moreover, it felt that exact model parameters and information on interactions can be better obtained / calculated from mathematical model of undesired relay responses for MIMO systems. Biased relay tests are carried out on 2 x 2 MIMO processes and undesirable responses are modeled.

#### **2.2 Process description**

The real time process considered for the sequential auto-tuning algorithm (**Figure 8**) is a multiple tank process in which two tanks are considered. The MIMO process considered here is the couple tank system, which consists of two cylindrical tanks connected in interacting fashion at two different heights. The system is configured as a MIMO system with two input and two output variables.

The couple tank system, which consists of two cylindrical tanks connected in interacting fashion at two different heights. To measure the interaction packers' valve are being used in the connecting lines between the two tanks. Water is pumped into the two tanks from the reservoir using the pumps. The levels of water in the two tanks are measured using differential pressure transmitters (DPT). Control valves are provided at the inflow lines of the two tanks in order to regulate the flow of water into the

tanks. To measure the interaction, the valve is being used. Two I/P converters are used to change stem position of control valve and to measure the level in the tanks two DPT are used. All these sensors and actuators are connected through NI-6008 interface DAQ card; it has 8 analog inputs, 2 analog outputs, 4 digital inputs and 4 digital outputs compatible with LabVIEW software.

#### **2.3 Biased relay feedback test**

LabVIEW (short for Laboratory Virtual Instrumentation Engineering Workbench) is a system design platform and development environment for a programming language from National Instruments. LabVIEW is commonly used for Instrument control, and industrial automation. LabVIEW provides three key elements. They are Data acquisition tools, Data analysis tools and Data visualization tools. Data acquisition is the process of gathering or generating information in an automated fashion from analog and digital measurement sources such as sensors and devices under test. Here DAQ 6008 is used for Data Acquisition.

A key benefit of Lab VIEW over other development environments is the extensive support for accessing instrumentation hardware. Lab VIEW provides tight softwarehardware integration. It has the ability to solve and execute complex algorithms in real time.

#### *2.3.1 Sequential autotuning*

The sequential identification and modified Ziegler-Nichols controller design method form the basic structure for the multivariable auto tuner. This method considers undesired response (output: for example, y1 in case of 2x2 system) and input (u1) for system identification and analysis, they give exact analytical expressions for the relay responses that may help in analyzing the interaction behavior between input/output and may provide information regarding closed-loop parameters (PID using model based tuning rules) of the MIMO system.

*Perspective Chapter: Modeling and Identification of Two Tank System Using Relay Feedback… DOI: http://dx.doi.org/10.5772/intechopen.106568*

#### *2.3.2 Sequential autotuning—step 1*

Based on the concept of sequential auto tuning method each controller is designed in sequence. Consider a 2 x 2 MIMO system with a known pairing (y1–u1**)** and (y2u2) under decentralized PI control. Relay is placed between y1 and u1, while loop 2 is on manual mode. Following the relay-feedback test, a controller can be designed from the ultimate gain and ultimate frequency. The block diagram shows in **Figure 9** gives the sequential auto tuning step 1, the DAQ Assistant is used to acquire analog signals from the level transmitters which measure the level of the tank 1 and tank2. DAQ Assistant2 is used to generate analog signals to control the final control elements. The block also includes MATLAB script which performs the relay operation.

Relay is placed between y1 and u1, while loop 2 is on manual mode and a controller can be designed from the ultimate gain and ultimate frequency of the relay response (**Figure 10**).

#### *2.3.3 Sequential autotuning—step 2*

• Perform the relay-feedback test between y2 and u2 while loop 1 is on automatic mode. A controller can also be designed for loop 2 following the relay-feedback test.

The block diagram shows in **Figure 11** gives the sequential auto tuning step 2, the DAQ Assistant is used to acquire analog signals from the level transmitters which

**Figure 9.** *Block diagram of biased relay feedback test (step - 1).*

**Figure 10.** *Biased relay feedback response (step 1) obtained from two tank experimental setup.*

measure the level of the tank 1 and tank2. DAQ Assistant2 is used to generate analog signals to control the final control elements. The block also includes MATLAB script which performs the relay operation. Relay is placed between y2 and u2, while loop 1 is on automatic mode. A controller is designed from the ultimate gain and ultimate frequency (**Figure 12**) for the loop 2.

#### *2.3.4 Sequential autotunng—step 3*


The block diagram shows in **Figure 13** gives the sequential auto tuning step 3, the DAQ Assistant is used to acquire analog signals from the level transmitters which measure the level of the tank 1 and tank 2. DAQ Assistant 2 is used to generate analog signals to control the final control element. The block also includes MATLAB script which performs the relay operation. Another relay-feedback experiment will be

*Perspective Chapter: Modeling and Identification of Two Tank System Using Relay Feedback… DOI: http://dx.doi.org/10.5772/intechopen.106568*

**Figure 11.** *Block diagram of biased relay feedback test (step 2).*

**Figure 12.** *Biased relay feedback response (step 2) obtained from two tank experimental setup.*

performed between y1 and u1 while loop 2 is on automatic mode and a new set of tuning constants will be found (**Figure 14**) for the controller in loop 1. The controller's parameters are converged in 3–4 relay feedback tests [10–12].

**Figure 13.** *Block diagram of biased relay feedback test (step 3).*

**Figure 14.** *Biased relay feedback response (step 3) obtained from two tank experimental setup.*

Sequential auto tuning is conducted under decentralized PI control, with known pairing (y1–u1) and (y2–u2) 2x2 MIMO systems is employed. First relay is placed between y1 and u1, while loop 2 is on manual mode. Following the relay-feedback test, a controller can be designed from the ultimate gain and ultimate frequency. A relay

*Perspective Chapter: Modeling and Identification of Two Tank System Using Relay Feedback… DOI: http://dx.doi.org/10.5772/intechopen.106568*

feedback test is performed between y2 and u2, while loop 1 is on automatic mode. A controller is designed from the ultimate gain and ultimate frequency for the loop 2. Afterwards relay feedback test is performed between y1 and u1, while loop 2 is on automatic mode. Generally, new set of tuning constants are founded for the controller in loop1. Based on the concept of sequential auto tuning method each controller is designed in sequence. The controller parameters are converged in 3–4 relay feedback tests is shown in the **Table 7**.

#### *2.3.5 Derivation of analytical expression for two tank process*


Analytical expression is derived for the relay response obtained from the relay feedback test when biased relay is used.

The generalized analytical expression for the undesired relay response of tank 1 (2x2 MIMO process) is given by

$$\begin{split} y\_n &= \left[ k\_{12} \mu\_+^- \frac{-t}{2k12\mu\_0 e^{\tau\_{12}}} \left[ \frac{\mathbf{1} - e^{\frac{-\mu^{\tau\_{12}}}{\tau\_{12}}}}{\mathbf{1} + e^{\frac{-\mu^{\tau\_{12}}}{\tau\_{12}}}} \right] \right] - \left[ \frac{k\_{11} k\_{\epsilon 1}}{\tau\_{11}} (\mu - t\_1 + \mu t) - \frac{k\_{11} k\_{\epsilon 1}}{\tau\_{i1}} (\tau\_{i1} + \tau\_{i1}) \right] \\ &\times \left[ \mu\_+ - \frac{-t}{2\mu e^{\tau\_{11}}} \right] \left[ \frac{\mathbf{1} - e^{\frac{-\mu^{\tau\_{11}}}{\tau\_{11}}}}{\mathbf{1} + e^{\frac{-\mu^{\tau\_{11}}}{\tau\_{11}}}} \right] y\_n(t - \mathbf{1}) \end{split} \tag{17}$$

Similarly, the generalized analytical expression for the undesired relay response of tank 2 (2x2 MIMO process) is given by


**Table 7.**

*Controller Parameters computed from biased the relay feedback response of step 2 and step 3 by sequential auto tuning.*

$$\begin{split} y\_n &= \left[ k\_{12} \mu\_+^- \frac{-t}{2k2\mu\_0 e^{\tau\_{12}}} \left[ \frac{1 - e^{\frac{-\mu r^{\tau\_{12}}}{\tau\_{\tau 1}}}}{1 + e^{\frac{-\mu r}{\tau\_{\tau 1}}}} \right] \right] - \left[ \frac{k\_{22} k\_{\epsilon 2}}{\tau\_{\tau 2}} (\mu - t\_1 + \mu t) - \frac{k\_{22} k\_{\epsilon 2}}{\tau\_{\tau 2}} (\tau\_{\tau 1} + \tau\_{\tau 1}) \right] \\ &\times \left[ \mu\_+ - \frac{-t}{2\mu e^{\tau\_{11}}} \right] \left[ \frac{1 - e^{\frac{-\mu r^{\tau 2}}{\tau\_{\tau 2}}}}{1 + e^{\frac{-\mu r}{\tau 2}}} \right] y\_n (t - 1) \end{split} \tag{18}$$

The term ynð Þ <sup>t</sup>‐<sup>1</sup> in the above Eqs. (17) and (18) is one step ahead prediction of yn(t).

#### **3. Conclusion**

Relay feedback test is conducted on both simulation and real time two tank experimental setup using biased relay and the responses are obtained. Analytical expressions are derived for the undesired relay responses in time domain using biased relay feedback test. The model parameters are estimated for the two tank system (simulation). The simulation results indicate that the estimation of Kp (Process gain) is more accurate when biased relay is used. Real time implementation of two tank experimental set up using biased relay feedback test is done.

#### **Author details**

Devarapalli Kishore<sup>1</sup> \* and Vaska Lokesh<sup>2</sup>

1 Department of EIE, University College of Engineering, Adikavi Nannaya Univerity, Rajamahendravaram, Andhrapradesh, India

2 Department of Geology, University College of Science and Technology, Adikavi Nannaya University, Rajamahendravaram, Andhrapradesh, India

\*Address all correspondence to: kishorenit.trichy@gmail.com

© 2022 The Author(s). Licensee IntechOpen. This chapter is distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

*Perspective Chapter: Modeling and Identification of Two Tank System Using Relay Feedback… DOI: http://dx.doi.org/10.5772/intechopen.106568*

#### **References**

[1] Åström KJ, Hägglund T. Automatic tuning of simple regulators with specification on phase angle and amplitude margins. Automatica. 1984; **20**:645-651

[2] Yu CC. Auto tuning of PID Controllers: A Relay Feedback Approach. 2nd ed. London: Springer-Verlag; 2006

[3] Tao-Liu, Qing-Guowang, Huang H-p. A tutorial review on process identification from step or relay feedback test. Journal of Process Control. 2013;**23**(10):1597-1623

[4] Prasad S, Jaya Y. Process Identification, Identification and control of Relay-Stabilizable Processes. Germany: VDM Verlag Publishers; 2009. pp. 1-84

[5] Sujatha V, Panda RC. Estimation of parameter for disturbance model of 3 by-3 MIMO process using relay feedback test. Journal of Modeling, Simulation, Identification and Control. 2013;**1**:27-40

[6] Panda R, Sujatha V. Identification and control of Multivariable systems-Role of Relay Feedback. In: Introduction to PID Controller-Theory, Tuning and Applications to Frontier Areas. 2012

[7] Kalpana D, Thygarajan T, Thenral R. Improved identification and control of 2-by-2 MIMO system using relay feedback. IEEE. 2015;**17**:23-32

[8] Chidambaram M, Sree RP. Control of Unstable and Single and Multi-Variable Systems. New Delhi: Narosa Publishing house Pvt. Ltd; 2017

[9] Chidambaram M, Viveksathe. Relay Auto Tuning for Identification and Control. 1st ed. UK: Cambridge Press; 2014

[10] Oskar Vivero & Jesus Liceaga. MIMO Tool Box for MATLAB, 978-1- 4244-2156. 2008

[11] Wang Q-G, Chieh C, Zou B. Multivariable process identification and control from decentralized relay feedback. IEEE. 2000;**20**:4

[12] Shen Y, Cai WJ, Li S. Normalized Decoupling Control for High Dimensional MIMO Process for Application in room Temperature Control HVAC Systems. Elsevier; 2010. pp. 652-654

#### **Chapter 3**

## PID Cascade Controller Design for an Unstable System

*Mustafa Saad*

#### **Abstract**

In control engineering, the control of an unstable system is very concerning. An example of an unstable system that uses control principles is a ball and beam balancer system. This system consists of a long beam attached to a motor at its midpoint. A steel ball moves on the top of the beam with an acceleration proportional to the beam angle. If the system is uncontrolled well, the steel ball may fall from the beam. This paper presents an approach to modeling and controlling the ball position for the ball and beam balancer system. Starting with a mathematical equation for the nonlinear system, the model of the system is produced. Then, it demonstrates the design of the PID cascade controller system to stabilize the system and regulate the ball to its reference position. The performance of the system was evaluated and tested for setpoint tracking signal and disturbance rejection test. The simulation studies were done using Matlab Simulink and the results indicated that the proposed approach yields robust closed-loop performance.

**Keywords:** unstable system, ball and beam balancer, modeling, controller design, cascade, PID, tuning, disturbance rejection

#### **1. Introduction**

Most industrial processes can be considered nonlinear and uncertain processes. The control of these processes is considered dangerous if they were unstable. An unstable system such as a ball and beam balancer exists in some control laboratories. It is used by control engineering students to design and apply some of the controller techniques to control the ball position [1]. It is easy to understand that many classical and modern controller design techniques still can be applied to this system. Besides that, due to its simple construction, it can be classified as a feedback control system and it is used to learn how to stabilize the unstable system. It is normally related to real control problems such as horizontally stabilizing an airplane during landing and in turbulent airflow [2, 3].

The ball and beam system is seen as a standard control engineering system whose fundamental theory can be used to stabilize problems for various systems. For example, the balancing problem in moving robots that use to carry and spacecraft position control systems in aerospace engineering [4].

The ball and beam system is considered a nonlinear system [5]. Its nonlinearity is due to the dead zone and saturation characteristic, DC motor, and pulley drive nonlinearity and discontinuity of position measurement. The ball and beam system can show a general nonlinear control object. The model of the ball and beam balancer system can be used as a control item such as balancing mobile robots in goods carrying, controlling spacecraft, control of nonlinear actuators, and position control of space vehicles in aerospace engineering can all use the model, the ball, and the beam system as the control object [6].

A steel ball is positioned on the beam as shown in **Figure 1**. When the current control signal is applied to the motor, the beam is titled about its center axis. This causes the rolling of the ball on the beam. The control objective is to regulate the ball position on the beam automatically by changing the manipulated beam angle. This is a difficult mission due to the rolling of the ball on the length of the beam. This movement has an acceleration proportional to the angle of the beam [7]. The system is bounded input unbounded output because bounded beam angle gives unbounded ball position. Hence, the system must be controlled to position the ball at the desired position.

Most real systems are nonlinear and closed-loop controllers that can be a valuable approach to ensure sufficient performance. Many nonlinear control systems have been documented for academic education to learn about feedback control in control engineering courses. As a standard, the ball and beam is a classic example of such a system [8].

The ball and beam system has been studied and controlled using several control techniques by many researchers. In the previous, a good variety of methods is applied to this system [9–13].

#### **2. System mathematical modeling**

The ball and beam balancer system is a multi-loop system that contains two parts. The modeling of DC motor and the model derivation ball and beam. The system freebody of the ball and beam balancer system as shown in **Figure 2** illustrates two DOF systems. One is moving the ball up and down on the beam and the other one is the beam rotating around its center. The demonstration of the whole system dynamics is

**Figure 1.** *Ball beam balancer system.*

*PID Cascade Controller Design for an Unstable System DOI: http://dx.doi.org/10.5772/intechopen.105669*

**Figure 2.** *System free body [14].*

**Figure 3.** *Ball and beam balancer motion.*

complex therefore a simple mathematical model derivation is done to design the system controller.

The ball and beam balancer motion is shown in **Figure 3**. It can be seen that three forces acting on the ball are the rolling constraint force, the sliding component force due to gravity, which depends on the beam angle ∅, and the reacting force.

The balancing forces are;

$$
\sum F\_b = m \text{g} \sin \bigotimes -F\_r = m \ddot{\mathbf{x}} \tag{1}
$$

*Fr* is the ball rolling force and *x* is the ball position on the beam. By geometry, the ball position is rewritten as

$$
\mathfrak{x} = a \mathfrak{a}'\tag{2}
$$

*α* = angular displacement of the ball.*a*´ = distance between the axis of the ball and the point of contact of the ball with the beam

The ball balancing torque *τ<sup>b</sup>* is:

$$
\sum \mathbf{r}\_b = F\_r \mathbf{\tilde{a}} = J\_b \mathbf{\tilde{a}} \tag{3}
$$

$$J\_b = \frac{2}{5}ma^2\tag{4}$$

where*Jb*= ball moment inertia.*a*= ball radius.

The beam and motor torque can be derived. Since the beam bears the load of the ball and the input motor torque. The torque balance equation is:

$$
\sum \tau\_{bm} = \tau\_{\dot{m}} = J\_{bm} \ddot{\mathcal{Q}} \tag{5}
$$

where *bm* indicates the beam and motor and *τin* represents the torque produced by the motor.

$$
\tau\_{\rm in} = k\_t I\_{\rm in} \tag{6}
$$

Then, the obtained equations after simplification and substitution are;

$$
\left[1 + \frac{2}{5} \left(\frac{a}{\check{a}}\right)^2\right] \check{x}^\circ = \text{gsin}\mathfrak{D} \tag{7}
$$

and

$$J\_{bm}\ddot{\mathcal{Q}} = k\_t I\_{in} \tag{8}$$

For linearization, it is assumed that the system operates at about 0° beam angle, and for small angles approximation sin ∅ ¼ ∅. Then Eq. (7) can be rewritten as

$$\left[\mathbf{1} + \frac{2}{5} \left(\frac{a}{\acute{a}}\right)^2\right] \ddot{\mathbf{x}} = \mathbf{g} \mathfrak{Q} \tag{9}$$

$$\left(\frac{a}{\acute{a}}\right) \cong \mathbf{1}$$

$$\left[\mathbf{1} + \frac{2}{5}\right] \ddot{\mathbf{x}} = \mathbf{g} \mathfrak{Q} \tag{10}$$

Two transfer functions can be obtained, the first transfer function for the motor system relates to beam angles ∅ and the current input*Iin*. This transfer function is the inner loop transfer function. By taking Laplace transform of Eq. (8)

$$\frac{\mathcal{Q}(s)}{I\_{in}(s)} = \frac{k\_t}{J\_{bm}s^2} \tag{11}$$

The second transfer function is for the ball and beam system that relates to the ball position *x* and beam angle ∅is obtained by taking Laplace to transform of Eq. (10). The parameters of the system that have been used in this paper are given in **Table 1**.

$$\frac{\mathbf{x}(s)}{\mathbf{x}(s)} = \frac{\mathbf{5g}}{7} \mathbf{\*} \frac{\mathbf{1}}{s^2} \tag{12}$$

*PID Cascade Controller Design for an Unstable System DOI: http://dx.doi.org/10.5772/intechopen.105669*


#### **Table 1.**

*Parameters of the system.*

#### **Figure 4.**

*Feedback control system with PID controller.*

The overall open-loop system transfer function related to ball position to the current input can be written as.

$$\frac{\varkappa(s)}{I\_{in}(s)} = \frac{\mathsf{5g}k\_t}{\mathsf{7f}\_{bm}} \bullet \frac{\mathbf{1}}{s^4} \tag{13}$$

#### **3. PID controller**

A proportional–integral–derivative controller (PID controller) is a generic feedback control mechanism as shown in **Figure 4**. PID controller is the most commonly used feedback controller. The controller is used to minimize the difference between the reference input and the controlled variable by manipulating the manipulated variable. PID controllers are the top first choice in the unknown process. However, the PID controller gains should be tuned carefully to obtain the greatest performance.

The PID controller design includes three parameters: proportional, integral, and derivative gains, denoted as P, I, and D, respectively [15, 16]. The proportional gain determines the reaction based on the present error, the integral gain calculates the reaction based on the sum of recent errors, and the derivative gain determines the reaction based on the rate of change of error. The PID controller can operate in different control mods since some processes need to use one mode or two mods to obtain the desired control. This is done by eliminating undesired control action and

putting its gain to zero. PI, PD, P, or I controller mods is derived from standard PID controller in the nonappearance of the relevant control actions. PI controllers are the most commonly used in industrial applications since derivative action is sensitive to noise. On the other hand, the presence of an integral action makes the system output reaches its setpoint.

The PID controller operates on the error signal *e(t)*, which is the difference between the desired setpoint *r(t)* and the measured output *y(t)* to yield a control signal *u(t).* PID controller has the general form.

$$u(t) = K\_p e(t) + K\_i \int\_0^t e(t)dt + K\_d \frac{d}{dt} e(t) \tag{14}$$

The desired closed-loop system characteristics are achieved by correcting the controller parameters *Kp*, *Ki*, and *Kd*, often by "tuning." The proportional term may only achieve stability. The integral term guarantees the disturbance rejection. The derivative term is responsible for response shaping. Even though, the PID controllers are the most widely used class of control systems. It cannot be used in MIMO systems due to system complicity.

Applying Laplace transformation of Eq. (14) results in the transformed PID controller equation

$$\frac{U(s)}{E(s)} = K\_p + \frac{K\_i}{s} + K\_d s \tag{15}$$

where the controller parameters are:

Proportional gain, *Kp.*

A large value gives a faster response and a very large value may cause an oscillatory and unstable system.

Integral gain, *Ki.*

A large value eliminates the steady-state error more quickly and increases overshoot. Very large values invite instability, integrator windup, or actuator saturation.

Derivative gain, *Kd* with a.

With a larger value, the response reaches the desired response faster, decreases overshoot, slows the transient response, and may cause instability.

The desired characteristic specification involves a good tuning control using adjusting the controller parameters to their optimum values. Sometimes PID controllers often offer satisfactory control; even their parameters did not tune. However, carefully tuning can improve the system performance, and poor tuning gives an unacceptable performance.

PID tuning is a challenging problem, the complex performance criteria should be achieved even with the limitations of PID control [17]. Different methods for loop tuning and techniques are more sophisticated and subjected to patents. In manual methods for loop tuning, if the system is online, the tuning method requires setting integral and derivative gains to zero and increasing the proportional gain until the output response becomes oscillatory. According to the "quarter amplitude decay" response, this gain is set to half of that value. Next, increase the integral gain until zero error is obtained for the process. Finally, increase the derivative gain, if necessary, until the output response reaches its steady state in presence of disturbance. However, too much high derivative gain will cause the output response overshoots the desired input and some systems cannot accept overshoot. **Table 2** describes the effectiveness


**Table 2.**

*Effect of increasing Kp, Ki, and Kd parameters in the closed-loop system.*


#### **Table 3.**

*PID controller tuning methods.*

of increasing controller gains on the system response. The table is used as a reference and only help in calculating the values of controller parameters. Because these parameters are dependent on each other, changing one of these parameters can change the effect of the other two.

There are different control techniques for tuning a PID controller. Many of these techniques require the transfer function of the system, then selecting controller gains *Kp, Ki,* and *Kd* depending on the system dynamics. The PID controller tuning method should be carefully chosen, for example, in practice, some loops take minutes or longer to reach steady state. Therefore, manual tuning could be an inefficient choice. In addition, tuning method selection depends on the type of tuning online or offline method. In the case of offline systems, some tuning methods require applying a step change to the system and measuring the output response. Then, use this reopens to calculate the controller gains. **Table 3** illustrates the comparison between some tuning methods.

#### **4. PID cascade controller design for ball beam balancer system**

The suggested block diagram of the ball and beam balancer system consists of twoloop as shown in **Figure 5**. The inner loop is the motor control loop, and the outer loop, which is the ball beam control loop. In this paper, the design approach has to stabilize the inner loop first, and then the outer loop is controlled. The motor loop control is in series with ball and beam loop control.

In this research, the design of the PID cascade controller [18–21] includes two control loops, an inner loop with a primary PID controller to control the beam angle, and an outer loop with a secondary PID controller to control the ball position.

**Figure 5.**

*Block diagram of ball and beam balancer system.*

The secondary controller is designed to offset the effect of disturbances before it significantly affects the output of the controlled system. While the output of the first controller provides the reference for the second loop. This reduces any unexpected changes from the inner loop. The secondary controller adjusts the motor. There is a relation between motor angle and beam angle, any change in the motor angle causes a change in the beam angle. Therefore, the ball position is controlled by changing the beam angle.

The PID controller parameters are hard to tune when system parameters besides control action change while the system operates. So, in this research, the controller gains are tuned using a tuning tool in Simulink control design based on the performance and robustness of the system. The obtained PID controller parameters to control the system are illustrated in **Table 4**. The simulation of the PID cascade diagram of the ball and beam balancer system is shown in **Figure 6**.

In this design, the ITAE criterion is used to compute the best PID controller gains. The ITAE index is the best selectivity of the performance indices because it measures the integration of error for a specific time as the system parameters are varied [7]. The ITAE index formula is:


#### **Table 4.**

*The best tuned PID controller parameters.*

**Figure 6.** *Simulation diagram of ball and beam balancer system.* *PID Cascade Controller Design for an Unstable System DOI: http://dx.doi.org/10.5772/intechopen.105669*

$$ITAE = \int\_{0}^{t} t|e(t)|dt$$

#### **5. Simulation results and discussion**

The system is an unstable system with poles at the origin. So a controller is needed to control the system by looking at the response of the system that specifies some criteria. One of the first things that must be done is to decide upon a criterion for measuring how good the response is. For example, it does not matter how the response reaches a steady state but the steady state itself either has a large error or not. However, transient behavior is also important in determining the best response. For that, some criterion is introduced such as settling time, rise time, overshoot, and steady-state error.

PID cascade controller was designed and built using Matlab Simulink as shown in the previous section. To verify its performance in regulating the ball position of the ball and beam balancer system for different positions, such as step signal, step-change tracking signal, and sinusoidal signals were used as input test types for ball position in the simulation model. To test the robustness of the controlled system, the input disturbance step signal was used.

Firstly, the inner loop was controlled by tuning the controller parameters of the first PID until getting the best performance of the inner loop. After that, the second PID in the outer loop was tuned to get a stable response from the system. **Figure 7** shows the step response of ball position using PID controllers, where the best parameters that give this response were tableted in **Table 4**.

The designed controller stabilized the system with zero steady-state error. The system has a fast response and reached the steady-state value in short time. However, the system has an overshoot; but it is still an acceptable overshoot because the ball was still on the specified beam length. The performance specifications of the system using the PID cascade controller are summarized in **Table 5**.

**Figure 7.** *Step response of ball position using PID controller.*


#### **Table 5.**

*performance specifications of the system using PID controller.*

**Figure 8.** *Step tracking response of the system.*

The setpoint change was performed, starting with the ball at the middle of the beam (origin). At 1 second, the ball position changed to 0.1-meter position of the beam center, at 2 seconds the desired ball position changed to the left side of the beam at position 0.1 meter, the ball position was changed again to 0.15 meter at 4 seconds on the left side, and return to the middle of the beam at 5.5 seconds.

The setpoint-tracking signal contained changing the set point during the operation as shown in **Figure 8**. Using the PID controller, the system output perfectly tracked the setpoint changes of ball position.

In this research, the ball and beam balancer system was subjected to input disturbance with a force acting on a ball with 0.05 step at a time of 2 sec. The output response under effective of this disturbance is shown in **Figure 9**.

The result showed that the disturbance affected the system at a time between 2 and 2.5 seconds. The designed controller can overcome this disturbance and reposition the ball to its setpoint in approximately 0.5 seconds.

The controlled system is tested for sinusoidal input and the result of the ball position is shown in **Figure 10**. As shown, the output response reached the desired input in a fast time, without overshoot and almost zero steady-state error.

*PID Cascade Controller Design for an Unstable System DOI: http://dx.doi.org/10.5772/intechopen.105669*

**Figure 9.** *Output response of the system with input disturbance.*

**Figure 10.** *Output response for sinusoidal input.*

#### **6. Conclusion**

This research was successfully elaborated and the PID cascade controllers were magnificently been designed. A model for a ball and beam system is well derived such that the free-rolling ball of the ball and beam system can be positioned at any desired location on the beam without falling. For a ball and beam balancer system, the most required criterion is that the system has a small or no overshoot and zero steady-state error. From the results and discussion in the previous section. The simulation results have shown that the suggested approach can stabilize the system efficiently. Furthermore, the performance during the transient period of the system is good where a small overshoot was obtained. In addition, the PID cascade controller provided a very small settling time and zero steady-state error. Moreover, the proposed controller has successfully tracked the step tracking signal and sinusoidal signal. For the disturbance, the designed controller approved its ability to reject the effect of disturbance. Therefore, it can be concluded that the robustness of PID cascade controllers was achieved.

### **Author details**

Mustafa Saad College of Electronic Technology, Bani Walid, Libya

\*Address all correspondence to: mustafasaad9@yahoo.com

© 2022 The Author(s). Licensee IntechOpen. This chapter is distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

*PID Cascade Controller Design for an Unstable System DOI: http://dx.doi.org/10.5772/intechopen.105669*

#### **References**

[1] Márton L, Lantos B. Stable adaptive ball and beam control. In: 2006 IEEE International Conference on Mechatronics. 2006. pp. 507-512

[2] Rahmat MF, Wahid H, Wahab NA. Application of intelligent controller in a ball and beam control system. International Journal on Smart Sensing and Intelligent Systems. 2010;**3**(1): 45-60. DOI: 10.21307/ijssis-2017-378

[3] Ali AT, Taha OA. Design and implementation of ball and beam system using PID controller. Automatic Control and Information Sciences. 2017;**3**(1):1-4

[4] Oh S-K, Jang H-J, Pedrycz W. The design of a fuzzy cascade controller for ball and beam system: A study in optimization with the use of parallel genetic algorithms. Engineering Applications of Artificial Intelligence. 2009;**22**(2):261-271

[5] Wang L-X. Stable and optimal fuzzy control of linear systems. IEEE Transactions on Fuzzy Systems. 1998; **6**(1):137-143. DOI: 1063-6706(98) 00848-0

[6] Meiling Ding BL, Wang L. Position control for ball and beam system based on active disturbance rejection control. Systems Science & Control Engineering. 2019;**7**(1):97-108. DOI: 10.1080/ 21642583.2019.1575297

[7] Amjad KIM, Abdullah SS, Shareef Z. A simplified intelligent controller for ball and beam system. In: Presented at the International Conference on Education Technology and Computer. 2010

[8] Al-dujaili Q, Humaidi AJ, Pereira DA, Ibraheem IK. Adaptive backstepping control design for ball and beam system. International Review of Applied Sciences and Engineering. Jul 2021;**12**(3):211-221. DOI: 10.1556/1848.2021.00193

[9] Hauser J, Sastry S, Kokotovic P. Nonlinear control via approximate input-output linearization: The ball and beam example. In: IEEE Transactions on Automatic Control. Mar 1992;**37**(3):392- 398. DOI: 10.1109/9.119645

[10] Nguyen HNPCX, Lukianov AD, Pham TD, Nguyen TT, Nguyen HT. Design of an embedded control system based on the quasi-time optimal control law when limiting the control signal for the ball and beam system. In: XV International Scientific-Technical Conference "Dynamic of Technical Systems". 2018

[11] Mustafa Saad MK. Design and implementation of an embedded ballbeam controller using PID algorithm. Universal Journal of Control and Automation. 2017;**5**(4):63-70. DOI: 10.13189/ujca.2017.050402

[12] Meenakshipriya B, Kalpana K. Modelling and control of ball and beam system using coefficient diagram method (CDM) based PID controller. IFAC Proceedings Volumes. 2014;**47**(1): 620-626

[13] Aguilar-Ibañez C, Suarez-Castanon MS, Rubio JDJ. Stabilization of the ball on the beam system by means of the inverse Lyapunov approach. Mathematical Problems in Engineering. 2012;**2012**: 810597. DOI: 10.1155/2012/810597

[14] Rosales EA. A Ball-on-Beam Project Kit. MA, USA: Massachusetts Institute of Technology; 2004

[15] Lloyds Raja G, Ali A. New PI-PD controller design strategy for industrial unstable and integrating processes with dead time and inverse response. Journal of Controlled Automative Electrical System. 2021;**32**:266-280. DOI: 10.1007/ s40313-020-00679-5

[16] Kumari S, Aryan P, Raja GL. Design and simulation of a novel FOIMC-PD/P double-loop control structure for CSTRs and bioreactors. International Journal of Chemical Reactor Engineering. 2021; **19**(12):1287-1303. DOI: 10.1515/ijcre-2021-0140

[17] Lawrence B. Tuning of A PID Controller for Optimal Performance of Ball and Beam System. IJERT; 2020

[18] Ali GLR. Design of Cascade control structure for stable processes using method of moments. In: 2nd International Conference on Power and Embedded Drive Control. 2019

[19] Raja AAGL. Modified series cascade control strategy for integrating processes. In: Indian Control Conference. 2018

[20] Raja AAGL. Enhanced series Cascade control strategies for unstable processes. In: 4th International Conference on Electrical Energy Systems. 2018

[21] Lloyds Raja AAG. Enhanced tuning of Smith predictor based series cascaded control structure for integrating processes. ISA Transactions. 2021;**14**: 191-205. DOI: 10.1016/j. isatra.2020.12.045

#### **Chapter 4**

## MIMO PID Control Retuning Using the Closed-Loop Frequency Response

*Anna Paula V. de A. Aguiar, George Acioli Júnior and Péricles R. Barros*

#### **Abstract**

Proportional integral derivative (PID) controllers are the most used in practice for regulatory control. This is due to the good performance achieved by this controller for a variety of processes. However, about 60% have performance issue. This problem can be more evident in multivariable processes, due to the coupling between the loops. One way to improve performance is to retune the parameters. Thus, in this article, a PID control retune technique is presented. Frequency domain data are used to compute gain increments. Identification of the parametric model of the process is not necessary. The method can be applied to multivariable processes with time delay, integrative dynamics, and nonminimum phase zeros. Simulation results and the effectiveness of the method are discussed.

**Keywords:** PID controller, multivariable systems, frequency domain, retuning, process control

#### **1. Introduction**

The proportional integral derivative (PID) controllers are the most used in the industry [1]. Numerous tuning methods for PID parameters are found in the literature. Despite this, about 60% controllers do not reach the desired performance [2]. This occurs due to actuator wear, process changes, and mainly poorly due to tuned controllers.

Several methods of evaluation of control loops and controller tuning are found in the literature. Many of these controller evaluation and retuning methods are based on process models. In [3], process data are used to identify a model. Then, the PID controller is designed based on obtained model and desired performance is computed by output prediction. If the performance index is adequate, the PID controller is retuned using the model; otherwise, another process model must be identified.

However, identifying the model may not be an easy task. An alternative is to use data-based methods. In these methods, knowledge of the system parametric model is not necessary and time or frequency domain data are used.

In [2], a PID controller performance assessment and retuning method is proposed. The closed-loop step response data are used, and the controller is retuned so that it

approximates a closed-loop reference model. The reference model is a second-order plus time delay transfer function. In [4], frequency domain constraints are inserted into the problem proposed in [2] to improve robustness and ensure stability.

Only frequency domain data are used in the method presented in [5] to retune the PID controller. In the last two, the reference model is a first order plus time delay transfer function, of which the time constant is defined according to the desired stability margins.

However, these methodologies are applied to single-input single-output (SISO) processes. On the other hand, the processes are increasingly complex due the demand for high-quality and energy-saving products is ever-increasing. Multivariable (MIMO) processes are common in the industry. PID control structures for MIMO processes are classified as decentralized control and centralized control. The decentralized control is a diagonal matrix, in which each nonzero element is a PID controller. The centralized control consists of a full matrix.

Ensuring a desired performance of the PID controller in MIMO processes is an even more difficult task, due to the coupling between the loops. The coupling represents the interactions between input and output variables. Thus, MIMO PID controller performance evaluation and retuning methods are necessary.

One way to use methods developed for SISO processes in MIMO processes is to apply them sequentially and iteratively. In [6], the method presented in [4] was used to evaluate and retune the decentralized PID controller. In [7], the PI controller retuning method presented in [5] has been extended to MIMO processes.

In this paper, the retuning method presented in [7] is reviewed and extended to PID controllers. The increments of the initial MIMO PID controllers gains are computed using only frequency domain data. The process parametric model is not required. The objective is to retune the controller so that the new closed loop is as close as possible to a given reference model. Simulation examples show that the method can be applied to multivariable process with time delay, integrative dynamics, and nonminimum phase zero.

The paper is organized as follows. In Section 2, the problem statement is presented. The frequency domain retuning method for MIMO process is proposed in Section 3. The simulation and experimental results are presented in sections 4. The conclusion is presented in Section 5.

#### **2. Problem statement**

Consider a multivariable process **<sup>G</sup>**ð Þ*<sup>s</sup>* <sup>∈</sup> *<sup>n</sup>*�*<sup>n</sup>* with *<sup>n</sup>* inputs and *<sup>n</sup>* outputs and a centralized or descentralized PI/PID controller **<sup>C</sup>**ð Þ*<sup>s</sup>* <sup>∈</sup> *<sup>n</sup>*�*<sup>n</sup>*, respectively, given by:

$$\mathbf{C}(s) = \begin{bmatrix} \mathbf{C}\_{11}(s) & \mathbf{C}\_{12}(s) & \cdots & \mathbf{C}\_{1n}(s) \\ \mathbf{C}\_{21}(s) & \mathbf{C}\_{22}(s) & \cdots & \mathbf{C}\_{2n}(s) \\ \vdots & \vdots & \ddots & \vdots \\ \mathbf{C}\_{n1}(s) & \mathbf{C}\_{n2}(s) & \cdots & \mathbf{C}\_{nn}(s) \end{bmatrix},\tag{1}$$
 
$$\mathbf{C}(s) = \begin{bmatrix} \mathbf{C}\_{11}(s) & \mathbf{0} & \cdots & \mathbf{0} \\ \mathbf{0} & \mathbf{C}\_{22}(s) & \cdots & \mathbf{0} \\ \vdots & \vdots & \ddots & \vdots \\ \mathbf{0} & \mathbf{0} & \cdots & \mathbf{C}\_{nn}(s) \end{bmatrix},\tag{2}$$

where each non-null element is given by:

$$\mathbf{C}\_{\vec{\eta}}(\mathfrak{s}) = \mathbf{K}p\_{\vec{\eta}} \ \ + \ \ \ \ \frac{\mathbf{K}i\_{\vec{\eta}}}{\mathfrak{s}} \ \ + \ \ s\mathbf{K}d\_{\vec{\eta}},\tag{3}$$

and *Kpij*, *Kiij*, and *Kdij* are the proportional, integrative, and derivative gains, respectively, and *i*, *j* ¼ 1, 2, … , *n*.

The closed loop is given by:

$$\mathbf{T(s)} = \begin{pmatrix} \mathbf{I} \ \mathbf{ \ } \end{pmatrix} + \begin{array}{c} \mathbf{G(s)} \mathbf{C(s)} \end{array} \begin{pmatrix} \mathbf{\hat{z}(s)} \mathbf{C(s)} \end{pmatrix} , \tag{4}$$

where **<sup>T</sup>**ð Þ*<sup>s</sup>* <sup>∈</sup> *n*�*<sup>n</sup>* must be stable, **<sup>I</sup>** <sup>∈</sup> *<sup>n</sup>*�*<sup>n</sup>* is an identity matrix. The sensitivity function is given by:

$$\mathbf{S}(\mathfrak{s}) = \begin{pmatrix} \mathbf{I} \ & + & \mathbf{G}(\mathfrak{s})\mathbf{C}(\mathfrak{s}) \end{pmatrix}^{-1}. \tag{5}$$

Given a closed-loop reference model **T***r*ð Þ*s* and the initial closed-loop frequency response data **T***i*ð Þ *jω* on a finite frequency set Ω ¼ ½ � *ω*1, *ω*2, ⋯, *ω<sup>N</sup>* , with *ω*<sup>1</sup> >0, the problem statement is: obtain a new PID controller so that the designed closed loop is close to the desired one, without the knowledge of the parametric model of the process.

#### **3. Frequency domain retuning**

Consider an initial closed-loop **<sup>T</sup>***i*ðÞ¼ *<sup>s</sup>* ð Þ **<sup>I</sup>** <sup>þ</sup> **<sup>G</sup>**ð Þ*<sup>s</sup>* **<sup>C</sup>***i*ð Þ*<sup>s</sup>* �<sup>1</sup> **G**ð Þ*s* **C***i*ð Þ*s* with a MIMO PI/PID controller Eq. (1) and (2). The retuned controller **C**ð Þ*s* is given by:

$$
\overline{\mathbf{C}}(\mathfrak{s}) = \mathbf{C}(\mathfrak{s}) \; \; \; + \; \; \mathbf{C}^{\Delta}(\mathfrak{s}), \tag{6}
$$

where the **<sup>C</sup>**<sup>Δ</sup>ð Þ*<sup>s</sup>* parameters are again increments of the initial controller gain.

The goal of the retune is to compute the **<sup>C</sup>**<sup>Δ</sup>ð Þ*<sup>s</sup>* parameters, so that the new closed loop with the new controller **<sup>C</sup>**ð Þ*<sup>s</sup>* is close to the desired one. To compute the **<sup>C</sup>**<sup>Δ</sup>ð Þ*<sup>s</sup>* parameters, first the frequency response **<sup>C</sup>**<sup>Δ</sup>ð Þ *<sup>j</sup><sup>ω</sup>* is computed considering the iterations between the loops, as shown in lemma 1.1. The frequency response of the initial closed loop, the reference model, and the initial controller is the algorithm input data.

Lemma 1.1 *Given a desired reference* **<sup>T</sup>***r*ð Þ*<sup>s</sup> and an initial* **<sup>T</sup>***i*ð Þ*<sup>s</sup> closed loop, the* **<sup>C</sup>**<sup>Δ</sup>ð Þ *<sup>j</sup><sup>ω</sup> can be computed as:*

$$\mathbf{C}^{\Delta}(j\rho) = \mathbf{C}(j\rho)\mathbf{T}\_i(j\rho)^{-1}(\mathbf{T}\_r(j\rho) - \mathbf{T}\_i(j\rho))\mathbf{S}\_r^{-1}(j\rho),\tag{7}$$

where **S***r*ðÞ¼ *s* **I** � **T***r*ð Þ*s* is the reference model sensitivity function.

**Proof of Lemma 1.1:** The proof is found in [7].

Once the **<sup>C</sup>**<sup>Δ</sup>ð Þ*<sup>s</sup>* frequency response is computed, the gain increments can be obtained. By definition each element of **<sup>C</sup>**<sup>Δ</sup>ð Þ *<sup>j</sup><sup>ω</sup>* <sup>∈</sup> *<sup>n</sup>*�*<sup>n</sup>* is of the form:

$$\mathcal{C}^{\Delta}\_{\vec{\eta}}(\mathfrak{s}) = \mathcal{K}p^{\Delta}\_{\vec{\eta}} \ \ + \ \ \ \frac{\mathcal{K}^{\Delta}\_{\vec{\eta}}}{\mathfrak{s}} \ \ + \ \ \mathcal{K}d^{\Delta}\_{\vec{\eta}}\mathfrak{s}.\tag{8}$$

From lemma 1.1 and matrix equality, the parameters *Kp*<sup>Δ</sup> *ij* , *Ki*<sup>Δ</sup> *ij* and *Kd*<sup>Δ</sup> *ij* of each element of *C*<sup>Δ</sup> *ij*ð Þ*s* are computed as shown in lemma 1.2.

Lemma 1.2 Given the *C*<sup>Δ</sup> *ij*ð Þ*<sup>s</sup>* frequency response, parameters *Kp*<sup>Δ</sup> *ij* , *Ki*<sup>Δ</sup> *ij* , and *Kd*<sup>Δ</sup> *ij* of the **<sup>C</sup>**<sup>Δ</sup>ð Þ *<sup>j</sup><sup>ω</sup>* are computed by:

$$K p\_{\vec{ij}}^{\Delta} = \left(\boldsymbol{\Phi}\_{r\_{\vec{ij}}}^{T} \boldsymbol{\Phi}\_{r\_{\vec{ij}}}\right)^{-1} \boldsymbol{\Phi}\_{r\_{\vec{ij}}}^{T} \boldsymbol{\Omega}\_{r\_{\vec{ij}}} \tag{9}$$

$$
\begin{bmatrix} K\dot{i}\_{\dot{\boldsymbol{y}}}^{\Delta} \\ K d\_{\dot{\boldsymbol{y}}}^{\Delta} \end{bmatrix} = \left( \boldsymbol{\Phi}\_{\dot{i}\_{\dot{\boldsymbol{y}}}}^{T} \boldsymbol{\Phi}\_{\dot{i}\_{\dot{\boldsymbol{y}}}} \right)^{-1} \boldsymbol{\Phi}\_{\dot{i}\_{\dot{\boldsymbol{y}}}}^{T} \boldsymbol{\Omega}\_{\dot{i}\_{\dot{\boldsymbol{y}}}}, \tag{10}
$$

where

$$\Phi\_{r\_{\vec{\eta}}} = \begin{bmatrix} \mathbf{1} \\ \mathbf{1} \\ \vdots \\ \mathbf{1} \end{bmatrix}, \tag{11}$$

$$\Omega\_{r\_{\vec{\eta}}} = \begin{bmatrix} \Re \left( C\_{\vec{\eta}}^{\Delta} (\
jo\_1) \right) \\ \Re \left( C\_{\vec{\eta}}^{\Delta} (\
jo\_2) \right) \\ \vdots \\ \Re \left( C\_{\vec{\eta}}^{\Delta} (\
jo\_N) \right) \end{bmatrix}, \tag{12}$$

$$\boldsymbol{\Phi}\_{i\_{\vec{\eta}}} = \begin{bmatrix} -\mathbf{1} & \begin{pmatrix} \boldsymbol{o}\boldsymbol{o}\_{1} & \boldsymbol{o}\boldsymbol{o}\_{1} \\ -\mathbf{1} & \begin{pmatrix} \boldsymbol{o}\boldsymbol{o}\_{2} & \boldsymbol{o}\boldsymbol{o}\_{2} \\ \vdots & \vdots \\ -\mathbf{1} & \begin{pmatrix} \boldsymbol{o}\boldsymbol{o}\_{N} & \boldsymbol{o}\boldsymbol{o}\_{N} \end{pmatrix} \end{bmatrix}, \tag{13}$$
 
$$\begin{bmatrix} \boldsymbol{\sigma} \left( \boldsymbol{C}^{\Delta} \boldsymbol{(i}\_{\vec{\eta}} \boldsymbol{o}\_{N} \right) \end{bmatrix} \end{bmatrix}$$

$$\Omega\_{i\_{\vec{\imath}}} = \begin{bmatrix} \mathfrak{S}\left(C\_{\vec{\imath}\vec{\jmath}}^{\Delta}(\ jo\_1)\right) \\ \mathfrak{S}\left(C\_{\vec{\imath}\vec{\jmath}}^{\Delta}(\ jo\_2)\right) \\ \vdots \\ \mathfrak{S}\left(C\_{\vec{\imath}\vec{\jmath}}^{\Delta}(\ jo\_N)\right) \end{bmatrix},\tag{14}$$

*C*<sup>Δ</sup> *ij*ð Þ *jω* is computed as shown in lemma 1.1, ℜð Þ*:* and ℑð Þ*:* are the real and imaginary parts, respectively, Ω ¼ ½ � *ω*1,*ω*2, ⋯, *ω<sup>N</sup>* is finite frequency set, and *ω*<sup>1</sup> >0.

**Proof of Lemma 1.2:** As **<sup>C</sup>**<sup>Δ</sup>ð Þ *<sup>j</sup><sup>ω</sup>* is a complex matrix so each element is given by:

$$\mathbf{C}\_{\vec{\eta}}^{\Delta}(j\rho) = \mathfrak{a}\_{\vec{\eta}} \; \; \; \; \; jb\_{\vec{\eta}}, \tag{15}$$

where *aij* <sup>¼</sup> <sup>ℜ</sup> **<sup>C</sup>**<sup>Δ</sup> *ij*ð Þ *jω* � � and *bij* <sup>¼</sup> <sup>ℑ</sup> **<sup>C</sup>**<sup>Δ</sup> *ij*ð Þ *jω* � �. For each frequency point, we have:

$$\mathbf{C}^{\Delta}(j\omega) = \begin{bmatrix} a\_{11} & + \ j b\_{11} & a\_{12} & + \ j b\_{12} & \cdots & a\_{1n} & + \ j b\_{1n} \\ a\_{21} & + \ j b\_{21} & a\_{22} & + \ j b\_{22} & \cdots & a\_{2n} & + \ j b\_{2n} \\ & \vdots & & \vdots & \ddots & & \vdots \\ a\_{n1} & + \ j b\_{n1} & a\_{n2} & + \ j b\_{n2} & \cdots & a\_{nn} & + \ j b\_{nn} \end{bmatrix},\tag{16}$$

*MIMO PID Control Retuning Using the Closed-Loop Frequency Response DOI: http://dx.doi.org/10.5772/intechopen.105015*

**K***p* þ **K***i <sup>j</sup><sup>ω</sup>* <sup>þ</sup> *<sup>j</sup>ω***K***<sup>d</sup>* <sup>¼</sup> *a*<sup>11</sup> þ *jb*<sup>11</sup> *a*<sup>12</sup> þ *jb*<sup>12</sup> ⋯ *a*1*<sup>n</sup>* þ *jb*1*<sup>n</sup> a*<sup>21</sup> þ *jb*<sup>21</sup> *a*<sup>22</sup> þ *jb*<sup>22</sup> ⋯ *a*2*<sup>n</sup>* þ *jb*2*<sup>n</sup>* ⋮ ⋮⋱⋮ *an*<sup>1</sup> þ *jbn*<sup>1</sup> *an*<sup>2</sup> þ *jbn*<sup>2</sup> ⋯ *ann* þ *jbnn* 2 6 6 6 4 3 7 7 7 5*:* (17)

Thus, for each element:

$$K p\_{\vec{\eta}}^{\Delta} \quad + \quad \frac{K \dot{t}\_{\vec{\eta}}^{\Delta}}{j o} \quad + \quad j o \mathcal{K} d\_{\vec{\eta}}^{\Delta} = \mathfrak{a}\_{\vec{\eta}} \quad + \quad j b\_{\vec{\eta}}, \tag{18}$$

$$K p\_{\vec{ij}}^{\Delta} \, \left. \, \frac{K \dot{\mathbf{r}}\_{\vec{ij}}^{\Delta}}{j \alpha} \right. \left. \, \left. \, \left. j \alpha \text{K} d\_{\vec{ij}}^{\Delta} \right| \, \left. \, \text{F} \left( j \alpha \right) \right) \right. \left. \, \left. \, \left. j \mathfrak{T} \right| \mathcal{C}\_{\vec{ij}}^{\Delta} (j \alpha) \right) \right. \tag{19}$$

$$\text{Kp}\_{\vec{\text{ij}}}^{\Delta} - j \frac{\text{Ki}\_{\vec{\text{ij}}}^{\Delta}}{\alpha} \; \left. + \; j\alpha \text{Kd}\_{\vec{\text{ij}}}^{\Delta} = \Re \left( \mathbf{C}\_{\vec{\text{ij}}}^{\Delta} (j\alpha) \right) \; \left. + \; j\Re \left( \mathbf{C}\_{\vec{\text{ij}}}^{\Delta} (j\alpha) \right) \right) . \tag{20}$$

Separating the real and imaginary terms of each element, we have:

$$K p\_{ij}^{\Delta} = \Re \left( \mathbf{C}\_{ij}^{\Delta} (j\rho) \right), \tag{21}$$

$$-\frac{K\dot{i}\_{\vec{\eta}}^{\Delta}}{\alpha} \quad + \quad \alpha \text{K}d\_{\vec{\eta}}^{\Delta} = \mathfrak{F}\Big(\mathbf{C}\_{\vec{\eta}}^{\Delta}(j\omega)\Big). \tag{22}$$

Considering a set with *N* points frequencies points:

$$
\begin{bmatrix} 1 \\ 1 \\ \vdots \\ 1 \end{bmatrix} K p\_{\vec{y}}^{\Delta} = \begin{bmatrix} \Re \left( \mathbf{C}\_{\vec{y}}^{\Delta} / jo\_{1} \right) \\ \Re \left( \mathbf{C}\_{\vec{y}}^{\Delta} / jo\_{2} \right) \\ \vdots \\ \Re \left( \mathbf{C}\_{\vec{y}}^{\Delta} / jo\_{N} \right) \end{bmatrix}, \tag{23}
$$

$$
\begin{bmatrix} -\mathbf{1} & \prime & o\_{1} & o\_{1} \\ -\mathbf{1} & \prime o\_{2} & o\_{2} & o\_{2} \\ \vdots & \vdots & \vdots \\ -\mathbf{1} & \prime o\_{N} & o\_{N} \end{bmatrix} \begin{bmatrix} K\_{\vec{y}}^{\Delta} \ K d\_{\vec{y}}^{\Delta} \\ \Re \left( \mathbf{C}\_{\vec{y}}^{\Delta} / jo\_{2} \right) \\ \Re \left( \mathbf{C}\_{\vec{y}}^{\Delta} / jo\_{2} \right) \\ \vdots \\ \Re \left( \mathbf{C}\_{\vec{y}}^{\Delta} / jo\_{N} \right) \end{bmatrix}. \tag{24}
$$

A particular case, presented in [7] is given by lemma 1.3, where a PI controller is considered:

$$\mathbf{C}\_{\vec{\eta}}^{\Delta}(\mathfrak{s}) = K p\_{\vec{\eta}}^{\Delta} \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \;$$

Lemma 1.3 Given the *C*<sup>Δ</sup> *ij*ð Þ*<sup>s</sup>* frequency response, parameters *Kp*<sup>Δ</sup> *ij* and *Ki*<sup>Δ</sup> *ij* of the **<sup>C</sup>**<sup>Δ</sup>ð Þ *<sup>j</sup><sup>ω</sup>* are computed by:

$$K p\_{\vec{ij}}^{\Delta} = \left(\Phi\_{r\_{\vec{ij}}}^T \Phi\_{r\_{\vec{ij}}}\right)^{-1} \Phi\_{r\_{\vec{ij}}}^T \Omega\_{r\_{\vec{ij}}} \tag{26}$$

$$K\dot{\mathbf{i}}\_{\dot{\boldsymbol{ij}}}^{\Delta} = \left(\boldsymbol{\Phi}\_{i\_{\dot{\boldsymbol{ij}}}}^{T}\boldsymbol{\Phi}\_{i\_{\dot{\boldsymbol{ij}}}}\right)^{-1}\boldsymbol{\Phi}\_{i\_{\dot{\boldsymbol{ij}}}}^{T}\boldsymbol{\Omega}\_{i\_{\dot{\boldsymbol{ij}}}},\tag{27}$$

where

$$\Phi\_{r\_{\vec{\imath}}} = \begin{bmatrix} \mathbf{1} \\ \mathbf{1} \\ \vdots \\ \mathbf{1} \end{bmatrix}, \tag{28}$$

$$\Omega\_{\overline{\gamma}} = \begin{bmatrix} \Re\left(C\_{\overline{\gamma}}^{\Delta}(\
jo\_1)\right) \\ \Re\left(C\_{\overline{\gamma}}^{\Delta}(\
jo\_2)\right) \\ \vdots \\ \Re\left(C\_{\overline{\gamma}}^{\Delta}(\
jo\_N)\right) \end{bmatrix},\tag{29}$$

$$\Phi\_{\overline{\gamma}} = \begin{bmatrix} -1 \ / \ o\_1 \\ -1 \ / \ o\_2 \end{bmatrix}$$

$$\Phi\_{i\bar{\eta}} = \left\lfloor \begin{array}{c} \cdot \\ \vdots \\ -\mathbf{1} \end{array} \right\rfloor,\tag{30}$$

$$\left\lceil \mathfrak{S} \left( C\_{\bar{\eta}}^{\Delta} (\ j o\_1) \right) \right\rceil$$

$$\Omega\_{i\_{\vec{\eta}}} = \begin{bmatrix} \mathbf{C}^{\Delta}(\cdot)^{\Delta - 1} \\ \mathfrak{S} \left( \mathbf{C}^{\Delta}\_{\vec{\eta}}(\cdot jo\_2) \right) \\ \vdots \\ \mathfrak{S} \left( \mathbf{C}^{\Delta}\_{\vec{\eta}}(\cdot jo\_N) \right) \end{bmatrix}, \tag{31}$$

*C*<sup>Δ</sup> *ij*ð Þ *jω* is given by Eq. (7), ℜð Þ*:* and ℑð Þ*:* are the real and imaginary parts, respectively, Ω ¼ ½ � *ω*1, *ω*2, ⋯, *ω<sup>N</sup>* is finite frequency set and *ω*<sup>1</sup> >0.

**Proof of Lemma 1.3:** The proof is similar to that of lemma 1.2 and can be found in [7].

#### **4. Simulation results**

In this section, the effectiveness of the proposed PID controller retune method through simulated examples. For this, the Wood-Berry, the integrating process of the distillation column and drum boiler are considered. The initial controller can be centralized or decentralized PI or PID type. The maximum singular value of the sensitivity function (**S**ð Þ*s* ) is used to show the closed-loop robustness property. The maximum singular value must be less than 2 to guarantee greater stability margin and robustness [8].

#### **4.1 Example 1**

Consider the Wood-Berry binary distillation column [9] given by 2 � 2 matrix, of which each element is a first-order plus time delay transfer function:

*MIMO PID Control Retuning Using the Closed-Loop Frequency Response DOI: http://dx.doi.org/10.5772/intechopen.105015*

$$\mathbf{G}\_{\rm ex1}(s) = \begin{bmatrix} \frac{12.8e^{-s}}{16.7s} + \frac{-18.9e^{-3s}}{21s} \\ \frac{6.6e^{-7s}}{10.9s} + \frac{-19.4e^{-3s}}{14.4s} \end{bmatrix} \tag{32}$$

and the decentralized PID controller:

$$\mathbf{C}\_{\text{ex1}}(\mathbf{s}) = \begin{bmatrix} \mathbf{0.34} & + \ \frac{\mathbf{0.0198}}{\mathfrak{s}} & + \ \mathbf{0.167s} & & \mathbf{0} \\\\ & \mathbf{0} & -\mathbf{0.14} - \frac{\mathbf{0.0081}}{\mathfrak{s}} - \mathbf{0.19s} \end{bmatrix}.\tag{33}$$

To obtain the closed loop as decoupled as possible, the reference model is given by the diagonal matrix given by:

$$\mathbf{T}\_{\rm rxx1}(s) = \begin{bmatrix} 1 \\ \frac{1}{5.9s} + \frac{1}{1}e^{-3s} & 0 \\ 0 & \frac{1}{5.8s} + \frac{1}{1}e^{-3s} \end{bmatrix}.\tag{34}$$

The retuned controller using the lemma 1.2 is given by:

$$\overline{\mathbf{C}}\_{\text{ex1}}(\mathbf{s}) = \begin{bmatrix} 0.172 & + \frac{0.0299}{\mathfrak{s}} & + & 0.204\mathfrak{s} & -0.0666 - \frac{0.016}{\mathfrak{s}} & + & 0.046\mathfrak{s} \\ -0.009 & + & \frac{0.0097}{\mathfrak{s}} - 0.0548\mathfrak{s} & -0.0512 - \frac{0.012}{\mathfrak{s}} & + & 0.106\mathfrak{s} \end{bmatrix} . \tag{35}$$

Note in **Figure 1** that the retuned closed loop is more robust than the initial one. This occurs because the peak of maximum singular value curve of the initial sensitivity function (**S**ð Þ*s* ), Eq. (5), is above the peak of the curve corresponding to the reprojected.

In **Figure 2**, the step responses of the initial and retuned closed loop are shown. The retuned closed-loop response approaches the desired. In addition, coupling reduction is observed when the retuned controller is used. This occurs because the coupling is compensated by the off-diagonal elements of the controller matrix. Consequently, the variation of the control signal increased, as can be observed in **Figure 3**.

In this example, the decentralized PID controller was retuned so that the new closed-loop close to a diagonal reference model. As result of the proposed method, a centralized PID controller was obtained.

#### **4.2 Example 2**

Consider the integrating process of the distillation column process [10]:

$$\mathbf{G}\_{\rm ex2}(s) = \begin{bmatrix} \frac{3.04}{s} & \frac{-278.28}{s(30s^2 + 1)(6s^2 + 1)}\\ \frac{0.052}{s} & \frac{319.47}{s(30s^2 + 1)(6s^2 + 1)} \end{bmatrix} \tag{36}$$

and the decentralized PI controller:

**Figure 1.** *Maximum singular value of sensitivity function (***S**ð Þ*s ) - example 1.*

**Figure 2.** *Closed-loop step response - example 1.*

$$\mathbf{C}\_{\text{cc2}}(\mathbf{s}) = \begin{bmatrix} \mathbf{16.181} & + & \frac{202.265}{\mathfrak{s}} & & & \mathbf{0} \\\\ \mathbf{0} & & \mathbf{23.614} & + & \frac{64.926}{\mathfrak{s}} & + & 2.147 \mathfrak{s} \end{bmatrix} . \tag{37}$$

The reference model is given by:

$$\mathbf{T}\_{\text{rx2}}(\mathbf{s}) = \begin{bmatrix} \mathbf{1} & & & \mathbf{0} \\ \hline \mathbf{0.05s} & + & \mathbf{1} & & \\ & \mathbf{0} & & \mathbf{1} \\ & & \mathbf{0.2s} + & \mathbf{1} \end{bmatrix} . \tag{38}$$

*MIMO PID Control Retuning Using the Closed-Loop Frequency Response DOI: http://dx.doi.org/10.5772/intechopen.105015*

The initial controller is decentralized and has PI-type and PID-type elements. Thus, the lemma 1.2 was used to retune the PID and for the lemma 1.3 was used to the other elements. The integral gain obtained for elements *C*11, *C*12, and *C*<sup>22</sup> was approximately zero. The retuned controller is given by:

$$
\overline{\mathbf{C}}\_{\text{ex2}}(s) = \begin{bmatrix} 5.73 & 1.25 \\ -1.65 & + \end{bmatrix}\_{\text{st}} \begin{array}{c} 1.25 \\ 2.4 \end{array} + \begin{array}{c} 1.25 \\ 0.49s \end{array} \tag{39}
$$

In **Figure 4**, the maximum singular value curve of the sensitivity function (**S**ð Þ*s* ) is presented. Observe that the curve peak of the proposed loop is smaller when compared with the initial loop.

The closed-loop step response is show in **Figure 5**. The retuned loop outputs were close to the desired and slower than the initial loop. Consequently, the overshoot was reduced. Also, the control signal became smoother as shown in **Figure 6**.

**Figure 4.** *Maximum singular value of sensitivity function (***S**ð Þ*s ) - example 2.*

**Figure 5.** *Closed-loop step response - example 2.*

**Figure 6.** *Control signal - example 2.*

#### **4.3 Example 3**

Consider the drum boiler integrating process [11]:

$$\mathbf{G}\_{\rm ex3}(s) = \begin{bmatrix} \frac{0.349}{38.19s + 1} & \frac{-0.1587}{203.9s + 1} \\ \frac{-0.0059}{s} \frac{1 - 401.2s}{1 + 21.15s} & \frac{0.01033}{s} \end{bmatrix} \tag{40}$$

and the decentralized PI controller:

*MIMO PID Control Retuning Using the Closed-Loop Frequency Response DOI: http://dx.doi.org/10.5772/intechopen.105015*

$$\mathbf{C}\_{\rm ex3}(s) = \begin{bmatrix} 2.5 & + & \frac{0.05}{s} & & & 0\\ & & & & & \\ & \mathbf{0} & & & \mathbf{1.25} & + & \frac{0.025}{s} \end{bmatrix}.\tag{41}$$

Note the presence of a zero in the right-half plane (RHP) that affects output 2. A zero in the RHP must appear in the element of the reference model transfer matrix referring to output 2. The reference model is given by:

$$\mathbf{T}\_{\text{rx3}}(\mathbf{s}) = \begin{bmatrix} \mathbf{1} & & & \mathbf{0} \\ \overline{\mathbf{3}\mathbf{0}\mathbf{s}} + \mathbf{1} & & & \\ \mathbf{0} & & \overline{\mathbf{1} - \mathbf{3}\mathbf{2}} & \\ \mathbf{0} & & \overline{\mathbf{7}\mathbf{0}\mathbf{s}} + \mathbf{1} \end{bmatrix}.\tag{42}$$

In this case, as the process is integrative, the retune was performed using lemma 1.3 so that the retuned controller was of the PI type:

$$\overline{\mathbf{C}}\_{\text{cc3}}(\varsigma) = \begin{bmatrix} 1.89 & + & \frac{0.0751}{\mathfrak{s}} & 0.198 & + & \frac{0.011}{\mathfrak{s}} \\ -17.3 & + & \frac{0.201}{\mathfrak{s}} & 2.69 & + & \frac{0.117}{\mathfrak{s}} \end{bmatrix}. \tag{43}$$

In **Figure 7**, the maximum singular value curve of the sensitivity function (**S**ð Þ*s* ) is presented. Note that the retuned loop is significantly more robust than the initial.

The closed-loop step response is show in **Figure 8**. With the retuned controller, the interaction of input 1 with output 2 has been reduced. However, the variation of the control signal increased, as shown in **Figure 9**.

In this example, the centralized PI controller has been retuned from closed loop with a decentralized PI controller. It was possible to obtain a decoupled closed loop with a greater gain margin.

**Figure 7.** *Maximum singular value of sensitivity function (***S**ð Þ*s ) - example 3.*

**Figure 8.** *Closed-loop step response - example 3.*

**Figure 9.** *Control signal - example 3.*

#### **5. Conclusions**

In this article, a PID controller retuning method for multivariable processes was presented. This method is an extension of the one presented in [7]. The controller is retuned so that the closed loop approximates the desired one. The method is based on closed-loop frequency domain data. Knowledge of the parametric model of the process is not necessary.

Controller gain increments are computed from the closed-loop reference model, the initial controller, and closed-loop frequency domain data. The initial controller can be centralized or decentralized, PI or PID type.

In the simulation examples, it can be seen that the retuned controller is centralized. In example 2, it is shown that the P/PD controller can be obtained from the initial PI/ PID controller. In this case, the integral gain of some controllers was approximately

zero. With the redesign, it was possible to reduce the coupling between the loops and improve the robustness properties of the closed loop.

As future work, there is the application of the method to unstable processes and a methodology to define the reference model.

#### **Acknowledgements**

This work is supported by CNPq (Conselho Nacional de Desenvolvimento Científico e Tecnológico) and COPELE (Coord. de Pós-graduação em Eng. Elétrica da UFCG).

#### **Author details**

Anna Paula V. de A. Aguiar†, George Acioli Júnior\*† and Péricles R. Barros† Electrical Engineering Department—DEE, Federal University of Campina Grande—UFCG, Campina Grande, Brazil

\*Address all correspondence to: georgeacioli@dee.ufcg.edu.br

† These authors contributed equally.

© 2022 The Author(s). Licensee IntechOpen. This chapter is distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

### **References**

[1] Nisi K, Nagaraj B, Agalya A. Tuning of a PID controller using evolutionary multi objective optimization methodologies and application to the pulp and paper industry. International Journal of Machine Learning and Cybernetics. 2019;**10**:2015-2025

[2] Gao X, Yang F, Shang C, Huang D. A novel data-driven method for simultaneous performance assessment and retuning of PID controllers. Industrial & Engineering Chemistry Research. 2017;**56**:2127-2139

[3] Yu S, Li X. Proportional–integral– derivative controller performance assessment and retuning based on general process response data. ACS Omega. 2021;**6**:10207-10223

[4] Moreira LJDS, Acioli Júnior G, Barros PR. IMC PI control loops frequency and time domains performance assessment and retuning. IFAC-PapersOnLine. 2018;**51**:148-153

[5] Moreira LJDS, Acioli Júnior G, Barros PR. Closed-loop frequency datadriven PID retuning. In: Proceedings of the 2018 IEEE Conference on Control Technology and Applications (CCTA). Copenhagen, Denmark: IEEE; 2018. pp. 1352-1357

[6] Moreira LJDS, Aguiar APVDA, Acioli Jãnior G, Barros PR. Time and frequency data-driven PID retuning applied in MIMO process. IFAC-PapersOnLine. 2021;**54**:469-474 16th IFAC Symposium on Advanced Control of Chemical Processes ADCHEM 2021

[7] Aguiar APVDA, Acioli Júnior G, Barros PR, Perkusich A. A closed-loop frequency domain PI retuning technique for multivariable systems. IFAC-PapersOnLine. 2021;**54**:463-468 16th

IFAC Symposium on Advanced Control of Chemical Processes ADCHEM. 2021

[8] Åström KJ, Hägglund T, Astrom KJ. Advanced PID Control. Vol. 461. NC, United States of America: ISA-The Instrumentation, Systems, and Automation Society Research Triangle Park; 2006

[9] Wood R, Berry M. Terminal composition control of a binary distillation column. Chemical Engineering Science. 1973;**28**:1707-1717

[10] Hu W, Cai WJ, Xiao G. Decentralized control system design for MIMO processes with integrators/ differentiators. Industrial & Engineering Chemistry Research. 2010;**49**: 12521-12528

[11] Balko P, Rosinová D. Robust decentralized control of nonlinear drum boiler. IFAC-PapersOnLine. 2015;**48**: 432-437 8th IFAC Symposium on Robust Control Design ROCOND 2015

## **Chapter 5** PID Control for Nonlinear Processes

### *Taieb Adel, Kanzari Bilel and Chaari Abdelkader*

#### **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.

**Keywords:** PID controller, Takagi-Sugeno, nonlinear system, MIMO system, particle swarm optimization-cuckoo search

#### **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:

$$\mathbf{y}\_{i}(k+1) = f\_{i}(\mathbf{x}(k)), i = \mathbf{1}, 2, \dots, n \tag{1}$$

with the regression vector given by:

$$\mathbf{x}\_{i}(k) = \begin{bmatrix} y\_{i}(k)\_{0}^{n}, \ u\_{1}(k - d\_{1i})\_{0}^{n}, \ \dots, \ u\_{n}(k - d\_{ni})\_{0}^{n} \end{bmatrix} \tag{2}$$

here *k* denotes the discrete time sample, *n* is an integer relative to the order of the system, and *dij* 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:

$$\mathbf{y}(k+1) = A\mathbf{y}(k) + Bu(k) + a \tag{3}$$

where *α* is a constant. A fuzzy TS model makes an attempt to approximate the unknown function *f*ð Þ*:* . 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: \operatorname{if } \mathfrak{x}(k) \text{ is } \Omega \text{ then } \operatorname{y}^l(k+1) = A \operatorname{y}(k) + B \operatorname{ul}(k) + a\_l \mathfrak{l} = \mathfrak{1}, \dots \mathfrak{r} \tag{4}$$

where Ω*<sup>l</sup>* is the antecedent fuzzy set of the *l th* rule, *Al* <sup>¼</sup> *Al*<sup>1</sup> ½ � , *Al*2, … , *Aln* and *Bl* ¼ *Bl*<sup>1</sup> ½ � , *Bl*2, … , *Bln* are the vectors of the two polynomials *Al* and *Bl*, and *r* is the number of rules.

Rule (4) can then be written as:

$$\begin{aligned} \mathcal{R}\_l: &\text{if } \propto\_1(k) \text{ is } \Omega\_{l1} \text{ and } \propto\_p(k) \text{ is } \Omega\_{lp} \text{ then} \\ \mathcal{Y}^l(k+1) &= A \mathcal{Y}(k) + B \mu(k) + a\_l \mathcal{l} = \mathbf{1}, \dots, r \end{aligned} \tag{5}$$

with *<sup>p</sup>* <sup>¼</sup> <sup>2</sup>*n*<sup>2</sup> <sup>þ</sup> <sup>1</sup>

The output of the TS model is then evaluated by:

$$y(k+1) = \frac{\sum\_{i=1}^{r} \mu\_i(\mathbf{x}(k)) y^i(k+1)}{\sum\_{i=1}^{r} \mu\_i(\mathbf{x}(k))}\tag{6}$$

By asking:

$$\Phi\_j(\mathbf{x}, c\_i, \sigma\_i) = \frac{\mu\_j(\mathbf{x}(k))}{\sum\_{i=1}^r \mu\_i(\mathbf{x}(k))} \tag{7}$$

with Φ*<sup>j</sup> x*, *ci* ð Þ , *σ<sup>i</sup>* is the validation function of the Gaussian function having as parameters the centers *ci* and the variances *σi*.

$$\mu\_i(\mathbf{x}(k)) = \exp\left(\frac{-(\mathbf{x}\_i - c\_{i1})^2}{2\sigma\_{i1}^2}\right) \dots \exp\left(\frac{-\left(\mathbf{x}\_i - c\_{ip}\right)^2}{2\sigma\_{ip}^2}\right) \tag{8}$$

The formula (6) becomes

$$\boldsymbol{y}(k+1) = \sum\_{i=1}^{r} \boldsymbol{y}^{i}(k+1)\boldsymbol{\Phi}\_{j}(\mathbf{x}, \mathbf{c}\_{i}, \sigma\_{i}) \tag{9}$$

A Gaussian function's center and variance (*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 = \begin{bmatrix} \cdots \\ \varkappa(k) \\ \cdots \\ \varkappa(N-1) \end{bmatrix}, Y = \begin{bmatrix} \cdots \\ \varkappa(k) \\ \cdots \\ \cdots \\ \varkappa(N-1) \end{bmatrix} \\ Z^T = [X \ Y] \tag{10}$$

where *N* is the number of observations.

The data set *Z* will be divided into *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 *<sup>U</sup>* <sup>¼</sup> *<sup>μ</sup>ik* ½ �*Nc*�*<sup>N</sup>* with *<sup>μ</sup>ik* <sup>∈</sup>½ � 0 1 represents the degree to which observation *xk* belongs to group *i*. Each group is characterized by a center *ci*, where *C* ¼ *c*1, … , *CNc* ½ � is the vector of the centers. And a covariance matrix *F* ¼ *F*1, … , *FNc* ½ � describes the variance of the data in the group *Fi*.

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

$$\Omega\_{\vec{\eta}}(\mathbf{x}\_{\vec{\jmath}}(k)) = \exp\left(\frac{-\mathbf{1}}{2}\frac{(\mathbf{x}\_{\vec{\jmath}}) - c\_{\vec{\imath}\vec{\jmath}}^2}{\sigma\_{\vec{\imath}\vec{\jmath}}^2}\right) \tag{11}$$

The parameters of the consequents *θ<sup>i</sup>* ¼ *Ai*, *Bi*, *Ci* ½ � are estimated separately by the recursive least squares algorithm by minimizing the following objective function:

$$\min \, \_{\theta\_i} \frac{1}{N} [Y - \xi \theta\_i] \, ^T Q\_i [Y - \xi \theta\_i] \tag{12}$$

with *ξ* ¼ ½ � *X* 1 is the regression matrix augmented by adding a unit column vector and *Qi* is a matrix containing the values of the validity functions Φ*<sup>i</sup>* of the *i th* local linear model of each data group:

$$Q = \begin{bmatrix} \Phi\_{-}(\mathfrak{x}(1), \ c\_{i}, \ \sigma\_{i}) & 0 & \dots & 0 \\ 0 & \Phi\_{-}(\mathfrak{x}(2), \ c\_{i}, \ \sigma\_{i}) & \dots & 0 \\ & \vdots & \vdots & \vdots & \vdots \\ & 0 & 0 & \dots & \Phi\_{-}(\mathfrak{x}(N), \ c\_{i}, \ \sigma\_{i}) \end{bmatrix},\tag{13}$$

The least squares estimate of the consequent parameters, (*θ<sup>i</sup>* ¼ *ai*<sup>1</sup> … *aim*, *bi*<sup>1</sup> … *bim*, *ci* ½ �) is given by:

$$\theta\_i = \left[\xi^T Q\_i \xi\right]^{-1} \xi^T Q\_i Y \tag{14}$$

#### **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]:

$$\phi\_i(k+1) = \left[\mu\_i \wp(k-1), \dots \mu\_i \wp(k-m), \mu\_i u(k-1), \dots \mu\_i u(k-m), 1\right]$$

$$\theta\_i = \left[a\_{i1}\dots a\_{im}, b\_{i1}\dots b\_{im}, c\_i\right]$$

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

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

Obtaining the estimated parameter vector ^*θ* from the system parameter vector *θ* is given by the following least squares algorithm:

$$\hat{\theta}(k) = \hat{\theta}(k-1) + \frac{p(k-1)\overline{\phi}(k-1)\overline{e}\_1(k-1)}{\lambda(k) + \overline{\phi}^T(k-1)p(k-1)\overline{\phi}(k-1)}\tag{15}$$

With:

$$\bullet \ \overline{e}\_1(k) = \overline{\jmath}(k) - \overline{\phi}^T(k-1)\theta(k-1)$$

$$\bullet \ p(k) = p(k-1) - \frac{p(k-1)\overline{\phi}(k)\overline{\phi}^\flat(k)p^\flat(k-1)}{\overline{\phi}^\flat(k) + \overline{\phi}^\Gamma(k-1)p^\flat(k-1)\overline{\phi}(k-1)}$$

$$\bullet \text{ } \lambda(k) = \lambda\_0 \lambda(k-1) + (1-\lambda\_0), \text{ with } \lambda\_0 = 0.95, p(0) = \varepsilon I, 0 < \varepsilon < \infty$$

• *y k*ð Þ¼ *y k*ð Þ *<sup>η</sup>*ð Þ *<sup>k</sup>*�<sup>1</sup> , *<sup>ϕ</sup>*ð Þ¼ *<sup>k</sup> <sup>ϕ</sup>*ð Þ *<sup>k</sup>*�<sup>1</sup> *<sup>η</sup>*ð Þ *<sup>k</sup>*�<sup>1</sup> , *<sup>η</sup>*ð Þ¼ *<sup>k</sup> max* 1,∥*ϕ*ð Þ*<sup>k</sup>* <sup>∥</sup>, *<sup>η</sup>*ð Þ*<sup>k</sup>* 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:

$$
\mu\_p(t) = K\_p e(t) \tag{16}
$$

where *ε*ðÞ¼ *t yc*ðÞ�*t y t*ð Þ represents the error signal with *yc*ð Þ*t* the set point and *y t*ð Þ the output of the system to be controlled.

• Integral action (I):

$$u\_I(t) = \frac{K\_p}{T\_i} \int\_0^t \varepsilon(\tau) \tag{17}$$

where *Ti* denotes the integration constant.

• Derivative action (D) is as follows:

$$
\mu\_D(t) = K\_p T\_d \frac{d\varepsilon(t)}{dt} \tag{18}
$$

where *Td* denotes the derivative constant. Hence, the command signal is as follows:

$$u(t) = u\_p(t) + u\_i(t) + u\_d(t)\tag{19}$$

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

$$R(p) = K\_p \left( 1 + \frac{1}{T\_i p} + t\_d p \right) \tag{20}$$

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}} \text{ with } N\_f \ge 5 \tag{21}$$

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

$$R(p) = K\_p \left( 1 + \frac{1}{T\_i p} + \frac{T\_d p}{1 + \frac{T\_d}{N\_f p}} p \right) \tag{22}$$

#### **3.2 Discrete version of the PID controller**

To obtain the discrete version, we replace:

• the integral by a discrete sum:

$$\mu\_I(t) = \frac{K\_p}{T\_i} \int\_0^t \varepsilon(\tau) d\tau \to \mu\_I(k) = \frac{K\_p T\_\epsilon}{T\_i} \sum\_{i=0}^k \varepsilon(k) \tag{23}$$

where *Te* represents the sampling period and k is an integer.

• The derivative by a difference:

$$
\mu\_D(t) = K\_p T\_d \frac{d\epsilon(t)}{dt} \to \mu\_D(k) = K\_p T\_d \frac{\epsilon(k) - \epsilon(k-1)}{T\_\epsilon} \tag{24}
$$

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

$$u\_{\rm PID}(k) = u\_p(k) + u\_i(k) + u\_d(k) = K\_p \left( \varepsilon(k) + \frac{T\_\varepsilon}{T\_i} \sum\_{i=0}^k \varepsilon(k) + \frac{T\_d}{T\_\varepsilon} \left( \varepsilon(k) - \varepsilon(k-1) \right) \right) \tag{25}$$

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\_\epsilon}{T\_i} \frac{1}{1 - z^{-1}} = \frac{T\_\epsilon}{T\_i} \frac{z}{z - 1}$$

• For the derived term:

$$T\_d p \to \frac{T\_d}{T\_\epsilon} \left( 1 - z^{-1} \right) = \frac{T\_d}{T\_\epsilon} \frac{z - 1}{z}$$

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

$$R\left(\mathbf{z}^{-1}\right) = K\_p \left(\mathbf{1} + \frac{T\_\epsilon}{T\_i} \frac{\mathbf{1}}{\mathbf{1} - \mathbf{z}^{-1}} + \frac{N\_f T\_d}{\left(N\_f T\_\epsilon + T\_d\right)} \frac{\mathbf{1} - \mathbf{z}^{-1}}{\mathbf{1} - \frac{T\_d}{N\_f T\_\epsilon + T\_d} \mathbf{z}^{-1}}\right)^{-1}$$

#### **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\_{\rm Par}\left(\mathbf{z}^{-1}\right) = K\_p + K\_i \frac{\mathbf{1}}{\mathbf{1} - \mathbf{z}^{-1}} + K\_d \frac{\mathbf{1} - \mathbf{z}^{-1}}{\mathbf{1} - \mathbf{z}^{-1}} \tag{26}$$

with *Ki* ¼ *Kp Te Ti* , *Kd* ¼ *Kp Nf Td Nf Te*þ*Td* and *<sup>ι</sup>* <sup>¼</sup> *Td Nf Te*þ*Td*

Note 1: if there is no filter on the derivative, that is to say that *Nf* ! ∞, then we will have

$$
\iota \to 0 \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\_{\rm Min}(z^{-1}) = k\_p \left( \mathbf{1} + k\_i \frac{\mathbf{1}}{\mathbf{1} - z^{-1}} + k\_d \frac{\mathbf{1} - z^{-1}}{\mathbf{1} - u z^{-1}} \right) \tag{27}$$

with *kp* <sup>¼</sup> *Kp*, *ki* <sup>¼</sup> *Te Ti* and *dd* <sup>¼</sup> *Nf Td Nf Te*þ*Td*

#### *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\_{Ser}(z^{-1}) = \frac{r\_0 + r\_1 z^{-1} + r\_2 z^{-2}}{(1 - z^{-1})(1 + s\_1 z^{-1})} \tag{28}$$

with

$$r\_0 = K\_p \left( \mathbf{1} + \frac{T\_\epsilon}{T\_i} - N\_{f^\sharp 1} \right)$$

$$r\_1 = K\_p \left( s\_1 \left( \mathbf{1} + \frac{T\_\epsilon}{T\_i} + 2N\_f \right) - \mathbf{1} \right)$$

$$r\_2 = -K\_p s\_1 \left( \mathbf{1} + N\_f \right)$$

$$s\_1 = \frac{-T\_d}{N\_f T\_\epsilon + td}$$

Note 2: the continuous equivalent does not always exist! The existing condition requires that:

$$-1 \le s\_1 \le 1 \text{ that is tooby } \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**):

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

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

The transfer function in a closed loop is given by:

$$H\_{Bf} \left( z^{-1} \right) = \frac{B(z^{-1})}{A(z^{-1})S(z^{-1}) + B(z^{-1})R(z^{-1})} \tag{29}$$

again:

$$H\_{\mathcal{B}'}\left(z^{-1}\right) = \frac{B(z^{-1})}{P'(z^{-1})}\tag{30}$$

with:

**Figure 1.** *Standard form of RST.*

*PID Control for Nonlinear Processes DOI: http://dx.doi.org/10.5772/intechopen.106820*

**Figure 2.** *RST structure.*

$$\begin{aligned} e(k) &= r(k) - \jmath(k) \\ B\big(z^{-1}\big) &= b\_1 z^{-1} + b\_2 z^{-2} \\ A\big(z^{-1}\big) &= 1 + a\_1 z^{-1} + a\_2 z^{-2} \\ R\big(z^{-1}\big) &= r\_0 + r\_1 z^{-1} + r\_2 z^{-2} \\ S\big(z^{-1}\big) &= \left(1 - z^{-1}\right)\left(1 + s\_1 z^{-1}\right) \end{aligned}$$

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

$$S(z^{-1})\mu(z^{-1}) = R(z^{-1})e(z^{-1})\tag{31}$$

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

$$u(k) = (\mathbf{1} - \varepsilon\_1)u(k - \mathbf{1}) + \varepsilon\_1 u(k - \mathbf{2}) + r\_0 \varepsilon(k) + r\_1 \varepsilon(k - \mathbf{1}) + r\_2 \varepsilon(k - \mathbf{2})\tag{32}$$

#### **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 :


#### **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 *Hbf* possesses the following features:


For the choice of the sampling period *Te*, 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*<sup>0</sup> 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\_\epsilon < 1,5$$

$$0,7 < m < 1$$

Because, the two settings *w*<sup>0</sup> and *m* directly affect the behavior of the system, rise time *tm*, settling time *ts*, and the maximal overshoot *D*1. The sampling period is smaller than *m* is small.

Given the desired characteristic polynomial system:

$$P(\mathbf{z}^{-1}) = \mathbf{1} + p\_1 \mathbf{z}^{-1} + p\_2 \mathbf{z}^{-2} \tag{33}$$

where *p*<sup>1</sup> and *p*<sup>2</sup> depend on the imposed specifications (rise time, overshoot, and damping), where:

$$\begin{aligned} p\_1 &= -2e^{-m w\_0 T\_\epsilon} \cos \left( w\_0 T\_\epsilon \text{sqrt} \left( 1 - m^2 \right) \right), \\ p\_2 &= e^{-m w\_0 T\_\epsilon} \end{aligned}$$

where equation Diophantine is written as:

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

such as *A z*�<sup>1</sup> ð Þ, *B z*�<sup>1</sup> ð Þ, and *P z*�<sup>1</sup> ð Þ are the known polynomials, and *S z*�<sup>1</sup> ð Þ and *R z*�<sup>1</sup> ð Þ are the wanted polynomials. We recall that:

$$\frac{B(z^{-1})}{A(z^{-1})} = \frac{b - \mathbf{1}z^{-1} + b\_2 z^{-2}}{\mathbf{1} + a\_1 z^{-1} + a\_2 z^{-2}}\tag{35}$$

whence (34) becomes

$$\begin{split} P(\boldsymbol{z}^{-1}) &= \left(\mathbf{1} + a\_1 \mathbf{z}^{-1} + a\_2 \mathbf{z}^{-2}\right) \left(\mathbf{1} - \mathbf{z}^{-1}\right) \left(\mathbf{1} + s\_1 \mathbf{z}^{-1}\right) \\ &+ \left(b - \mathbf{1} \mathbf{z}^{-1} + b\_2 \mathbf{z}^{-2}\right) \left(r\_0 + r\_1 \mathbf{z}^{-1} + r\_2 \mathbf{z}^{-2}\right) \end{split} \tag{36}$$

By identification, we can write

$$\begin{cases} a\_1 - \mathbf{1} + \mathbf{s}\_1 + r\_0 b\_1 = p\_1 \\ (a\_1 - \mathbf{1})\mathbf{s}\_1 - a\_1 + a\_2 + r\_0 b\_2 + r\_1 b\_1 = p\_2 \\ (a\_2 - a\_1)\mathbf{s}\_1 - a\_2 + r\_1 b\_2 + r\_2 b\_1 = \mathbf{0} \\ -a\_2 \mathbf{s}\_1 + r\_2 b\_2 = \mathbf{0} \end{cases} \tag{37}$$

where *<sup>p</sup>*<sup>1</sup> and *<sup>p</sup>*<sup>2</sup> are the coefficients of the characteristic polynomial *P z*�<sup>1</sup> ð Þ:

$$\begin{cases} +s\_1 + r\_0 b\_1 = p\_1 - a\_1 + \mathbf{1} \\ (a\_1 - \mathbf{1})s\_1 + r\_0 b\_2 + r\_1 b\_1 = p\_2 - a\_2 + a\_1 \\ (a\_2 - a\_1)s\_1 + r\_1 b\_2 + r\_2 b\_1 = a\_2 \\ -a\_2 s\_1 + r\_2 b\_2 = \mathbf{0} \end{cases} \tag{38}$$

We can write equations in a matrix form, we then have *MQ* ¼ *P* with:

$$M = \begin{bmatrix} 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{bmatrix} \\ Q = \begin{bmatrix} s\_1 \\ r\_0 \\ r\_1 \\ r\_2 \end{bmatrix} \\ P = \begin{bmatrix} p\_1 - a\_1 + 1 \\ p\_2 - a\_2 + a\_1 \\ a\_2 \\ 0 \end{bmatrix},\tag{39}$$

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\tag{40}$$

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{aligned} K\_p &= \frac{\left(r\_0 s\_1 - r\_1 - (2 + s\_1)r\_2\right)}{\left(1 + s\_1\right)^2} \\ T\_i &= T\_\epsilon \frac{k\_p \left(1 + s\_1\right)}{r\_0 + r\_1 + r\_2} \\ T\_d &= T\_\epsilon \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\_\epsilon}{1 + s\_1} \end{aligned} \tag{41}$$

#### **4.2 PID controller based on PSOCS**

#### *4.2.1 PSO*

This approach was introduced in 1995 [11]. The swarm *Xp*,*p* ¼ 1,2, … ,*Np* which has *Np* particles is considered in the standard PSO. In D-dimensional space, the *pth* particle *Xp* ¼ *xp*1, … ,*xpD* is a potential solution to the researched problem. The velocity of the *pth* particle is *Vp* <sup>¼</sup> *vp*1, … ,*vpD*. For the *<sup>p</sup>th* particle, the best position in the *<sup>p</sup>th* step is expressed as *pbestp*ðÞ¼ *t pbestp*1ð Þ*t* , … ,*pbestpD*ð Þ*t* . Meanwhile, the best position of the entire swarm is defined as *Gbest t*ðÞ¼ *Gbest*1ð Þ*t* , … ,*GbestD*ð Þ*t* . Thus, at the ð Þ *<sup>t</sup>* <sup>þ</sup> <sup>1</sup> *th* step, the new position of the *pth* particle is calculated as follows:

$$\boldsymbol{\varkappa}\_{\vec{p}\vec{j}}(t+\mathbf{1}) = \boldsymbol{\varkappa}\_{\vec{p}\vec{j}}(t) + \boldsymbol{\nu}\_{\vec{p}\vec{j}}(t+\mathbf{1})\tag{42}$$

$$v\_{\overline{p}\dagger}(t+1) = \alpha v\_{\overline{p}\dagger}(t) + c\_1 \times rand\_1 \times \left[ \text{Pbest}\_{\overline{p}\dagger}(t) - \text{x}\_{\overline{p}\dagger}(t) \right] + c\_2 \times rand\_2 \times \left[ \text{Gbest}\_p(t) - \text{x}\_{\overline{p}\dagger}(t) \right] \tag{43}$$

where *xpj*ð Þ*<sup>t</sup>* and *vpj*ð Þ*<sup>t</sup>* represent the position and velocity of the *<sup>p</sup>th* particle with respect to the *j th* dimension, *ω* is the inertia weight factor, *c*<sup>1</sup> and *c*<sup>2</sup> are two positive constants called acceleration coefficients, and *rand*<sup>1</sup> and *rand*<sup>2</sup> are two uniformly distributed random values in the range 0, 1 ½ �.

#### *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 *pa* of the *np* 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 x*ð Þ, *x* ¼ *x*1, … ,*xD* 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 *Np* host nests *xp*,*p* ¼ 1, … ,*Np*. Each nest *xp* represents a potential answer to the optimization problem that needs to be solved. The representations are the same as the particle *xp* in PSO. Finding the new population of nests *xp*ð Þ *t* þ 1 is one of the CS's crucial phases. Furthermore, Lévy flying is used to obtain the new nests:

$$
\varkappa\_{pj}(t+1) = \varkappa\_{pj}(t) + a \oplus Levy(\lambda) \tag{44}
$$

where *α* is the step size, ⊕ represents the entry-wise multiplication operation, and *λ* 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(\mathfrak{s}, \mathfrak{y}, \mathfrak{\nu}) = \begin{cases} \sqrt{\frac{\chi}{2\Pi}} \exp\left[ -\frac{\chi}{2(\mathfrak{s}-\nu)} \frac{\mathbf{1}}{(\mathfrak{s}-\nu)^{3/2}} \right], \mathfrak{0} < \mathfrak{\nu} < \mathfrak{s} < \mathbf{1} \\ \mathbf{0} & \text{Otherwise} \end{cases} \tag{45}$$

where *ν*>0 is a minimum step and *γ* is a scale parameter. In Mantegna's algorithm, the step length *s* can be calculated by:

$$s = \frac{u}{|v|^{1/\beta}}\tag{46}$$

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

$$\begin{cases} \quad u \sim N(\mathbf{0}, \ \sigma\_u^2), \\ \quad v \sim N(\mathbf{0}, \ \sigma\_v^2), \\ \quad \sigma\_u = \frac{\tau (1+\beta) (\sin \left( \pi \beta / 2 \right))^{1/\beta}}{\tau [(1+\beta) / 2] \beta^{2(1+\beta)/2}} \\ \qquad \sigma\_v = \mathbf{1} \end{cases} \tag{47}$$

where *<sup>τ</sup>*ð Þ*<sup>z</sup>* is the Gamma function *<sup>τ</sup>*ð Þ¼ *<sup>z</sup>* <sup>Ð</sup> <sup>∞</sup> 0 *t <sup>z</sup>*�1*e*�*<sup>t</sup> dt*. In the case when *z* ¼ *np* is an integer, *τ np* � � <sup>¼</sup> *np* � <sup>1</sup> � �!

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

Although the entry-wise product ⊕ 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 *rand*<sup>1</sup> and *rand*2, when we set the parameters *ω*, *c*1, and *c*<sup>2</sup> 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:


*while t*< *Maxgeneratio* or other stop criterion.

• Step 4. The fitness function of each particle *f Xp* � � <sup>¼</sup> *<sup>f</sup> <sup>p</sup>* is evaluated. The best position *Pbest p*ð Þð Þ*t* and the best position of the whole swarm *Gbest t*ð Þ are found (similar to PSO):

$$Pbest\_p(t+1) = \begin{pmatrix} \propto\_p(t+1) \text{ if } f\left(\mathbf{x}\_p(t+1)\right) \lhd f\left(\mathbf{x}\_p(t)\right) \\\quad Pbest\_p \quad \text{2cm} \\\ \text{Gbest}(t) = \min Pbest\_1(t), Pbest\_2(t), \dots, Pbest\_{N\_p}(t) \end{pmatrix} \tag{48}$$


$$
\rho = \alpha\_{\text{max}} - \frac{\alpha\_{\text{max}} - \alpha\_{\text{min}}}{T} \times t \tag{49}
$$

The velocity *vp* and position *Xp* of each particle are updated according to Eqs. (42) and (43). Parameters *rand*<sup>1</sup> and *rand*2vary with Lévy flight pattern following (45) to (47),which is different from the former PSO (hybrid of PSO and CS):

$$\begin{split} \upsilon\_{\vec{p}\vec{j}}(\mathfrak{t}+\mathfrak{1}) &= \alpha \upsilon\_{\vec{p}\vec{j}}(\mathfrak{t}) + (c\_1 \oplus \mathit{Levy}(\lambda)) \times \left[ \mathit{Pbest}\_{\vec{p}\vec{j}}(\mathfrak{t}) - \mathfrak{x}\_{\vec{p}\vec{j}}(\mathfrak{t}) \right] \\ &+ (c\_2 \oplus \mathit{Levy}(\lambda)) \times \left[ \mathit{Gbest}\_{\vec{j}}(\mathfrak{t}) - \mathfrak{x}\_{\vec{p}\vec{j}}(\mathfrak{t}) \right] \end{split} \tag{50}$$

$$
\boldsymbol{\kappa}\_{p\circ}(t+\mathbf{1}) = \boldsymbol{\kappa}\_{p\circ}(t) + \boldsymbol{\upsilon}\_{p\circ}(t+\mathbf{1})\tag{51}
$$


The problem's dimension determines the size of the swarm *Np*. 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, *<sup>τ</sup>*ð Þ¼ *<sup>z</sup>* <sup>Ð</sup> <sup>∞</sup> 0 *t <sup>z</sup>*�<sup>1</sup>*e*�*<sup>t</sup> dt* is realized by Gamma function *τ* ¼ *γ*ð Þ*z* using MATLAB. The inertia weight changes from *ω min* to *ω max* .

The friction *pa* 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:

$$J = \left(\mathbf{1} - \mathbf{e}^{(-\beta)}\right) \left(D + E\_{\rm s}\right) + \mathbf{e}^{-\beta} \left(\mathbf{t}\_{\rm s} - \mathbf{t}\_{m}\right)$$

$$\begin{aligned} \operatorname{minimize} I(G)\_{G = r\_{0}, r\_{1}, r\_{2}, \mathbf{s}\_{1}} \\ G\_{\rm min} < G < G\_{\rm max} \end{aligned} \tag{52}$$

Within the initial search space (*R min* ,*R max* ), the controller parameters (*r*0,*r*1,*r*2,*s*1) will be optimized using the optimization algorithms. The graphic below shows the

*PID Control for Nonlinear Processes DOI: http://dx.doi.org/10.5772/intechopen.106820*

**Figure 3.** *Block diagram of a PID-PSOCS algorithm.*

regulator block diagram for a nonlinear system using a PID controller and the PSOCS optimization algorithm (**Figure 3**):

The PID-PSOCS algorithm is summarized by the following steps

#### **Algorithm 1 :** Proposed AFMPC-PSOCS

	- 1. Step 1. Set parameters *Np*, *c*1, *c*2. Initialize the position and velocity of each particle.
		- Repeat for h1=1,2,...
	- 2. Step 2. For *t* ¼ 1 ! *t max* do **For each particle do:**
	- 3. Step 3. Calculate the control law *ui*ð Þ*k* .
	- 4. Step 4. Determine the parameters *r*0, *r*1, *r*<sup>2</sup> and *s*<sup>1</sup>
	- 5. Step 5. Calculate the value of the fitness function from (52).
	- 6. Step 6. Update the value of *Pbest* and *Gbest*.
	- 7. Step 7. Update the position and velocity of each particle from Eqs. (??). **End for**

#### **End for**


#### **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]:

$$\begin{cases} \mathbf{C}\_{a}(k+1) = \mathbf{C}\_{a}(k) + T\_{\epsilon} \left( \frac{1}{\nu} q(k) \right) \left( \mathbf{C}\_{A\circ} - \mathbf{C}\_{a}(k) - k\_{0} \mathbf{C}\_{a}(k) e^{\overline{\nu T(k)}} \right) \\\\ T(k+1) = T(k) + T\_{\epsilon} \left( \frac{1}{\nu} q(k) \right) \left( T\_{f} - T(k) \right) + k\_{1} \mathbf{C}\_{a}(k) e^{\overline{\nu T(k)}} \\\\ \qquad + k\_{2} q\_{\epsilon}(k) \left( \mathbf{1} - e^{-\left(\frac{k\_{0}}{\alpha^{\boxplus}}\right)} \left( T\_{\epsilon f} - T(k) \right) \right) \end{cases} \tag{53}$$

with *<sup>k</sup>*<sup>1</sup> ¼ � <sup>Δ</sup>*Hk*<sup>0</sup> *<sup>ρ</sup>Cp* , *<sup>k</sup>*<sup>2</sup> <sup>¼</sup> *<sup>ρ</sup>cCpc <sup>ρ</sup>Cp<sup>υ</sup>* and *<sup>k</sup>*<sup>3</sup> <sup>¼</sup> *ha <sup>ρ</sup>cCpc* 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 *qc*. The nominal system parameters are listed in **Table 1**. We applied the PID-PSOCS algorithm to optimize the PID controller parameters for the *Ca* and *T* outputs.

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 *Ca k*ð Þ � <sup>1</sup> *Ca k*ð Þ � <sup>2</sup> *qc*ð Þ *<sup>k</sup>* � <sup>1</sup> *qc*ð Þ *<sup>k</sup>* � <sup>2</sup> *q k*ð Þ � <sup>1</sup> � � and *T k*ð Þ � <sup>1</sup> *T k*ð Þ � <sup>2</sup> *q k*ð Þ � <sup>1</sup> *q k*ð Þ � <sup>2</sup> *qc*ð Þ *<sup>k</sup>* � <sup>1</sup> � � 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: *tmax* ¼ 150, *Np* ¼ 10, *w max* ¼ 0*:*9, *w min* ¼ 0*:*4, *c*<sup>1</sup> ¼ 1*:*5, and *c*<sup>2</sup> ¼ 2*:*5. The three algorithms' optimal parameters for the two outputs (*Ca* and *T*) are provided in **Tables 2** and **3**, respectively.


**Table 1.** *CSTR system parameters.*

#### *PID Control for Nonlinear Processes DOI: http://dx.doi.org/10.5772/intechopen.106820*


#### **Table 2.**

*Different values of PID parameters for Ca output.*


#### **Table 3.**

*Different values of PID parameters for T output.*


#### **Table 4.**

*Performance of different methods for Ca output.*


#### **Table 5.**

*Performance of different methods for T output.*

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

**Figures 4** and **5** shows the evolution of the two system outputs *Ca* and *T* at their reference setpoints. The PID-PSOCS algorithm superimposed the other two algorithms in terms of the other performance indices namely *tm*, *ts*, and *D*ð Þ % .

**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

*PID Control for Nonlinear Processes DOI: http://dx.doi.org/10.5772/intechopen.106820*

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.

### **Author details**

Taieb Adel\*, Kanzari Bilel and Chaari Abdelkader Laboratory of Engineering of Industrial Systems and Renewable Energy, National High School of Engineers of Tunis (ENSIT), Tunis, Tunisia

\*Address all correspondence to: adel2taieb@gmail.com

© 2022 The Author(s). Licensee IntechOpen. This chapter is distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

### **References**

[1] Ghoshal S. Optimizations of pid gains byparticle swarm optimizations in fuzzy based automaticgeneration control. Electric Power Systems Research. 2004; **70**(3):203-212

[2] Jalili A, Shayeghi H, Shayanfar HA. T-sfuzzy parallel distribution compensation controller forpower system stabilizer. In: 5th International Conferenceon Technical and Physical Problems of Engineering (ICTPE-2009). Spain: Bilbao; 2009. pp. 180-184

[3] Guerra T, Vermeiren L. Lmi-based relaxednonquadratic stabilization conditions for nonlinearsystems in the takagi-sugeno's form. Automatica. 2004; **40**(500):823-829

[4] Khodabakhshian A, Edrisi M. A new robustpid load frequency controller. Control Engineering Practice. 2008; **16**(500):1069-1080

[5] Tan W. Tuning of pid load frequency controllerfor power systems. Energy Conversion and Management. 2009; **50**(500):1465-1472

[6] Aryan P, Raja GL. A novel equilibrium optimized double-loop control scheme for unstable and integrating chemical processes involving dead time. International Journal of Chemical Reactor Engineering. 2022;**20**(6)

[7] Anandh A, Aryan P, Kumari N, Raja G. Type-2 fuzzy-based branched controller tuned using arithmetic optimizer for load frequency control. Energy Sources, Part A: Recovery, Utilization, and Environmental Effects. 2022;**44**(2):4575-4596

[8] Shayeghi H, Shayanfar HA, Jalili A. Multistage fuzzy load frequency control using pso. EnergyConversion and Management. 2008;**49**(500):2570-2580

[9] Abonyi J, Babuška R, Szeifert F. Modified gath-geva fuzzy clustering for identification of takagi-sugeno fuzzy models. IEEE Transactions on Systems, Man and Cybernetics. 2002;**32**(5): 612-621

[10] Praly L. Robustness of indirect adaptive control based on pole placement design. In: Proceedings of the IFAC Workshop, San Francisco, USA, 20–22 June 1983. **1984**:55-60

[11] Kennedy J, Eberhart R. Particle swarm optimization. In: IEEE International Conferences on Neural Networks. Proceedings of ICNN'95 - International Conference on Neural Networks, Perth, Australia; 1995. pp. 1942-1948

[12] Yang XS, Deb S. Multi-objective cuckoo search for design optimization. Computers and Operations Research. 2013;**40**(500):1616-1624

[13] Payne RB, Sorenson MD, Klitz K. The Cuckoos. 1st edition. Oxford University Press; September 15, 2015

[14] Yang XS. Nature-Inspired Metaheuristic Algorithms. 2nd ed. UK: Luniver Press; 2010

[15] Chiroma H, Herawan T, Fister I. Bio-inspired computation:recent development on the modifications of the cuckoosearch algorithm. Applied Soft Computing. 2017;**61**(500):149-173

[16] Ouaarab A, Ahiod B, Yang X-S. Discrete cuckoo searchalgorithm for job shop scheduling problem. In: Yang X-S, editor. Proceedings of the International Symposium on Intelligent Control (ISIC). France; 2014. pp. 1872-1876

*PID Control for Nonlinear Processes DOI: http://dx.doi.org/10.5772/intechopen.106820*

[17] Zhao F, Liu Y, Zhang C. A selfadaptive harmony psosearch algorithm and its performance analysis. Expert Systemswith Applications. 2015;**42**(21): 7436-7455

[18] Alfi A. PSO with adaptive mutation and inertia weight and itsapplication in parameter estimation of dynamic systems. Acta Automatica Sinica. 2011; **37**(5):500-549

[19] Li X, Yin M. Modified cuckoo search algorithm with selfadaptive parameter method. Information Sciences. 2015; **298**(500):80-97

[20] Ikravesh M. Dynamic Neural Network Control. Columbia, SC: Universityof South Carolina; 1994

#### **Chapter 6**

## A Study on Numerical Solution of Fractional Order microRNA in Lung Cancer

*Mohammed Baba Abdullahi and Amiru Sule*

#### **Abstract**

The foremost cause of death resulting from cancer is lung cancer. From the statistics, 2.09 million new cases and 1.7 million deaths from lung cancer were estimated. In this chapter, the analytical solution of the concerned model was studied with help of the Laplace-Adomian Decomposition Method. To obtain the model's numerical scheme of fractional differential equations, the Caputo fractional derivative operator of order *α*∈ ð � 0, 1 is used. To find an approximate solution to a system of nonlinear fractional differential equations, the Laplace-Adomian Decomposition Method is used. Numerical simulations are presented to show the method's reliability and simplicity.

**Keywords:** fractional order, microRNA, numerical solution, cancer-related deaths, lung cancer

#### **1. Introduction**

The largest cause of cancer-related deaths worldwide is lung cancer. In the United States, a projected 236,740 people will receive a lung cancer diagnosis in 2022, making it the 16th most common cancer overall (1 in 15 males and 1 in 17 women). Smoking causes 80% of lung cancer fatalities and is the main risk factor for the disease. Twenty percent of lung cancer deaths occur in people who have never smoked. The second most important risk factor for lung cancer is radon gas exposure [1, 2]. Depending on the average radon level and the incidence of smoking in a nation, radon contributes to anywhere between 3 and 14% of lung cancer cases. Smokers are 25 times more likely than non-smokers to develop lung cancer from radon than non-smokers are likely to develop the cancer [3].

Early detection of high-risk lung cancer cases can reduce the chance of death by up to 20%. If you smoke now or have in the past, ask your doctor if lung cancer screening may be right for you. Approximately 8 million Americans are at high risk for lung cancer and could benefit from a lung cancer screening and yet only 5.7% actually get screened [4]. The dismal statistics associated with lung cancer are a result of both a lack of early detection and a lack of effective target therapy. Therefore, it is likely that developments in both of these areas will end in better results.

MicroRNAs (miRNAs) are a class of short nonprotein-coding RNAs (20–25 nucleotides in length) that predominantly inhibit the expression of target messenger RNAs (mRNAs) by directly interacting with their 30-untranslated regions (30UTRs) [5]. Numerous biological processes, from organismic development to tumor progression, depend on microRNAs in one way or another. These microRNAs have a crucial regulatory function in the pathogenesis of cancer in oncology, which forms the basis for investigating the influences on clinical characteristics using transcriptome data [6]. The seed match architecture between the mRNA seeding and miRNA binding regions determines the fate of the target mRNA. Perfect miRNA complementarity with the seeding sequence induces mRNA degradation, but imperfect or partial complementarities decrease protein translation [5].

Given the fact that a single miRNA may regulate tens to hundreds of genes, understanding the importance of an individual miRNA in cancer biology can be challenging. This is further complicated by observations that the dysregulation of several miRNAs is often required to cause a given phenotype [7]. To date, few models exist to elucidate the mechanisms by which multiple miRNAs contribute both individually and in tandem to promote tumor initiation and progression [8]. However, applying mathematical modeling to miRNA biology provides an opportunity to understand these complex relationships. In the work of Bersimbaev et al. [8], they developed for the first time a mathematical model focusing on miRNAs (miR-9 and let-7) in the context of lung cancer as a mathematical model system.

The organization of this chapter is as follows. In Section 2, some mathematical preliminaries of fractional calculus are needed to demonstrate the main results. The formulation of the Laplace-Adomian Decomposition Method (LADM) and Differential Transform Method (DTM) are given in Section 3. In Section 4, the numerical simulations are presented. In Section 5, the conclusions are given.

For the details of the integer mathematical model see Ref. [8], and the model is given below.

$$\begin{aligned} \frac{d}{dt}S(t) &= \mu\_{\rm S} \frac{S\_{\rm tot} - S}{S\_{\rm tot} - S + K\_{\rm S}} - \delta\_{\rm S} \text{Ek} \frac{S}{S + K\_{\rm S}}\\ \frac{d}{dt}R(t) &= \mu\_{\rm R} \frac{R\_{\rm tot} - R}{R\_{\rm tot} - R + K\_{\rm R}} \frac{K\_{\rm R2}}{L + K\_{\rm R2}} - \delta\_{\rm R} \frac{R}{R + K\_{\rm R3}}\\ \frac{d}{dt}\text{Ek}(t) &= \mu\_{\rm Ek} \frac{Ek\_{\rm tot} - Ek}{Ek\_{\rm tot} - Ek + K\_{\rm Ek1}} - \delta\_{\rm Ek} \frac{Ek}{Ek + K\_{\rm Ek2}}\\ \frac{d}{dt}\text{C}(t) &= \mu\_{\rm C} \text{Ek} - \delta\_{\rm C} \text{C}\\ \frac{d}{dt}\text{M}(t) &= \mu\_{\rm M} \frac{C^{4}}{C^{4} + K\_{\rm M}} - \delta\_{\rm M} \text{M} \\ \frac{d}{dt}\text{L}(t) &= \mu\_{\rm L} \frac{K\_{\rm L}}{C + K\_{\rm L}} - \delta\_{\rm L} \text{L} \\ \frac{d}{dt}\text{H}(t) &= \mu\_{\rm H} \text{L} \frac{K\_{\rm H}}{M + K\_{\rm H}} - \delta\_{\rm H} \text{H} \\ \frac{d}{dt}\text{P}(t) &= \mu\_{\rm P} - \delta\_{\rm P} \frac{H}{H + K\_{\rm P}} \end{aligned} \tag{1}$$

With the given initial condition.

*A Study on Numerical Solution of Fractional Order microRNA in Lung Cancer DOI: http://dx.doi.org/10.5772/intechopen.111387*

*S*ð Þ¼ 0 *n*1, *R*ð Þ¼ 0 *n*2, *Ek*ð Þ¼ 0 *n*3,*C*ð Þ¼ 0 *n*4, *M*ð Þ¼ 0 *n*5, *L*ð Þ¼ 0 *n*6, *H*ð Þ¼ 0 *n*7, *P*ð Þ¼ 0 *n*8, where tables give the descriptions of the state variables and parameters.

$$\begin{aligned} \,^cD^m S(t) &= \mu\_S E \frac{S\_{tot} - S}{S\_{tot} - S + K\_{S1}} - \delta\_S E k \frac{S}{S + K\_{S2}} \\ \,^cD^m R(t) &= \mu\_R S \frac{R\_{tot} - R}{R\_{tot} - R + K\_{R1}} \frac{K\_{R2}}{L + K\_{R2}} - \delta\_R \frac{R}{R + K\_{R3}} \\ \,^cD^m E k(t) &= \mu\_{Ek} R \frac{E k\_{tot} - E k}{E k\_{tot} - E k + K\_{Ek1}} - \delta\_{Ek} \frac{E k}{E k + K\_{Ek2}} \\ \,^cD^m S(t) &= \mu\_C E k - \delta\_C C \\ \,^cD^m M(t) &= \mu\_M \frac{C^4}{C^4 + K\_M} - \delta\_M M \\ \,^cD^m L(t) &= \mu\_L \frac{K\_L}{C + K\_L} - \delta\_L L \\ \,^cD^m H(t) &= \mu\_H L \frac{K\_H}{M + K\_H} - \delta\_H H \end{aligned} \tag{2}$$

With given initial condition.

*S*ð Þ¼ 0 *n*1, *R*ð Þ¼ 0 *n*2, *Ek*ð Þ¼ 0 *n*3,*C*ð Þ¼ 0 *n*4, *M*ð Þ¼ 0 *n*5, *L*ð Þ¼ 0 *n*6, *H*ð Þ¼ 0 *<sup>n</sup>*7, *<sup>P</sup>*ð Þ¼ <sup>0</sup> *<sup>n</sup>*8, where. *<sup>c</sup>*

*<sup>D</sup><sup>α</sup>* <sup>0</sup><*xi* <sup>≤</sup><sup>1</sup> *for i* <sup>¼</sup> 0, 1, 2 is the Caputo's derivate of fractional order and *<sup>x</sup>* shows fractional time derivative.

In model 2, the initial conditions are independent of each other and satisfy the relation.

*N*ð Þ¼ 0 *S t*ðÞþ *R t*ðÞþ *Ek t*ð Þþ *C t*ðÞþ *M t*ðÞþ *L t*ðÞþ *H t*ðÞþ *P t*ð Þ, where*N t*ð Þis the total population.

$$S(\mathbf{0}) = n\_1, R(\mathbf{0}) = n\_2, Ek(\mathbf{0}) = n\_3, C(\mathbf{0}) = n\_4, M(\mathbf{0}) = n\_5, L(\mathbf{0}) = n\_6, H(\mathbf{0}) = n\_7, P(\mathbf{0}) = n\_8 \tag{3}$$

#### **2. Preliminaries**

This section focuses on some basic definitions and outcomes from fractional calculus. For more in-depth, detailed research [9–11].

Definition 2.1. The fractional integral of Riemann-Liouville type of order *α*∈ ð Þ 0, 1 of a function *f* ∈*L*<sup>1</sup> ð Þ ½ � 0, *T* , ℜ is defined as:

The Caputo fractional order derivative of a function *f* on the interval at 0, ½ � *T* is defined by the following:

$${}^{c}D\_{0+}^{a}f(t) = \frac{1}{\Gamma(a)} \int\_{0}^{1} (t-s)^{n-a-1} f^{(n)}(s)ds,\tag{4}$$

when *n* ¼ j j *x* þ 1 *and x*j j represents the integer part of *x*. In particularity, for 0<*x*<1, Caputo derivative becomes

$$^cD\_{0+}^a f(t) = \frac{1}{\Gamma(a)} \int\_0^1 \frac{f(s)}{(1-s)} ds. \tag{5}$$

Lemma 2.1. The next outcome holds for fractional differential equations.

$$I^a(^cD^a h)(t) = h(t) + \sum\_{i=0}^{n-1} \frac{h^i(0)}{i!} t^i. \tag{6}$$

*for  arbitrary x*>0, *<sup>i</sup>* <sup>¼</sup> 0, 1, 2, … , *<sup>n</sup>* � 1, *when n* <sup>¼</sup> j j *<sup>x</sup>* <sup>þ</sup> <sup>1</sup>*and x*j j*represents*the integer part of x

Definition 2.2. We recall the definition of Laplace transform of Caputo derivative as:

$$\mathcal{E}\left\{{}^{c}D^{a}\boldsymbol{y}(t)\right\} = \mathfrak{s}^{a}h(\mathfrak{s}) - \sum\_{k=0}^{n-1} \mathfrak{s}^{a-i-1} \mathfrak{y}^{(k)}(\mathbf{0}), \ n - 1 < a < n, \ n \in N. \tag{7}$$

*for  arbitrary x* <sup>&</sup>gt;0, *<sup>i</sup>* <sup>¼</sup> 0, 1, 2, … , *<sup>n</sup>* � 1, *when n* <sup>¼</sup> j j *<sup>x</sup>* <sup>þ</sup> <sup>1</sup>*and x*j j*represents*the integer part of x*:*

#### **2.1 The Laplace-Adomian decomposition method**

This section focuses on model (3)'s overall operation under specified initial conditions. When both sides of the model are transformed using the Caputo fractional derivative system (3), the following results are obtained:

$$\begin{aligned} \begin{aligned} \{^C D^m S(t)\} &= \mu\_S E \frac{S\_{tot} - S}{S\_{tot} - S + K\_{S1}} - \delta\_S Ek \frac{S}{S + K\_{S2}} \\ L\{^C D^m R(t)\} &= \mu\_R S \frac{R\_{tot} - R}{R\_{tot} - R + K\_{R1}} \frac{K\_{R2}}{L + K\_{R2}} - \delta\_R \frac{R}{R + K\_{R3}} \\ L\{^C D^m Ek(t)\} &= \mu\_{Ek} R \frac{Ek\_{tot} - Ek}{Ek\_{tot} - Ek + K\_{Ek1}} - \delta\_{Ek} \frac{Ek}{Ek + K\_{Ek2}} \\ L\{^C D^m S(t)\} &= \mu\_C Ek - \delta\_C C \end{aligned} \tag{8} \\ \begin{aligned} \Delta\{^C D^m M(t)\} &= \mu\_M \frac{C^4}{C^4 + K\_M} - \delta\_M M \\ L\{^C D^m L(t)\} &= \mu\_L \frac{K\_L}{C + K\_L} - \delta\_L L \end{aligned} \tag{9} \\ L\{^C D^m H(t)\} &= \mu\_H L \frac{K\_H}{M + K\_H} - \delta\_H H \end{aligned} \tag{9}$$

*A Study on Numerical Solution of Fractional Order microRNA in Lung Cancer DOI: http://dx.doi.org/10.5772/intechopen.111387*

thus indicates

*s <sup>α</sup>*1*L <sup>c</sup> <sup>D</sup><sup>α</sup>*<sup>1</sup> f g *S t*ð Þ � *<sup>s</sup> α*1�1 *<sup>S</sup>*ð Þ¼ <sup>0</sup> *<sup>L</sup> <sup>μ</sup>SE Stot* � *<sup>S</sup> Stot* � *S* þ *KS*<sup>1</sup> � *<sup>δ</sup>SEk <sup>S</sup> S* þ *KS*<sup>2</sup> *s <sup>α</sup>*2*L <sup>c</sup> <sup>D</sup>*f g *<sup>α</sup>*2*R t*ð Þ � *<sup>s</sup> α*1�1 *<sup>R</sup>*ð Þ¼ <sup>0</sup> *<sup>L</sup> <sup>μ</sup>RS Rtot* � *<sup>R</sup> Rtot* � *R* þ *KR*<sup>1</sup> *KR*<sup>2</sup> *L* þ *KR*<sup>2</sup> � *δ<sup>R</sup> R R* þ *KR*<sup>3</sup> *s <sup>α</sup>*3*L <sup>c</sup> <sup>D</sup><sup>α</sup>*<sup>3</sup> f g *Ek t*ð Þ � *<sup>s</sup> α*1�1 *Ek*ð Þ¼ <sup>0</sup> *<sup>L</sup> <sup>μ</sup>EkR Ektot* � *Ek Ektot* � *Ek* þ *KEk*<sup>1</sup> � *δEk Ek Ek* þ *KEk*<sup>2</sup> *s <sup>α</sup>*4*L <sup>c</sup> <sup>D</sup>*f g *<sup>α</sup>*4*C t*ð Þ � *<sup>s</sup> α*1�1 *C*ð Þ¼ 0 *L*f g *μCEk* � *δcC s <sup>α</sup>*5*L <sup>c</sup> <sup>D</sup>α*5*M t*ð Þ � *<sup>s</sup> α*1�1 *M*ð Þ¼ 0 *L μ<sup>M</sup> C*4 *<sup>C</sup>*<sup>4</sup> <sup>þ</sup> *KM* � *δMM s <sup>α</sup>*6*L <sup>c</sup> <sup>D</sup>α*<sup>6</sup> f g *L t*ð Þ � *<sup>s</sup> α*1�1 *L*ð Þ¼ 0 *L μ<sup>L</sup> KL C* þ *KL* � *δLL s <sup>α</sup>*7*L <sup>c</sup> <sup>D</sup>*f g� *<sup>α</sup>*7*H t*ð Þ *<sup>s</sup> α*1�1 *<sup>H</sup>*ð Þ¼ <sup>0</sup> *<sup>L</sup> <sup>μ</sup>HL KH M* þ *KH* � *δHH s <sup>α</sup>*8*L <sup>c</sup> <sup>D</sup>α*<sup>8</sup> f g *P t*ð Þ � *<sup>s</sup> α*1�1 *<sup>P</sup>*ð Þ¼ <sup>0</sup> *<sup>L</sup> <sup>μ</sup><sup>P</sup>* � *<sup>δ</sup>PP <sup>H</sup> H* þ *KP* (9)

Using the initial conditions and taking inverse Laplace transform to system (5), we have:

$$\begin{aligned} S(t) &= S\_0 + L^{-1} \left\{ \mu\_S E \frac{S\_{\rm tot} - S}{S\_{\rm hot} - K + K\_{S1}} - \delta\_S E k \frac{S}{S + K\_{S2}} \right\} \\ R(t) &= R\_0 + L^{-1} \left\{ \mu\_R S \frac{R\_{\rm tot} - R}{R\_{\rm hot} - R + K\_{R1} L + K\_{R2}} - \delta\_R \frac{R}{R + K\_{R3}} \right\} \\ Ek(t) &= Ek\_0 + L^{-1} \left\{ \mu\_{Ek} R \frac{Ek\_{\rm tot} - Ek}{Ek\_{\rm hot} - Ek + K\_{Ek1}} - \delta\_{Ek} \frac{Ek}{Ek + K\_{Ek2}} \right\} \\ C(t) &= C\_0 + L^{-1} \{\mu\_C Ek - \delta\_C C\} \\ M(t) &= M\_0 + L^{-1} \left\{ \mu\_M \frac{C^4}{C^4 + K\_M} - \delta\_M M \right\} \\ L(t) &= L\_0 + L^{-1} \left\{ \mu\_L \frac{K\_L}{C + K\_L} - \delta\_L L \right\} \\ H(t) &= H\_0 + L^{-1} \left\{ \mu\_H L \frac{K\_H}{M + K\_H} - \delta\_H H \right\} \end{aligned} \tag{10}$$

$$P(t) = P\_0 + L^{-1} \left\{ \mu\_P - \delta\_P P \frac{H}{H + K\_P} \right\}$$

Using the values of the initial condition in Eq. (6), we get:

$$\begin{aligned} S(t) &= n\_1 + L^{-1} \left\{ \mu\_S E \frac{S\_{\rm tot} - S}{S\_{\rm tot} - S + K\_{S1}} - \delta\_S E k \frac{S}{S + K\_{S2}} \right\} \\ R(t) &= n\_2 + L^{-1} \left\{ \mu\_R S \frac{R\_{\rm tot} - R}{R\_{\rm tot} - R + K\_{R1}} \frac{K\_{R2}}{L + K\_{R2}} - \delta\_R \frac{R}{R + K\_{R3}} \right\} \\ Ek(t) &= n\_3 + L^{-1} \left\{ \mu\_{Ek} R \frac{Ek\_{\rm tot} - Ek}{Ek\_{\rm tot} - Ek + K\_{Ek1}} - \delta\_{Ek} \frac{Ek}{Ek + K\_{Ek2}} \right\} \\ C(t) &= n\_4 + L^{-1} \{\mu\_C Ek - \delta\_C C\} \\ M(t) &= n\_5 + L^{-1} \left\{ \mu\_M \frac{C^4}{C^4 + K\_M} - \delta\_{\rm bf} M \right\} \\ L(t) &= n\_6 + L^{-1} \left\{ \mu\_L \frac{K\_L}{C + K\_L} - \delta\_{\rm bf} L \right\} \\ H(t) &= n\_7 + L^{-1} \left\{ \mu\_H L \frac{K\_H}{M + K\_H} - \delta\_H H \right\} \\ P(t) &= n\_8 + L^{-1} \left\{ \mu\_p - \delta\_P H \frac{H}{H + K\_P} \right\} \end{aligned} \tag{11}$$

Assume that the solutions, *S t*ð Þ, *R t*ð Þ, *Ek t*ð Þ,*C t*ð Þ, *M t*ð Þ, *L t*ð Þ, *H t*ð Þ, *P t*ð Þ in the form of infinite series, are given by:

$$\begin{aligned} S(t) &= \sum\_{n=0}^{\infty} S\_n(t), R(t) = \sum\_{n=0}^{\infty} R\_n(t), \; Ek(t) = \sum\_{n=0}^{\infty} Ek\_n(t). \; \mathcal{C}(t) = \sum\_{n=0}^{\infty} C\_n(t) \\\ M(t) &= \sum\_{n=0}^{\infty} M\_n(t), L(t) = \sum\_{n=0}^{\infty} L\_n(t), \; H(t) = \sum\_{n=0}^{\infty} H\_n(t). \; P(t) = \sum\_{n=0}^{\infty} P\_n(t) \end{aligned} \tag{12}$$

While the nonlinear term involved in the model is *S t*ð Þ*Ek t*ð Þ, *S t*ð Þ*R t*ð Þ, *F t*ð Þ*R t*ð Þ, *P t*ð Þ*H t*ð Þ and is decomposed as follows, where *Xn*, *Yn and Zn* are the Adomian polynomials defined as are:

$$\begin{cases} X\_n = \frac{1}{\Gamma(n+1)} \frac{d^n}{dt^n} \left[ \sum\_{k=0}^\infty \lambda^k S\_k \sum\_{k=0}^\infty \lambda^k E k\_k \right] |\lambda=0 \\\\ Y\_n = \frac{1}{\Gamma(n+1)} \frac{d^n}{dt^n} \left[ \sum\_{k=0}^\infty \lambda^k S\_k \sum\_{k=0}^\infty \lambda^k R\_k \right] |\lambda=0 \\\\ Z\_n = \frac{1}{\Gamma(n+1)} \frac{d^n}{dt^n} \left[ \sum\_{k=0}^\infty \lambda^k E k\_k \sum\_{k=0}^\infty \lambda^k R\_k \right] |\lambda=0 \\\\ W\_n = \frac{1}{\Gamma(n+1)} \frac{d^n}{dt^n} \left[ \sum\_{k=0}^\infty \lambda^k P\_k \sum\_{k=0}^\infty \lambda^k H\_k \right] |\lambda=0 \end{cases} \tag{13}$$

**82**

*A Study on Numerical Solution of Fractional Order microRNA in Lung Cancer DOI: http://dx.doi.org/10.5772/intechopen.111387*

The first three polynomials are given by:

$$\begin{cases} X\_0 = S\_0(t)Ek\_0(t), \\ X\_1 = S\_0(t)Ek\_1(t) + S\_1(t)Ek\_0(t) \\ X\_2 = 2S\_0(t)Ek\_2(t) + 2S\_1(t)Ek\_1(t) + 2S\_2(t)Ek\_0(t) \\ Y\_0 = S\_0(t)R\_0(t), \\ Y\_1 = S\_0(t)R\_1(t) + S\_1(t)R\_0(t) \\ Y\_2 = 2S\_0(t)R\_2(t) + 2S\_1(t)R\_1(t) + 2S\_2(t)R\_0(t) \\ Z\_0 = Ek\_0(t)R\_0(t), \\ Z\_1 = Ek\_0(t)R\_1(t) + Ek\_1(t)R\_0(t) \\ Z\_2 = 2Ek\_0(t)R\_2(t) + 2Ek\_1(t)R\_1(t) + 2Ek\_2(t)R\_0(t) \\ W\_0 = P\_0(t)H\_0(t), \\ W\_1 = P\_0(t)H\_1(t) + P\_1(t)H\_0(t) \\ W\_2 = 2P\_0(t)H\_2(t) + 2P\_1(t)H\_1(t) + 2P\_2(t)H\_0 \end{cases} (14)$$

Using Eqs. (8) and (10) in model (6), yields

*<sup>L</sup>* <sup>X</sup><sup>∞</sup> *n*¼0 *Sk*ð Þ*t* ( ) ¼ *S*0 *s* þ 1 *<sup>s</sup><sup>α</sup> <sup>L</sup> <sup>μ</sup>SE Stot* � *<sup>S</sup> Stot* � *S* þ *KS*<sup>1</sup> � *<sup>δ</sup>SEk <sup>S</sup> S* þ *KS*<sup>2</sup> � � � � *<sup>L</sup>* <sup>X</sup><sup>∞</sup> *n*¼0 *Rk*ð Þ*t* ( ) <sup>¼</sup> *<sup>R</sup>*<sup>0</sup> *s* þ 1 *<sup>s</sup><sup>α</sup> <sup>L</sup> <sup>μ</sup>RS Rtot* � *<sup>R</sup> Rtot* � *R* þ *KR*<sup>1</sup> *KR*<sup>2</sup> *L* þ *KR*<sup>2</sup> � *δ<sup>R</sup> R R* þ *KR*<sup>3</sup> � � � � *<sup>L</sup>* <sup>X</sup><sup>∞</sup> *n*¼0 *Ekk*ð Þ*t* ( ) <sup>¼</sup> *Ek*<sup>0</sup> *s* þ 1 *<sup>s</sup><sup>α</sup> <sup>L</sup> <sup>μ</sup>EkR Ektot* � *Ek Ektot* � *Ek* þ *KEk*<sup>1</sup> � *δEk Ek Ek* þ *KEk*<sup>2</sup> � � � � *<sup>L</sup>* <sup>X</sup><sup>∞</sup> *n*¼0 *Ck*ð Þ*t* ( ) <sup>¼</sup> *<sup>C</sup>*<sup>0</sup> *s* þ 1 *<sup>s</sup><sup>α</sup> <sup>L</sup>*f g *<sup>μ</sup>CEk* � *<sup>δ</sup>cC* � � *<sup>L</sup>* <sup>X</sup><sup>∞</sup> *n*¼0 *Mk*ð Þ*t* ( ) <sup>¼</sup> *<sup>M</sup>*<sup>0</sup> *s* þ 1 *<sup>s</sup><sup>α</sup> <sup>L</sup> <sup>μ</sup><sup>M</sup> C*4 *<sup>C</sup>*<sup>4</sup> <sup>þ</sup> *KM* � *δMM* � � � � *<sup>L</sup>* <sup>X</sup><sup>∞</sup> *n*¼0 *Lk*ð Þ*t* ( ) <sup>¼</sup> *<sup>L</sup>*<sup>0</sup> *s* þ 1 *<sup>s</sup><sup>α</sup> <sup>L</sup> <sup>μ</sup><sup>L</sup> KL C* þ *KL* � *δLL* � � � � *<sup>L</sup>* <sup>X</sup><sup>∞</sup> *n*¼0 *Hk*ð Þ*t* ( ) <sup>¼</sup> *<sup>H</sup>*<sup>0</sup> *s* þ 1 *<sup>s</sup><sup>α</sup> <sup>L</sup> <sup>μ</sup>HL KH M* þ *KH* � *δHH* � � � � *<sup>L</sup>* <sup>X</sup><sup>∞</sup> *n*¼0 *Pk*ð Þ*t* ( ) <sup>¼</sup> *<sup>P</sup>*<sup>0</sup> *s* þ 1 *<sup>s</sup><sup>α</sup> <sup>L</sup> <sup>μ</sup><sup>P</sup>* � *<sup>δ</sup>PP <sup>H</sup> H* þ *KP* � � � � (15)

Now, comparing like terms on both sides, yields

*L S*½ �¼ <sup>0</sup>ð Þ*<sup>t</sup> <sup>n</sup>*<sup>1</sup> *<sup>s</sup>* , *L R*½ �¼ <sup>0</sup>ð Þ*<sup>t</sup> <sup>n</sup>*<sup>2</sup> *<sup>s</sup>* , *L Ek* ½ �¼ <sup>0</sup>ð Þ*<sup>t</sup> <sup>n</sup>*<sup>3</sup> *<sup>s</sup>* , *L C*½ �¼ <sup>0</sup>ð Þ*<sup>t</sup> <sup>n</sup>*<sup>4</sup> *s* , *L M*½ �¼ <sup>0</sup>ð Þ*<sup>t</sup> <sup>n</sup>*<sup>5</sup> *<sup>s</sup>* , *L L*½ �¼ <sup>0</sup>ð Þ*<sup>t</sup> <sup>n</sup>*<sup>6</sup> *<sup>s</sup>* , *L H*½ �¼ <sup>0</sup>ð Þ*<sup>t</sup> <sup>n</sup>*<sup>7</sup> *<sup>s</sup>* , *L P*½ �¼ <sup>0</sup>ð Þ*<sup>t</sup> <sup>n</sup>*<sup>8</sup> *s* , *L S*ð Þ¼ <sup>1</sup> *μSE sα Stot* � *S*<sup>0</sup> *Stot* � *S*<sup>0</sup> þ *KS*<sup>1</sup> � *δS sα X*<sup>0</sup> *S*<sup>0</sup> þ *KS*<sup>2</sup> � � 1 *<sup>s</sup><sup>α</sup>*1þ<sup>1</sup> , *L R*ð Þ¼ <sup>1</sup> *μRS*<sup>0</sup> *Rtot Rtot* � *R*<sup>0</sup> þ *KR*<sup>1</sup> *KR*<sup>2</sup> *L*<sup>0</sup> þ *KR*<sup>2</sup> � *μ<sup>R</sup> Y*<sup>0</sup> *Rtot* � *R*<sup>0</sup> þ *KR*<sup>1</sup> *KR*<sup>2</sup> *L*<sup>0</sup> þ *KR*<sup>2</sup> � *δ<sup>R</sup> R*0 *R*<sup>0</sup> þ *KR*<sup>3</sup> � � 1 *<sup>s</sup><sup>α</sup>*2þ<sup>1</sup> , *L Ek* ð Þ¼ <sup>1</sup> *μEkR*<sup>0</sup> *Ektot Ektot* � *Ek*<sup>0</sup> þ *KEk*<sup>1</sup> � *μEk Z*0 *Ektot* � *Ek*<sup>0</sup> þ *KEk*<sup>1</sup> � *δEk Ek*<sup>0</sup> *Ek*<sup>0</sup> þ *KEk*<sup>2</sup> � � 1 *<sup>s</sup><sup>α</sup>*3þ<sup>1</sup> , *L C*ð Þ¼ <sup>1</sup> ð Þ *μCEk*<sup>0</sup> � *δcC*<sup>0</sup> 1 *<sup>s</sup>α*4þ<sup>1</sup> , *L M*ð Þ¼ <sup>1</sup> *μ<sup>M</sup> C*<sup>0</sup> 4 *C*<sup>0</sup> <sup>4</sup> <sup>þ</sup> *KM* � *δMM*<sup>0</sup> � � 1 *<sup>s</sup>α*5þ<sup>1</sup> , *L L*ð Þ¼ <sup>1</sup> *μ<sup>L</sup> KL C*<sup>0</sup> þ *KL* � *δLL*<sup>0</sup> � � 1 *<sup>s</sup>α*6þ<sup>1</sup> , *L H*ð Þ¼ <sup>1</sup> *μHL*<sup>0</sup> *KH M*<sup>0</sup> þ *KH* � *δHH*<sup>0</sup> � � 1 *<sup>s</sup>α*7þ<sup>1</sup> , *L P*ð Þ¼ <sup>1</sup> *μ<sup>P</sup>* � *δ<sup>P</sup> W*<sup>0</sup> *H*<sup>0</sup> þ *KP* � � 1 *<sup>s</sup>α*8þ<sup>1</sup> , … *L S*ð Þ¼ *<sup>n</sup>*þ<sup>1</sup> *<sup>μ</sup>SE Stot* � *Sn Stot* � *Sn* þ *KS*<sup>1</sup> � *δ<sup>S</sup> Xn Sn* þ *KS*<sup>2</sup> � � 1 *<sup>s</sup>α*1þ<sup>1</sup> , *L R*ð Þ¼ *<sup>n</sup>*þ<sup>1</sup> *μRSn Rtot Rtot* � *Rn* þ *KR*<sup>1</sup> *KR*<sup>2</sup> *Ln* þ *KR*<sup>2</sup> � *μ<sup>R</sup> Yn Rtot* � *Rn* þ *KR*<sup>1</sup> *KR*<sup>2</sup> *Ln* þ *KR*<sup>2</sup> � *δ<sup>R</sup> Rn Rn* þ *KR*<sup>3</sup> � � 1 *<sup>s</sup>α*2þ<sup>1</sup> , *L Ek* ð Þ¼ *<sup>n</sup>*þ<sup>1</sup> *μEkRn Ektot Ektot* � *Ekn* þ *KEk*<sup>1</sup> � *μEk Zn Ektot* � *Ekn* þ *KEk*<sup>1</sup> � *δEk Ekn Ekn* þ *KEk*<sup>2</sup> � � 1 *<sup>s</sup>α*3þ<sup>1</sup> , *L C*ð Þ¼ *<sup>n</sup>*þ<sup>1</sup> ð Þ *μCEkn* � *δCCn* 1 *<sup>s</sup>α*4þ<sup>1</sup> , *L M*ð Þ¼ *<sup>n</sup>*þ<sup>1</sup> *μ<sup>M</sup> Cn* 4 *Cn* <sup>4</sup> <sup>þ</sup> *KM* � *δMMn* � � 1 *<sup>s</sup>α*5þ<sup>1</sup> , *L L*ð Þ¼ *<sup>n</sup>*þ<sup>1</sup> *μ<sup>L</sup> KL Cn* þ *KL* � *δLLn* � � 1 *<sup>s</sup>α*6þ<sup>1</sup> , *L H*ð Þ¼ *<sup>n</sup>*þ<sup>1</sup> *μHLn KH Mn* þ *KH* � *δHHn* � � 1 *<sup>s</sup>α*7þ<sup>1</sup> , *L P*ð Þ¼ *<sup>n</sup>*þ<sup>1</sup> *μ<sup>P</sup>* � *δ<sup>P</sup> Wn Hn* þ *KP* � � 1 *<sup>s</sup>α*8þ<sup>1</sup> *:* 8 >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>< >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>: (16)

Taking Laplace inverse of (11) and considering the first two terms at different values of*α* ¼ 1, 0*:*95, 0*:*85 and 0*:*75: and using the following values in **Tables 1** and **2**.


*A Study on Numerical Solution of Fractional Order microRNA in Lung Cancer DOI: http://dx.doi.org/10.5772/intechopen.111387*


#### **Table 1.**

*The state variables of the model.*


**Table 2.** *The parameters of the model.*

From *α* ¼ 1, 12 ð Þ *obtained*

$$\begin{cases} S(t) = 394.5868 + 21.03377825t + 39.51795700t^2 \\ R(t) = 0.0053 + 8874.951815t + 6221.509715t^2 \\ Ek(t) = 0.2488 - 0.00877591439t + 62646.35055t^2 \\ C(t) = 0.2189 - 0.00047867t + 0.0000862605909t^2 \\ M(t) = 0.000018t + 0.00002672931626t - 1.924510749x10^{-7}t^2 \\ L(t) = 0.0023 + 0.000006684617295t^2 \\ H(t) = 0.1 - 0.00023333t + 0.0000029641441673t^2 \\ P(t) = 1.157x10^{-13} + 3.4x10^{-20}t + 4.941947609x10^{-17}t^2 \end{cases} \tag{17}$$

From *α* ¼ 0*:*95, 12 ð Þ *obtained*

$$\begin{cases} S(t) = 394.5868 + 21.46565321t^{0.95} + 41.36344349t^{1.90} \\ R(t) = 0.0053 + 9057.176303t^{0.95} + 6512.053888t^{1.90} \\ Ek(t) = 0.2488 - 0.0047867t^{0.95} + 65571.93179t^{1.90} \\ C(t) = 0.2189 - 0.0004884982670t^{0.95} + 0.00009029571799t^{1.90} \\ M(t) = 0.000018 + 0.00002727813456t^{0.95} - 2.014385299x10^{-7}t^{1.90} \\ L(t) = 0.0023 + 0.000006996788568t^{1.90} \\ H(t) = 0.1 - 0.000238120836t^{0.95} + 0.000003102566932t^{1.90} \\ P(t) = 1.157x10^{-13} + 3.469810324x10^{-20}t^{0.95} + 5.17227363x10^{-17}t^{1.90} \end{cases} \tag{18}$$

From *α* ¼ 0*:*85, 10 ð Þ � 13 *obtained*

$$\begin{cases} S(t) = 394.5868 + 22.2435804t^{0.85} + 45.17936838t^{1.70} \\ R(t) = 0.0053 + 935.413409t^{0.85} + 712.814038t^{1.70} \\ Ek(t) = 0.2488 - 0.009280679638t^{0.85} + 71621.71589t^{1.70} \\ C(t) = 0.2189 - 0.000506201758t^{0.85} + 0.0000986258194t^{1.70} \\ M(t) = 0.000018 + 0.0000226670933t^{0.85} - 2.200219512x10^{-7}t^{1.70} \\ L(t) = 0.0023 + 0.000007642267217t^{1.70} \\ H(t) = 0.1 - 0.00024675046770^{0.85} + 0.000003388789774t^{1.70} \\ P(t) = 1.157x10^{-13} + 3.5955818010^{-20}t^{0.85} + 5.649939635x10^{-17}t^{1.70} \end{cases} \tag{19}$$

From *α* ¼ 0*:*75, 12 ð Þ *obtained*

$$\begin{cases} S(t) = 394.5868 + 22.88612324t^{0.75} + 49.1470382t^{1.50} \\ R(t) = 0.0053 + 9656.526685t^{0.75} + 7736.466899t^{1.50} \\ E(t) = 0.2488 - 0.00954876704t^{0.75} + 77900.93395t^{1.50} \\ C(t) = 0.2189 - 0.0005208241943t^{0.75} + 0.000107273391t^{1.50} \\ M(t) = 0.000018 + 0.00002908324024t^{0.75} - 2.393135170 \text{x} 10^{-7}t^{1.50} \\ L(t) = 0.0023 + 0.00000831234631t^{1.50} \\ H(t) = 0.1 - 0.0002538782653^{0.75} + 0.0000003685919493t^{1.50} \\ P(t) = 1.157 \text{x} 10^{-13} + 3.699421857 \text{x} 10^{-20}t^{0.75} + 6.145327396 \text{x} 10^{-17}t^{1.50} \end{cases} \tag{20}$$

*A Study on Numerical Solution of Fractional Order microRNA in Lung Cancer DOI: http://dx.doi.org/10.5772/intechopen.111387*

#### **2.2 Differential transform method**

The following recurrence relation to the system (2) with respect to time ð Þ*t* is obtained.

*S k*ð Þ¼ þ 1 1 *<sup>k</sup>* <sup>þ</sup> <sup>1</sup> *<sup>μ</sup>SE Stot* � *S k*ð Þ *Stot* � *S k*ð Þþ *KS*<sup>1</sup> *<sup>∂</sup>*ð Þ� *<sup>k</sup> <sup>δ</sup><sup>S</sup>* P*n i*¼0 *S l*ð Þ*Ek k*ð Þ � *l S k*ð Þþ *KS*<sup>2</sup> 2 6 6 4 3 7 7 5, *R k*ð Þ¼ þ 1 1 *<sup>k</sup>* <sup>þ</sup> <sup>1</sup> *<sup>μ</sup>RS k*ð Þ *Rtot Rtot* � *R k*ð Þþ *KR*<sup>1</sup> *KR*<sup>2</sup> *L k*ð Þþ *KR*<sup>2</sup> � *μ<sup>R</sup>* P*n i*¼0 *S l*ð Þ*R k*ð Þ � *l Rtot* � *R k*ð Þþ *KR*<sup>1</sup> *KR*<sup>2</sup> *L k*ð Þþ *KR*<sup>2</sup> � *δ<sup>R</sup> R k*ð Þ *R k*ð Þþ *KR*<sup>3</sup> 2 6 6 4 3 7 7 5, *Ek k*ð Þ¼ þ 1 1 *<sup>k</sup>* <sup>þ</sup> <sup>1</sup> *<sup>μ</sup>EkR k*ð Þ *Ektot Ektot* � *Ek k*ð Þþ *KEk*<sup>1</sup> � *μEk* P*n i*¼0 *Ek l*ð Þ*R k*ð Þ � *l Ektot* � *Ek k*ð Þþ *KEk*<sup>1</sup> � *δEk Ek k*ð Þ *Ek k*ð Þþ *KEk*<sup>2</sup> 2 6 6 4 3 7 7 5, *C k*ð Þ¼ þ 1 1 *k* þ 1 ½ � *μCEk k*ð Þ� *δcC k*ð Þ , *M k*ð Þ¼ þ 1 1 *<sup>k</sup>* <sup>þ</sup> <sup>1</sup> *<sup>μ</sup><sup>M</sup> <sup>C</sup>*4ð Þ *<sup>k</sup> <sup>C</sup>*4ð Þþ *<sup>k</sup> KM* � *<sup>δ</sup>MM k*ð Þ � �, *L k*ð Þ¼ þ 1 1 *<sup>k</sup>* <sup>þ</sup> <sup>1</sup> *<sup>μ</sup><sup>L</sup> KL C*<sup>0</sup> þ *KL* � *<sup>δ</sup>LL k*ð Þ � �, *H k*ð Þ¼ þ 1 1 *<sup>k</sup>* <sup>þ</sup> <sup>1</sup> *<sup>μ</sup>HL k*ð Þ *KH M k*ð Þþ *KH* � *<sup>δ</sup>HH k*ð Þ � �, *P k*ð Þ¼ þ 1 1 *<sup>k</sup>* <sup>þ</sup> <sup>1</sup> *<sup>μ</sup><sup>P</sup>* � *<sup>δ</sup><sup>P</sup>* P*n i*¼0 *P l*ð Þ*H k*ð Þ � *l H k*ð Þþ *KP* 2 6 6 4 3 7 7 5, (21)

The inverse differential transform of *S k*ð Þis defined as: When *t*0is taken as zero, the given function *y x*ð Þ is declared by a finite series and the above equation can be written in the form *S t*ðÞ¼ <sup>P</sup><sup>2</sup> *<sup>i</sup>*¼<sup>0</sup>*S k*ð Þ*<sup>i</sup> k* .

By solving the above equation for

$$\{S(k+1), R(k+1), Ek(k+1), C(k+1), M(k+1), L(k+1), H(k+1) and P(k+1) \quad \text{(22)}\}$$

up to order 2, we get the function.

*S k*ð Þ, *R k*ð Þ, *Ek k*ð Þ*C k*ð Þ, *M k*ð Þ, *L k*ð Þ, *H k*ð Þ*and P k*ð Þ*S k*ð Þ, *E k*ð Þ,*I k*ð Þ *and R k*ð Þ of respectively

$$S(t) = \sum\_{i=0}^{2} S(k)\dot{r}^k, R(t) = \sum\_{i=0}^{2} R(k)\dot{r}^k, \; Ek(t) = \sum\_{n=0}^{\infty} Ek(k)\dot{r}^k, C(t) = \sum\_{n=0}^{\infty} C(k)\dot{r}^k \tag{23}$$

$$M(t) = \sum\_{n=0}^{\infty} M(k)\dot{r}^k, L(t) = \sum\_{n=0}^{\infty} L(k)\dot{r}^k, \; H(t) = \sum\_{n=0}^{\infty} H(k)\dot{r}^k, \; P(t) = \sum\_{n=0}^{\infty} P(k)\dot{r}^k$$

$$\begin{cases} S(t) = 394.5868 + 21.03377825t + 40.55133050t^2 \\\\ R(t) = 0.0053 + 8874.951815t + 6221.509905t^2 \\\\ E(t) = 0.2488 - 0.00877591439t + 62646.35060t^2 \\\\ C(t) = 0.2189 - 0.00047867t + 0.0000862605090t^2 \\\\ M(t) = 0.000018 + 0.00002672931626t - 1.924510768x10^{-7}t^2 \\\\ L(t) = 0.0023 + 0.000006684617295t^2 \\\\ H(t) = 0.1 - 0.00023333t + 0.00000029641441658t^2 \\\\ P(t) = 1.157x10^{-13} + 3.4x10^{-20}t + 4.941947609x10^{-17}t^2 \end{cases} \tag{24}$$

#### **3. Numerical results**

The plots below show the population of each compartment for different values of *αi*ð Þ *i* ¼ 1, 2, 3, 4 (**Figures 1**–**8**).

#### **3.1 The comparison plots of the LADM and DTM of different compartments**

*A Study on Numerical Solution of Fractional Order microRNA in Lung Cancer DOI: http://dx.doi.org/10.5772/intechopen.111387*

*The plot shows the population of active Ras concentration for αi*,ð Þ *i* ¼ 1, 2, 3 *.*

**Figure 3.** *The plot shows the population of active ERK concentration for αi*,ð Þ *i* ¼ 1, 2, 3 *.*

**Figure 4.** *The plot shows the population of active MYC protein concentration for αi*,ð Þ *i* ¼ 1, 2, 3 *.*

**Figure 5.** *The plot shows the population of miR-9 concentration for αi*,ð Þ *i* ¼ 1, 2, 3 *.*

**Figure 6.** *The plot shows the population of let* � 7 *concentration for αi*,ð Þ *i* ¼ 1, 2, 3 *.*

**Figure 7.** *The plot shows the population of E-cadherin concentration for αi*,ð Þ *i* ¼ 1, 2, 3 *.*

*A Study on Numerical Solution of Fractional Order microRNA in Lung Cancer DOI: http://dx.doi.org/10.5772/intechopen.111387*

**Figure 8.** *The plot shows the population of MMP in RNA concentration for αi*,ð Þ *i* ¼ 1, 2, 3 *.*

#### **4. Conclusion**

In this chapter, a fractional order differential equation model is considered. The model was investigated and a scheme for the numerical solution for the fractional order differential equation microRNA in lung cancer using LADM (**Figures 9**–**16**). The LADM is an effective technique to solve nonlinear mathematical models and

is extensively applied in engineering and applied mathematics. Applying Laplace-Adomian Decomposition Method to obtain the series solution of fractional the model and comparing the results of the model at *α* ∈ð � 0, 1 with the classical Differential

*A Study on Numerical Solution of Fractional Order microRNA in Lung Cancer DOI: http://dx.doi.org/10.5772/intechopen.111387*

**Figure 12.**

*The comparison plots of the dynamics of MYC protein concentration using LADM and DTM.*

**Figure 13.**

*The comparison plots of the dynamics of miR-9 concentration using LADM and DTM.*

Transform Method is the main contribution of the work. The solution obtained through this method strongly agrees with DTM as shown in **Figures 1**–**16**. The effect of fractional parameters on our obtained solutions is presented through graphs.

**Figure 15.**

*The comparison plots of the dynamics of E-cadherin concentration using LADM and DTM.*

*A Study on Numerical Solution of Fractional Order microRNA in Lung Cancer DOI: http://dx.doi.org/10.5772/intechopen.111387*

**Figure 16.** *The comparison plots of the dynamics of MMP mRNA concentration using LADM and DTM.*

#### **Author details**

Mohammed Baba Abdullahi<sup>1</sup> \* and Amiru Sule<sup>2</sup>

1 Department of Mathematics/Statistics, School of Mathematics and Computing, Kampala International University, Kampala, Uganda

2 Faculty of Science, Department of Mathematical Science, Federal University Gusau, Nigeria

\*Address all correspondence to: kutigibaba@gmail.com

© 2023 The Author(s). Licensee IntechOpen. This chapter is distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

### **References**

[1] Ferlay J, Ervik M, Lam F, Colombet M, Mery L, Piñeros M, Soerjomataram IBF. International Agency for Research on Cancer. Global Cancer Observatory Cancer Today. 2020;**419**:1-2

[2] American Cancer Society. Cancer Facts and Figures 2022. Atlanta: American Cancer Society; 2022

[3] Radon and Health. The Fact Sheets, World Health Organization (WHO). Switzerland: WHO; 2021. p. 1211

[4] American Lung Association. State of Lung Cancer. 2022. Available from: https://www.lungcancerresearchfounda tion.org/wp-content/uploads/FactSheet-2023.01.pdf

[5] Wu SG, Chang TH, Liu YN, Shih JY. MicroRNA in lung cancer metastasis. Cancers. 2019;**11**(2):265

[6] Kong D, Wang K, Zhang QN, Bing ZT. Systematic analysis reveals key microRNAs as diagnostic and prognostic factors in progressive stages of lung cancer. 2022. DOI: 10.48550/arXiv. 2201.05408

[7] Kang HW, Crawford M, Fabbri M, Nuovo G, Garofalo M, Nana-Sinkam SP, et al. A mathematical model for microRNA in lung cancer. PLoS One. 2013;**8**(1):e53663

[8] Bersimbaev R, Pulliero A, Bulgakova O, Kussainova A, Aripova A, Izzotti A. Radon biomonitoring and microRNA in lung cancer. International Journal of Molecular Sciences. 2020; **21**(6):2154

[9] Yakubu AA, Abdullah FA, Abdullahi MB. Analysis and numerical solution of fractional order control of COVID-19 using Laplace Adomian decomposition

method expansion. In: AIP Conference Proceedings. AIP Publishing LLC; 2021; **2423**(1):020017. DOI: 10.1063/ 5.0075649

[10] Haq F, Shah K, Rahman G, Shahzad M. Numerical solution of fractional order smoking model via Laplace Adomian decomposition method. Alexandria Engineering Journal. 2018;**57**(2):1061-1069

[11] Anjam YN, Shafqat R, Sarris IE, Ur Rahman M, Touseef S, Arshad M. A fractional order investigation of smoking model using Caputo-Fabrizio differential operator. Fractal and Fractional. 2022;**6**(11):623

#### **Chapter 7**

## Performance Comparison of the Ball and Beam System using Linear Quadratic Regulator Controller

*Abubakar Umar, Muhammed B. Mu'azu, Zaharuddeen Haruna, Ore-Ofe Ajayi, Nafisa S. Usman, Onoshoho J. Oghenetega and Abdulfatai D. Adekale*

#### **Abstract**

This paper proposes the performance comparison of a linear quadratic regulator (LQR) controller for the ball and beam system (BBS). The BBS is a standard benchmark control system, which has two degree-of-freedom (2 DOF). It is an open loop and a highly nonlinear unstable system. This makes its parameter difficult to be estimated accurately, hence designing a controller for it is a challenging task. MatheThe BBS was modelled using Euler–Lagrange modeling technique, while the LQR controller was used for the stabilization of the ball on the beam. Simulation was done in MATLAB/Simulink 2022b environment, and the results simulated showed that for the two weighting matrices (*Q and R*), the state weighting matrix had a higher penalty on the ball displacement, ball velocity, beam angle, and beam angular velocity at lower values of *Q*. For the state weighting matrix had a better effect of penalty performance on the BBS with lower values. Also, as the diagonal element of the state weighting matrix *Q* increases from 0.1 to 20, the values of the optimal controller *K* increase, the reduced Ricatti matrix *P* increases, and the estimated eigenvalues *E* reduce. This implies that the ball displacement, ball velocity, beam angle, and beam angular velocity are better at lower values of *Q*.

**Keywords:** ball and beam system (BBS), linear quadratic regulator (LQR), weighting matrices, optimal controller

#### **1. Introduction**

Nonlinear systems play an important role in the field of control engineering. This is because suitable control techniques are used to improve the system performance [1]. The ball and beam system (BBS) is a highly nonlinear benchmark control problem in the field of control engineering. This is similar to practical control problems like balance control, position control, and tracking control problems [2]. BBS consists of a rigid beam that rotates freely in a vertical plane around the axis, while the ball rolls along the beam. The system can be categorized into two configurations, the first

configuration is the ball and beam balancer, in which the beam is supported in the middle, and it rotates against its central axis. The second configuration is built with the beam supported by two-level arms on both sides. One of the level arms acts as the pivot, while the other is coupled to the motor output gear [3, 4]. The purpose of the BBS is to hold the ball in a desired position on the beam while controlling the ball position by adjusting the angle of the beam [5, 6].

The BBS is used to implement and analyze the results of different modern control algorithms [7]. The control structure of this system is used for many different schemes in practical applications. It is used for demonstrating control applications like aircraft roll-yaw applications [8]. It is widely used due to its nonlinearity, simplicity, and open-loop instability. The control objective is the stabilization of the ball on the beam while tracking the reference trajectory [1, 9].

The BBS comprises of the base, the ball, a beam, support block, gear, and motor. This is shown in **Figure 1**.

The beam consists of two parts; the first part of the beam consists of a rigid shaft, while the other part rotates up and down on which the ball moves freely on it [10].

However, there are some research done on the system in applying different control algorithms to stabilize and perform trajectory tracking of the ball on the beam. Rahmat et al. [11] investigated the performance of some control techniques that consist of a proportional-integral-derivative (PID) controller, linear quadratic regulator (LQR) controller, and neural network (NN) designed in terms of stabilization and trajectory tracking. It showed that the LQR controller had a better satisfactory result. In Ref. [12], the particle swarm optimization (PSO) algorithm was used to tune the gains of the PID controller for the BBS. The optimized PSO-PID controller was compared with fuzzy logic controller (FLC) and integral of time multiplied by absolute error (ITAE), which the optimized PID outperformed the other two techniques. Kazemi et al. [13] designed cascade PD and fuzzy cascade controllers for stabilization of the BBS. The gains of the PD controller were optimized using the asexual reproduction optimization (ARO) algorithm. The results of PD-optimized ARO were compared with the fuzzy-cascaded controller in which the tuned PD-ARO outperformed the fuzzy-cascaded PD. Also, Ezzabi et al. [14] demonstrated a nonlinear backstepping technique for controlling the ball position of the BBS. The results were compared with LQR controller, it showed that the nonlinear technique required less input magnitude to achieve a better performance than the LQR controller. In Ref. [15], control strategies that were based on optimal control synthesis were presented. LQR and H2 controllers were used to control the ball on the beam. The control systems were implemented on a real BBS with a data acquisition card of DSP F28335. A new control strategy was proposed by Ref. [16] to control the stabilization of the BBS by the use of

**Figure 1.** *Ball and beam system [10].*

#### *Performance Comparison of the Ball and Beam System using Linear Quadratic Regulator… DOI: http://dx.doi.org/10.5772/intechopen.110513*

active disturbance rejection control (ADRC) on the system. The simulated results were compared with the proportional integration differentiation controller in which ADRC had a better performance than the integration differentiation controller. While Howimanporn et al. [17] developed a nonlinear discrete optimal control technique for the regulation of all the state variables in the discrete mode of the BBS. The proposed controller showed passivity, stability, and optimality properties during closed-loop operation. In Ref. [18], the BBS was designed using pole placement and LQR. The ball was able to be stabilized on the beam, and the results showed that LQR performed better than the pole-placement method. While an adaptive control was implemented in Ref. [19] for the BBS. Linear state-feedback model reference adaptive control (MRAC) was used in synchronizing the states of the BBS with the given reference model. Results showed that the error convergence was improved for different sets of the sinusoidal reference signal for the MRAC with modified feedback gains.

The main contribution of this article is the investigation of the performance effect of LQR controller on the BBS. This will be done by varying the *Q* and *R* matrix on the system, and observing which of the weighting matrix has a penalty effect on the minimization of the performance index of the LQR controller. However, the simulation was implemented in MATLAB/Simulink 2022b environment by adopting Lagrange's equation for modeling the system.

The rest of the paper is organized as follows: Section one introduces the BBS, while section two presents the mathematical model of the system. Section three discusses the controller design of the system and section four gives the simulated results. Finally, the conclusion is given in section five.

#### **2. Mathematical model of the ball and beam system**

The BBS is a two degree-of-freedom (DOF) system, the lateral movement of the ball is represented by its position on the horizontal axis, while the vertical movement of the beam is represented by the angle with the horizontal axis [16, 18]. The ball position is given by a sensor allocated at one end of the beam. The angle of the beam is adjusted by a toque provided by an actuator placed at the other end, where there is a connected axis [20]. The BBS can be simplified and linearized using the following assumptions [21]:


The equation which describes the dynamics of the system is obtained by using Euler Lagrange equation based on the energy balance of the system as follows [22, 23]:

$$\left(\frac{I\_b}{R^2} + m\right)\ddot{r} + mg\sin\theta - mr\dot{\theta}^2 = 0\tag{1}$$

$$\left(mr^2 + I + I\_b\right)\ddot{\theta} + 2mr\dot{r}\dot{\theta} + mgr\cos\theta = \mu\tag{2}$$

where *m* is the mass of the ball, *g* is the acceleration due to gravity, *I* is the beam moment of inertia, *Ib* is the ball moment of inertia, *r* is the position of the ball, *R* is the radius of the ball, *θ* is the angle of the beam, and *u* is the torque applied to the beam.

The model of the system can be described by the following state variables as *x*<sup>1</sup> represents the ball position along the beam, *x*<sup>2</sup> is the velocity of the ball, *x*<sup>3</sup> is the beam angle, and *x*<sup>4</sup> is the beam angular velocity. The generalized coordinate is given as [24]:

$$\begin{bmatrix} \mathbf{x}\_1 \ \mathbf{x}\_2 \ \mathbf{x}\_3 \ \mathbf{x}\_4 \end{bmatrix} = \begin{bmatrix} r \\ \dot{r} \\ \theta \\ \dot{\theta} \end{bmatrix} \tag{3}$$

The state space representation of the system is given by [25]:

$$
\dot{\mathfrak{x}}\_1 = \mathfrak{x}\_2 \tag{4}
$$

$$
\dot{x}\_2 = -a\mathbf{x}\_3 \tag{5}
$$

$$
\dot{\mathfrak{x}}\_3 = \mathfrak{x}\_4 \tag{6}
$$

$$\mathbf{x}\_4 = -\beta \mathbf{x}\_3 + \chi \mathbf{x}\_4 \tag{7}$$

where

$$a = \frac{\mathcal{M}\_{\text{g}}}{\frac{I\_{\text{b}}}{R\_{\text{b}}} + \mathcal{M}} \tag{8}$$

$$\beta = \frac{\mathbf{M\_g}}{J + J\_b} \tag{9}$$

$$\gamma = \frac{1}{J + f\_b} \tag{10}$$

The parameters of the BBS used for this article are given in **Table 1**.


**Table 1.** *Ball and beam system parameters.* *Performance Comparison of the Ball and Beam System using Linear Quadratic Regulator… DOI: http://dx.doi.org/10.5772/intechopen.110513*

#### **3. Linear quadratic regulator (LQR) controller**

The LQR controller takes the state equation of the system as feedback and generates a feedback error signal, as shown in **Figure 2**.

Given the system dynamics, the optimization procedure is to find an optimal control law that minimizes the performance index *J*, which is given as [26]:

$$J = \int\_{0}^{t\_1} (\mathbf{x}^T Q \mathbf{x} + \mathbf{u}^T R \mathbf{u}) dt \tag{11}$$

The optimal control law and the optimal controller are given as [26]:

$$
\mu\_{opt} = -\mathbf{R}^{-1} \mathbf{B}^T \mathbf{P} \mathbf{x} \tag{12}
$$

$$K = \mathbb{R}^{-1} \mathbb{B}^{T} P \tag{13}$$

Substituting the value of the Eq. (13), into (12), the optimal control law is given as [23]:

$$
\mu\_{opt} = -\mathbf{K}\mathbf{x} \tag{14}
$$

From Eq. (11), the *P* matrix must satisfy the reduced matrix equation, which is given as [26]:

$$A^T P + PA - PBR^{-1}B^T P + Q = \mathbf{0} \tag{15}$$

#### **3.1 Selection of** *Q* **and** *R***matrices**

Using Eq. (11), matrices *Q* and *R* penalizes the performance of the states, which control the response of the system and the cost of the energy consumed by the system. To improve on the system performance, the *Q* matrix is been considered, while to improve on the cost, *R* matrix is been focused on. The *Q* and *R* matrices are chosen to be a diagonal matrix while considering the effect of increasing or decreasing their values.

The *Q* and *R* matrices used for this research are:

$$Q\_0 = \text{diag}([\mathbf{0.1} \quad \mathbf{0.1} \quad \mathbf{0.1} \quad \mathbf{0.1}]) \tag{16}$$

**Figure 2.** *LQR control structure.*

$$Q\_1 = \text{diag}([\mathbf{1} \quad \mathbf{1} \quad \mathbf{1} \quad \mathbf{1}]) \tag{17}$$

$$Q\_2 = \text{diag}([\mathbf{10} \quad \mathbf{10} \quad \mathbf{10} \quad \mathbf{10}]) \tag{18}$$

*Q*<sup>3</sup> ¼ *diag*ð Þ ½ � 20 20 20 20 (19)

$$R\_0 = \text{diag}([1])\tag{20}$$

#### **4. Results and discussion**

The results of the BBS were generated from MATAB 2022b environment. The system has been proven to be controllable and observable, which gives way to further perform analysis of stabilizing the ball on the beam. Various values of the states weighting matrices *Q* were simulated against the control weighting matrix *R*, and the various values of the optimal controllers *K*, reduced Ricatti matrix *P*, and the estimated eigenvalues *E* were extracted.

For *Q*0, *R*0, the values of the *K*0, *P*0, *E*<sup>0</sup> extracted are:

$$K\_0 = \begin{bmatrix} -1.0741 & -0.6445 & 1.8936 & 0.4192 \end{bmatrix} \tag{21}$$

$$\begin{bmatrix} 0.3761 & 0.1577 & -0.2447 & -0.0215 \\ 0.1577 & 0.1393 & -0.2487 & -0.0129 \end{bmatrix} \tag{22}$$

$$P\_0 = \begin{bmatrix} 0.1577 & 0.1393 & -0.2487 & -0.0129 \\ -0.2447 & -0.2487 & 0.7035 & 0.0379 \\ -0.0215 & -0.0129 & 0.0379 & 0.0084 \end{bmatrix} \tag{22}$$

$$E\_0 = \begin{bmatrix} -1.6935 + 0.0000i \\ -1.7190 + 2.1616i \\ -1.7190 - 2.1616i \\ -15.8277 + 0.0000i \end{bmatrix} \tag{23}$$

The effect of the ball displacement, ball velocity, beam angle, and beam angular velocity is shown in **Figure 3**.

From **Figure 3**, it is seen that the ball stabilizes at about 2.8 secs on the beam on the x-axis, while the ball's velocity was stabilized at about 5.05 secs. Also, the beam angle was stabilized at about 4.7 secs, while the beam angular velocity was stabilized at about 5.1 secs. This shows that the state weighting matrix *Q* has a penalty on the ball position and also on the beam angle.

Also, for *Q*1, *R*0, the values of the *K*1, *P*1, *E*<sup>1</sup> extracted are:

$$K\_1 = \begin{bmatrix} -1.6043 & -1.6123 & 5.0314 & 1.0960 \end{bmatrix} \tag{24}$$

$$P\_1 = \begin{bmatrix} 1.7958 & 0.7998 & -1.2208 & -0.0321 \\ 0.7998 & 0.9835 & -1.7350 & -0.0322 \\ -1.2208 & -1.7350 & 5.2886 & 0.1006 \\ -0.0321 & -0.0322 & 0.1006 & 0.0219 \end{bmatrix} \tag{25}$$

$$E\_1 = \begin{bmatrix} -1.1087 + 0.0000i \\ -1.8503 + 1.9018i \\ -1.8503 - 1.9018i \\ -49.9866 + 0.0000i \end{bmatrix} \tag{26}$$

*Performance Comparison of the Ball and Beam System using Linear Quadratic Regulator… DOI: http://dx.doi.org/10.5772/intechopen.110513*

**Figure 3.** *Q*0, *R*<sup>0</sup> *weighting matrices.*

The effect of the ball displacement, ball velocity, beam angle, and beam angular velocity is shown in **Figure 4**.

From **Figure 4**, it is seen that the ball stabilizes at about 5.5 secs on the beam on the x-axis, while the ball's velocity was stabilized at about 7.5 secs. Also, the beam angle was stabilized at about 5.2 secs, while the beam angular velocity was stabilized at about 5.8 secs. This shows that the state weighting matrix *Q* has a penalty on the ball position and also on the beam angle.

Also, for *Q*2, *R*0, the values of the *K*2, *P*2, *E*<sup>2</sup> extracted are:

$$K\_2 = \begin{bmatrix} -3.6906 & -4.8907 & 15.2386 & 3.2572 \end{bmatrix} \tag{27}$$

$$P\_2 = \begin{bmatrix} 15.6506 & 6.9593 & -10.4235 & -0.0738 \\ 6.9593 & 9.1483 & -15.8562 & -0.0978 \\ -10.4235 & -15.8562 & 48.9502 & 0.3048 \\ -0.0738 & -0.0987 & 0.3048 & 0.0652 \end{bmatrix} \tag{28}$$

$$E\_2 = 10^2 \begin{bmatrix} -0.0101 + 0.0000i \\ -0.0187 + 0.0187i \\ -0.0187 - 0.0187i \\ -1.5809 + 0.0000i \end{bmatrix} \tag{29}$$

The effect of the ball displacement, ball velocity, beam angle, and beam angular velocity is shown in **Figure 5**.

**Figure 4.** *Q*1, *R*<sup>0</sup> *weighting matrices.*

From **Figure 5**, it is seen that the ball stabilizes at about 7.1 secs on the beam on the x-axis, while the ball's velocity was stabilized at about 8.2 secs. Also, the beam angle was stabilized at about 6.1 secs, while the beam angular velocity was stabilized at about 6.6 secs. This shows that the state weighting matrix *Q* has a penalty on the ball position and also on the beam angle.

However, for *Q*3, *R*0, the values of the *K*3, *P*3, *E*<sup>3</sup> extracted are:

$$K\_3 = \begin{bmatrix} -4.9895 & -6.8948 & 21.4451 & 4.5670 \end{bmatrix} \tag{30}$$

$$P\_3 = \begin{bmatrix} \mathbf{31.0192} & \mathbf{13.7689} & -\mathbf{20.5469} & -\mathbf{0.0998} \\ \mathbf{13.7689} & \mathbf{18.1689} & -\mathbf{31.3889} & -\mathbf{0.1379} \\ -\mathbf{20.5469} & -\mathbf{31.3889} & \mathbf{96.9744} & \mathbf{0.4289} \\ -\mathbf{0.0998} & -\mathbf{0.1379} & \mathbf{0.4289} & \mathbf{0.0914} \end{bmatrix} \tag{31}$$

$$E\_3 = 10^2 \begin{bmatrix} -0.0101 + 0.0000i \\ -0.0187 + 0.0187i \\ -0.0187 - 0.0187i \\ -2.2358 + 0.0000i \end{bmatrix} \tag{32}$$

The effect of the ball displacement, ball velocity, beam angle, and beam angular velocity is shown in **Figure 6**.

*Performance Comparison of the Ball and Beam System using Linear Quadratic Regulator… DOI: http://dx.doi.org/10.5772/intechopen.110513*

**Figure 5.** *Q*2, *R*<sup>0</sup> *weighting matrices.*

**Figure 6.** *Q*3, *R*<sup>0</sup> *weighting matrices.*

From **Figure 6**, it is seen that the ball stabilizes at about 7.6 secs on the beam on the x-axis, while the ball's velocity was stabilized at about 8.7 secs. Also, the beam angle was stabilized at about 6.5 secs, while the beam angular velocity was stabilized at about 7.2 secs. This shows that the state weighting matrix *Q* has a penalty on the ball position and also on the beam angle.

From **Figures 3**–**6**, it can be deduced that as the leading element of the state weighting matrix *Q* increases from 0.1 to 20, the values of optimal controller *K* increases, the reduced Ricatti matrix *P* also increases, and the estimated eigenvalues also reduces. This implies that the ball displacement, ball velocity, beam angle, and beam angular velocity are better at lower values of *Q*. This shows that the state weighting matrix has an effect on penalty performance on the BBS within lower values of *Q*.

#### **5. Conclusion**

An analysis on the effect of the state and control weighting matrices ð Þ *Q andR* on a benchmark control problem, the ball and beam system (BBS) has been studied. The system is a 2 DOF, an open loop, and a highly nonlinear unstable system, which makes estimating its parameters difficult, hence designing a controller was a difficult task. The BBS is an underactuated system with multiple input multiple output characteristics. Lagrange modeling technique was used for modeling the system, and LQR controller was used for the stabilization of the ball on the beam. A simulation was done in MATLAB/Simulink 2022b environment, and it showed that the state weighting matrix had a higher penalty on the ball displacement, ball velocity, beam angle, and beam angular velocity at lower values of *Q*. Also, it showed that as *Q* increases from 0.1 to 20, the values of the optimal controller *K* increase, the reduced Ricatti matrix *P* increases, and the estimated eigenvalues *E* also reduce. This implies that the ball displacement, ball velocity, beam angle, and beam angular velocity are better at lower values of *Q*. Future research will consider the effect of varying the control weighting matrix *R* over some range of values on the BBS system, which can be further applied to the field of autonomous vehicles.

#### **Author details**

Abubakar Umar\*, Muhammed B. Mu'azu, Zaharuddeen Haruna, Ore-Ofe Ajayi, Nafisa S. Usman, Onoshoho J. Oghenetega and Abdulfatai D. Adekale Computer Engineering Department, Ahmadu Bello University, Zaria, Nigeria

\*Address all correspondence to: abubakaru061010@gmail.com

© 2023 The Author(s). Licensee IntechOpen. This chapter is distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

*Performance Comparison of the Ball and Beam System using Linear Quadratic Regulator… DOI: http://dx.doi.org/10.5772/intechopen.110513*

### **References**

[1] Mehedi IM, Al-Saggaf UM, Mansouri R, Bettayeb M. Two degrees of freedom fractional controller design: Application to the ball and beam system. Measurement. 2019;**135**:13-22

[2] Srivastava V, Srivastava S. Hybrid optimization based PID control of ball and beam system. Journal of Intelligent Fuzzy Systems. 2022;**42**(2): 919-928

[3] Rosa MR, Romdlony MZ, Trilaksono BR. The ball and beam system: Cascaded LQR-FLC design and implementation. International Journal of Control, Automation Systems Science Control Engineering. 2023;**21**(1):201-207

[4] Umar A, Muazu MB, Aliyu UD, Musa U, Haruna Z, Oyibo PO. Position and trajectory tracking control for the ball and plate system using mixed H∞ sensitivity problem. Covenant Journal of Informatics and Communication Technology, (CJICT). 2018;**6**(1):1-15

[5] Nguyen C, Phan HN, Hoang L, Tran HN, editors. The design of a quasitime optimal cascade controller for ball and beam system. In: IOP Conference Series: Materials Science and Engineering. IOP Publishing; 2021

[6] Mehedi IM, Al-Saggaf UM, Mansouri R, Bettayeb M. Two degrees of freedom fractional controller design: Application to the ball and beam system. Measurement. 2019;**135**:13-22

[7] Howimanporn S, Chookaew S, Silawatchananai C. Monitoring and controlling of a real-time ball beam fuzzy predicting based on PLC network and information technologies. Journal of Advances in Information Technology. 2022;**13**

[8] Kim Y, Kim S-K, Ahn CK. Variable cut-off frequency observer-based positioning for ball-beam systems without velocity and current feedback considering actuator dynamics. IEEE Transactions on Circuits Systems I: Regular Papers. 2020;**68**(1):396-405

[9] Umar A, Muazu MB, Aliyu UD, Musa U, Haruna Z, Oyibo PO. Position and trajectory tracking control for the ball and plate system using mixed H∞ sensitivity problem. Covenant Journal of Informatics and Communication Technology, (CJICT). 2018;**6**(1):1-15

[10] Borgohain N, Buragohain M. Comparative study of optimal controller application on nonlinear systems. In: Modeling, Simulation and Optimization. Springer; 2021. pp. 417-428

[11] Rahmat MF, Wahid H, Wahab NA. Application of intelligent controller in a ball and beam control system. International Journal on Smart Sensing Intelligent Systems. 2010;**3**(1): 45-60

[12] Rana MA, Usman Z, Shareef Z, editors. Automatic control of ball and beam system using particle swarm optimization. In: IEEE 12th International Symposium on Computational Intelligence and Informatics (CINTI), 2011. IEEE; 2011

[13] Kazemi M, Najafi J, Menhaj MB. Fuzzy PD, Cascade controller design for ball and beam system based on an improved ARO technique. Journal of Computer Robotics. 2012;**5**(1):1-6

[14] Ezzabi AA, Cheok KC, Alazabi FA, editors. A nonlinear backstepping control design for ball and beam system. In: 2013 IEEE 56th International Midwest Symposium on Circuits and

Systems (MWSCAS). Columbus, OH, USA: IEEE; 2013

[15] Hung B, You S, Kim H, Lim T. Embedded controller building for ball and beam system using optimal control synthesis. Journal of Engineering Science and Technology. 2017;**12**(6): 1460-1474

[16] Ding M, Liu B, Wang L. Position control for ball and beam system based on active disturbance rejection control. Systems Science Control Engineering. 2019;**7**(1):97-108

[17] Danilo Montoya O, Gil-González W, Ramírez-Vanegas C. Discrete-inverse optimal control applied to the ball and beam dynamical system: A passivitybased control approach. Symmetry. 2020;**12**(8):1359

[18] Rai A, Suman SK, Kumar A, Yadav S, editors. Impact of control stability using LQR and pole-placement for ball and beam system. In: 2021 5th International Conference on Intelligent Computing and Control Systems (ICICCS). Madurai, India: IEEE; 2021

[19] Romdlony MZ, Rosa MR, Syamsudin EMP, Trilaksono BR, ASJJoM W. Electrical Power,. Design and application of models reference adaptive control (MRAC) on ball and beam. Journal of Mechatronics, Electrical Power, Vehicular Technology. 2022; **13**(1):15-23

[20] Edutech QI. Ball and Beam Laboratory Manual. Quanser Inc.; 2012. pp. 1-21

[21] Hajipour S, Pourhashem H, Chegini SN, Bagheri A. Optimized neuro observer-based sliding mode control for a nonlinear system using fuzzy static sliding surface. Applied Soft Computing. 2022;**124**:108904

[22] Saleem MK, Shahid MLUR, Nouman A, Zaki H, Tariq MAUR. Design and implementation of adaptive neurofuzzy inference system for the control of an uncertain ball and beam apparatus. Mehran University Research Journal of Engineering Technology. 2022;**41**(2): 178-184

[23] Do Van D, editor Researching and applying sliding control method for ball and beam system. In: 2022 International Conference on Electrical, Computer, Communications and Mechatronics Engineering (ICECCME). Maldives, Maldives: IEEE; 2022

[24] Khan R, Malik FM, Raza A, Mazhar N, Ullah H, Umair M, editors. Robust nonlinear control design and disturbance estimation for ball and beam system. In: 2020 3rd International Conference on Computing, Mathematics and Engineering Technologies (iCoMET). Sukkur, Pakistan: IEEE; 2020

[25] Soni R, editor. Optimal Control of a Ball and Beam System through LQR and LQG. In: 2018 2nd International Conference on Inventive Systems and Control (ICISC). Coimbatore, India: IEEE; 2018

[26] Umar A, Haruna Z, Mu'azu MB, Shilintang SG, Usman NS, Adekale AD, editors. Performance analysis of a ballon-sphere system using linear quadratic regulator controller. In: 2022 IEEE Nigeria 4th International Conference on Disruptive Technologies for Sustainable Development (NIGERCON). Lagos, Nigeria: IEEE; 2022

### *Edited by Mohammad Shamsuzzoha and G. Lloyds Raja*

A proportional-integral-derivative (PID) controller is an instrument used in industrial control applications to regulate process variables such as temperature, pressure, flow, etc. PID controllers are still very much preferred in the industry due to their simplicity and ability to yield reasonable closed-loop performance. About 90% of industrial controllers are of the PID type. To meet the continuously evolving challenges in industrial process control, it is essential to formulate PID-based control strategies which can yield improved performance. The contents of this book will serve as a good introduction to PID controllers and equip readers to design such controllers for various industrial applications.

Published in London, UK © 2023 IntechOpen © struvictory / iStock

PID Control for Linear and Nonlinear Industrial Processes

PID Control for Linear

and Nonlinear Industrial

Processes

*Edited by Mohammad Shamsuzzoha* 

*and G. Lloyds Raja*