An internal model control method for wave rotor refrigeration process based on Harris hawks optimization

Abstract

A novel internal model control (IMC) method for wave rotor refrigeration process (WRRP) based on Harris hawks optimization (HHO) is presented in this paper. Firstly, an identification model is established by the test data of the WRRP. In order to reduce the complexity of the model, the Routh order reduction method was used to reduce the order of the identification model. This procedure yields an approximate second order model with time delay for the WRRP. Secondly, an IMC method based on HHO is proposed after analyzing the time delay characteristics of the WRRP model. This method improves the traditional IMC structure, and introduces a three degree of freedom IMC structure. The parameters of the controller are optimized by HHO to make the comprehensive performance index of the system reach the best. Simulation results demonstrate that the proposed method can achieve a good performance of tracking and disturbance rejection.

Introduction

Wave rotor refrigerator process (WRRP) is a gas expansion refrigeration device which uses the pressure wave moving in the oscillating tube to realize the separation of cold and heat. Its refrigeration principle mainly uses the expansion cooling of gas expansion wave and shock wave pressurization to achieve the energy transfer and transformation between gases. WRRP has the advantages of simple structure, low production cost, suitable for liquid containing work and strong adaptability to the working environment. At present, WRRP has been successfully applied to the pressure energy recovery of natural gas, light hydrocarbon recovery, low-temperature wind tunnel and other related engineering fields. ¹ The proportional-integral-derivative (PID) control scheme is usually adopted for the temperature control of wave rotor refrigeration process. Due to the unique structure and operation mode of the WRRP, the control system exists large time delay, which will increase the difficulty of PID parameter tuning and the performance of control system is not good enough. How to feed back the control results to the WRRP timely and reliably to improve the refrigeration efficiency is an urgent problem to be solved.

In recent years, many scholars have studied the control methods of time delay systems and proposed different types of PID controller design methods. ² Anwar and Pan ³ proposed a PID controller design method based on frequency response matching for the integral processes with time delay. Ajmeri and Ali ⁴ designed a controller of the modified Smith predictor for the pure integration, integration plus first-order and double integration processes with large time delay, and the direct synthesis method is used to adjust the set point tracking controller. In Panda ⁵ and Ribić and Mataušek ⁶ parameters of PID controller are set based on the principle of dead-time compensation. These methods can enhance the disturbance rejection performance of the system and have the advantages of simple structure, but the PID controller can only be applied to the system with small time delay. The ratio of time delay to time constant of this kind of system model must be limited to a certain range. Garcia and Morari ⁷ proposed internal model control (IMC) strategy based on identification model. According to the internal correspondence between IMC and traditional PID, many scholars designed IMC-PID controller, which achieved good control effect.^8–14 Some PID controller tuning methods based on IMC for stable and unstable time-delay processes were proposed, and achieved expected control effect. After analyzing the analytical relationship between the robustness index of the control system and the controller parameters based on the maximum sensitivity, Jin et al. ¹² proposed a simple IMC-PID tuning method, which converted the IMC controller to PID form in terms of the time domain rather than the frequency domain. For the second-order process with time delay, an IMC-PID controller cascaded with a lead-lag compensator is designed in Shamsuzzoha and Lee, ¹³ the control effect is selected by adjusting the parameters. In Chia and Lefkowitz, ¹⁴ it modeled the integration process as a first-order lag with a very large time constant, and designed a PID controller to assure zero steady-state error for setpoint changes and process disturbance inputs. Singh et al. ¹⁵ proposed a two degree of freedom IMC combined with the model reduction method, which was applied to the power system, and achieved desired dynamic response and robustness against load disturbance. Although the IMC parameters tuning rules are given in Refs.,^8–15 there are still many parameters to be adjusted, and the control performance is not necessarily the best. Especially when the reference model is inaccurate, the capability of target tracking and disturbance rejection as well as robustness of the control system must be compromised. At this time, it is more difficult to meet the design requirements to adjust the parameters by experience. Therefore, Refs.^16–20 applied some optimization algorithms to adjust the parameters, such as particle swarm optimization algorithm (PSO),^16,17 genetic algorithm (GA)^18,19 and so on. In Zeng et al, ¹⁶ an IMC-PID controller based on PSO algorithm for core power control of molten salt breeder reactor was designed. Wang et al. ¹⁹ proposed a multiple-model based IMC method for power control of small pressurized water reactor, and the IMC filter time constant is determined by multi-objective optimization using the non-dominated sorting genetic algorithm-II. Nisi et al. ²⁰ proposed a multi-objective optimization method to tune the PID controller parameters, improved the performance of process in terms of performance index, set point tracking and also provided stability. In this paper, as the WRRP system exists large time delay, which will increase the difficulty of PID parameter tuning and the performance of control system is not good enough. So, we want to research an advanced control method based on model identification method and optimization method for the WRRP to improve the performance of control system. Harris hawks optimization (HHO) is a naturally inspired algorithm proposed by Heidari et al. ²⁴ The HHO is able to efficiently avoid local optimum and immature convergence drawbacks, HHO algorithm provides very promising and competitive results compared to other intelligence optimization algorithms. For these reasons, we present a three degree of freedom IMC control structure based on HHO algorithm for the temperature control of wave rotor refrigeration process.

We firstly collecting the test data and get an approximate second order model with time delay for the WRRP. The IMC method based on HHO is proposed after analyzing the characteristics of the identified model. This modified method improves the traditional IMC structure, and introduces a three degree of freedom IMC control structure. The parameters of the controller are optimized by HHO to make the comprehensive performance index of the system reach the best. Simulation results demonstrate that the proposed method can enable the WRRP to obtain good performance of set point tracking and disturbance rejection as well as robustness simultaneously.

The organization of this paper is arranged as follows. Section 2 presents the dynamical model of WRRP. Design of improved IMC is presented in section 3. The HHO algorithm is introduced in section 4. In section 5, simulation results are performed to validate, to show the effectiveness of the designed controllers and are discussed to emphasize the performance of the proposed control technique in this work, and finally, the conclusion is given in section 6.

Modeling of WRRP

Description of WRRP

A schematic diagram of the WRRP is shown in Figure 1, which mainly including air compressor, dryer, heat exchanger and wave rotor refrigerator (WRR). There are four ports on the WRR, high-pressure (HP) port, low-temperature (LT) port, medium pressure (MP) port and low-pressure (LP) port. Furthermore, the HP and LP ports are air inlets, while the MP and LT ports are air outlets. Each port is equipped with a valve to control the gas pressure. The gas is compressed by the compressor and dried by the adsorption dryer, and then HP gas flows into the WRR from the HP port, the wave rotor rotates to refrigerate. After energy exchange, the gas that flows out of WRR through the MP port flows back to the wave rotor through the LP port after dumping heat to the cooling water, finally, the cold gas flows out of the LT port.^1,21

Figure 1.

Schematic diagram of the WRRP.

In the WRRP, the temperature and pressure of the output gas are directly related to the working efficiency of the WRR, which are important indexes to measure the performance of the WRR. This paper mainly discusses the relationship between the temperature TT103 of the output gas and the inlet flow TV101 of the WRRP. For the existence of large time delay process dynamic characteristics of the WRRP, a three degree of freedom IMC control structure based on HHO algorithm is proposed in this paper. Firstly, the WRRP model is obtained by system identification and model reduction method based on the process test data, the specific identification algorithm and high-order model reduction algorithm are described in the next subsections.

Model identification algorithm of WRRP

Assuming that the discrete system model of WRRP is as follows:

$$G(z^{-1})=\frac{b_1+b_2{z}^{-1}+\cdots +b_m{z}^{-m+1}}{1+a_1{z}^{-1}+a_2{z}^{-2}+\cdots +a_n{z}^{-n}}{z}^{-d}$$

(1)

Its corresponding auto-regressive exogenous (ARX) model difference equation model is as follows:

$$\begin{array}{c}y(t)+a_{1}y(t-1)+a_{2}y(t-2)+\cdots +a_{n}y(t-n)\hfill \\ =b_{1}u(t-d)+b_{2}u(t-d-1)+\cdots \hfill \\ +b_{m}u(t-d-m+1)+\epsilon (t)\hfill \end{array}$$

(2)

where $\epsilon (t)$ is the residual signal, $y(t)$ is the output signal, and $y(t-1)$ is used to represent the function value of the output signal $y(t)$ at the previous sampling period. Assuming that a set of input signals $u=[u(1),u(2),\dots ,u(M){]}^T$ and a set of output signals $y=[y(1),y(2),\dots ,y(M){]}^T$ have been measured, the relationship between them can be described by Eq. (3):

$$\begin{array}{c}y(1)=-a_{1}y(0)-\cdots -a_{n}y(1-n)+b_{1}u(1-d)\hfill \\ +\cdots +b_{m}u(2-m-d)+\epsilon (1)\hfill \\ y(2)=-a_{1}y(1)-\cdots -a_{n}y(2-n)+b_{1}u(2-d)\hfill \\ +\cdots +b_{m}u(3-m-d)+\epsilon (2)\hfill \\ ⋮\hfill \\ y(M)=-a_{1}y(M-1)-\cdots -a_{n}y(M-n)+b_{1}u(M-d)\hfill \\ +\cdots +b_{m}u(M-m-d)+\epsilon (M)\hfill \end{array}$$

(3)

The initial values of $y(t)$ and $u(t)$ are assumed to be zero when $t\le 0$ . The Eq. (3) is written in matrix form as:

$$y=\mathrm{\Phi }\theta +\epsilon$$

(4)

where

$$\mathrm{\Phi }=|\begin{array}{cccccc}y(0)\hfill & \cdots \hfill & y(1-n)\hfill & u(1-d)\hfill & \cdots \hfill & u(2-m-d)\hfill \\ \multicolumn{1}{c}{y(1)}& \cdots \hfill & y(2-n)\hfill & u(2-d)\hfill & \cdots \hfill & u(3-m-d)\hfill \\ \multicolumn{1}{c}{⋮}& \hfill & ⋮\hfill & ⋮\hfill & \hfill & ⋮\hfill \\ \multicolumn{1}{c}{y(M-1)}& \cdots \hfill & y(M-n)\hfill & u(M-d)\hfill & \cdots \hfill & u(M+1-m-d)\hfill \end{array}|,$$

$$\begin{array}{c}{\theta }^T=[-a_1,-a_2,\cdots ,-a_n,b_1,b_2\cdots ,b_m],\hfill \\ {\epsilon }^T=[\epsilon (1),\epsilon (2),\cdots \epsilon (M)]\hfill \end{array}$$

To minimize the sum of squares of residuals, namely $\underset{\theta }{min}\sum _{i=1}^M{\epsilon }^2(i)$ , the optimal estimation of undetermined parameters can be determined by the least square method as:

$$\theta ={[{\mathrm{\Phi }}^T\mathrm{\Phi }]}^{-1}{\mathrm{\Phi }}^{T}y$$

(5)

Routh order reduction method

Nevertheless, analyzing and designing the high-order model has always been a tough task, it is necessary to select the appropriate model reduction method which can approximate a high-order model to a low-order model keeping the inherent property of the original system. The Routh approximation method characterizes simple calculation, easy acquisition and ensuring the stability of the model. Therefore, the Routh approximation method is used to reduce the order of the model in this paper.^22,23

Assuming the transfer function of a higher-order system is in the following form:

$$H(s)=\frac{N_0+N_{1}s+N_2{s}^2+\cdots +N_{n-1}{s}^{n-1}}{d_0+d_{1}s+d_2{s}^2+\cdots +d_{n-1}{s}^{n-1}+d_n{s}^n}$$

(6)

where $N_i(i=0,1,\dots n-1)$ , $d_i(i=0,1,\dots n)$ are constant coefficients, $d_n=1$ . The Routh kth-order reduced model is as following form:

$${R_r}^{(k)}(s)=\frac{W_0+W_{1}s+W_2{s}^2+\cdots +W_{k-1}{s}^{k-1}}{z_0+z_{1}s+z_2{s}^2+\cdots +z_{k-1}{s}^{k-1}+z_k{s}^k}$$

(7)

where $W_i(i=0,1,\dots k-1)$ , $z_i(i=0,1,\dots k)$ are constant coefficients, $z_k=1$ and $k\le n-1$ , the ${R_r}^{(k)}(s)$ can be obtained by defining the kth-order Routh convergent, that is an approximation of $H(s)$ . Suppose that

$${{\tilde{R}}_r}^{(k)}(s)=\frac{B_k(s)}{A_k(s)}$$

(8)

where

$$\{\begin{array}{c}{A}_k(s)={\alpha }_{k}s{A}_{k-1}(s)+A_{k-2}(s)\hfill \\ \multicolumn{1}{c}{B_k(s)={\alpha }_{k}s{B}_{k-1}(s)+B_{k-2}(s)+{\beta }_k}\end{array}k\ge 1$$

(9)

and . $A_{-1}(s)=0$ ., $A_0(s)=1$ , $A_1(s)={\alpha }_{1}s+1$ , $B_{-1}(s)=0$ , $B_0(s)=1$ , the coefficients ${\alpha }_k$ and ${\beta }_k$ can be obtained from the Tables 1 and 2 respectively, which are calculated by Eq. (10–15) as follows:

$${\alpha }_i=\frac{{d_0}^{i-1}}{{d_0}^i}i=1,2,3,\dots$$

(10)

$${d_j}^{i-1}=d_{i+j-1}i=1,2;j=0,2,4,\dots$$

(11)

$${d_j}^{i+1}=d_{j+2}^{i-1}-{\alpha }_i{d}_{j+2}^{i}i=1,2,3,\dots ;j=0,2,4,\dots$$

(12)

$${\beta }_i=\frac{N_0^i}{d_0^i}i=1,2,3,\dots$$

(13)

$$N_j^i=N_{i+j-1}i=1,2;j=0,2,4,\dots$$

(14)

$$N_j^{i+2}=N_{j+2}^i-{\beta }_i{d}_{j+2}^{i}i=1,2,3,\dots ;j=0,2,4,\dots$$

(15)

where the $N_i$ and $d_i$ are the coefficients of Eq. (6).

Table 1.

Computation of ${\alpha }_k$ Routh parameters.

	$\begin{array}{c}{d_0}^0=d_0\hfill \\ {d_0}^1=d_1\hfill \end{array}$	$\begin{array}{c}{d_2}^0=d_2\hfill \\ {d_2}^1=d_3\hfill \end{array}$	$\begin{array}{c}{d_4}^0=d_4\hfill \\ {d_4}^1=d_5\hfill \end{array}$	……
$\begin{array}{c}{\alpha }_1={d_0}^0/{d_0}^1\hfill \\ {\alpha }_2={d_0}^1/{d_0}^2\hfill \\ {\alpha }_3={d_0}^2/{d_0}^3\hfill \end{array}$	$\begin{array}{c}{d_0}^2={d_2}^0-{\alpha }_1{d_2}^1\hfill \\ {d_0}^3={d_2}^1-{\alpha }_2{d_2}^2\hfill \\ {d_0}^4={d_2}^2-{\alpha }_3{d_3}^2\hfill \end{array}$	$\begin{array}{c}{d_2}^2={d_4}^0-{\alpha }_1{d_4}^1\hfill \\ {d_2}^3={d_4}^1-{\alpha }_2{d_4}^2\hfill \\ {d_2}^4={d_4}^2-{\alpha }_3{d_4}^2\hfill \end{array}$	$\begin{array}{c}{d_4}^2={d_6}^0-{\alpha }_1{d_6}^1\hfill \\ {d_4}^3={d_6}^1-{\alpha }_2{d_6}^2\hfill \\ {d_4}^4={d_6}^2-{\alpha }_3{d_6}^3\hfill \end{array}$	… … …
⋮	⋮	⋮	⋮	⋮

Table 2.

Computation of ${\beta }_k$ Routh parameters.

	$\begin{array}{c}{N_0}^1=N_0\hfill \\ {N_0}^2=N_1\hfill \end{array}$	$\begin{array}{c}{N_2}^1=N_2\hfill \\ {N_2}^2=N_3\hfill \end{array}$	$\begin{array}{c}{N_4}^1=N_4\hfill \\ {N_4}^2=N_5\hfill \end{array}$	……
$\begin{array}{c}{\beta }_1={N_0}^1/{d_0}^1\hfill \\ {\beta }_2={N_0}^2/{d_0}^2\hfill \\ {\beta }_3={N_0}^3/{d_0}^3\hfill \end{array}$	$\begin{array}{c}{N_0}^3={N_2}^1-{\beta }_1{d_2}^1\hfill \\ {N_0}^4={N_2}^2-{\beta }_2{d_2}^2\hfill \\ {N_0}^5={N_2}^3-{\beta }_3{d_3}^2\hfill \end{array}$	$\begin{array}{c}{N_2}^3={N_4}^1-{\beta }_1{d_4}^1\hfill \\ {N_2}^4={N_4}^2-{\beta }_2{d_4}^2\hfill \\ {N_2}^5={N_4}^3-{\beta }_3{d_4}^2\hfill \end{array}$	$\begin{array}{c}{N_4}^3={N_6}^1-{\beta }_1{d_6}^1\hfill \\ {N_4}^4={N_6}^2-{\beta }_2{d_6}^2\hfill \\ {N_4}^5={N_6}^3-{\beta }_3{d_6}^2\hfill \end{array}$	… … …
⋮	⋮	⋮	⋮	⋮

Finally, the low-frequency kth-order Routh approximant of the $H(s)$ is acquired by the reciprocal transformation:

$${R_r}^{(k)}(s)=\frac{1}{s}{{\tilde{R}}_r}^{(k)}\left(\frac{1}{s}\right)$$

(16)

Modeling results and analysis of WRRP

The experimental sample data of output gas temperature TT103 and inlet flow TV101 were recorded in Figure 1. The relationship between output gas temperature TT103 and inlet flow TV101 is modeled, and the identification model is obtained as:

$$\begin{array}{c}G(s)=\hfill \\ \frac{-0.0003122s^4-0.001407s^3-0.00498s^2-0.007774s-0.05839}{s^5+2.897s^4+15.77s^3+19.4s^2+33.76s-0.1195}{e}^{-25s}\hfill \end{array}$$

(17)

Verify the identification model with the original data, the validation result is as shown in Figure 2, the original data is basically consistent with the data obtained by the ARX model, and the fit index reaches 94.03%. Therefore, it can be considered that the model is very accurate.

Figure 2.

Validation results of the attained ARX model of WRRP.

Since the ARX model is a high-order model, to facilitate the IMC controller design, the Routh approximation method is used to reduce the order and the reduced order model is as follows:

$$G_{p1}(s)=\frac{-0.002231s-0.02723}{s^2+15.77s-0.05573}{e}^{-25s}$$

(18)

It is necessary to verify $G_{p1}(s)$ and $G(s)$ , Figures 3 and 4 show the bode diagram and step response curve of $G_{p1}(s)$ and $G(s)$ respectively. From Figure 3, the frequency characteristics of amplitude and phase of $G_{p1}(s)$ are good consistent with $G(s)$ . Figure 4 shows that their step responses are basically consistent. Hence, it can be considered that the reduced order model $G_{p1}(s)$ is accurate and reliable.

Figure 3.

Bode diagram of $G_{p1}(s)$ and $G(s)$ .

Figure 4.

Step response of $G_{p1}(s)$ and $G(s)$ .

There is an unstable pole in the original model $G(s)$ , the unstable pole will still exist after model reduction. Eq. (19) and (20) are the zero pole forms of the original model $G(s)$ and the approximate model $G_{p1}(s)$ respectively.

$$\begin{array}{c}G(s)=\hfill \\ \frac{-0.00031218(s^2+2.312s+2.14)(s^2+2.194s+8.742)}{(s-0.003533)(s^2+1.398s+3.242)(s^2+1.503s+10.43)}{e}^{-25s}\hfill \end{array}$$

(19)

$$G_{p1}(s)=\frac{-0.0003995(s+0.751)}{(\mathrm{s}+1.738)(s-0.003533)}{e}^{-25s}$$

(20)

The controlled object is unstable with large time delay, and the expected control effect is hard to achieve by using PID control scheme. Therefore, an improved three degree of freedom IMC strategy with optimization algorithm is proposed for the model $G_{p1}(s)$ obtained by the system identification and model reduction method.

IMC controller design

Structure of improved IMC

An improved IMC structure is proposed in this paper for the controlled object $G_{p1}(s)$ as shown in Figure 5.

Figure 5.

Block diagram of the improved IMC structure.

In Figure 5, $r$ and $y$ are the system set value input and output respectively, $d_1$ and $d_2$ are the disturbance signals, $G_{p1}(s)$ is the controlled object, $Q_1(s)$ is the controller to stabilize $G_{p1}(s)$ , $G_p(s)$ is the generalized controlled object, $G_{pm}(s)$ is the process model, $Q_2(s)$ and $Q_3(s)$ are the IMC controller to be designed. The $G_p(s)$ can be formulated from Figure 5 as:

$$G_p\left(s\right)=\frac{G_{p1}\left(s\right)}{1+G_{p1}\left(s\right)Q_1\left(s\right)}$$

(21)

The transfer function of the system is formulated as:

$$\begin{array}{c}{G}_{\mathrm{yr}}\left(s\right)=\hfill \\ \frac{G_{p1}\left(s\right)Q_3\left(s\right)}{1+Q_2\left(s\right)[G_{p1}\left(s\right)-G_{\mathrm{pm}}\left(s\right)]+G_{p1}\left(s\right)Q_1\left(s\right)[1-G_{\mathrm{pm}}\left(s\right)Q_2\left(s\right)]}\hfill \end{array}$$

(22)

$$\begin{array}{l}{G}_{yd1}\left(s\right)=\\ \frac{G_{p1}\left(s\right)[1-G_{pm}\left(s\right)Q_2\left(s\right)]}{1+Q_2\left(s\right)[G_{p1}\left(s\right)-G_{pm}\left(s\right)]+G_{p1}\left(s\right)Q_1\left(s\right)[1-G_{pm}\left(s\right)Q_2\left(s\right)]}\end{array}$$

(23)

$$\begin{array}{c}{G}_{\mathrm{yd}2}\left(s\right)=\hfill \\ \frac{1+G_{p1}\left(s\right)Q_1\left(s\right)-G_{\mathrm{pm}}\left(s\right)Q_2\left(s\right)-G_{\mathrm{pm}}\left(s\right)Q_2\left(s\right)G_{p1}\left(s\right)Q_1\left(s\right)}{1+Q_2\left(s\right)[G_{p1}\left(s\right)-G_{\mathrm{pm}}\left(s\right)]+G_{p1}\left(s\right)Q_1\left(s\right)[1-G_{\mathrm{pm}}\left(s\right)Q_2\left(s\right)]}\hfill \end{array}$$

(24)

Substituting Eq. (21) into Eq. (22- 24) gives following forms:

$$G_{\mathrm{yr}}\left(s\right)=\frac{G_p\left(s\right)Q_3\left(s\right)}{1+Q_2\left(s\right)[G_p\left(s\right)-G_{\mathrm{pm}}\left(s\right)]}$$

(25)

$$G_{\mathrm{yd}1}\left(s\right)=\frac{G_p\left(s\right)[1-G_{\mathrm{pm}}\left(s\right)Q_2\left(s\right)]}{1+Q_2\left(s\right)[G_p\left(s\right)-G_{\mathrm{pm}}\left(s\right)]}$$

(26)

$$G_{\mathrm{yd}2}\left(s\right)=\frac{1-G_{\mathrm{pm}}\left(s\right)Q_2\left(s\right)}{1+Q_2\left(s\right)[G_p\left(s\right)-G_{\mathrm{pm}}\left(s\right)]}$$

(27)

Assuming $G_p\left(s\right)=G_{\mathrm{pm}}\left(s\right)$ , the Eq. (25–27) can be expressed as:

$$G_{\mathrm{yr}}\left(s\right)=G_p\left(s\right)Q_3\left(s\right)$$

(28)

$$G_{\mathrm{yd}1}\left(s\right)=G_p\left(s\right)[1-G_{\mathrm{pm}}\left(s\right)Q_2\left(s\right)]$$

(29)

$$G_{\mathrm{yd}2}\left(s\right)=1-G_{\mathrm{pm}}\left(s\right)Q_2\left(s\right)$$

(30)

The output of the system can be expressed as:

$$\begin{array}{c}\multicolumn{1}{c}{Y\left(s\right)=G_p\left(s\right)Q_3\left(s\right)R\left(s\right)+G_p\left(s\right)[1-G_{\mathrm{pm}}\left(s\right)Q_2\left(s\right)]D_1\left(s\right)}\\ +[1-G_{\mathrm{pm}}\left(s\right)Q_2\left(s\right)]D_2\left(s\right)\hfill \end{array}$$

(31)

According to section 2, the model of WRRP is a second-order unstable time-delay model, and its format is as follows:

$$G_{p1}(s)=\frac{K(T_{1}s+1)}{(T_{2}s+1)(T_{3}s-1)}{e}^{-\mathrm{Ls}}$$

(32)

where $K=-0.0483$ , $T_1=1.3316$ , $T_2=0.5754$ , $T_3=283.0456$ , $L=25$ . In the subsequent design, considering that the process of heating and cooling are symmetric and consistent, the minus sign in Eq. (32) is ignored, which has no impact on the performance index analysis of the control system, therefore we change the value of $K$ to a positive number, namely $K=0.0483$ .

Design the secondary loop controller Q₁(s)

To make the controlled object stable, the controller $Q_1(s)$ is designed first. Since the structure of proportional controller is relatively simple and easy to implement, $Q_1(s)$ is selected as proportional controller, whose value is a proportional gain $K_1$ . According to Eq. (21) and (32), the $G_p(s)$ can be formulated as:

$$G_p(s)=\frac{K(T_{1}s+1)e^{-\mathrm{Ls}}}{(T_{2}s+1)(T_{3}s-1)+K_{1}K(T_{1}s+1)e^{-\mathrm{Ls}}}$$

(33)

The time delay term in the denominator of the Eq. (33) is approximated by Taylor expansion: $e^{-Ls}\approx 1-Ls$ , we can get:

$$\begin{array}{c}{G}_p(s)=\hfill \\ \frac{K(T_{1}s+1)e^{-\mathrm{Ls}}}{(T_2{T}_3-T_{1}L{K}_{1}K)s^2+(T_2+T_3+K_{1}K{T}_1-K_1\mathrm{KL})s+K_{1}K-1}\hfill \end{array}$$

(34)

In order to make the system stable, according to the characteristic equation in Eq. (34), the value scope of $K_1$ that can be determined by the Hurwitz stability criterion as:

$$\frac{1}{K}<K_1<\min (\frac{T_2{T}_3}{T_1\mathrm{LK}},\frac{T_2+T_3}{(L-T_1)K})$$

(35)

By adjusting the value of $K_1$ in the above range, it can be ensured that the system has good closed-loop response characteristics.

Design the primary loop controller $Q_2(s)$

According to Eq. (31), $Q_2(s)$ is designed to reject disturbance, based on the IMC principle, $Q_2(s)$ can be designed as:

$$Q_2\left(s\right)=G_{\mathrm{pm}-}^{-1}\left(s\right)f_2(s)$$

(36)

where $G_{\mathrm{pm}}(s)$ is process model and $G_{\mathrm{pm}-}$ is selected as inverse of minimum phase component of $G_{\mathrm{pm}}(s)$ , $f_2(s)$ is a low pass filter which is the in following form:

$$f_2(s)=\frac{1}{{\lambda }_{2}s+1}$$

(37)

where ${\lambda }_2$ is the time constant of $f_2(s)$ that can be tuned.

From Eq. (34) and (36), the $Q_2(s)$ can be determined as:

$$\begin{array}{c}{Q}_2(s)=\hfill \\ \frac{(T_2{T}_3-T_{1}L{K}_{1}K)s^2+(T_2+T_3+K_{1}K{T}_1-K_1\mathrm{KL})s+K_{1}K-1}{K(T_{1}s+1)({\lambda }_{2}s+1)}\hfill \end{array}$$

(38)

Design the setpoint filter controller $Q_3(s)$

From Eq. (31), the controller $Q_3(s)$ is related to the setpoint following characteristics of the system. According to the IMC principle, $Q_3(s)$ can be designed as:

$$Q_3(s)=G_{p-}^{-1}\left(s\right)f_3\left(s\right)$$

(39)

where $G_{p-}$ is the inverse of minimum phase component of $G_p(s)$ , the filter $f_3(s)$ can be formulated as:

$$f_3(s)=\frac{1}{{\lambda }_{3}s+1}$$

(40)

where ${\lambda }_3$ is a tuning parameter.

According to Eq. (34) and (39), the $Q_2(s)$ can be determined as:

$$\begin{array}{c}{Q}_3(s)=\hfill \\ \frac{(T_2{T}_3-T_{1}L{K}_{1}K)s^2+(T_2+T_3+K_{1}K{T}_1-K_1\mathrm{KL})s+K_{1}K-1}{K(T_{1}s+1)({\lambda }_{3}s+1)}\hfill \end{array}$$

(41)

HHO algorithm for optimized tuning

In this paper, Harris hawks optimization algorithm is adopted for tuning controller parameters. The HHO is a naturally inspired algorithm proposed by Heidari et al. ²⁴ in 2019. It is derived from observing the behavior of Harris hawks when they prey on rabbits and using mathematical formulas to simulate the strategy of Harris hawks in catching prey under different mechanisms. In HHO, Harris hawks are the candidate solutions, and the best candidate solution in each step is considered as the intended prey or nearly the optimum. HHO includes two parts: exploration phase and exploitation phase.

Exploration phase

During initialization in this part Harris hawks randomly perch at a certain position in the search space (LB, UB), and wait to detect a prey based on two strategies, they selectively update their positions with probability $q$ according to the position of other individuals and prey, which is expressed as:

$$\begin{array}{c}X(t+1)=\hfill \\ \{\begin{array}{cc}{X}_{\mathrm{rand}}(t)-r_1|X_{\mathrm{rand}}(t)-2r_{2}X(t)|\hfill & \mathrm{q}\ge 0.5\hfill \\ \multicolumn{1}{c}{X_{\mathrm{rabbit}}(t)-X_m(t)-r_3(\mathrm{LB}+r_4(\mathrm{UB}-\mathrm{LB}))}& \mathrm{q}<0.5\hfill \end{array}\hfill \end{array}$$

(42)

where, $X(t+1)$ is the position of hawks in the next iteration, $X_{\mathrm{rabbit}}(t)$ is the position of the rabbit, $X_{\mathrm{rand}}(t)$ is the position of a hawk randomly selected in the current population, and $r_1$ , $r_2$ , $r_3$ , $r_4$ , $q$ are the random number between (0,1), and $X_m=\frac{1}{N}\sum _{i=1}^N{X}_i(t)$ is the average position of N individuals in the current population.

According to the escape energy of prey, HHO algorithm transforms from exploration phase to exploitation phase. The energy change of prey is expressed as:

$$E=2E_0\left(1-\frac{t}{T}\right)$$

(43)

where $E$ is the escaping energy of the prey, $E_0$ is the initial energy of the prey between (−1,1), $T$ is the maximum number of iterations, and $t$ is the current iteration. When $\left|E\right|\ge 1$ , HHO will perform exploration phase, and when $\left|E\right|<1$ , HHO will perform exploitation phase.

Exploitation phase

In this phase, HHO designs four ways to simulate the hunting strategy of hawks, $\theta$ indicates the chance of the prey in successfully escaping ( $\theta <0.5$ ) or not successfully escaping ( $\theta \ge 0.5$ ). At the same time, hawks will choose different besiege strategies by the energy of their prey, when $\left|E\right|\ge 0.5$ , they execute soft besiege, and when $\left|E\right|<0.5$ , the hard besiege is executed.

Soft besiege

When $\theta \ge 0.5$ and $\left|E\right|\ge 0.5$ , the prey has enough energy and attempt to escape by some random misleading jumps, but it fails to escape at last. The Harris hawks will gradually consume the energy of the prey and choose the best position to catch prey. The position is updated as following formula:

$$X(t+1)=\mathrm{\Delta }X(t)-E|J{X}_{\mathrm{rabbit}}(t)-X(t)|$$

(44)

$$\mathrm{\Delta }X(t)=X_{\mathrm{rabbit}}(t)-X(t)$$

(45)

where, $\mathrm{\Delta }X(t)$ represents the position deviation between the hawk and the prey in the t iteration, $J=2(1-r_5)$ is the random jump strength when the prey escapes, and $r_5$ is a random number between (0,1).

Hard besiege

When $\theta \ge 0.5$ and $\left|E\right|<0.5$ , the prey is exhausted. Harris hawks surround the prey and make a raid. In this case, the position of the hawks is updated as:

$$X(t+1)=X_{\mathrm{rabbit}}(t)-E|\mathrm{\Delta }X(t)|$$

(46)

Soft besiege with progressive rapid dives

When $\theta <0.5$ and $\left|E\right|\ge 0.5$ , the prey is full of energy, and the Harris hawks will choose to conduct a soft besiege before the raid. In this situation, the expression of the updated position is as:

$$\begin{array}{c}X(t+1)=\hfill \\ \{\begin{array}{ccc}Y:X(t)-E|J{X}_{\mathrm{rabbit}}(t)-X(t)|\hfill & \mathrm{if}\hfill & F(Y)<F(X(t))\hfill \\ \multicolumn{1}{c}{Z:Y+S\times \mathrm{LF}(D)}& \mathrm{if}\hfill & F(Z)<F(X(t))\hfill \end{array}\hfill \end{array}$$

(47)

where $D$ is the dimension of Harris hawks position vector, $S$ is the random vector of $1\times D$ dimension, and $\mathrm{LF}$ is the Levy flight function which is formulated as :

$$\mathrm{LF}(x)=0.01\times \frac{u\times \sigma }{{\left|v\right|}^{\frac{1}{\beta }}},\sigma =\left(\frac{\Gamma (1+\beta )\times \sin \left(\frac{\pi \beta }{2}\right)}{\Gamma (\frac{1+\beta }{2})\times \beta \times 2^{(\frac{\beta -1}{2})}}\right)$$

(48)

where $u$ , $v$ are random values within (0,1), $\beta$ is a default constant set to 1.5.

Hard besiege with progressive rapid dives

When $\theta <0.5$ and $\left|E\right|<0.5$ , the prey does not have enough energy to escape. Before raiding and killing the prey, Harris hawks will organize a hard besiege to reduce the distance between their average position and the prey position, and the position is updated as follows:

$$\begin{array}{c}X(t+1)=\hfill \\ \{\begin{array}{ccc}Y:X_{\mathrm{rabbit}}(t)-E|J{X}_{\mathrm{rabbit}}(t)-X_m(t)|\hfill & \mathrm{if}\hfill & F(Y)<F(X(t))\hfill \\ \multicolumn{1}{c}{Z:Y+S\times \mathrm{LF}(D)}& \mathrm{if}\hfill & F(Z)<F(X(t))\hfill \end{array}\hfill \end{array}$$

(49)

The flowchart of the HHO algorithm is shown in Figure 6.

Figure 6.

The flowchart of the HHO algorithm.

Results and discussion

Multi-objective optimization for WRRP

As there are several conflicting design objectives need to be simultaneously achieved in the control design problems. If these synthesis objectives are analytically represented as a set of design objective functions subject to the existing constraints, the synthesis problem could be formulated as a multi-objective optimization problem.^20,21 In section 3, an improved three degree of freedom IMC structure for the WRRP model is proposed. The controller parameters $K_1$ , ${\lambda }_2$ and ${\lambda }_3$ should be set appropriately to improve the performance and robustness of the control system. The objective function for evaluating the set-point tracking performance and actuator work is designed as follow:

$$J=w_1\mathrm{ISE}+w_2\mathrm{Cons}+w_3\mathrm{Ts}+w_4\mathrm{Tr}+w_5\mathrm{Mp}$$

(50)

where $w_1$ , $w_2$ , $w_3$ , $w_4$ , $w_5$ are the weights corresponding to each objective, $\mathrm{Ts}$ is the settling time, $\mathrm{Tr}$ is the rise time, and $\mathrm{Mp}$ is the peak overshoot percentage, $\mathrm{ISE}$ is the integral square error which is defined as:

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

(51)

where $e$ is the error. $\mathrm{Cons}$ is obtained by integrating the output of the controllers, it is defined as:

$$\mathrm{Cons}=\sum _{j=1}^3{w}_{2,j}{\int }_0^{\infty }{u}_j(t)\mathrm{dt}$$

(52)

where $u_j$ is $j$ th controller output, $w_{2,1}$ , $w_{2,2}$ , $w_{2,3}$ are the preference coefficients on each controller output, and they are set to 0.01 by default in subsequent simulations.

In this paper, different $w_i(i=1,2,\dots ,5)$ are selected and five objective functions are considered as follows:

$$J_1=0.20\mathrm{ISE}+0.20\mathrm{Ts}+0.80\mathrm{Mp}+0.20\mathrm{Tr}+0.60\mathrm{Cons}$$

(53)

$$J_2=0.60\mathrm{ISE}+0.40\mathrm{Ts}+0.50\mathrm{Mp}+0.00\mathrm{Tr}+0.40\mathrm{Cons}$$

(54)

$$J_3=0.40\mathrm{ISE}+0.80\mathrm{Ts}+0.50\mathrm{Mp}+0.40\mathrm{Tr}+0.10\mathrm{Cons}$$

(55)

$$J_4=0.20\mathrm{ISE}+0.60\mathrm{Ts}+0.20\mathrm{Mp}+0.80\mathrm{Tr}+0.80\mathrm{Cons}$$

(56)

$$J_5=0.80\mathrm{ISE}+0.10\mathrm{Ts}+0.80\mathrm{Mp}+0.80\mathrm{Tr}+0.20\mathrm{Cons}$$

(57)

Compare different optimization algorithms

In this paper, the controller parameters $K_1$ , ${\lambda }_2$ and ${\lambda }_3$ are tuned by HHO algorithm, the parameters of HHO algorithm and the boundaries of controller parameters are given in Table 3. As a comparison, genetic algorithm (GA),^18,19 particle swarm optimization (PSO),^16,17 Henry gas solubility algorithm (HGSO) ²⁵ and water cycle algorithm (WCA) ²⁶ are employed to tune controller parameters, and their related parameters are given in Table 4. For a fair comparison, the population size, the maximum iteration number, the boundaries of controller parameters are consistent with HHO algorithm in Table 3.

Table 3.

The parameters of HHO algorithm and the boundaries for controller parameters.

Parameter	Value
Population size	40
Maximum iteration number	30
Dimension of optimization problem	3
Lower bound of [ $K_1$ ${\lambda }_2$ ${\lambda }_3$ ]	[20.7039 0.01 0.01]
Upper bound of [ $K_1$ ${\lambda }_2$ ${\lambda }_3$ ]	[100.2897 100 100]

Table 4.

Compared algorithms and related parameter values.

Algorithm	Parameters settings
GA	pcorss = 0.4, pmutation = 0.2
PSO	c1 = 1.3, c2 = 1.7, vmax = 1, vmin = −1,
HGSO	l1 = 5e-03, l2 = 100, l3 = 1e-02, alpha = 1, beta = 1, M1 = 0.1, M2 = 0.2
WCA	Nsr = 6, dmax = 1e-5

Simulation results for step input

1. Select $J_1$ as the objective function, the simulation result is given in Table 5, the step response of the system optimized by different algorithms is shown in Figure 7.

2. Select $J_2$ as the objective function, the simulation result is given in Table 6, the step response of the system optimized by different algorithms is shown in Figure 8.

3. Select $J_3$ as the objective function, the simulation result is given in Table 7, the step response of the system optimized by different algorithms is shown in Figure 9.

4. Select $J_4$ as the objective function, the simulation result is given in Table 8, the step response of the system optimized by different algorithms is shown in Figure 10.

5. Select $J_5$ as the objective function, the simulation result is given in Table 9, the step response of the system optimized by different algorithms is shown in Figure 11.

Table 5.

Comparison of optimization results with $J_1$ .

Algorithm	Controller parameters and performance indices
	K1	λ2	λ3	ISE	Ts	Mp(%)	Tr	Cons
GA	21.8215	10.9743	22.3833	36.9838	86.30	0.0000	75.18	157.6522
PSO	20.7039	59.8118	6.4554	28.7601	48.86	0.7028	43.26	150.8196
HGSO	21.1108	56.1740	27.2999	39.5492	105.93	0.4049	88.75	150.1597
WCA	22.9124	2.7846	2.7846	26.8108	40.88	4.9980	34.57	172.4467
HHO	20.7039	50.2535	0.06381	25.1288	35.89	0.8431	26.15	152.1609

Figure 7.

Step response of the system with $J_1$ .

Table 6.

Comparison of optimization results with $J_2$ .

Algorithm	Controller parameters and performance indices
	K1	λ2	λ3	ISE	Ts	Mp(%)	Tr	Cons
GA	24.6867	80.4306	7.1300	29.2468	51.02	0.8112	45.42	113.5687
PSO	20.7039	69.5223	12.9515	32.2062	66.08	0.5784	57.16	149.4785
HGSO	23.1810	79.3967	7.5838	29.4463	51.90	0.6889	46.17	173.7835
WCA	22.9087	2.7800	2.7800	26.8081	40.87	4.9961	34.56	172.4130
HHO	20.7039	100	0.0649	25.1264	36.07	0.3714	26.17	152.1435

Figure 8.

Step response of the system with $J_2$ .

Table 7.

Comparison of optimization results with $J_3$ .

Algorithm	Controller parameters and performance indices
	K1	λ2	λ3	ISE	Ts	Mp(%)	Tr	Cons
GA	25.4152	36.2094	0.3930	25.4004	38.95	1.6686	27.72	196.8294
PSO	34.2639	73.1962	0.0100	25.2617	41.42	2.0100	32.44	282.4737
HGSO	45.6764	84.3527	2.9948	27.5082	46.39	3.5106	40.34	387.4251
WCA	22.9102	2.7819	2.7819	26.8092	40.87	4.9992	34.56	172.4268
HHO	20.7039	81.3905	0.0353	25.1075	35.95	0.4982	26.13	152.2502

Figure 9.

Step response of the system with $J_3$ .

Table 8.

Comparison of optimization results with $J_4$ .

Algorithm	Controller parameters and performance indices
	K1	λ2	λ3	ISE	Ts	Mp(%)	Tr	Cons
GA	24.7062	81.2106	2.4960	26.6413	41.12	0.8120	34.15	189.5329
PSO	20.7039	40.3318	14.1783	32.8405	69.08	0.8848	59.73	149.2286
HGSO	42.0086	59.4100	12.4988	32.8371	65.01	3.2751	57.68	346.7392
WCA	22.9104	2.7821	2.7821	26.8093	40.87	4.9986	34.57	172.4285
HHO	20.7039	33.5428	0.0543	25.1255	35.88	1.2199	26.14	152.1777

Figure 10.

Step response of the system with $J_4$ .

Table 9.

Comparison of optimization results with $J_5$ .

Algorithm	Controller parameters and performance indices
	K1	λ2	λ3	ISE	Ts	Mp(%)	Tr	Cons
GA	33.6594	92.4341	2.8740	27.0733	44.26	1.7194	37.25	274.1146
PSO	21.4940	32.8836	0.4524	25.3865	36.83	1.3238	27.14	159.5617
HGSO	22.8593	8.0661	1.5006	26.0536	38.76	3.6963	30.77	172.2692
WCA	22.9102	2.7819	2.7819	26.8092	40.87	4.9992	34.56	172.4268
HHO	20.7039	46.1860	0.0352	25.1104	35.86	0.9135	26.12	152.2577

Figure 11.

Step response of the system with $J_5$ .

As shown in Tables 5–9, it can be clearly found that HHO algorithm has significant advantages in $\mathrm{ISE}$ $\mathrm{Cons}$ , $\mathrm{Ts}$ and $\mathrm{Tr}$ among the five algorithms under the same condition. As for $\mathrm{Mp}$ , the mean and root mean square error (RMSE) of $\mathrm{Mp}$ of five algorithms in the above simulation are summarized in Table 10, the best mean of $\mathrm{Mp}$ is obtained by HHO and the RMSE of $\mathrm{Mp}$ obtained by HHO is ranked second, hence, HHO is considered to have the best performance in $\mathrm{Mp}$ . Figures 7–11 show the step response of the system optimized by different algorithms, which further demonstrate that the controlled outputs optimized by HHO track their desired values in a very satisfying manner. The value of the objective function obtained in the above five simulation are listed in Table 11. Compared with the other algorithms, the value of objective function obtained by HHO optimization is the smallest, the system optimized by HHO has better dynamic characteristics.

Table 10.

The mean and root mean square error of $\mathrm{Mp}$ of five algorithms.

Algorithm	Mean(%)	RMSE
GA	1.0022	1.4267
PSO	1.1000	1.1638
HGSO	2.3152	3.2484
WCA	4.9982	0.0026
HHO	0.7692	0.6790

Table 11.

Objective function values optimized by five algorithms in simulation.

Algorithm	The value of objective function J
	J1	J2	J3	J4	J5
GA	134.2841	113.5687	72.9254	209.1090	112.0830
PSO	115.2300	105.8363	85.4690	215.3600	78.6755
HGSO	137.2656	108.2856	104.7491	369.7638	86.7458
WCA	196.5031	103.8961	76.9860	196.4824	91.6671
HHO	108.7370	90.5469	64.7292	169.2097	75.0292

Simulation results for step input with disturbance

To test the disturbance rejection performance of the proposed method in this paper. Select $J_4$ as the objective function, and the parameter optimization results in Table 8 are used, the following simulation experiments are conducted.

1.  Set the simulation time to 1000 s, a unit step change in set point at time t = 0 and a disturbance $d_1$ with magnitude of 0.3 is added when t = 600 s and keep it for 40 s, namely $d_1=-0.3[e(t-600)-e(t-640)]$ , the simulation result is given in Table 12, and the system response is shown in Figure 12.

2.  Set the simulation time to 1000 s, a unit step change in set point at time t = 0 and a disturbance $d_2$ with magnitude of 0.1 is added when t = 600 s, namely $d_2=-0.1e(t-600)$ , the simulation result is given in Table 13, and the system response is shown in Figure 13.

Table 12.

Comparison of optimization results with $J_4$ when $d_1$ is added.

Algorithm	Performance indices and the value of $J_4$
	ISE	Ts	Mp(%)	Tr	Cons	$J_4$
GA	26.6423	41.12	0.8120	34.15	332.6395	323.5945
PSO	32.8419	69.08	0.8848	59.73	252.4283	297.9200
HGSO	32.8393	65.01	3.2751	57.68	662.4536	622.3358
WCA	26.8099	40.87	4.9986	34.57	297.6758	296.6803
HHO	25.1256	35.88	1.2199	26.14	225.3566	227.9944

Figure 12.

Step response of the system with $J_4$ when $d_1$ is added.

Table 13.

Comparison of optimization results with $J_4$ when $d_2$ is added.

Algorithm	Performance indices and the value of $J_4$
	ISE	Ts	Mp(%)	Tr	Cons	$J_4$
GA	27.7724	41.12	0.8120	34.15	332.0851	323.3770
PSO	33.5931	69.08	0.8848	59.73	253.2613	298.7366
HGSO	33.7629	65.01	3.2751	57.68	663.6567	623.4830
WCA	27.0894	40.87	4.9986	34.57	300.0726	298.6537
HHO	25.7917	35.88	1.2199	26.14	256.4654	261.0146

Figure 13.

Step response of the system with $J_4$ when $d_2$ is added.

As is shown in Figures 12 and 13, when disturbance is introduced, the control system optimized by HHO still has better response characteristics than GA, HGSO and PSO. Although WCA has faster response than HHO, overshoot occurs in the response process. According to the values of the objective function $J_4$ in Tables 12 and 13, we hold that HHO method outperforms the other methods in disturbance rejection performance.

Compare the proposed method with other control methods

In this simulation case, we compare the proposed method with the other methods proposed by Wang et al. ⁹ and Shamsuzzoha and Lee. ¹³ A new objective function is designed as follow:

$$J=\sum _{i=1}^N(w_1|y_i-r_i|+w_2{u}_{i,j}+w_3|u_{i,j}-u_{i-1,j}|)+w_4\mathrm{Mp}j=1,2,3$$

(58)

where, $N$ is the number of sampling points, $y$ , $r$ represents the output and the input respectively, $u_{i,j}$ is the $i$ th sampled value of the $j$ th controller output, $\mathrm{Mp}$ is the peak overshoot percentage, $w_1$ , $w_2$ , $w_3$ , $w_4$ are the weights corresponding to each objective and set as $w_1=5$ , $w_2=0.05$ , $w_3=0.05$ , $w_4=1$ . The proposed method is based the objective function $J$ in Eq. (58) to optimize controller parameters, and the parameters of HHO algorithm are given in Table 3. Apply the methods proposed by Wang et al. ⁹ and Shamsuzzoha and Lee ¹³ to WRRP, the related controller setting parameters are listed in Table 14. The unit step change in set point at time t = 0, a disturbance with negative magnitude of 0.2 is introduced when t=600s and keep it for 20 s, the simulation result is given in Table 15. Figure 14 shows the output response of the proposed controller and other controllers. It can be clearly observed that the proposed method shows faster tracking speed, smaller overshoot, as well as better disturbance rejection performance. Table 15 and Figure 14 demonstrate the superiority of the proposed control method over the other control methods.

Table 14.

The parameters of related controller.

Method	Controller Parameters
Wang et al. 9	KC = 0.2161, TI = 0.5854, TD = 0.0021
Shamsuzzoha and Lee 13	α2 = 8.517, α1 = 283, λ = 3.23, r = 2, γ = 0.8
Proposed	K1 = 20.3079, λ2 = 100, λ3 = 0.01232

Table 15.

Comparison result of proposed method and other methods.

Method	Performance Indices
	Tr	Ts	ISE	Mp(%)
Wang et al. 9	58.08	121.60	38.3584	19.30
Shamsuzzoha and Lee 13	75.11	82.80	31.5668	2.08
Proposed	26.12	36.08	25.0962	1.22

Figure 14.

System responses of proposed method and other methods.

Conclusions

In this paper, a three degree of freedom IMC method based on HHO algorithm and identification model is proposed for WRRP. The dynamic model of WRRP is determined by least square identification and Routh approximation method. Based on the model, a three degree of freedom IMC structure is designed. In the improved structure, three controllers are intended for stabilizing the controlled object, set point tracking as well as disturbance rejection separately. The multi-objective optimization of controller parameters based on HHO algorithm is conducted so that the comprehensive performance of the control system can achieve desired results. In the simulation, the recently proposed control methods and other optimization algorithms are employed as the comparison. The simulation results demonstrate the effectiveness and better performance produced by the proposed method for WRRP. As a future recommendation, the accuracy of the process model is a key to achieve satisfactory performance, the time delay term is approximated by Taylor expansion in this paper. We will continue to explore other approaches to obtain accurate process model such as data driven and the more complex responses such as sine response will be considered.