Design of Optimal PID Controller with 𝜺-Routh Stability for Different Processes

Abstract

This paper presents a design method of the optimal proportional-integral-derivative (PID) controller with 𝜀-Routh stability for different processes through Lyapunov approach. The optimal PID controller could be acquired by minimizing an augmented integral squared error (AISE) performance index which contains control error and at least first-order error derivative, or even may contain nth-order error derivative. The optimal control problem could be transformed into a nonlinear constraint optimization (NLCO) problem via Lyapunov theorems. Therefore, optimal PID controller could be obtained by solving NLCO problem through interior method or other optimization methods. The proposed method can be applied for different processes, and optimal PID controllers under various control weight matrices and 𝜀-Routh stability are presented for different processes. Control weight matrix and 𝜀-Routh stability’s effects on system performances are studied, and different tuning methods’ system performances are also discussed. 𝜀-Routh stability’s effects on disturbance rejection ability are investigated, and different tuning methods’ disturbances rejection ability is studied. To further illustrate the proposed method, experimental results of coupled water tank system (CWTS) under different set points are presented. Both simulation results and experiment results show the effectiveness and usefulness of the proposed method.

1. Introduction

PID controllers have been widely used in various industrial processes due to simple structure, few parameters, and easy implementation. In [1–43], many tuning methods such as empirical methods, robust methods, frequency methods, numerical methods, and intelligent methods had been presented to design PID controllers. In [3, 4], Ziegler-Nichols (Z-N) methods include continuous cycling (C-C) method and process reaction curve (PRC) method, and they can be called as first Z-N method and second Z-N method, respectively. The PRC method is especialy proposed for processes with time delay. In [5], the Cohen-Coon method is presented for first-order plus time delay (FOPTD) processes. In [6], the integral absolute error (IAE) performance index is used to design the optimal PID controller for FOPTD and second-order plus time delay (SOPTD) plants. In [7–13], the immune algorithm, modified genetic algorithm, modified ant colony algorithm, particle swarm algorithm (PSA), and modified PSA are used to design optimal PID controller for different systems. In [14], multiple tabu search algorithm (MTSA) is applied to design the optimal fuzzy-logic PID controller for load frequency control system. In [15], some formulas are developed to directly compute the optimal PID controller for certain processes. In [16], zero-pole assignment method is directly applied to compute optimal PID controller for continuous-time processes. In [17], direct search algorithm is used to design the optimal PID controller. In [18], the optimal PID controller is presented for FOPTD plants via dimensional analysis and numerical optimization techniques. In [19], a design method is presented to design the optimal PID controller via minimizing load disturbances and solving a constrained optimization problem in frequency domain. In [20, 21], design methods are proposed to develop optimal PI/PID controller via solving a nonconvex optimization problem in frequency domain. In [22], hybrid Taguchi genetic algorithm (HTGA) is used to design optimal PID controller for pulse width modulation (PWM) feedback systems through solving a constrained optimization problem. In [23–25], linear quadratic regulator (LQR) method is proposed to design the optimal PID controller for different systems and processes. In [26], the loop transfer recovery technique is applied to devise multivariable robust optimal PID controller for a nonminimum phase boiler system. In [27], optimal PID controller with 𝐻₂/𝐻_∞ norm constraint is presented for the superheated steam temperature system via PSA. In [28], a parameter space approach is proposed to design the minimum variance PID controller for linear time-invariant (LTI) plants. In [29], a design method is proposed to design PID controller by considering a transient performance item. In [30], a design method of multivariable fixed-structure PID controller is proposed for multivariable plants by solving an optimization problem with linear matrix inequality (LMI) constraint. In [31], iterative feedback tuning (IFT) method is used to design the optimal PID controller for achieving fast response to set point variations. In [32], a two-layer online autotuning algorithm is presented to design the optimal PID controller via general predictive control method. In [33], the LQR method is presented to design the optimal PID controller for robot arms. In [39], a design method of PID controller is proposed for SOPTD processes via internal model control (IMC) framework. In [40], analytic rules of the PID controller are presented for FOPTD and SOPTD processes via the IMC method and model reduction method. In [41], a performance assessment method of PID controller is proposed for integral processes via focusing on step responses of the set point change and load disturbance. In [42], a design method of the PID controller with feedforward compensator is proposed for the FOPTD process. In [43], a design method of the multiloop PI/PID controller is presented for square multiple-input and multiple-output (MIMO) process which is decomposed into a series of equivalent effective open-loop process (EOP) with reasonable gain margins and phase margins. In [44], an on-line relay feedback approach is used to identify, assess, and tune the PID controller. In [45], the pole placement method is used to design the PID controller. In [46], the multiobjective genetic algorithm is presented to design the optimal fractional-order PID controller. In [47], a design method is presented for optimal PID controller with dynamic performance constrained.

Design methods of the PID controller in [1–43] can be roughly divided into the following cases: the empirical methods, intelligent algorithm methods, numerical parametric methods, classic LQR/LQG methods, robust design methods, IMC methods, and pole placement methods. Many design methods of the PID controller are proposed for first-order and second-order processes. In this paper, a novel systematic design method is presented to design optimal PID controller with ε-Routh stability for different processes. The optimal PID controller with ε-Routh stability is proposed by minimizing augmented integral squared error (AISE) performance index. The proposed method is unlike classic optimal control design methods which need to solve the Riccati equation (RE). The proposed optimal control problem is equivalently transformed into the NLCO problem via Lyapunov theorems. Therefore, optimal controller parameters can be obtained by solving a NLCO problem. The proposed method is used to design optimal PID controller for a third-order, SOPTD, and FOPTD processes, respectively. The proposed method is used to control the liquid level of the CWTS under different set points. The simulation results and experiment results show effectiveness and usefulness of the proposed method. The proposed method is convenient to design the optimal PID controller which can provide high performance for control system.

This paper is arranged as follows: in Section 2, preliminary and optimal control problem statements are described. In Section 3, CWTS model is established, and systematic design procedures of the optimal PID controller are presented in detail. Characteristics and robustness of 𝜀-Routh stability are analyzed, and the interior method is recommended to solve the NLCO problem of the proposed optimal PID controller. In Section 4, optimal PID controllers and simulation results are presented for different processes. The control weight matrix and ε-Routh stability’s effects on system performances are studied. The system performances and disturbance rejection ability of different tuning methods are investigated. 𝜀-Routh stability’s effects on disturbance rejection ability are discussed as well. To further validate the proposed method, the proposed method is used to control the liquid level of the CWTS under different set points. In Section 5, study contents are reviewed, and some conclusions are made.

2. Preliminary and Optimal Control Problem Statements

Preliminary. Based on Lyapunov theorems, a systematic tuning method is proposed to design the optimal PID controller with ε-Routh stability. The following theorems and ε-Routh stability definition are utilized to design the optimal PID controller.

Theorem 1 (see [1, 2]). Sufficient and necessary conditions for linear time-invariant system 𝑥̇=𝐴𝑥 asymptotically stable in a large scope are that: for any given matrix 𝑄=𝑄^𝑇>0, there is a matrix 𝑃=𝑃^𝑇>0 that satisfies the Lyapunov algebraic equation (LAE): 𝐴^𝑇𝑃+𝑃𝐴=−𝑄.

Theorem 2 (see [1, 2]). The linear time-invariant (LTI) system 𝑥̇=𝐴𝑥, and 𝑥(0)=𝑥₀ is asymptotically stable, for any given matrix 𝑄=𝑄^𝑇>0, then the performance index 𝐽 will have equivalent relationships: 𝐽=^∞₀𝑥^𝑇𝑄𝑥𝑑𝑡=𝑥^𝑇(0)𝑃𝑥(0),𝐴^𝑇𝑃+𝑃𝐴=−𝑄. (1)

Definition of 𝜀-Routh stability. 𝐷(𝑠) is characteristic polynomial of the linear time-invariant control system. For any nonnegative real 𝜀, if all elements of the Routh array’s first column of characteristic polynomial 𝐷(𝑠)=𝑎₀𝑠^𝑛+𝑎₁𝑠^𝑛−1+⋯+𝑎_𝑛−1𝑠+𝑎_𝑛 can satisfy the following inequality relationships, then the linear time-invariant control system is considered to possess the 𝜀-Routh stability: 𝑅_1,1𝑎₀,𝑎₁,…,𝑎_𝑛>𝜀,𝑅_2,1𝑎₀,𝑎₁,…,𝑎_𝑛𝑅>𝜀,…,_𝑛,1𝑎₀,𝑎₁,…,𝑎_𝑛>𝜀,𝑅_𝑛+1,1𝑎₀,𝑎₁,…,𝑎_𝑛(),>𝜀𝜀≥0 (2)

where 𝑅_1,1,𝑅_2,1,𝑅_3,1,…𝑅_𝑛,1,  and  𝑅_𝑛+1,1 are the elements of the Routh array’s first column. It calls Routh array’s first column as the Routh column, and the general computation formulas of the Routh column can be founded in [48, 49]. It could be known that the ε-Routh stability degenerates into the Routh-Hurwitz stable when the 𝜀 is zero (𝜀=0). The motivation for the definition of 𝜀-Routh stability comes from stability analysis of linear systems. As known, if any one element in Routh column of polynomial 𝐷(𝑠) is zero or negative real, then linear control systems may bring oscillation behaviors or unstable behaviors. Each element in the Routh column of the polynomial 𝐷(𝑠) is much closer to the zero 0⁺, at least it is more possible to bring oscillation behaviors. Consequently, the definition of 𝜀-Routh stability is motivated to describe how far away each Routh column from the zero 0⁺. It can be known that 𝜀-Routh stability will imply system robustness in some degree.

Optimal Control Problem Statements. We assume that the processes can be described by general transfer function 𝐺_𝑝(𝑠): 𝐺_𝑝(𝑠)=𝑏₁𝑠^𝑤+𝑏₂𝑠^𝑤−1+⋯+𝑏_𝑤𝑠+𝑏_𝑤+1𝑠^𝑛+𝑎₁𝑠^𝑛−1+⋯+𝑎_𝑛−1𝑠+𝑎_𝑛=∑^𝑤+1_𝑗=1𝑏_𝑗𝑠^{𝑤+1−𝑗}𝑠^𝑛+∑^𝑛_𝑖=1𝑎_𝑖𝑠^𝑛−𝑖=𝑁(𝑠)𝑀(𝑠), (3) where 𝑏₁,…,𝑏_𝑤+1,𝑎₁,…,  and  𝑎_𝑛 are system parameters, and 𝑀(𝑠) and 𝑁(𝑠) are the denominator and numerator, respectively. The denominator and numerator’s orders are the deg{𝑀(𝑠)}=𝑛 and deg{𝑁(𝑠)}=𝑤, respectively. For most practical processes, the number of poles 𝑛 is greater than the number of zeros 𝑤(𝑛>𝑤), and 𝑛 and 𝑤 are the positive integers. It could be known that for the all-pole plants 𝑏₁,𝑏₂,…,𝑏_𝑤=0, additional for the type I plants 𝑎_𝑛=0, and for the type II plants 𝑎_𝑛−1=𝑎_𝑛=0, and so on. The diagram of the control system is shown in Figure 1, and the control system can be described by the following equation: 𝑢(𝑠)=𝐺_𝑐(𝑠)𝑒(𝑠),𝐺_𝑐(𝑠)=𝑘_𝑝𝑠+𝑘_𝑖+𝑘_𝑑𝑠²𝑠,𝑟(𝑠)−𝑦(𝑠)≜𝑒(𝑠),𝑦(𝑠)=𝐺_𝑝(𝑠)𝑢(𝑠), (4) where 𝑒(𝑡) is control error, 𝑟(𝑡) and 𝑦(𝑡) are reference command and system output, respectively, 𝐺_𝑐(𝑠) is PID controller, and 𝑘_𝑝, 𝑘_𝑖, and 𝑘_𝑑 are controller parameters. Equation (4) can be equally expressed as 𝑟(𝑠)(𝑠)(𝑠)(𝑠)=𝑁(𝑠)𝑢(𝑠)=𝑦+𝑒,𝑦𝑀(𝑠),𝑢(𝑠)=𝑘_𝑝𝑠+𝑘_𝑖+𝑘_𝑑𝑠²𝑒(𝑠)𝑠. (5) Therefore, error transfer function can be obtained as 𝑒(𝑠)𝑟(𝑠)=11+𝐺_𝑝(𝑠)𝐺_𝑐(𝑠(𝑠)𝑒(𝑠)(𝑠)𝑟(𝑠),𝐷(𝑠)(𝑠)(𝑠)𝑘)⟺𝐷=𝑠𝑀=𝑠𝑀+𝑁_𝑝𝑠+𝑘_𝑖+𝑘_𝑑𝑠²,deg{𝐷(𝑠)}=𝑛+1, (6) where 𝐷(𝑠) is the control system’s characteristic polynomial. The error differential equation can be obtained by inversing Laplace transform of (6): 𝛿𝑀(𝛿)+𝑁(𝛿)𝑘_𝑝𝛿+𝑘_𝑖+𝑘_𝑑𝛿²𝑒(𝑡)=𝛿𝑀(𝛿)𝑟(𝑡), (7) where 𝛿 is the differential operator. When reference command 𝑟(𝑡) is constant signal or zero signal (𝑟(𝑡)=0), or even the piecewise constant signal, the error differential equation will yield 𝛿^𝑛+1+^𝑛_𝑖=1𝑎_𝑖𝛿^{𝑛+1−𝑖}+𝑘_𝑝𝛿+𝑘_𝑖+𝑘_𝑑𝛿^2𝑤+1_𝑗=1𝑏_𝑗𝛿^{𝑤+1−𝑗}𝑒(𝑡)=0. (8) Considering the polynomial differential operator 𝑀(𝛿) and 𝑁(𝛿), error differential equation (8) could be expressed as 𝛿^𝑛+1+^𝑛_𝑖=1𝑎_𝑖𝛿^{𝑛+1−𝑖}+^𝑤+1_𝑗=1𝑏_𝑗𝑘_𝑝𝛿^{𝑤+2−𝑗}+^𝑤+1_𝑗=1𝑏_𝑗𝑘_𝑖𝛿^{𝑤+1−𝑗}+^𝑤+1_𝑗=1𝑏_𝑗𝑘_𝑑𝛿^{𝑤+3−𝑗}𝑒(𝑡)=0. (9) Therefore, (9) can be equally expressed as 𝑒(𝑡)^(𝑛+1)+𝑎₁𝑒(𝑡)^(𝑛)+𝑎₂𝑒(𝑡)^(𝑛−1)+⋯+𝑎_{𝑛−𝑤−2}𝑒(𝑡)^(𝑤+3)+𝑎_{𝑛−𝑤−1}+𝑏₁𝑘_𝑑𝑒(𝑡)^(𝑤+2)+𝑎_𝑛−𝑤+𝑏₁𝑘_𝑝+𝑏₂𝑘_𝑑𝑒(𝑡)^(𝑤+1)+^𝑤−1_𝑗=1𝑎_{𝑛−𝑤+𝑗}+𝑏_𝑗+1𝑘_𝑝+𝑏_𝑗𝑘_𝑖+𝑏_𝑗+2𝑘_𝑑𝑒(𝑡)^{(𝑤+1−𝑗)}+𝑎_𝑛+𝑏_𝑤+1𝑘_𝑝+𝑏_𝑤𝑘_𝑖𝑒̇(𝑡)+𝑏_𝑤+1𝑘_𝑖𝑒(𝑡)=0. (10) For the symbol simplicity, the error differential equation (10) is equally expressed as 𝑒(𝑡)^(𝑛+1)+𝑎₁𝑒(𝑡)^(𝑛)+𝑎₂𝑒(𝑡)^(𝑛−1)+⋯+𝑎_{𝑛−𝑤−2}𝑒(𝑡)^(𝑤+3)+𝑚₁𝑒(𝑡)^(𝑤+2)+𝑚₂𝑒(𝑡)^(𝑤+1)+^𝑤−1_𝑗=1𝑚_𝑗+2𝑒(𝑡)^{(𝑤+1−𝑗)}+𝑚_𝑤+2𝑒̇(𝑡)+𝑚_𝑤+3𝑒(𝑡)=0, (11) where 𝑚_𝑖(𝑖=1,2,…,𝑤+3) is the decision parameter. Parameter 𝑚_𝑖(𝑖=1,2,…,𝑤+3) is obtained as 𝑚₁=𝑎_{𝑛−𝑤−1}+𝑘_𝑑𝑏₁,𝑚₂=𝑎_𝑛−𝑤+𝑏₁𝑘_𝑝+𝑏₂𝑘_𝑑𝑚_𝑗+2=𝑎_{𝑛−𝑤+𝑗}+𝑏_𝑗+1𝑘_𝑝+𝑏_𝑗𝑘_𝑖+𝑏_𝑗+2𝑘_𝑑𝑚𝑗=1,2,…,𝑤−1_𝑤+2=𝑎_𝑛+𝑏_𝑤+1𝑘_𝑝+𝑏_𝑤𝑘_𝑖,𝑚_𝑤+3=𝑏_𝑤+1𝑘_𝑖. (12) Thus, the corresponding state-space model of the differential equation (11) can be obtained: 𝑥̇=𝐴𝑥, (13) where 𝑥 is a state vector, 𝐴 is a state matrix, 𝐶_1×𝑛 is a parameter vector, 𝑂_𝑛×1 is a zero vector, and 𝐼_𝑛×𝑛 is an unit matrix. They are obtained as 𝑂𝐴=^𝑛×1𝐼_𝑛×𝑛−𝑚_𝑤+3𝐶_1×𝑛,𝑥=𝑒(𝑡)𝑒̇(𝑡)𝑒̈(𝑡)⋯𝑒(𝑡)^(𝑛)𝑇,𝐶_1×𝑛=−𝑚_𝑤+2⋯−𝑚₁−𝑎_{𝑛−𝑤−2}⋯−𝑎₂−𝑎₁. (14) To obtain the optimal PID controller, it needs to select suitable optimized performance indices. The performance indices such as the integral squared error (ISE), integral absolute error (IAE), integral time absolute error (ITAE), and integral time squared error (ITSE) are often used for criteria when designing optimal PID controller. The precision and steady-state property of control system are directly reflected by control error. In this paper, augmented integral squared error (AISE) in [34, 36] is used to design an optimal PID controller. Hence, optimal PID controller not only stabilizes the control system with 𝜀-Routh stability, but also minimizes the performance index. Therefore, the proposed optimal control problem can be stated as (i) the optimal PID controller stabilizes the control system with 𝜀-Routh stability; (ii) the optimal PID controller minimizes the given performance index. Eventually, the proposed optimal control problem can be depicted by the following equation: 𝐽^∗𝑘^∗_𝑝,𝑘^∗_𝑖,𝑘^∗_𝑑=min_{𝐺∗𝑐(𝑠)}𝐽𝑚₁,…,𝑚_𝑤+3=min^∞₀𝑥^𝑇(𝑡)𝑄𝑥(𝑡)𝑑𝑡=min^𝑛_𝑖=0^∞₀𝑞_{𝑖+1,𝑖+1}𝑒(𝑡)^(𝑖)2𝑑𝑡subjectto(s.t)∶(1)̇(2)𝑥=𝐴𝑥;Reeig(𝐴){()}(3)<0,Re𝑠∣det𝑠𝐼−𝐴=0<0𝜀≥0,𝑅_1,1𝑚₁,…,𝑚_𝑤+3𝑅>𝜀,_2,1𝑚₁,…,𝑚_𝑤+3𝑅>𝜀,…,_𝑛+1,1𝑚₁,…,𝑚_𝑤+3>𝜀,𝑅_𝑛+2,1𝑚₁,…,𝑚_𝑤+3(4)𝑚>𝜀;₄=𝑓₁𝑚₁,𝑚₂,𝑚₃,…,𝑚_𝑤+3=𝑓_𝑤𝑚₁,𝑚₂,𝑚₃;(5)0<𝑘_𝑝≤𝜏⁰_𝑝,0<𝑘_𝑖≤𝜏⁰_𝑖,0<𝑘_𝑑≤𝜏⁰_𝑑, (15) where 𝑄 is a control weight matrix, and at least a positive semidefinite real symmetric matrix, and 𝑓₁,…,𝑓_𝑤 is the linear function of the controller parameters 𝑘_𝑝, 𝑘_𝑖, and 𝑘_𝑑, and 𝑘^∗_𝑝, 𝑘^∗_𝑖, and 𝑘^∗_𝑑 are optimal controller parameters. Constraint (1) is the control system’s state constraint, constraint (2) guarantees the control system stable in a large scope, and constraint (3) guarantees the control system with the 𝜀-Routh stability. However, there are three independent decision variables among decision variables 𝑚_𝑖(𝑖=1,2,…,𝑤+3); the rest decision variables are interrelated variables. Constraint (4) reflects decision variables’ inner relationships. Constraint (5) can assure the existence of feasible solutions and can reduce the searching space. Constraint (5) can also assure PID controller to be physically achieved. In some actual control system, PID controller parameters are limited in certain range. Thus, the controller parameters cannot be too large. The parameters 𝜏⁰_𝑝, 𝜏⁰_𝑖, and 𝜏⁰_𝑑 are the scope of controller parameters.

Figure 1 

Feedback control system with optimal PID controller.

3. CWTS Model and Proposed Design Procedures

The proposed design method is used to design optimal PID controller which is employed to control the liquid level of the CWTS. The experiment setup is shown in Figure 2. It can use the experimental setup to construct different physical processes which are shown in Figure 3.

Figure 2 

Experimental setup of CWTS.

Figure 3 

The liquid level process of the CWTS.

The physical process in Figure 3 contains the interacting behavior because the liquid level ℎ₂ of tank-2 depends on the liquid level ℎ₁ of tank-1. In this paper, it does not consider nonlinearities of the coupled water tank system. Thus, the water tank system is considered as the linear process. Based on the volume balance principle, the liquid level equation of tank-1 can be obtained: 𝑄_𝑖𝑐−𝑄_1𝑐−𝑄₁₂=𝐶₁𝑑ℎ₁,𝑑𝑡 (16) where 𝑄_𝑖𝑐 is the flow of control valve, 𝑄_1𝑐 is the flow of valve-1, 𝑄₁₂ is the flow of valve-3, 𝐶₁ is the hydraulic capacity of tank-1, and ℎ₁ is the liquid level of tank-1. Likewise, the liquid level equation of tank-2 can be obtained: 𝑄₁₂−𝑄_2𝑐=𝐶₂𝑑ℎ₂,𝑑𝑡 (17) where 𝑄_2𝑐 is the flow of valve-2, 𝐶₂ is the hydraulic capacity of tank-2, and ℎ₂ is the liquid level of tank-2. The water flows equation of valve-1, valve-2, and valve-3 can be obtained: 𝑄_1𝑐=ℎ₁𝑅₁,𝑄_2𝑐=ℎ₂𝑅₂,𝑄₁₂=ℎ₁−ℎ₂𝑅₁₂, (18) where 𝑅₁ and 𝑅₂ are the liquid resistance of valve-1 and valve-2, respectively, and 𝑅₁₂ is the liquid resistance of valve-3. Thus, the liquid level model of water tank-1 is acquired: 𝐺_ℎ(𝑠)=𝐻₁(𝑠)𝑄_𝑖𝑐(𝑠)=𝑏₁𝑠+𝑏₂𝑎₀𝑠²+𝑎₁𝑠+𝑎₂=𝑏₁𝑠/𝑎₀+𝑏₂/𝑎₀𝑠²+𝑎₁𝑠/𝑎₀+𝑎₂/𝑎₀=𝑁₁(𝑠)𝑀₁(𝑠), (19) where 𝐺_ℎ(𝑠) is the liquid level transfer function of tank-1. The transfer function 𝐺_ℎ(𝑠) shows that physical process is the second-order process with a negative zero. System parameters are obtained as 𝑏₁=𝑅₁𝑅₁₂𝑅₂𝐶₂,𝑏₂=𝑅₁𝑅₁₂+𝑅₂,𝑎₀=𝑅₁𝐶₁𝑅₁₂𝑅₂𝐶₂,𝑎₁=𝑅₁𝐶₁𝑅₂+𝑅₁𝐶₁𝑅₁₂+𝑅₁𝑅₂𝐶₂+𝑅₁₂𝑅₂𝐶₂,𝑎₂=𝑅₁+𝑅₁₂+𝑅₂. (20)(21)

3.1. Proposed Design Procedures

Without losing generalities, a third-order controlled process is considered in the unity feedback system. The process (𝑛=3,𝑤=2) and PID controller have forms, respectively, 𝐺_𝑝(𝑠)=𝑏₁𝑠²+𝑏₂𝑠+𝑏₃𝑎₀𝑠³+𝑎₁𝑠²+𝑎₂𝑠+𝑎₃=𝑁(𝑠)𝑀(𝑠),𝐺_𝑐(𝑠)=𝑘_𝑝+𝑘_𝑖𝑠+𝑘_𝑑𝑠. (22) Thus, open-looped transfer function of control system yields 𝐺(𝑠)=𝐺_𝑐(𝑠)𝐺_𝑝(𝑠)=𝑘_𝑝𝑠+𝑘_𝑖+𝑘_𝑑𝑠²𝑁(𝑠)(𝑠).𝑠𝑀 (23) Then, the error equation and of the control system are obtained: 𝑒(𝑠)𝑟(𝑠)=11+𝐺_𝑝(𝑠)𝐺_𝑐(𝑠(𝑠)𝑒(𝑠)(𝑠)𝑟(𝑠),𝐷(𝑠)=(𝑠)(𝑠)𝑘)⟺𝐷=𝑠𝑀𝑠𝑀+𝑁_𝑝𝑠+𝑘_𝑖+𝑘_𝑑𝑠²,deg{𝐷(𝑠)}=4, (24) where 𝐷(𝑠) is the control system’s characteristic polynomial. The error differential equation can be obtained by inversing Laplace transform of the error equation: 𝛿𝑀(𝛿)+𝑁(𝛿)𝑘_𝑝𝛿+𝑘_𝑖+𝑘_𝑑𝛿²𝑒(𝑡)=𝛿𝑀(𝛿)𝑟(𝑡), (25) where δ is a differential operator. When reference command 𝑟(𝑡) is the constant signal, or zero signal (𝑟(𝑡)=0), or even piecewise constant signal, then the error differential equation will yield 𝑎₀+𝑘_𝑑𝑏₁𝑒(𝑡)⁽⁴⁾+𝑎₁+𝑘_𝑝𝑏₁+𝑘_𝑑𝑏₂𝑒⃛(𝑡)+𝑎₂+𝑘_𝑖𝑏₁+𝑘_𝑝𝑏₂+𝑘_𝑑𝑏₃𝑒̈(𝑡)+𝑎₃+𝑘_𝑖𝑏₂+𝑘_𝑝𝑏₃𝑒̇(𝑡)+𝑘_𝑖𝑏₃𝑒(𝑡)=0. (26) For symbol simplicity, error differential equation is equally expressed as 𝑚₁𝑒(𝑡)⁽⁴⁾+𝑚₂𝑒⃛(𝑡)+𝑚₃𝑒̈(𝑡)+𝑚₄𝑒̇(𝑡)+𝑚₅𝑒(𝑡)=0, (27) where 𝑚_𝑖(𝑖=1,…,5) is the decision parameter. The decision parameter yields 𝑚₁=𝑎₀+𝑘_𝑑𝑏₁,𝑚₂=𝑎₁+𝑘_𝑝𝑏₁+𝑘_𝑑𝑏₂,𝑚₅=𝑘_𝑖𝑏₃,𝑚₃=𝑎₂+𝑘_𝑖𝑏₁+𝑘_𝑝𝑏₂+𝑘_𝑑𝑏₃,𝑚₄=𝑎₃+𝑘_𝑖𝑏₂+𝑘_𝑝𝑏₃. (28) Thus, the corresponding state-space model of error differential equation (27) can be obtained: ̇(𝑒(𝑡)𝑒̇(𝑡)𝑒̈(𝑡)𝑒⃛(𝑡𝑥=𝐴𝑥𝑥=))^𝑇𝑂,𝐴=^3×1𝐼_3×3−𝑚¹₄𝐶_1×3,𝐶_1×3=−𝑚¹₃−𝑚¹₂−𝑚¹₁, (29) where 𝑥 is state vector, 𝐴 is state matrix, 𝑚¹_𝑖 is the middle decision variables, and the 𝑚¹_𝑖=𝑚_𝑖+1/𝑚₁(𝑖=1,2,3,4). Therefore, the optimal control problem is formulated as 𝐽^∗𝑘^∗_𝑝,𝑘^∗_𝑖,𝑘^∗_𝑑=min_𝐺∗𝑐𝐽𝑚₁,…,𝑚₅=min^∞₀𝑥^𝑇(𝑡)𝑄𝑥(𝑡)𝑑𝑡=min³_𝑖=0^∞₀𝑞_{𝑖+1,𝑖+1}𝑒(𝑡)^(𝑖)2𝑑𝑡s.t∶(1)̇(2)𝑥=𝐴𝑥;Reeig(𝐴)<0,or{()}(3)Re𝑠∣det𝑠𝐼−𝐴=0<0;𝜀≥0,𝑅_1,1𝑚₁,…,𝑚_𝑤+3𝑅>𝜀,_2,1𝑚₁,…,𝑚_𝑤+3𝑅>𝜀,_3,1𝑚₁,…,𝑚_𝑤+3>𝜀,𝑅_4,1𝑚₁,…,𝑚_𝑤+3𝑅>𝜀,_3,1𝑚₁,…,𝑚_𝑤+3>𝜀,𝑅_4,1𝑚₁,…,𝑚_𝑤+3𝑅>𝜀,_5,1𝑚₁,…,𝑚_𝑤+3(4)𝑚>𝜀;₄=𝑓₁𝑚₁,𝑚₂,𝑚₃,𝑚₅=𝑓₂𝑚₁,𝑚₂,𝑚₃;(5)0<𝑘_𝑝≤𝜏⁰_𝑝,0<𝑘_𝑖≤𝜏⁰_𝑖,0<𝑘_𝑑≤𝜏⁰_𝑑. (30) If optimal control problem (30) is solved, then the proposed optimal PID controller could be obtained. The constraint (2) assures that the control systems asymptotically stable in a large scope. The determinant is det(𝑠𝐼−𝐴)=𝑠⁴+𝑚¹₁𝑠³+𝑚¹₂𝑠²+𝑚¹₃𝑠+𝑚¹₄. Based on LTI stablility theory, constraint (2) is equivalently transformed into the following inequality:Reeig(𝐴)<0,orRe{𝑠∣det(𝑠𝐼−𝐴)=0}<0⟺𝑚¹₁>0,𝑚¹₁𝑚¹₂−𝑚¹₃𝑚>0,₁₁𝑚¹₂𝑚¹₃−𝑚¹₃²−𝑚¹₄𝑚¹₁²>0,𝑚¹₄>0. (31) The constraint (3) assures that the control system possesses the 𝜀-Routh stability. The characteristic polynomial is 𝐷(𝑠)=𝑚₁𝑠⁴+𝑚₂𝑠³+𝑚₃𝑠²+𝑚₄𝑠+𝑚₅. According to 𝜀-Routh stability definition, the constraint (3) is formulated by the following inequalities: 𝜀≥0,𝑅_1,1𝑚₁,…,𝑚₅>𝜀,𝑅_2,1𝑚₁,…,𝑚₅𝑅>𝜀,_3,1𝑚₁,…,𝑚₅𝑅>𝜀,_4,1𝑚₁,…,𝑚₅>𝜀,𝑅_5,1𝑚₁,…,𝑚₅>𝜀⟺𝑚₁>𝜀,𝑚₂>𝜀,𝑚₅𝑚>𝜀,₁𝑚₄−𝑚₂𝑚₃>𝜀𝑚₂,𝑚₂²𝑚₅−𝑚₁𝑚₄²−𝑚₂𝑚₃𝑚₄>𝜀𝑚₁𝑚₄−𝑚₂𝑚₃. (32) The inequalities (31) and (32) assure control system asymptotic stability and 𝜀-Routh stability, respectively. Therefore, the control system is asymptotically stable in a large scope; based on the Lyapunov Theorems 1 and 2, performance index 𝐽 is equivalently transformed into the performance index 𝐽^∗=min^∞₀𝑥^𝑇𝑄𝑥𝑑𝑡=min𝑥^𝑇(0)(0),𝐴𝑃𝑥_𝑇𝑃+𝑃𝐴=−𝑄, (33) where 𝑃 is a positive definite real symmetric matrix, and matrix 𝑃 meets the LAE: 𝐴^𝑇𝑃+𝑃𝐴=−𝑄. The performance index 𝐽 is determined by the matrix 𝑃 and initial state. Therefore, the proposed optimal control problem (30) can be equivalently transformed into the following NLCO problem: 𝐽^∗𝑘^∗_𝑝,𝑘^∗_𝑖,𝑘^∗_𝑑=min_{𝐺∗𝑐(𝑠)}𝐽𝑚₁,…,𝑚₅=min𝑥(0)^𝑇𝑃𝑥(0)s.t∶(1)𝐴^𝑇(2)𝑚𝑃+𝑃𝐴=−𝑄;₁₁>0,𝑚¹₄>0,𝑚¹₁𝑚¹₂−𝑚¹₃𝑚>0,₁₁𝑚¹₂𝑚¹₃−𝑚¹₃²−𝑚¹₄𝑚¹₁²(3)𝑚>0;₁>𝜀,𝑚₂>𝜀,𝑚₅>𝜀,𝑚₁𝑚₄−𝑚₂𝑚₃>𝜀𝑚₂,𝑚₂²𝑚₅−𝑚₁𝑚₄²−𝑚₂𝑚₃𝑚₄>𝜀𝑚₁𝑚₄−𝑚₂𝑚₃;(4)𝑚₄=𝑓₁𝑚₁,𝑚₂,𝑚₃,𝑚₅=𝑓₂𝑚₁,𝑚₂,𝑚₃;(5)0<𝑘_𝑝≤𝜏⁰_𝑝,0<𝑘_𝑖≤𝜏⁰_𝑖,0<𝑘_𝑑≤𝜏⁰_𝑑. (34) The proposed optimal control problem (30) is equivalently transformed into a NLCO problem via applying the Lyapunov theorems. The optimal PID controller with 𝜀-Routh stability can be obtained by solving NLCO problem (34). The design problem of optimal PID controller with 𝜀-Routh stability is deduced into the issue that: for the given control weight matrix and initial states, we pursue a suitable matrix 𝑃 to minimize performance index J. We assume that the control system is static at the beginning, then a unit step command 𝑟(𝑡)=1(𝑡) is inputted into the system, and the initial states can be obtained: 𝑥₁(0)=1,𝑥₂(0)=0,𝑥₃(0)=0, and 𝑥₄(0)=0. The control weight matrix 𝑄 is a diagonal matrix with 𝑄=diag(𝑞₁𝑞₂𝑞₃𝑞₄). Apparently, NLCO problem (34) has five middle decision variables. In fact, it only has three independent variables. Different optimization methods such as the Newton method, quasi-Newton method, Lagrange method, conjugate gradient method, interior point method, linear programming method, particle swarm algorithm, genetic algorithm, evolution algorithm, and other intelligent algorithms in [50–63] have been well established to solve the NLCO problem. Therefore, optimal parameters 𝑚^∗_𝑖 and proposed PID controller are gotten: 𝐺^∗_𝑐(𝑠)=𝑘^∗_𝑝+𝑘^∗_𝑖𝑠+𝑘^∗_𝑑𝑠,𝑘^∗_𝑝=𝑚^∗₄−𝑎₃−𝑏₂𝑘^∗_𝑖𝑏₃,𝑘^∗_𝑖=𝑚^∗₅𝑏₃,𝑘^∗_𝑑=𝑚^∗₁−𝑎₀𝑏₁. (35)

It is noted that (i) the proposed method is also suitable for first-order, second-order, and high-order processes; (ii) the proposed method is also suitable for nonunit feedback system. The non-unit feedback system can be transformed into unit feedback system through transfer function equivalent principle; (iii) the proposed method is also suitable for processes with time delay which can be approximated by transfer function; (iv) the proposed method can be used for different processes.

3.2. Characteristic Analysis of 𝜀-Routh Stability

The 𝜀-Routh stability constraint will imply control system robustness in some degree. The 𝜀-Routh stability constraint could avoid too conservative system performances and may help to improve the robustness of control systems. The closed-loop characteristic polynomial 𝐷(𝑠) is 𝑠^𝑛+1+𝑎₁𝑠^𝑛+⋯+𝑎_{𝑛−𝑤−2}𝑠^𝑤+3+𝑚₁𝑠^𝑤+2+𝑚₂𝑠^𝑤+1+⋯+𝑚_𝑤+2𝑠+𝑚_𝑤+3. The polynomial coefficient 𝑚_𝑘(𝑘=1,2,…,𝑤+3) is a linear mapping function of the PID controller parameters 𝑘_𝑝, 𝑘_𝑖, and 𝑘_𝑑. The parameters 𝑘_𝑝, 𝑘_𝑖, and 𝑘_𝑑 are constrained in intervals (0,𝜏⁰_𝑝], (0,𝜏⁰_𝑖], and (0,𝜏⁰_𝑑] respectively. Therefore, the PID controller parameters can be seen as the uncertain bounded parameters in some degree. In this sense, polynomial 𝐷(𝑠) can be seen as interval polynomial with three uncertain parameters 𝑘_𝑝, 𝑘_𝑖, and 𝑘_𝑑. To analyze the robustness of 𝜀-Routh stability, without losing generality, interval polynomial 𝐹(𝑠) has the general form: 𝐹(𝑠)=𝑠^𝑛+𝑝₁𝑠^𝑛−1+𝑝₂𝑠^𝑛−2+⋯+𝑝_𝑛−1𝑠+𝑝_𝑛, (36) where 𝑝₁,𝑝₂,…,and  𝑝_𝑛 are the uncertain coefficients of the interval polynomial 𝐹(𝑠). The uncertain coefficients 𝑝₁,𝑝₂,…,  and  𝑝_𝑛 are continuous mapping functions of uncertain bounded parameters 𝑞₁,𝑞₂,…,  and  𝑞_𝑚: 𝑝₁=𝑓₁𝑞₁,…,𝑞_𝑚,…,𝑝_𝑛=𝑓_𝑛𝑞₁,…,𝑞_𝑚Ω_𝑞=𝑞₁,…,𝑞_𝑚𝑞⁻₁,𝑞⁺₁×⋯×𝑞⁻_𝑚,𝑞⁺_𝑚𝑞₁∈𝑞⁻₁,𝑞⁺₁,…,𝑞_𝑚∈𝑞⁻_𝑚,𝑞⁺_𝑚, (37) where Ω_𝑞 is the uncertain parameter space, 𝑞⁻_𝑖 and 𝑞⁺_𝑖 are the lower bound and upper bound of uncertain parameter 𝑞_𝑖, and 𝑓_𝑖(⋅) is continuous mapping function of polynomial 𝐹(𝑠) for uncertain parameters. Based on 𝜀-Routh stability definition for interval polynomial 𝐹(𝑠), the following inequalities can be obtained: 𝑅_1,1𝑞₁,…,𝑞_𝑚=𝑓_𝑟,1𝑝₁,…,𝑝_𝑛𝑅>𝜀,…,_𝑛+1,1𝑞₁,…,𝑞_𝑚=𝑓_𝑟,𝑛+1𝑝₁,…,𝑝_𝑛Ω>𝜀_𝑟,𝜀=𝑞₁,…,𝑞_𝑚∣𝑓_𝑟,1𝑝₁,…,𝑝_𝑛𝑓>𝜀,…,_𝑟,𝑛+1𝑝₁,…,𝑝_𝑛>𝜀,(𝜀≥0), (38) where 𝑓_𝑟,𝑖(⋅) is the element of first column of Routh array. According to the interval polynomial stable criterion in [36, 64], if the disjoint set Ω_𝑓 of feasible set Ω_𝑟,𝜀 and uncertain parameters space Ω_𝑞 is not a empty set, then the interval polynomial 𝐹(𝑠) will satisfy the interval polynomial stable criterion: min_{𝑞1,…,𝑞𝑚∈Ω𝑓}𝑓_𝑟,1𝑝₁,…,𝑝_𝑛,…,𝑓_𝑟,𝑛+1𝑝₁,…,𝑝_𝑛Ω>0,_𝑓=Ω_𝑟,𝜀Ω_𝑞≠∅. (39) Therefore, the interval polynomial 𝐹(𝑠) is Routh-Hurwitz stable in the feasible set Ω_𝑓, and the 𝜀-Routh stability could guarantee the robustness of linear systems in a certain degree. For any positive real 𝜀₂>𝜀₁>0, two inequalities for 𝜀-Routh stability are obtained: I∶Ω₁=𝑞₁,…,𝑞_𝑚∣𝑓_𝑟,1𝑝₁,…,𝑝_𝑛>𝜀₁𝑓,…,_𝑟,𝑛+1𝑝₁,…,𝑝_𝑛>𝜀₁,II∶Ω₂=𝑞₁,…,𝑞_𝑚∣𝑓_𝑟,1𝑝₁,…,𝑝_𝑛>𝜀₂𝑓,…,_𝑟,𝑛+1𝑝₁,…,𝑝_𝑛>𝜀₂. (40) The inequality (I) and inequality (II) have the feasible spaces Ω₁ and Ω₂, and the first inequality of inequalities (I) and (II) has the feasible set Ω¹₁=𝑞₁,…,𝑞_𝑚∣𝑓_𝑟,1𝑝₁,…,𝑝_𝑛>𝜀₁,Ω¹₂=𝑞₁,…,𝑞_𝑚∣𝑓_𝑟,1𝑝₁,…,𝑝_𝑛>𝜀₂. (41) Then, any point 𝑀^∗₁ in the feasible set Ω¹₂ will have the relationship ∃𝑀^∗₁=𝑞^∗₁₁,𝑞^∗₂₁,…,𝑞^∗_𝑚,1,𝑓_𝑟,1𝑀^∗₁=𝜀^∗₂>𝜀₂. (42) Likewise, any point 𝑀^∗₂ in the feasible set Ω¹₁ will have the relationship ∃𝑀^∗₂=𝑞^∗₁₂,𝑞^∗₂₂,…,𝑞^∗_𝑚,2,𝑓_𝑟,1𝑀^∗₂=𝜀^∗₁>𝜀₁. (43) For the feasible set Ω¹₁, there at least a point 𝑁^∗₁ which will satisfy the relationshipCase1∶∃𝑁^∗₁=𝑞^∗₁₁,𝑞^∗₂₁,…,𝑞^∗_𝑚,1,𝜀₁<𝑓_𝑟,1𝑁^∗₁<𝜀₂⟺𝑓_𝑟,1𝑁^∗₁∈𝜀₁,𝜀₂⟹𝑁^∗₁∈Ω¹₁,𝑁^∗₁∉Ω¹₂⟹Ω¹₁⊃Ω¹₂.Case2∶𝑓_𝑟,1𝑁^∗₁∉𝜀₁,𝜀₂⟹Ω¹₁=Ω¹₂. (44) Therefore, the feasible sets Ω^𝑘₁ and Ω^𝑘₂ will have the following relationship: Case1∶Ω^𝑘₁⊃Ω^𝑘₂,(∀𝑘=1,2,…,𝑛,𝑛+1)Case2∶Ω^𝑘₁=Ω^𝑘₂,()Ω∀𝑘=1,2,…,𝑛,𝑛+1_𝑘1=𝑞₁,…,𝑞_𝑚∣𝑓_𝑟,𝑘𝑝₁,…,𝑝_𝑛>𝜀₁Ω^𝑘₂=𝑞₁,…,𝑞_𝑚∣𝑓_𝑟,𝑘𝑝₁,…,𝑝_𝑛>𝜀₂. (45) Then, the overall feasible spaces Ω₁ and Ω₂ will have following relationships: Ω₁⊇Ω₂⟺Ω₁⊃Ω₂∪Ω₁=Ω₂Ω₁=^𝑛+1_𝑗=1Ω^𝑗₁,Ω₂=^𝑛+1_𝑗=1Ω^𝑗₂. (46) The overall feasible set will be gradually condensed with the ε increasing. With 𝑎 real ε sequence 𝜀_𝑡+1>𝜀_𝑡>𝜀_𝑡−1>⋯>𝜀₂>𝜀₁, one will get following relationships: Ω_𝑡=𝑞₁,…,𝑞_𝑚∣𝑓_𝑟,1𝑝₁,…,𝑝_𝑛>𝜀_𝑡𝑓,…,_𝑟,𝑛+1𝑝₁,…,𝑝_𝑛>𝜀_𝑡⟺Ω_𝑡=Ω¹_𝑡∩⋯∩Ω_𝑡^𝑛+1=⋂^𝑛+1_𝑗=1Ω^𝑗_𝑡,𝑡=1,2,3,...,𝑡∈ℤ⁺Ω^𝑘_𝑡=𝑞₁,…,𝑞_𝑚∣𝑓_𝑟,𝑘𝑝₁,…,𝑝_𝑛>𝜀_𝑡,(𝑘=1,2,…,𝑛,𝑛+1). (47) Therefore, the feasible set Ω_𝑡 will has the following relationship: Ω_𝑡+1⊆Ω_𝑡⊆⋯⊆Ω₂⊆Ω₁⟺Ω_𝑡+1Δ_𝑡Ω_𝑡⊂⋯Δ₂Ω₂Δ₁Ω₁.Ω_𝑡+1=Ω_𝑡=⋯=Ω₂=Ω₁,∃𝑘,Δ_𝑘=“⊂”,𝑘∈{1,2,…,𝑡}. (48) When the feasible set Ω_𝑡 is strictly condensed with 𝜀_𝑡, then the limit of feasible space will yield the following relationships: Case1∶lim_{𝜀𝑡→+∞}𝑞₁,…,𝑞_𝑚∣𝑓_𝑟,1𝑝₁,…,𝑝_𝑛>𝜀_𝑡𝑓,…,_𝑟,𝑛+1𝑝₁,…,𝑝_𝑛>𝜀_𝑡=lim_{𝜀𝑡→+∞}Ω_𝑡=lim_{𝜀𝑡→+∞}Ω¹_𝑡∩Ω²_𝑡∩⋯∩Ω_𝑡^𝑛+1=lim_{𝜀𝑡→+∞}^𝑛+1_𝑗=1Ω^𝑗_𝑡=Ω^∗=𝑞^∗_1,0𝑞^∗_2,0,…,𝑞^∗_𝑚,0.Case2∶lim_{𝜀𝑡→+∞}𝑞₁,…,𝑞_𝑚∣𝑓_𝑟,1𝑝₁,…,𝑝_𝑛>𝜀_𝑡𝑓,…,_𝑟,𝑛+1𝑝₁,…,𝑝_𝑛>𝜀_𝑡=lim_{𝜀𝑡→+∞}Ω_𝑡=lim_{𝜀𝑡→+∞}Ω¹_𝑡∩Ω²_𝑡∩⋯∩Ω_𝑡^𝑛+1=lim_{𝜀𝑡→+∞}^𝑛+1_𝑗=1Ω^𝑗_𝑡=∅. (49) If the 𝜀_𝑡 tends to infinity, then the limit of the feasible set will be condensed into a point Ω^∗ or an empty set ∅. If the feasible set is condensed into a point Ω^∗, then the interval polynomial degenerates into a usual polynomial. It reflects that the point Ω^∗ is a fixed point in the stable region of linear control systems. If the feasible set tends to empty set ∅, therefore it has a real 𝜀^∗_𝑡  to satisfy the following relationships: 𝜀_𝑡⟶+∞,lim_{𝜀𝑡→+∞}Ω_𝑡=∅⟺∃𝜀^∗_𝑡,𝜀_𝑡>𝜀^∗_𝑡,Ω_𝑡=∅,0<𝜀_𝑡≤𝜀^∗_𝑡,Ω_𝑡𝜀≠∅,_∗𝑡>𝜀_𝑡>𝜀_𝑡−1>⋯>𝜀₂>𝜀₁Ω≥0,_∗𝑡⊂Ω_𝑡⊂Ω_𝑡−1⊂⋯⊂Ω₂⊂Ω₁ (50)when the feasible set has the relationship Ω_𝑡+1=Ω_𝑡=Ω_𝑡−1=⋯=Ω₂=Ω₁. However, this is a special event which can be depicted by the probability event 𝑃Ω_𝑡+1=Ω_𝑡=Ω_𝑡−1=⋯=Ω₂=Ω₁⟺𝑃Ω_𝑡+1=Ω_𝑡,Ω_𝑡=Ω_𝑡−1,…,Ω₂=Ω₁. (51) It assumes that all the probability events are independent and identically distributed. Hence, the probability can be written as 𝑃Ω_𝑡+1=Ω_𝑡,Ω_𝑡=Ω_𝑡−1,…,Ω₂=Ω₁=𝑃Ω_𝑡+1=Ω_𝑡×𝑃Ω_𝑡=Ω_𝑡−1×⋯×𝑃Ω₂=Ω₁=^𝑡+1_𝑗=2𝑃Ω_𝑗=Ω_𝑗−1=𝑃Ω₂=Ω₁^𝑡,𝑡=1,2,…,𝑡∈ℤ⁺. (52) Considering probability event {Ω₂=Ω₁} first, the probability will have 𝑃Ω₁=Ω₂=𝑃Ω¹₁∩⋯∩Ω₁^𝑛+1=Ω¹₂∩⋯∩Ω₂^𝑛+1=𝑃^𝑛+1_𝑗=1Ω^𝑗₁=^𝑛+1_𝑗=1Ω^𝑗₂⟺𝑃Ω¹₂∩⋯∩Ω₂^𝑛+1⊂Ω^𝑘₁,Ω¹₁∩⋯∩Ω₁^𝑛+1⊂Ω^𝑘₂,(∀𝑘=1,2,…,𝑛,𝑛+1). (53) The probability 𝑃{Ω₂=Ω₁} for the probability event {Ω₂=Ω₁} will yield 𝑃Ω¹₂∩Ω²₂∩⋯∩Ω^𝑛₂⊂Ω¹₁=𝜆₁𝑃Ω,…,₁₂∩Ω²₂∩⋯∩Ω^𝑛₂⊂Ω₁^𝑛+1=𝜆₁,𝑃Ω¹₁∩Ω²₁∩⋯∩Ω^𝑛₁⊂Ω¹₂=𝜆₂𝑃Ω,…,₁₁∩Ω²₁∩⋯∩Ω^𝑛₁⊂Ω₂^𝑛+1=𝜆₂,𝑃Ω₁=Ω₂=𝑃^𝑛+1_𝑗=1Ω^𝑗₂⊂Ω^𝑘₁,(∀𝑘=1,2,…,𝑛,𝑛+1)×𝑃^𝑛+1_𝑗=1Ω^𝑗₁⊂Ω^𝑘₂,()=∀𝑘=1,2,…,𝑛,𝑛+1_𝑛+1𝑘=1𝑃^𝑛+1_𝑗=1Ω^𝑗₂⊂Ω^𝑘₁×^𝑛+1_𝑘=1𝑃^𝑛+1_𝑗=1Ω^𝑗₁⊂Ω^𝑘₂=𝜆₁^𝑛+1𝜆₂^𝑛+1,0<𝜆₁<1,0<𝜆₂<1. (54) Therefore, the probability of the probability event yields 𝑃Ω_𝑡+1=Ω_𝑡,Ω_𝑡=Ω_𝑡−1,…,Ω₂=Ω₁=𝜆₁^(𝑛+1)𝑡𝜆₂^(𝑛+1)𝑡𝑃Ω_𝑡+1Δ_𝑡Ω_𝑡Δ_𝑡−1Ω_𝑡−1⋯Δ₂Ω₂Δ₁Ω₁=1−𝜆₁^(𝑛+1)𝑡𝜆₂^(𝑛+1)𝑡. (55)

3.3. Interior Method for NLCO Problem

Without losing generality, NLCO problem (34) of the optimal PID controller can be written in the following general forms: min𝑓(𝑥)s.t∶ℎ(𝑥)=0,𝑔(𝑥)≤0, (56) where 𝑓(𝑥)∶𝑅^𝑛→𝑅,ℎ(𝑥)∶𝑅^𝑛→𝑅^𝑚, and 𝑔(𝑥)∶𝑅^𝑛→𝑅^𝑞 are the smooth and differentiable functions, 𝑥 is the decision variable, and 𝑛, 𝑚, and 𝑞 denote the number of the decision variable, equality constraint, and inequality constraint, respectively. The various optimization methods in [50–63] are developed to solve optimization problems. The optimization methods such as the Newton methods, conjugate gradient methods, steepest descent methods, interior point methods, genetic algorithms, and particle swarm algorithms in [50–63], are well established to solve constraint optimization problems. Based on interior point methods in [56–61], the interior method is recommended to solve the NLCO problem. Then, the NLCO problem (56) could be transformed into the following form: min𝑓(𝑥)−𝜐^𝑞_𝑖=1ln𝛿_𝑖s.t∶ℎ(𝑥)=0,𝑔(𝑥)+𝛿=0, (57) where 𝜐>0 is the barrier parameter, the slack vector 𝛿=(𝛿₁,𝛿₂,…,𝛿_𝑞)^𝑇>0 is set to be positive, and 𝑔(𝑥) is an expanded constraint inequality. It introduces the Lagrange multipliers 𝑦 and 𝑧 for the barrier problem (57): 𝐿𝑥,𝑦,𝑧,𝛿=𝑓(𝑥)−𝜐^𝑞_𝑖=1ln𝛿_𝑖+𝑦^𝑇𝑔(𝑥)+𝛿+𝑧^𝑇ℎ(𝑥), (58) where 𝐿(𝑥,𝑦,𝑧,𝛿) is Lagrange function, 𝐲=(𝑦₁,𝑦₂,…,𝑦_𝑞)^𝑇 and 𝐳=(𝑧₁,𝑧₂,…,𝑧_𝑚)^𝑇 are Lagrange multipliers for constraints 𝑔(𝑥)+𝛿 and ℎ(𝑥), respectively. Based on the Karush-Kuhn-Tucker (KKT) optimality conditions [57–59], the optimality conditions for NLCO problem (56) can be expressed as ∇_𝑥𝐿𝑥,𝑦,𝑧,𝛿=∇𝑓(𝑥)+∇𝑔(𝑥)^𝑇𝑦+(∇ℎ(𝑥))^𝑇∇𝑧=0_𝛿𝐿𝑥,𝑦,𝑧,𝛿=−𝜐𝑆_𝛿⁻¹𝑒+𝑦=0⟺−𝜐𝑒+𝑆_𝛿∇𝑌𝑒=0_𝑦𝐿(𝑥)∇𝑥,𝑦,𝑧,𝛿=𝑔+𝛿=0_𝑧𝐿𝑥,𝑦,𝑧,𝛿=ℎ(𝑥)=0, (59) where 𝑆_𝛿=diag(𝛿₁,𝛿₂,…,𝛿_𝑞) is the diagonal matrix and its elements are the components of the vector 𝜹, and 𝐞 is a vector of all ones, ∇ℎ(𝑥) and ∇𝑔(𝑥) are the Jacobian matrices of the vectors ℎ(𝑥) and 𝑔(𝑥), respectively, ∇𝑓(𝑥) is the grand of function 𝑓(𝑥), and 𝑌=diag(𝑦₁,𝑦₂,…,𝑦_𝑞) is a diagonal matrix and its elements are the components of vector 𝐲. The system (59) is the KKT condition of the NLCO problem (56). When in the search approach, it should remain 𝛿, 𝑦>0. To obtain the iteration direction, it can make the point (𝑥+Δ_𝑥,𝛿+Δ𝛿,𝑧+Δ_𝑧,𝑦+Δ𝑦) satisfying the KKT conditions(59); then the following system will be obtained: ∇²𝑓(𝑥)+∇²𝑔(𝑥)^𝑇𝑦+∇²ℎ(𝑥)^𝑇𝑧Δ𝑥+∇𝑔(𝑥)^𝑇(𝑥)Δ𝑦+∇ℎ^𝑇Δ𝑧+∇𝑓(𝑥)+∇𝑔(𝑥)^𝑇𝑦+∇ℎ(𝑥)^𝑇𝑧=0−𝜐𝑒+𝑆_𝛿𝑌𝑒+𝑆_𝛿𝑔(𝑥)(𝑥)(𝑥)(𝑥)Δ𝑦+𝑌Δ𝛿=0+∇𝑔Δ𝑥+𝛿+Δ𝛿=0,∇ℎΔ𝑥+ℎ=0. (60) System (60) is obtained by ignoring higher-order incremental and replaces the nonlinear terms with linear approximation in system (59). The system (60) is written in the matrix form: 𝐻𝑥,𝑦,𝑧0∇ℎ(𝑥)^𝑇∇𝑔(𝑥)^𝑇0𝑆_𝛿⁻¹(𝑥)(𝑥)=(𝑥)(𝑥)𝑌0𝐼∇ℎ000∇𝑔𝐼00Δ𝑥Δ𝛿Δ𝑦Δ𝑧−∇𝑓−∇𝑔^𝑇𝑦−∇ℎ(𝑥)^𝑇𝑧𝜐𝑆_𝛿⁻¹(𝑥)(𝑥)𝑒−𝑦−ℎ−𝑔−𝛿𝐻𝑥,𝑦,𝑧=∇²𝑓(𝑥)+∇²𝑔(𝑥)^𝑇𝑦+∇²ℎ(𝑥)^𝑇𝑧, (61)(62) where 𝐻(𝑥,𝑦,𝑧) is the Hessian matrix in system. Finally, the new iterate direction is obtained via solving the system (61), which is the essential process in the interior point method. Thus, the new iteration point can be gotten in the next iteration: 𝑥,𝛿,𝑧,𝑦⟵𝑥,𝛿,𝑧,𝑦+𝜁₁Δ𝑥,Δ𝛿,Δ𝑧,Δ𝑦, (63) where 𝜁₁ is the step size. Choosing the step size 𝜁₁ holds the 𝛿, 𝑦>0 in search process. In this paper, the interior point method is recommended to solve the NLCO problem. The interior algorithm framework is depicted as follows.

Interior Algorithm Framework for NLCO Problem. One has the following.

Step 1. Choose an initial iteration point (𝑥^(0,),𝛿⁽⁰⁾,𝑧⁽⁰⁾,𝑦⁽⁰⁾) in the feasible region set ={𝑥∣ℎ(𝑥)=0,𝑔(𝑥)≤0}, and the 𝛿⁽⁰⁾>0,𝑦⁽⁰⁾>0,𝑘=0.

Step 2. Construct current iterate; we have the current iterate values 𝑥^(𝑘), 𝛿^(𝑘), 𝑧^(𝑘), and 𝑦^(𝑘) of the primal variable 𝑥, the of the slack variable 𝛿, the multipliers 𝑦 and 𝑧, respectively.

Step 3. Calculate the Hessian matrix 𝐻(𝑥,𝑦,𝑧) of the Lagrange system 𝐿(𝑥,𝑦,𝑧,𝛿), and the Jacobian matrices ∇ℎ(𝑥) and ∇𝑔(𝑥) are of the vectors ℎ(𝑥) and 𝑔(𝑥) in the current iterate (𝑥^(𝑘),𝛿^(𝑘),𝑧^(𝑘),𝑦^(𝑘)).

Step 4. Solve the linear system (61) and construct the iterate direction (Δ𝑥,Δ𝛿,Δ𝑧,Δ𝑦). Solving the linear matrix equation (61), we obtain the primal solution Δ𝑥, multipliers solution Δ𝑧, Δ𝑦, and also the slack variable solution Δ𝛿.

Step 5. Choosing the step size 𝜁₁ holds the 𝛿, 𝑦>0 in search process, 𝜁₁∈[0,1]. Update the iterate values: (𝑥^(𝑘+1),𝛿^(𝑘+1),𝑧^(𝑘+1),𝑦^(𝑘+1))←(𝑥^(𝑘),𝛿^(𝑘),𝑧^(𝑘),𝑦^(𝑘))+𝜁₁(Δ𝑥,Δ𝛿,Δ𝑧,Δ𝑦),𝑘←𝑘+1.

Step 6. Check the ending conditions in region 𝛀. If it does not satisfied, go to Step 3, else minimum 𝑓_min of NLCO problem is obtained in feasible region Ω.

Step 7. End.

4. Simulation and Experiment Results

The proposed method is used to design the optimal PID controller with 𝜀-Routh stability for processes 𝐺_𝑝1(𝑠) [7], 𝐺_𝑝2(𝑠) [15], and 𝐺_𝑝3(𝑠) [6], respectively. 𝐺_𝑝1(𝑠) is a third-order all-pole process, and 𝐺_𝑝2(𝑠) and 𝐺_𝑝3(𝑠) are the SOPTD and FOPTD processes, respectively, 𝐺_𝑝1(𝑠)=15𝑠²(),𝐺+0.9𝑠+5𝑠+3_𝑝2(𝑠)=𝑒^−0.5𝑠(𝑠+1)²,𝐺_𝑝3(𝑠)=4𝑒^−10𝑠.10𝑠+1 (64)

In order to use the proposed method, time delay of plant.2 and plant.3 is approximated as 𝑒^−0.5𝑠≈(1−0.25𝑠)/(1+0.25𝑠) and 𝑒^−10𝑠≈(𝑠²−0.6𝑠+0.12)/(𝑠²+0.6𝑠+0.12), respectively. Thus, transfer functions of 𝐺_𝑝2(𝑠) and 𝐺_𝑝3(𝑠) are the (−𝑠+4)/(𝑠³+6𝑠²+9𝑠+4) and (0.4𝑠²−0.24𝑠+0.048)/(𝑠³+0.7𝑠²+0.18𝑠+0.012), respectively. The proposed method is used to design the optimal PID controller with 𝜀-Routh stability for different processes, and optimal PID controllers under different control weight matrices and 𝜀-Routh stability are shown in Table 1.

Table 1 

Proposed optimal PID controllers of plants.1–3 under various control weight matrices.

For the proposed method, the plant.1’s step response and system performances under various control weight matrices are shown in Figure 4 and Table 2, respectively. The control weight matrix has effects on system performances. Figure 4 and Table 2 obviously show that system performances are affected by control weight factors 𝑞₁, 𝑞₂, 𝑞₃, and 𝑞₄. The results show that the response velocity will be heightened by increasing weight factor 𝑞₁ or reducing weight factors 𝑞₂, 𝑞₃, and 𝑞₄. When control weight factor 𝑞₁ (with 𝑞₂, 𝑞₃, and 𝑞₄ fixed) is increased, the overshoot is increased, and the rise time, delay time, peak time, and setting time are reduced. Likewise, the overshoot will be increased, and the rise time, delay time, and peak time will be reduced by reducing weight factors 𝑞₂, 𝑞₃, and 𝑞₄ (with 𝑞₁ fixed). Response results show that system performances are affected by control weight matrix.

Table 2 

Plant.1’s system performances of proposed method under various control weight matrices.

(a)

(b)

Figure 4 

Step response of the proposed method under various control weight matrices (for plant.1).

For 𝑛th-order processes, it can be inferred that overshoot will be heightened, and the rise time, delay time, and peak time will be reduced by increasing control weight factor 𝑞₁ (with other factors fixed) or reducing other weight factors (with 𝑞₁ fixed). On the contrary, the overshoot will be reduced, and the rise time, delay time, and peak time will be increased by reducing weight factor 𝑞₁ (with other factors fixed) or increasing other factors (with 𝑞₁ fixed). Therefore, system performances are affected by control weight factors 𝑞₁,𝑞₂,…,𝑞_𝑛+1. In all, control weight matrix has effects on system performances.

For plant.2, step response results and system performances under different 𝜀-Routh stability are shown in Figure 5 and Table 3, respectively. The response results show that system performances are affected by 𝜀-Routh stability (with control weight matrix fixed). The results show that response velocity will be heightened by increasing 𝜀. When the 𝜀 (with control weight matrix fixed) is increased, the overshoot is increased, and the rise time, delay time, and peak time are reduced. However, setting time is gradually increased and then is decreased by increasing 𝜀. Therefore, the rise time, delay time, and peak time will be increased, and the overshoots will be reduced by reducing the 𝜀. In all, system performances are affected by the 𝜀-Routh stability too.

Table 3 

System performances of proposed method under different ε values (for plant.2).

(a)  

(b)

Figure 5 

Step response of the proposed method under various 𝜀-stability (for plant.2).

4.1. System Performances of Different Tuning Methods

The step response and system performances of different tuning methods for plant.1, plant.2, and plant.3 are shown in Figures 6, 7, and 8 and Table 4, respectively. For plant.1, the results obviously show that the Z-N method, IA-ISE method [7], and IA-ITSE method [7] show oscillatory behaviors, but the proposed method shows smoothly response without any oscillation. The Z-N method has small rise time, delay time, and peak time compared with the proposed method, IA-ISE method, and IA-ITSE method. However, the Z-N method has the largest overshoot and setting time among these tuning methods. The overshoot and setting time of the Z-N method is about 11.08% and 6.41 sec, and the overshoot and setting time of the proposed method, IA-ISE method, and IA-ITSE method are about 4.08% and 3.54 sec, 4.24% and 3.13 sec, and 5.27% and 4.88 sec, respectively. Thus, these tuning methods own advantages for plant.1.

Table 4 

System performances of different tuning methods for plants.1–3.

(a)  

(b)

Figure 6 

Step response of plant.1 under different tuning methods.

(a)  

(b)

Figure 7 

Step response of plant.2 under different tuning methods.

(a)  

(b)

Figure 8 

Step response of plant.3 under different tuning methods.

The plant.2 is a SOPTD plant with dead time ratio 𝜏/𝑇=0.5. The results in Figure 7 show that dynamic performances of the Z-N method, A-H method [19], gain and phase method [15], and proposed method have different features. The Z-N method, A-H method [19], gain and phase method [15], and the proposed method have small rise time, delay time and peak time, setting time, and overshoot, respectively. The proposed method has the largest rise time, delay time, and peak time, and the A-H method and Z-N method have the largest setting time and overshoot, respectively. The longest rise time, delay time, peak time, and setting time are about 4.778 sec, 1.841 sec, 6.622 sec, and 4.43 sec, respectively, and the largest overshoot is about 43.09%. The gap between the largest and the smallest items of the rise time, delay time, peak time, and setting time is small, but only the gap between the largest and the smallest overshoot is about 13 times of the smallest ones. Thus, the comprehensive performances of Z-N method are poorer than rest methods.

The plant.3 is a FOPTD plant with large dead time ratio 𝜏/𝑇=1.0. The results in Figure 8 show that all the tuning methods show oscillation response. The overshoot of the Z-N method (C-C) and Z-N method (PRC) are about 18.18% and 21.17%, respectively. The results demonstrate that the proposed method has the smallest peak time 20.87 sec and the smallest overshoot 6.67%, and the Cohen-Coon method has the smallest rise time 17.252 sec and the smallest delay time 14.306 sec, and the Z-N method (C-C) has about the smallest setting time 51.60 sec. However, the Z-N method (PRC) has the largest rise time 20.286 sec, delay time 14.635 sec, and peak time 26.43 sec, and the Cohen-Coon method has the largest setting time 67.05 sec and largest overshoot 38.06%. For plant.3, the Cohen-Coon method shows poorer system performances than other tuning methods. By studying different tuning methods’ dynamic performances, it can be known that the proposed method could provide good system performances.

4.2. Disturbance Rejection Ability of Different Methods

Disturbance response results for plant.1, plant.2, and plant.3 are shown in Figures 9, 10, and 11, respectively. For plant.1, the results show that the proposed method, Z-N method, IA-ISE and IA-ITSE methods have good ability to reject the disturbances, and they show roughly the same disturbance response amplitude (DRA) and disturbance attenuation time (DAT). The proposed method has the best disturbance rejection ability among the methods. The proposed method does not show oscillatory behaviors, but Z-N method, IA-ISE method, and IA-ITSE method show oscillatory behaviors. For plant.2, the disturbance response results of the proposed method, Z-N method, A-H method, and gain and phase method are shown in Figure 10. The results show that these methods have good ability to reject disturbances and they have roughly the same DAT. The DRAs of the proposed method, Z-N method, A-H method, and gain and phase method are about 0.13, 0.10, 0.10, and 0.11, respectively. The proposed method has larger DRA than the rest methods, but the proposed method shows smooth disturbance response. For plant.3, disturbance response results in Figure 11 show that these tuning methods have roughly the same ability to reject disturbances. The proposed method, Z-N method (C-C), Z-N method (PRC), and Cohen-Coon method have roughly the same DRA and DAT. Among these tuning methods, the DRA of the Z-N method (PRC) is about 0.275, the proposed method, Z-N method (C-C), and Cohen-Coon method have the same DRA which is about 0.25. Overall, as known from disturbance response results, it could be that the proposed method, A-H method, and gain and phase method roughly have roughly the same ability to reject disturbances.

Figure 9 

Disturbance rejection ability of different tuning methods for plant.1.

Figure 10 

Disturbance rejection ability of different tuning methods for plant.2.

Figure 11 

Disturbance rejection ability of different tuning methods for plant.3.

To study 𝜀-Routh stability’s effects on disturbance rejection ability, the disturbance response results under various 𝜀 are shown in Figure 12. The disturbance response results in Figure 12 clearly show that the ε has effects on the DRA and DAT. The DAT will be continuously reduced when the 𝜀 increases from 0.5 to 5.0. However, the DAR will be gradually reduced and then be increased when the 𝜀 increases from 0.5 to 5.0. The DAT reflects how fast the disturbances will be eliminated, and DRA reveals how much disturbance amplitude of the control system. Thus, 𝜀-Routh stability has effects on the disturbance rejection ability.

Figure 12 

Disturbance rejection ability of the proposed method under various ε-stability (for plant.2).

4.3. Experiment Results for the CWTS

The proposed PID controller is utilized to control the liquid level of a CWTS with interactive behavior. The proposed PID controllers’ experimental results of the CWTS under various control weight matrix and different set points are shown in Figures 14, 15, 16, 17, 18, 19, 20, and 21. The theoretical model of the CWTS is already obtained in previous section; therefore the identification model can be acquired via open-loop experiment. The liquid level’s real-time data is stored in the disk of monitoring computer, and the date can be accessed on-line for process control (OPC) database. The identification model of CWTS is identified as ℎ₁(𝑠)/𝑄_𝑖𝑐(𝑠)=(0.4465𝑠+0.3064)/(𝑠²+1.340𝑠+0.01264) via least square method (sample time is 1 second) when open ratio of control valve, output valve-1, output valve-2, and connected valve-3 is set at 75%, 100%, 30%, and 100%, respectively. The steady-state liquid level is about 294.00 mm, and it will take about 12 minutes to reach this liquid level. The actual liquid level and the identification model’s liquid level are shown in Figure 13. The identification model is consistent with the actual model of the CWTS; thus, the identification model describes the actual model accurately.

Figure 13 

Actual liquid level and identified model’s liquid level (open loop).

Figure 14 

Experimental response of the CWTS with the PID controller 𝑘∗𝑝=10.000, 𝑘∗𝑖=0.32549, 𝑘∗𝑑=0.001, and 𝜀=0 in set point 300 mm.

Figure 15 

Experimental response of the CWTS with the PID controllers 𝑘∗𝑝=10.000, 𝑘∗𝑖=0.10118, 𝑘∗𝑑=0.001, and 𝜀=0 in set point 300 mm.

Figure 16 

Experimental response of the CWTS with the PID controllers 𝑘∗𝑝=10.000, 𝑘∗𝑖=0.22986, 𝑘∗𝑑=0.001, and 𝜀=0 in set point 300 mm.

Figure 17 

Experimental response of the CWTS with the PID controllers 𝑘∗𝑝=6.000, 𝑘∗𝑖=0.13817, 𝑘∗𝑑=0.001, and 𝜀=0 in set point 300 mm.

Figure 18 

Experimental response of the CWTS with the PID controllers 𝑘∗𝑝=10.000, 𝑘∗𝑖=0.32549, 𝑘∗𝑑=0.001, and 𝜀=0 in set point 160 mm.

Figure 19 

Experimental response of the CWTS with the PID controllers 𝑘∗𝑝=10.000, 𝑘∗𝑖=0.10118, 𝑘∗𝑑=0.001, and 𝜀=0 in set point 160 mm.

Figure 20 

Experimental response of the CWTS with the PID controllers 𝑘∗𝑝=10.000, 𝑘∗𝑖=0.22986, 𝑘∗𝑑=0.001, and 𝜀=0 in set point 160 mm.

Figure 21 

Experimental response of the CWTS with the PID controllers 𝑘∗𝑝=6.000, 𝑘∗𝑖=0.13817, 𝑘∗𝑑=0.001, and 𝜀=0 in set point 160 mm.

In the closed-loop experiment, the proposed method is used to control the liquid level of the CWTS under different set points. Proposed PID controllers under different control weight matrices and 𝜀-Routh stability are used in the experiments. The proposed PID controllers (𝑘_𝑝=10.000,𝑘_𝑖=0.32549,𝑘_𝑑=0.001,𝜀=0;𝑘_𝑝=10.000,𝑘_𝑖=0.10118,𝑘_𝑑=0.001,𝜀=0;𝑘_𝑝=10.000,𝑘_𝑖=0.22986,𝑘_𝑑=0.001,𝜀=0;𝑘_𝑝=6.000,𝑘_𝑖=0.13817,𝑘_𝑑=0.001,and𝜀=0) are used to control liquid level of the CWTS with interactive behaviors. The first group experiments’ setting point of the liquid level ℎ₁ is set at 300.00 mm and second group experiments’ setting point ℎ₁ is set at 160.00 mm. The corresponding experiment results of the first group and second group experiments are shown in Figures 14–21, respectively. In Figures 14–21, the solid red line is the dynamic response of the liquid level ℎ₁, the solid magenta line is the dynamic response of the liquid level ℎ₂, and the solid cyan line is the control signal. In the first group experiments, the dynamic response results of the proposed method shown in Figures 14–17 obviously reveal that the proposed method provides good system performances for the CWTS with interactive behaviors. In first group experiments, the setting time for the CWTS is less than 7 minutes, and the overshoot of the proposed method is less than 4% under various control weight matrices of the proposed method.

In the second group experiments, the dynamic response results of the proposed method shown in Figures 18–21 clearly reveal that the proposed method provides good system performances and robustness. The proposed PID controller is designed for the liquid level of set point at 294 mm. The actual model of the CWTS has been changed in the second group experiments because the liquid level set point has made a huge change. In the second group experiments, the setting time of the proposed method for the CWTS is less than 3.5 minutes, and the overshoot of the proposed method is less than 12%. Figures 18 and 20 show that the liquid level has a small fluctuation; however, overall system performances are still accepted. The system parameters of liquid model will show perturbations when the set point of liquid level makes a huge change. The proposed method can still provide good system performances when the set point of the liquid level changes from 300 mm to 160 mm. Furthermore, the dynamic response results in Figures 18–21 also reflect that the proposed method shows good robustness to rejection model perturbation, system parameters perturbation, and set point variation.

Simulation results and experiment results reveal that the proposed method can provide good system performances for different processes. The proposed method also demonstrates good disturbance rejection ability and robustness for different processes. Therefore, the proposed method is useful to design the optimal PID controller. In short, the proposed method is suitable for different processes and can provide good system performances. The simulation results and experiment results show the effectiveness and usefulness of the proposed method.

5. Conclusions

In this paper, a systematic method is presented to design the optimal PID controller with 𝜀-Routh stability for different processes via Lyapunov approach. The optimal PID controller is obtained by minimizing AISE performance index which contains the control error and at least first-order error derivative, or even may contain the 𝑛th-order error derivative. The proposed optimal control problem could be equivalently transformed into a NLCO problem via Lyapunov theorems. Therefore, the optimal PID controller can be obtained by solving the NLCO problem via the interior method or other optimization methods. Optimal PID controllers under different control weight matrices and 𝜀-Routh stability are presented for different processes. For the proposed method, control weight matrix and 𝜀-Routh stability’s effects on system performances are analyzed. Different tuning methods’ system performances are discussed.

The control weight matrix and 𝜀-Routh stability have large effects on system performances. For third-order processes, system performances and response velocity are affected by the control weight factors 𝑞₁, 𝑞₂, 𝑞₃, and 𝑞₄. The peak overshoot and response velocity will be heightened by increasing weight factor 𝑞₁ (with other factors fixed). Likewise, the overshoot and response velocity will be heightened by reducing weight factors 𝑞₂, 𝑞₃, and 𝑞₄ (with 𝑞₁ fixed). On the contrary, the overshoot and response velocity will be reduced by reducing weight factor 𝑞₁ (with other factors fixed) or increasing weight factors 𝑞₂, 𝑞₃, and 𝑞₄ (with 𝑞₁ fixed). Through the study, for 𝑛th-order processes, the above conclusions still hold; it can be inferred that system performances are affected by the control weight factors 𝑞₁,𝑞₂,…,𝑞_𝑛+1. The 𝜀-Routh stability (with control weight matrix fixed) has effects on system performances. The results show that the overshoots and response velocity will be heightened by increasing the 𝜀. Overshoots will be increased, and the rise time, delay time, and peak time will be reduced by increasing the 𝜀. On the contrary, overshoots will be reduced, and the rise time, delay time, and peak time will be increased by reducing the 𝜀.

Disturbance rejection ability of different tuning methods is also investigated, and simulation results reveal that the proposed method shows good disturbance rejection ability. The 𝜀-Routh stability’s effects on disturbance rejection ability are investigated. The disturbance response results clearly show that the 𝜀 has effects on the DRA and DAT. The DAT will be continuously reduced by increasing the 𝜀. However, the DAR will be gradually reduced and then be increased by increasing the 𝜀.

To further validate the proposed method, the proposed method is utilized to control the liquid level of CWTS with interactive behaviors. The experimental results under different set points and control weight matrices are presented in detail. Both simulation results and experimental results show the effectiveness and usefulness of the proposed method.

Acknowledgments

A partial content of this paper, as the previous work and earlier version paper, was presented at the 2009 International Conference on Information Engineering and Computer Science (ICIECS-2009), Wuhan, China December 19-20, 2009. This work is supported by the Research Projects no. 2007BAF09B01 and 2009A A04Z155. The authors would like to thank the reviewers for suggestions.