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

Authors: XianHong Li; HaiBin Yu; MingZhe Yuan
Journal: Mathematical Problems in Engineering
 (2013)

---

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


---

## Body

## 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 H2/H∞ 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 Section2, 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 systemx˙=Ax asymptotically stable in a large scope are that: for any given matrix Q=QT>0, there is a matrix P=PT>0 that satisfies the Lyapunov algebraic equation (LAE): ATP+PA=-Q.Theorem 2 (see [1, 2]).
The linear time-invariant (LTI) systemx˙=Ax, and x(0)=x0 is asymptotically stable, for any given matrix Q=QT>0, then the performance index J will have equivalent relationships:
(1)J=∫0∞xTQxdt=xT(0)Px(0),ATP+PA=-Q.Definition of
ε
-Routh stability. D(s) 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 D(s)=a0sn+a1sn-1+⋯+an-1s+an can satisfy the following inequality relationships, then the linear time-invariant control system is considered to possess the ε-Routh stability:
(2)R1,1(a0,a1,…,an)>ε,R2,1(a0,a1,…,an)>ε,…,Rn,1(a0,a1,…,an)>ε,Rn+1,1(a0,a1,…,an)>ε(ε≥0),whereR1,1,R2,1,R3,1,…Rn,1,  and  Rn+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 D(s) 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 D(s) 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 Gp(s):
(3)Gp(s)=b1sw+b2sw-1+⋯+bws+bw+1sn+a1sn-1+⋯+an-1s+an=∑j=1w+1bjsw+1-jsn+∑i=1naisn-i=N(s)M(s),
where b1,…,bw+1,a1,…,  and  an are system parameters, and M(s) and N(s) are the denominator and numerator, respectively. The denominator and numerator’s orders are the deg{M(s)}=n and deg{N(s)}=w, respectively. For most practical processes, the number of poles n is greater than the number of zeros w(n>w), and n and w are the positive integers. It could be known that for the all-pole plants b1,b2,…,bw=0, additional for the type I plants an=0, and for the type II plants an-1=an=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)u(s)=Gc(s)e(s),Gc(s)=(kps+ki+kds2)s,r(s)-y(s)≜e(s),y(s)=Gp(s)u(s),
where e(t) is control error, r(t) and y(t) are reference command and system output, respectively, Gc(s) is PID controller, and kp, ki, and kd are controller parameters. Equation (4) can be equally expressed as
(5)r(s)=y(s)+e(s),y(s)=N(s)u(s)M(s),u(s)=(kps+ki+kds2)e(s)s.
Therefore, error transfer function can be obtained as
(6)e(s)r(s)=1(1+Gp(s)Gc(s))⟺D(s)e(s)=sM(s)r(s),D(s)=sM(s)+N(s)(kps+ki+kds2),deg{D(s)}=n+1,
where D(s) is the control system’s characteristic polynomial. The error differential equation can be obtained by inversing Laplace transform of (6):
(7)(δM(δ)+N(δ)(kpδ+ki+kdδ2))e(t)=δM(δ)r(t),
where δ is the differential operator. When reference command r(t) is constant signal or zero signal (r(t)=0), or even the piecewise constant signal, the error differential equation will yield
(8){∑j=1w+1bjδw+1-j(δn+1+∑i=1naiδn+1-i)∑j=1w+1bjδw+1-j+(kpδ+ki+kdδ2)∑j=1w+1bjδw+1-j}e(t)=0.
Considering the polynomial differential operator M(δ) and N(δ), error differential equation (8) could be expressed as
(9){δn+1+∑i=1naiδn+1-i+∑j=1w+1bjkpδw+2-j+∑j=1w+1bjkiδw+1-j+∑j=1w+1bjkdδw+3-j}e(t)=0.
Therefore, (9) can be equally expressed as
(10)e(t)(n+1)+a1e(t)(n)+a2e(t)(n-1)+⋯+an-w-2e(t)(w+3)+(an-w-1+b1kd)e(t)(w+2)+(an-w+b1kp+b2kd)e(t)(w+1)+∑j=1w-1(an-w+j+bj+1kp+bjki+bj+2kd)e(t)(w+1-j)+(an+bw+1kp+bwki)e˙(t)+bw+1kie(t)=0.
For the symbol simplicity, the error differential equation (10) is equally expressed as
(11)e(t)(n+1)+a1e(t)(n)+a2e(t)(n-1)+⋯+an-w-2e(t)(w+3)+m1e(t)(w+2)+m2e(t)(w+1)+∑j=1w-1mj+2e(t)(w+1-j)+mw+2e˙(t)+mw+3e(t)=0,
where mi(i=1,2,…,w+3) is the decision parameter. Parameter mi(i=1,2,…,w+3) is obtained as
(12)m1=(an-w-1+kdb1),m2=(an-w+b1kp+b2kd)mj+2=(an-w+j+bj+1kp+bjki+bj+2kd)(j=1,2,…,w-1)mw+2=(an+bw+1kp+bwki),mw+3=bw+1ki.
Thus, the corresponding state-space model of the differential equation (11) can be obtained:
(13)x˙=Ax,
where x is a state vector, A is a state matrix, C1×n is a parameter vector, On×1 is a zero vector, and In×n is an unit matrix. They are obtained as
(14)A=(On×1In×n-mw+3C1×n),x=(e(t)e˙(t)e¨(t)⋯e(t)(n))T,C1×n=(-mw+2⋯-m1-an-w-2⋯-a2-a1).
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:
(15)J*(kp*,ki*,kd*)=minGc*(s)J(m1,…,mw+3)=min∫0∞xT(t)Qx(t)dt=min{∑i=0n∫0∞qi+1,i+1(e(t)(i))2dt}subjectto(s.t):(1)x˙=Ax;(2)Re(eig(A))<0,Re{s∣det(sI-A)=0}<0(3)ε≥0,R1,1(m1,…,mw+3)>ε,R2,1(m1,…,mw+3)>ε,…,Rn+1,1(m1,…,mw+3)>ε,Rn+2,1(m1,…,mw+3)>ε;(4)m4=f1(m1,m2,m3),…,mw+3=fw(m1,m2,m3);(5)0<kp≤τp0,0<ki≤τi0,0<kd≤τd0,
where Q is a control weight matrix, and at least a positive semidefinite real symmetric matrix, and f1,…,fw is the linear function of the controller parameters kp, ki, and kd, and kp*, ki*, and kd* 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 mi(i=1,2,…,w+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 τp0, τi0, and τd0 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 Figure2. 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 Figure3 contains the interacting behavior because the liquid level h2 of tank-2 depends on the liquid level h1 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:
(16)Qic-Q1c-Q12=C1dh1dt,
where Qic is the flow of control valve, Q1c is the flow of valve-1, Q12 is the flow of valve-3, C1 is the hydraulic capacity of tank-1, and h1 is the liquid level of tank-1. Likewise, the liquid level equation of tank-2 can be obtained:
(17)Q12-Q2c=C2dh2dt,
where Q2c is the flow of valve-2, C2 is the hydraulic capacity of tank-2, and h2 is the liquid level of tank-2. The water flows equation of valve-1, valve-2, and valve-3 can be obtained:
(18)Q1c=h1R1,Q2c=h2R2,Q12=h1-h2R12,
where R1 and R2 are the liquid resistance of valve-1 and valve-2, respectively, and R12 is the liquid resistance of valve-3. Thus, the liquid level model of water tank-1 is acquired:
(19)Gh(s)=H1(s)Qic(s)=b1s+b2a0s2+a1s+a2=b1s/a0+b2/a0s2+a1s/a0+a2/a0=N1(s)M1(s),
where Gh(s) is the liquid level transfer function of tank-1. The transfer function Gh(s) shows that physical process is the second-order process with a negative zero. System parameters are obtained as
(20)b1=R1R12R2C2,b2=R1(R12+R2),a0=R1C1R12R2C2,a1=R1C1R2+R1C1R12+R1R2C2+R12R2C2,(21)a2=R1+R12+R2.


### 3.1. Proposed Design Procedures

Without losing generalities, a third-order controlled process is considered in the unity feedback system. The process(n=3,w=2) and PID controller have forms, respectively,
(22)Gp(s)=b1s2+b2s+b3a0s3+a1s2+a2s+a3=N(s)M(s),Gc(s)=kp+kis+kds.
Thus, open-looped transfer function of control system yields
(23)G(s)=Gc(s)Gp(s)=(kps+ki+kds2)N(s)sM(s).
Then, the error equation and of the control system are obtained:
(24)e(s)r(s)=1(1+Gp(s)Gc(s))⟺D(s)e(s)=sM(s)r(s),D(s)=(sM(s)+N(s)(kps+ki+kds2)),deg{D(s)}=4,
where D(s) is the control system’s characteristic polynomial. The error differential equation can be obtained by inversing Laplace transform of the error equation:
(25)(δM(δ)+N(δ)(kpδ+ki+kdδ2))e(t)=δM(δ)r(t),
where δ is a differential operator. When reference command r(t) is the constant signal, or zero signal (r(t)=0), or even piecewise constant signal, then the error differential equation will yield
(26)(a0+kdb1)e(t)(4)+(a1+kpb1+kdb2)e⃛(t)+(a2+kib1+kpb2+kdb3)e¨(t)+(a3+kib2+kpb3)e˙(t)+kib3e(t)=0.
For symbol simplicity, error differential equation is equally expressed as
(27)m1e(t)(4)+m2e⃛(t)+m3e¨(t)+m4e˙(t)+m5e(t)=0,
where mi(i=1,…,5) is the decision parameter. The decision parameter yields
(28)m1=a0+kdb1,m2=a1+kpb1+kdb2,m5=kib3,m3=a2+kib1+kpb2+kdb3,m4=a3+kib2+kpb3.
Thus, the corresponding state-space model of error differential equation (27) can be obtained:
(29)x˙=Axx=(e(t)e˙(t)e¨(t)e⃛(t))T,A=(O3×1I3×3-m41C1×3),C1×3=(-m31-m21-m11),
where x is state vector, A is state matrix, mi1 is the middle decision variables, and the mi1=mi+1/m1(i=1,2,3,4). Therefore, the optimal control problem is formulated as
(30)J*(kp*,ki*,kd*)=minGc*J(m1,…,m5)=min∫0∞xT(t)Qx(t)dt=min{∑i=03∫0∞qi+1,i+1(e(t)(i))2dt}s.t:(1)x˙=Ax;(2)Re(eig(A))<0,orRe{s∣det(sI-A)=0}<0;(3)ε≥0,R1,1(m1,…,mw+3)>ε,R2,1(m1,…,mw+3)>ε,R3,1(m1,…,mw+3)>ε,R4,1(m1,…,mw+3)>ε,R3,1(m1,…,mw+3)>ε,R4,1(m1,…,mw+3)>ε,R5,1(m1,…,mw+3)>ε;(4)m4=f1(m1,m2,m3),m5=f2(m1,m2,m3);(5)0<kp≤τp0,0<ki≤τi0,0<kd≤τd0.
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(sI-A)=s4+m11s3+m21s2+m31s+m41. Based on LTI stablility theory, constraint (2) is equivalently transformed into the following inequality:(31)Re(eig(A))<0,orRe{s∣det(sI-A)=0}<0⟺m11>0,m11m21-m31>0,m11m21m31-(m31)2-m41(m11)2>0,m41>0.
The constraint (3) assures that the control system possesses the ε-Routh stability. The characteristic polynomial is D(s)=m1s4+m2s3+m3s2+m4s+m5. According to ε-Routh stability definition, the constraint (3) is formulated by the following inequalities:
(32)ε≥0,R1,1(m1,…,m5)>ε,R2,1(m1,…,m5)>ε,R3,1(m1,…,m5)>ε,R4,1(m1,…,m5)>ε,R5,1(m1,…,m5)>ε⟺m1>ε,m2>ε,m5>ε,m1m4-m2m3>εm2,(m2)2m5-m1(m4)2-m2m3m4>ε(m1m4-m2m3).
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 J is equivalently transformed into the performance index
(33)J*=min∫0∞xTQxdt=minxT(0)Px(0),ATP+PA=-Q,
where P is a positive definite real symmetric matrix, and matrix P meets the LAE: ATP+PA=-Q. The performance index J is determined by the matrix P and initial state. Therefore, the proposed optimal control problem (30) can be equivalently transformed into the following NLCO problem:
(34)J*(kp*,ki*,kd*)=minGc*(s)J(m1,…,m5)=minx(0)TPx(0)s.t:(1)ATP+PA=-Q;(2)m11>0,m41>0,m11m21-m31>0,m11m21m31-(m31)2-m41(m11)2>0;(3)m1>ε,m2>ε,m5>ε,m1m4-m2m3>εm2,(m2)2m5-m1(m4)2-m2m3m4>ε(m1m4-m2m3);(4)m4=f1(m1,m2,m3),m5=f2(m1,m2,m3);(5)0<kp≤τp0,0<ki≤τi0,0<kd≤τd0.
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 P to minimize performance index J. We assume that the control system is static at the beginning, then a unit step command r(t)=1(t) is inputted into the system, and the initial states can be obtained: x1(0)=1,x2(0)=0,x3(0)=0, and x4(0)=0. The control weight matrix Q is a diagonal matrix with Q=diag(q1q2q3q4). 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 mi* and proposed PID controller are gotten:
(35)Gc*(s)=kp*+ki*s+kd*s,kp*=m4*-a3-b2ki*b3,ki*=m5*b3,kd*=m1*-a0b1.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 D(s) is sn+1+a1sn+⋯+an-w-2sw+3+m1sw+2+m2sw+1+⋯+mw+2s+mw+3. The polynomial coefficient mk(k=1,2,…,w+3) is a linear mapping function of the PID controller parameters kp, ki, and kd. The parameters kp, ki, and kd are constrained in intervals (0,τp0], (0,τi0], and (0,τd0] respectively. Therefore, the PID controller parameters can be seen as the uncertain bounded parameters in some degree. In this sense, polynomial D(s) can be seen as interval polynomial with three uncertain parameters kp, ki, and kd. To analyze the robustness of ε-Routh stability, without losing generality, interval polynomial F(s) has the general form:
(36)F(s)=sn+p1sn-1+p2sn-2+⋯+pn-1s+pn,
where p1,p2,…,and  pn are the uncertain coefficients of the interval polynomial F(s). The uncertain coefficients p1,p2,…,  and  pn are continuous mapping functions of uncertain bounded parameters q1,q2,…,  and  qm:
(37)p1=f1(q1,…,qm),…,pn=fn(q1,…,qm)Ωq={(q1,…,qm)|[q1-,q1+]×⋯×[qm-,qm+]}q1∈[q1-,q1+],…,qm∈[qm-,qm+],
where Ωq is the uncertain parameter space, qi- and qi+ are the lower bound and upper bound of uncertain parameter qi, and fi(·) is continuous mapping function of polynomial F(s) for uncertain parameters. Based on ε-Routh stability definition for interval polynomial F(s), the following inequalities can be obtained:
(38)R1,1(q1,…,qm)=fr,1(p1,…,pn)>ε,…,Rn+1,1(q1,…,qm)=fr,n+1(p1,…,pn)>εΩr,ε={(q1,…,qm)∣fr,1(p1,…,pn)>ε,…,fr,n+1(p1,…,pn)>ε},(ε≥0),
where fr,i(·) is the element of first column of Routh array. According to the interval polynomial stable criterion in [36, 64], if the disjoint set Ωf of feasible set Ωr,ε and uncertain parameters space Ωq is not a empty set, then the interval polynomial F(s) will satisfy the interval polynomial stable criterion:
(39)min(q1,…,qm)∈Ωf{fr,1(p1,…,pn),…,fr,n+1(p1,…,pn)}>0,(Ωf=(Ωr,ε⋂Ωq)≠∅).
Therefore, the interval polynomial F(s) is Routh-Hurwitz stable in the feasible set Ωf, and the ε-Routh stability could guarantee the robustness of linear systems in a certain degree. For any positive real ε2>ε1>0, two inequalities for ε-Routh stability are obtained:
(40)I:Ω1={(q1,…,qm)∣fr,1(p1,…,pn)>ε1,…,fr,n+1(p1,…,pn)>ε1},II:Ω2={(q1,…,qm)∣fr,1(p1,…,pn)>ε2,…,fr,n+1(p1,…,pn)>ε2}.
The inequality (I) and inequality (II) have the feasible spaces Ω1 and Ω2, and the first inequality of inequalities (I) and (II) has the feasible set
(41)Ω11={(q1,…,qm)∣fr,1(p1,…,pn)>ε1},Ω21={(q1,…,qm)∣fr,1(p1,…,pn)>ε2}.
Then, any point M1* in the feasible set Ω21 will have the relationship
(42)∃M1*=(q11*,q21*,…,qm,1*),fr,1(M1*)=ε2*>ε2.
Likewise, any point M2* in the feasible set Ω11 will have the relationship
(43)∃M2*=(q12*,q22*,…,qm,2*),fr,1(M2*)=ε1*>ε1.
For the feasible set Ω11, there at least a point N1* which will satisfy the relationship(44)Case1:∃N1*=(q-11*,q-21*,…,q-m,1*),ε1<fr,1(N1*)<ε2⟺fr,1(N1*)∈(ε1,ε2)⟹N1*∈Ω11,N1*∉Ω21⟹Ω11⊃Ω21.Case2:fr,1(N1*)∉(ε1,ε2)⟹Ω11=Ω21.
Therefore, the feasible sets Ω1k and Ω2k will have the following relationship:
(45)Case1:Ω1k⊃Ω2k,(∀k=1,2,…,n,n+1)Case2:Ω1k=Ω2k,(∀k=1,2,…,n,n+1)Ω1k={(q1,…,qm)∣fr,k(p1,…,pn)>ε1}Ω2k={(q1,…,qm)∣fr,k(p1,…,pn)>ε2}.
Then, the overall feasible spaces Ω1 and Ω2 will have following relationships:
(46)Ω1⊇Ω2⟺(Ω1⊃Ω2)∪(Ω1=Ω2)Ω1=⋂j=1n+1Ω1j,Ω2=⋂j=1n+1Ω2j.
The overall feasible set will be gradually condensed with the ε increasing. With a real ε sequence εt+1>εt>εt-1>⋯>ε2>ε1, one will get following relationships:
(47)Ωt={(q1,…,qm)∣fr,1(p1,…,pn)>εt,…,fr,n+1(p1,…,pn)>εt}⟺Ωt=(Ωt1∩⋯∩Ωtn+1)=(⋂j=1n+1Ωtj),(t=1,2,3,...,t∈ℤ+)Ωtk={(q1,…,qm)∣fr,k(p1,…,pn)>εt},(k=1,2,…,n,n+1).
Therefore, the feasible set Ωt will has the following relationship:
(48)Ωt+1⊆Ωt⊆⋯⊆Ω2⊆Ω1⟺(Ωt+1ΔtΩt⊂⋯Δ2Ω2Δ1Ω1).⋃(Ωt+1=Ωt=⋯=Ω2=Ω1),∃k,Δk=“⊂”,k∈{1,2,…,t}.
When the feasible set Ωt is strictly condensed with εt, then the limit of feasible space will yield the following relationships:
(49)Case1:limεt→+∞{(q1,…,qm)∣fr,1(p1,…,pn)>εt,…,fr,n+1(p1,…,pn)>εt}=limεt→+∞Ωt=limεt→+∞(Ωt1∩Ωt2∩⋯∩Ωtn+1)=limεt→+∞(⋂j=1n+1Ωtj)=Ω*=(q1,0*q2,0*,…,qm,0*).Case2:limεt→+∞{(q1,…,qm)∣fr,1(p1,…,pn)>εt,…,fr,n+1(p1,…,pn)>εt}=limεt→+∞Ωt=limεt→+∞(Ωt1∩Ωt2∩⋯∩Ωtn+1)=limεt→+∞(⋂j=1n+1Ωtj)=∅.
If the εt 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 εt*  to satisfy the following relationships:
(50)εt⟶+∞,limεt→+∞Ωt=∅⟺∃εt*,εt>εt*,Ωt=∅,0<εt≤εt*,Ωt≠∅,εt*>εt>εt-1>⋯>ε2>ε1≥0,Ωt*⊂Ωt⊂Ωt-1⊂⋯⊂Ω2⊂Ω1when the feasible set has the relationship Ωt+1=Ωt=Ωt-1=⋯=Ω2=Ω1. However, this is a special event which can be depicted by the probability event
(51)P{Ωt+1=Ωt=Ωt-1=⋯=Ω2=Ω1}⟺P{Ωt+1=Ωt,Ωt=Ωt-1,…,Ω2=Ω1}.
It assumes that all the probability events are independent and identically distributed. Hence, the probability can be written as
(52)P{Ωt+1=Ωt,Ωt=Ωt-1,…,Ω2=Ω1}=P{Ωt+1=Ωt}×P{Ωt=Ωt-1}×⋯×P{Ω2=Ω1}=∏j=2t+1P{Ωj=Ωj-1}=(P{Ω2=Ω1})t,(t=1,2,…,t∈ℤ+).
Considering probability event {Ω2=Ω1} first, the probability will have
(53)P{Ω1=Ω2}=P{(Ω11∩⋯∩Ω1n+1)=(Ω21∩⋯∩Ω2n+1)}=P{(⋂j=1n+1Ω1j)=(⋂j=1n+1Ω2j)}⟺P{(Ω21∩⋯∩Ω2n+1)⊂Ω1k,(Ω11∩⋯∩Ω1n+1)⊂Ω2k,Ω1k(∀k=1,2,…,n,n+1)}.
The probability P{Ω2=Ω1} for the probability event {Ω2=Ω1} will yield
(54)P{(Ω21∩Ω22∩⋯∩Ω2n)⊂Ω11}=λ1,…,P{(Ω21∩Ω22∩⋯∩Ω2n)⊂Ω1n+1}=λ1,P{(Ω11∩Ω12∩⋯∩Ω1n)⊂Ω21}=λ2,…,P{(Ω11∩Ω12∩⋯∩Ω1n)⊂Ω2n+1}=λ2,P{Ω1=Ω2}=P{(⋂j=1n+1Ω2j)⊂Ω1k,(∀k=1,2,…,n,n+1)}×P{(⋂j=1n+1Ω1j)⊂Ω2k,(∀k=1,2,…,n,n+1)}=∏k=1n+1P{(⋂j=1n+1Ω2j)⊂Ω1k}×∏k=1n+1P{(⋂j=1n+1Ω1j)⊂Ω2k}=λ1n+1λ2n+1,(0<λ1<1,0<λ2<1).
Therefore, the probability of the probability event yields
(55)P{Ωt+1=Ωt,Ωt=Ωt-1,…,Ω2=Ω1}=λ1(n+1)tλ2(n+1)tP{Ωt+1ΔtΩtΔt-1Ωt-1⋯Δ2Ω2Δ1Ω1}=1-λ1(n+1)tλ2(n+1)t.


### 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:
(56)minf(x)s.t:h(x)=0,g(x)≤0,
where f(x):Rn→R,h(x):Rn→Rm, and g(x):Rn→Rq are the smooth and differentiable functions, x is the decision variable, and n, m, and q 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:
(57)minf(x)-υ∑i=1qlnδis.t:h(x)=0,g(x)+δ=0,
where υ>0 is the barrier parameter, the slack vector δ=(δ1,δ2,…,δq)T>0 is set to be positive, and g(x) is an expanded constraint inequality. It introduces the Lagrange multipliers y and z for the barrier problem (57):
(58)L(x,y,z,δ)=f(x)-υ∑i=1qlnδi+yT(g(x)+δ)+zTh(x),
where L(x,y,z,δ) is Lagrange function, y=(y1,y2,…,yq)T and z=(z1,z2,…,zm)T are Lagrange multipliers for constraints g(x)+δ and h(x), respectively.  Based on the Karush-Kuhn-Tucker (KKT) optimality conditions [57–59], the optimality conditions for NLCO problem (56) can be expressed as
(59)∇xL(x,y,z,δ)=∇f(x)+(∇g(x))Ty+(∇h(x))Tz=0∇δL(x,y,z,δ)=-υSδ-1e+y=0⟺-υe+SδYe=0∇yL(x,y,z,δ)=g(x)+δ=0∇zL(x,y,z,δ)=h(x)=0,
where Sδ=diag(δ1,δ2,…,δq) is the diagonal matrix and its elements are the components of the vector δ, and e is a vector of all ones, ∇h(x) and ∇g(x) are the Jacobian matrices of the vectors h(x) and g(x), respectively, ∇f(x) is the grand of function f(x), and Y=diag(y1,y2,…,yq) is a diagonal matrix and its elements are the components of vector y. The system (59) is the KKT condition of the NLCO problem (56). When in the search approach, it should remain δ, y>0. To obtain the iteration direction, it can make the point (x+Δx,δ+Δδ,z+Δz,y+Δy) satisfying the KKT conditions(59); then the following system will be obtained:
(60)(∇2f(x)+∇2g(x)Ty+∇2h(x)Tz)Δx+∇g(x)TΔy+∇h(x)TΔz+(∇f(x)+∇g(x)Ty+∇h(x)Tz)=0-υe+SδYe+SδΔy+YΔδ=0g(x)+∇g(x)Δx+δ+Δδ=0,∇h(x)Δx+h(x)=0.
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:
(61)(H(x,y,z)0∇h(x)T∇g(x)T0Sδ-1Y0I∇h(x)000∇g(x)I00)(ΔxΔδΔyΔz)=(-∇f(x)-∇g(x)Ty-∇h(x)TzυSδ-1e-y-h(x)-g(x)-δ)(62)H(x,y,z)=∇2f(x)+∇2g(x)Ty+∇2h(x)Tz,
where H(x,y,z) 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)(x,δ,z,y)⟵(x,δ,z,y)+ζ1(Δx,Δδ,Δz,Δy),
where ζ1 is the step size. Choosing the step size ζ1 holds the δ, y>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(x(0,),δ(0),z(0),y(0)) in the feasible region set ={x∣h(x)=0,g(x)≤0}, and the δ(0)>0,y(0)>0,k=0.Step 2.
Construct current iterate; we have the current iterate valuesx(k), δ(k), z(k), and y(k) of the primal variable x, the of the slack variable δ, the multipliers y and z, respectively.Step 3.
Calculate the Hessian matrixH(x,y,z) of the Lagrange system L(x,y,z,δ), and the Jacobian matrices ∇h(x) and ∇g(x) are of the vectors h(x) and g(x) in the current iterate (x(k),δ(k),z(k),y(k)).Step 4.
Solve the linear system (61) and construct the iterate direction (Δx,Δδ,Δz,Δy). Solving the linear matrix equation (61), we obtain the primal solution Δx, multipliers solution Δz, Δy, and also the slack variable solution Δδ.Step 5.
Choosing the step sizeζ1 holds the δ, y>0 in search process, ζ1∈[0,1]. Update the iterate values: (x(k+1),δ(k+1),z(k+1),y(k+1))←(x(k),δ(k),z(k),y(k))+ζ1(Δx,Δδ,Δz,Δy),k←k+1.Step 6.
Check the ending conditions in regionΩ. If it does not satisfied, go to Step 3, else minimum fmin of NLCO problem is obtained in feasible region Ω.Step 7.
End.


## 3.1. Proposed Design Procedures

Without losing generalities, a third-order controlled process is considered in the unity feedback system. The process(n=3,w=2) and PID controller have forms, respectively,
(22)Gp(s)=b1s2+b2s+b3a0s3+a1s2+a2s+a3=N(s)M(s),Gc(s)=kp+kis+kds.
Thus, open-looped transfer function of control system yields
(23)G(s)=Gc(s)Gp(s)=(kps+ki+kds2)N(s)sM(s).
Then, the error equation and of the control system are obtained:
(24)e(s)r(s)=1(1+Gp(s)Gc(s))⟺D(s)e(s)=sM(s)r(s),D(s)=(sM(s)+N(s)(kps+ki+kds2)),deg{D(s)}=4,
where D(s) is the control system’s characteristic polynomial. The error differential equation can be obtained by inversing Laplace transform of the error equation:
(25)(δM(δ)+N(δ)(kpδ+ki+kdδ2))e(t)=δM(δ)r(t),
where δ is a differential operator. When reference command r(t) is the constant signal, or zero signal (r(t)=0), or even piecewise constant signal, then the error differential equation will yield
(26)(a0+kdb1)e(t)(4)+(a1+kpb1+kdb2)e⃛(t)+(a2+kib1+kpb2+kdb3)e¨(t)+(a3+kib2+kpb3)e˙(t)+kib3e(t)=0.
For symbol simplicity, error differential equation is equally expressed as
(27)m1e(t)(4)+m2e⃛(t)+m3e¨(t)+m4e˙(t)+m5e(t)=0,
where mi(i=1,…,5) is the decision parameter. The decision parameter yields
(28)m1=a0+kdb1,m2=a1+kpb1+kdb2,m5=kib3,m3=a2+kib1+kpb2+kdb3,m4=a3+kib2+kpb3.
Thus, the corresponding state-space model of error differential equation (27) can be obtained:
(29)x˙=Axx=(e(t)e˙(t)e¨(t)e⃛(t))T,A=(O3×1I3×3-m41C1×3),C1×3=(-m31-m21-m11),
where x is state vector, A is state matrix, mi1 is the middle decision variables, and the mi1=mi+1/m1(i=1,2,3,4). Therefore, the optimal control problem is formulated as
(30)J*(kp*,ki*,kd*)=minGc*J(m1,…,m5)=min∫0∞xT(t)Qx(t)dt=min{∑i=03∫0∞qi+1,i+1(e(t)(i))2dt}s.t:(1)x˙=Ax;(2)Re(eig(A))<0,orRe{s∣det(sI-A)=0}<0;(3)ε≥0,R1,1(m1,…,mw+3)>ε,R2,1(m1,…,mw+3)>ε,R3,1(m1,…,mw+3)>ε,R4,1(m1,…,mw+3)>ε,R3,1(m1,…,mw+3)>ε,R4,1(m1,…,mw+3)>ε,R5,1(m1,…,mw+3)>ε;(4)m4=f1(m1,m2,m3),m5=f2(m1,m2,m3);(5)0<kp≤τp0,0<ki≤τi0,0<kd≤τd0.
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(sI-A)=s4+m11s3+m21s2+m31s+m41. Based on LTI stablility theory, constraint (2) is equivalently transformed into the following inequality:(31)Re(eig(A))<0,orRe{s∣det(sI-A)=0}<0⟺m11>0,m11m21-m31>0,m11m21m31-(m31)2-m41(m11)2>0,m41>0.
The constraint (3) assures that the control system possesses the ε-Routh stability. The characteristic polynomial is D(s)=m1s4+m2s3+m3s2+m4s+m5. According to ε-Routh stability definition, the constraint (3) is formulated by the following inequalities:
(32)ε≥0,R1,1(m1,…,m5)>ε,R2,1(m1,…,m5)>ε,R3,1(m1,…,m5)>ε,R4,1(m1,…,m5)>ε,R5,1(m1,…,m5)>ε⟺m1>ε,m2>ε,m5>ε,m1m4-m2m3>εm2,(m2)2m5-m1(m4)2-m2m3m4>ε(m1m4-m2m3).
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 J is equivalently transformed into the performance index
(33)J*=min∫0∞xTQxdt=minxT(0)Px(0),ATP+PA=-Q,
where P is a positive definite real symmetric matrix, and matrix P meets the LAE: ATP+PA=-Q. The performance index J is determined by the matrix P and initial state. Therefore, the proposed optimal control problem (30) can be equivalently transformed into the following NLCO problem:
(34)J*(kp*,ki*,kd*)=minGc*(s)J(m1,…,m5)=minx(0)TPx(0)s.t:(1)ATP+PA=-Q;(2)m11>0,m41>0,m11m21-m31>0,m11m21m31-(m31)2-m41(m11)2>0;(3)m1>ε,m2>ε,m5>ε,m1m4-m2m3>εm2,(m2)2m5-m1(m4)2-m2m3m4>ε(m1m4-m2m3);(4)m4=f1(m1,m2,m3),m5=f2(m1,m2,m3);(5)0<kp≤τp0,0<ki≤τi0,0<kd≤τd0.
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 P to minimize performance index J. We assume that the control system is static at the beginning, then a unit step command r(t)=1(t) is inputted into the system, and the initial states can be obtained: x1(0)=1,x2(0)=0,x3(0)=0, and x4(0)=0. The control weight matrix Q is a diagonal matrix with Q=diag(q1q2q3q4). 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 mi* and proposed PID controller are gotten:
(35)Gc*(s)=kp*+ki*s+kd*s,kp*=m4*-a3-b2ki*b3,ki*=m5*b3,kd*=m1*-a0b1.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 D(s) is sn+1+a1sn+⋯+an-w-2sw+3+m1sw+2+m2sw+1+⋯+mw+2s+mw+3. The polynomial coefficient mk(k=1,2,…,w+3) is a linear mapping function of the PID controller parameters kp, ki, and kd. The parameters kp, ki, and kd are constrained in intervals (0,τp0], (0,τi0], and (0,τd0] respectively. Therefore, the PID controller parameters can be seen as the uncertain bounded parameters in some degree. In this sense, polynomial D(s) can be seen as interval polynomial with three uncertain parameters kp, ki, and kd. To analyze the robustness of ε-Routh stability, without losing generality, interval polynomial F(s) has the general form:
(36)F(s)=sn+p1sn-1+p2sn-2+⋯+pn-1s+pn,
where p1,p2,…,and  pn are the uncertain coefficients of the interval polynomial F(s). The uncertain coefficients p1,p2,…,  and  pn are continuous mapping functions of uncertain bounded parameters q1,q2,…,  and  qm:
(37)p1=f1(q1,…,qm),…,pn=fn(q1,…,qm)Ωq={(q1,…,qm)|[q1-,q1+]×⋯×[qm-,qm+]}q1∈[q1-,q1+],…,qm∈[qm-,qm+],
where Ωq is the uncertain parameter space, qi- and qi+ are the lower bound and upper bound of uncertain parameter qi, and fi(·) is continuous mapping function of polynomial F(s) for uncertain parameters. Based on ε-Routh stability definition for interval polynomial F(s), the following inequalities can be obtained:
(38)R1,1(q1,…,qm)=fr,1(p1,…,pn)>ε,…,Rn+1,1(q1,…,qm)=fr,n+1(p1,…,pn)>εΩr,ε={(q1,…,qm)∣fr,1(p1,…,pn)>ε,…,fr,n+1(p1,…,pn)>ε},(ε≥0),
where fr,i(·) is the element of first column of Routh array. According to the interval polynomial stable criterion in [36, 64], if the disjoint set Ωf of feasible set Ωr,ε and uncertain parameters space Ωq is not a empty set, then the interval polynomial F(s) will satisfy the interval polynomial stable criterion:
(39)min(q1,…,qm)∈Ωf{fr,1(p1,…,pn),…,fr,n+1(p1,…,pn)}>0,(Ωf=(Ωr,ε⋂Ωq)≠∅).
Therefore, the interval polynomial F(s) is Routh-Hurwitz stable in the feasible set Ωf, and the ε-Routh stability could guarantee the robustness of linear systems in a certain degree. For any positive real ε2>ε1>0, two inequalities for ε-Routh stability are obtained:
(40)I:Ω1={(q1,…,qm)∣fr,1(p1,…,pn)>ε1,…,fr,n+1(p1,…,pn)>ε1},II:Ω2={(q1,…,qm)∣fr,1(p1,…,pn)>ε2,…,fr,n+1(p1,…,pn)>ε2}.
The inequality (I) and inequality (II) have the feasible spaces Ω1 and Ω2, and the first inequality of inequalities (I) and (II) has the feasible set
(41)Ω11={(q1,…,qm)∣fr,1(p1,…,pn)>ε1},Ω21={(q1,…,qm)∣fr,1(p1,…,pn)>ε2}.
Then, any point M1* in the feasible set Ω21 will have the relationship
(42)∃M1*=(q11*,q21*,…,qm,1*),fr,1(M1*)=ε2*>ε2.
Likewise, any point M2* in the feasible set Ω11 will have the relationship
(43)∃M2*=(q12*,q22*,…,qm,2*),fr,1(M2*)=ε1*>ε1.
For the feasible set Ω11, there at least a point N1* which will satisfy the relationship(44)Case1:∃N1*=(q-11*,q-21*,…,q-m,1*),ε1<fr,1(N1*)<ε2⟺fr,1(N1*)∈(ε1,ε2)⟹N1*∈Ω11,N1*∉Ω21⟹Ω11⊃Ω21.Case2:fr,1(N1*)∉(ε1,ε2)⟹Ω11=Ω21.
Therefore, the feasible sets Ω1k and Ω2k will have the following relationship:
(45)Case1:Ω1k⊃Ω2k,(∀k=1,2,…,n,n+1)Case2:Ω1k=Ω2k,(∀k=1,2,…,n,n+1)Ω1k={(q1,…,qm)∣fr,k(p1,…,pn)>ε1}Ω2k={(q1,…,qm)∣fr,k(p1,…,pn)>ε2}.
Then, the overall feasible spaces Ω1 and Ω2 will have following relationships:
(46)Ω1⊇Ω2⟺(Ω1⊃Ω2)∪(Ω1=Ω2)Ω1=⋂j=1n+1Ω1j,Ω2=⋂j=1n+1Ω2j.
The overall feasible set will be gradually condensed with the ε increasing. With a real ε sequence εt+1>εt>εt-1>⋯>ε2>ε1, one will get following relationships:
(47)Ωt={(q1,…,qm)∣fr,1(p1,…,pn)>εt,…,fr,n+1(p1,…,pn)>εt}⟺Ωt=(Ωt1∩⋯∩Ωtn+1)=(⋂j=1n+1Ωtj),(t=1,2,3,...,t∈ℤ+)Ωtk={(q1,…,qm)∣fr,k(p1,…,pn)>εt},(k=1,2,…,n,n+1).
Therefore, the feasible set Ωt will has the following relationship:
(48)Ωt+1⊆Ωt⊆⋯⊆Ω2⊆Ω1⟺(Ωt+1ΔtΩt⊂⋯Δ2Ω2Δ1Ω1).⋃(Ωt+1=Ωt=⋯=Ω2=Ω1),∃k,Δk=“⊂”,k∈{1,2,…,t}.
When the feasible set Ωt is strictly condensed with εt, then the limit of feasible space will yield the following relationships:
(49)Case1:limεt→+∞{(q1,…,qm)∣fr,1(p1,…,pn)>εt,…,fr,n+1(p1,…,pn)>εt}=limεt→+∞Ωt=limεt→+∞(Ωt1∩Ωt2∩⋯∩Ωtn+1)=limεt→+∞(⋂j=1n+1Ωtj)=Ω*=(q1,0*q2,0*,…,qm,0*).Case2:limεt→+∞{(q1,…,qm)∣fr,1(p1,…,pn)>εt,…,fr,n+1(p1,…,pn)>εt}=limεt→+∞Ωt=limεt→+∞(Ωt1∩Ωt2∩⋯∩Ωtn+1)=limεt→+∞(⋂j=1n+1Ωtj)=∅.
If the εt 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 εt*  to satisfy the following relationships:
(50)εt⟶+∞,limεt→+∞Ωt=∅⟺∃εt*,εt>εt*,Ωt=∅,0<εt≤εt*,Ωt≠∅,εt*>εt>εt-1>⋯>ε2>ε1≥0,Ωt*⊂Ωt⊂Ωt-1⊂⋯⊂Ω2⊂Ω1when the feasible set has the relationship Ωt+1=Ωt=Ωt-1=⋯=Ω2=Ω1. However, this is a special event which can be depicted by the probability event
(51)P{Ωt+1=Ωt=Ωt-1=⋯=Ω2=Ω1}⟺P{Ωt+1=Ωt,Ωt=Ωt-1,…,Ω2=Ω1}.
It assumes that all the probability events are independent and identically distributed. Hence, the probability can be written as
(52)P{Ωt+1=Ωt,Ωt=Ωt-1,…,Ω2=Ω1}=P{Ωt+1=Ωt}×P{Ωt=Ωt-1}×⋯×P{Ω2=Ω1}=∏j=2t+1P{Ωj=Ωj-1}=(P{Ω2=Ω1})t,(t=1,2,…,t∈ℤ+).
Considering probability event {Ω2=Ω1} first, the probability will have
(53)P{Ω1=Ω2}=P{(Ω11∩⋯∩Ω1n+1)=(Ω21∩⋯∩Ω2n+1)}=P{(⋂j=1n+1Ω1j)=(⋂j=1n+1Ω2j)}⟺P{(Ω21∩⋯∩Ω2n+1)⊂Ω1k,(Ω11∩⋯∩Ω1n+1)⊂Ω2k,Ω1k(∀k=1,2,…,n,n+1)}.
The probability P{Ω2=Ω1} for the probability event {Ω2=Ω1} will yield
(54)P{(Ω21∩Ω22∩⋯∩Ω2n)⊂Ω11}=λ1,…,P{(Ω21∩Ω22∩⋯∩Ω2n)⊂Ω1n+1}=λ1,P{(Ω11∩Ω12∩⋯∩Ω1n)⊂Ω21}=λ2,…,P{(Ω11∩Ω12∩⋯∩Ω1n)⊂Ω2n+1}=λ2,P{Ω1=Ω2}=P{(⋂j=1n+1Ω2j)⊂Ω1k,(∀k=1,2,…,n,n+1)}×P{(⋂j=1n+1Ω1j)⊂Ω2k,(∀k=1,2,…,n,n+1)}=∏k=1n+1P{(⋂j=1n+1Ω2j)⊂Ω1k}×∏k=1n+1P{(⋂j=1n+1Ω1j)⊂Ω2k}=λ1n+1λ2n+1,(0<λ1<1,0<λ2<1).
Therefore, the probability of the probability event yields
(55)P{Ωt+1=Ωt,Ωt=Ωt-1,…,Ω2=Ω1}=λ1(n+1)tλ2(n+1)tP{Ωt+1ΔtΩtΔt-1Ωt-1⋯Δ2Ω2Δ1Ω1}=1-λ1(n+1)tλ2(n+1)t.


## 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:
(56)minf(x)s.t:h(x)=0,g(x)≤0,
where f(x):Rn→R,h(x):Rn→Rm, and g(x):Rn→Rq are the smooth and differentiable functions, x is the decision variable, and n, m, and q 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:
(57)minf(x)-υ∑i=1qlnδis.t:h(x)=0,g(x)+δ=0,
where υ>0 is the barrier parameter, the slack vector δ=(δ1,δ2,…,δq)T>0 is set to be positive, and g(x) is an expanded constraint inequality. It introduces the Lagrange multipliers y and z for the barrier problem (57):
(58)L(x,y,z,δ)=f(x)-υ∑i=1qlnδi+yT(g(x)+δ)+zTh(x),
where L(x,y,z,δ) is Lagrange function, y=(y1,y2,…,yq)T and z=(z1,z2,…,zm)T are Lagrange multipliers for constraints g(x)+δ and h(x), respectively.  Based on the Karush-Kuhn-Tucker (KKT) optimality conditions [57–59], the optimality conditions for NLCO problem (56) can be expressed as
(59)∇xL(x,y,z,δ)=∇f(x)+(∇g(x))Ty+(∇h(x))Tz=0∇δL(x,y,z,δ)=-υSδ-1e+y=0⟺-υe+SδYe=0∇yL(x,y,z,δ)=g(x)+δ=0∇zL(x,y,z,δ)=h(x)=0,
where Sδ=diag(δ1,δ2,…,δq) is the diagonal matrix and its elements are the components of the vector δ, and e is a vector of all ones, ∇h(x) and ∇g(x) are the Jacobian matrices of the vectors h(x) and g(x), respectively, ∇f(x) is the grand of function f(x), and Y=diag(y1,y2,…,yq) is a diagonal matrix and its elements are the components of vector y. The system (59) is the KKT condition of the NLCO problem (56). When in the search approach, it should remain δ, y>0. To obtain the iteration direction, it can make the point (x+Δx,δ+Δδ,z+Δz,y+Δy) satisfying the KKT conditions(59); then the following system will be obtained:
(60)(∇2f(x)+∇2g(x)Ty+∇2h(x)Tz)Δx+∇g(x)TΔy+∇h(x)TΔz+(∇f(x)+∇g(x)Ty+∇h(x)Tz)=0-υe+SδYe+SδΔy+YΔδ=0g(x)+∇g(x)Δx+δ+Δδ=0,∇h(x)Δx+h(x)=0.
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:
(61)(H(x,y,z)0∇h(x)T∇g(x)T0Sδ-1Y0I∇h(x)000∇g(x)I00)(ΔxΔδΔyΔz)=(-∇f(x)-∇g(x)Ty-∇h(x)TzυSδ-1e-y-h(x)-g(x)-δ)(62)H(x,y,z)=∇2f(x)+∇2g(x)Ty+∇2h(x)Tz,
where H(x,y,z) 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)(x,δ,z,y)⟵(x,δ,z,y)+ζ1(Δx,Δδ,Δz,Δy),
where ζ1 is the step size. Choosing the step size ζ1 holds the δ, y>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(x(0,),δ(0),z(0),y(0)) in the feasible region set ={x∣h(x)=0,g(x)≤0}, and the δ(0)>0,y(0)>0,k=0.Step 2.
Construct current iterate; we have the current iterate valuesx(k), δ(k), z(k), and y(k) of the primal variable x, the of the slack variable δ, the multipliers y and z, respectively.Step 3.
Calculate the Hessian matrixH(x,y,z) of the Lagrange system L(x,y,z,δ), and the Jacobian matrices ∇h(x) and ∇g(x) are of the vectors h(x) and g(x) in the current iterate (x(k),δ(k),z(k),y(k)).Step 4.
Solve the linear system (61) and construct the iterate direction (Δx,Δδ,Δz,Δy). Solving the linear matrix equation (61), we obtain the primal solution Δx, multipliers solution Δz, Δy, and also the slack variable solution Δδ.Step 5.
Choosing the step sizeζ1 holds the δ, y>0 in search process, ζ1∈[0,1]. Update the iterate values: (x(k+1),δ(k+1),z(k+1),y(k+1))←(x(k),δ(k),z(k),y(k))+ζ1(Δx,Δδ,Δz,Δy),k←k+1.Step 6.
Check the ending conditions in regionΩ. If it does not satisfied, go to Step 3, else minimum fmin 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 Gp1(s) [7], Gp2(s) [15], and Gp3(s) [6], respectively. Gp1(s) is a third-order all-pole process, and Gp2(s) and Gp3(s) are the SOPTD and FOPTD processes, respectively,
(64)Gp1(s)=15(s2+0.9s+5)(s+3),Gp2(s)=e-0.5s(s+1)2,Gp3(s)=4e-10s10s+1.In order to use the proposed method, time delay of plant.2 and plant.3 is approximated ase-0.5s≈(1-0.25s)/(1+0.25s) and e-10s≈(s2-0.6s+0.12)/(s2+0.6s+0.12), respectively. Thus, transfer functions of Gp2(s) and Gp3(s) are the (-s+4)/(s3+6s2+9s+4) and (0.4s2-0.24s+0.048)/(s3+0.7s2+0.18s+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.
Plants
Control weight matrixQ = diag (q1q2q3q4)
Optimal PID controller
ε-Routh stability
k
p
*
k
i
*
k
d
*
Plant.1
diag (5 5 1 1)
0.13121
0.58138
0.58655
ε
=
0
diag (5 2 1 1)
0.01633
0.58138
0.51933
diag (10 2 1 1)
0.23263
0.82219
0.64449
diag (10 10 1 1)
0.47606
0.82219
0.77901
diag (10 10 4 2)
0.15570
0.58138
0.64101
Plant.2
diag (0.2 1 1 1)
1.15698
0.52966
0.86781
ε
=
0.5
diag (0.2 1 1 0.2)
2.16962
0.84500
1.82140
diag (0.2 2 2 0.5)
1.78668
0.58700
1.59290
diag (0.1 0.8 0.8 0.2)
1.85901
0.64942
1.62251
diag (0.4 0.5 5 1)
1.22806
0.62168
1.42095
Plant.3
diag (0.04 27 158 175)
0.32447
0.01537
1.45256
ε
=
0
diag (0.04 28 120 80)
0.33078
0.01543
1.48167
diag (0.02 25 150 175)
0.30719
0.01193
1.39141
diag (0.05 27 120 100)
0.33727
0.01693
1.50359
diag (0.06 30 158 175)
0.33603
0.01737
1.49509For the proposed method, the plant.1’s step response and system performances under various control weight matrices are shown in Figure4 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 q1, q2, q3, and q4. The results show that the response velocity will be heightened by increasing weight factor q1 or reducing weight factors q2, q3, and q4. When control weight factor q1 (with q2, q3, and q4 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 q2, q3, and q4 (with q1 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.
Control weight matrix
Plant.1
Rise time
Delay time
Peak time
Setting time
Peak overshoot
Q = diag (q1q2q3q4)
k
p
*
k
i
*
k
d
*
t
r (sec)
t
d (sec)
t
p (sec)
t
s (sec)
M
p (%)
diag (5 5 1 1)
0.13121
0.58138
0.58655
4.060
1.685
5.31
3.54
4.08
diag (5 2 1 1)
0.01633
0.58138
0.51933
3.681
1.801
5.00
6.07
6.69
diag (10 2 1 1)
0.23263
0.82219
0.64449
2.850
1.318
4.02
5.35
8.47
diag (10 10 1 1)
0.47606
0.82219
0.77901
3.161
1.140
4.38
2.756
4.58
diag (10 10 4 2)
0.15570
0.58138
0.64101
4.120
1.691
5.37
3.60
4.23Step response of the proposed method under various control weight matrices (for plant.1).
(a)
(b)Fornth-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 q1 (with other factors fixed) or reducing other weight factors (with q1 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 q1 (with other factors fixed) or increasing other factors (with q1 fixed). Therefore, system performances are affected by control weight factors q1,q2,…,qn+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).
Control matrixQ and ε value
k
p
k
i
k
d
Rise timetr (sec)
Delay timetd (sec)
Peak timetp (sec)
Setting timets (sec)
Peak overshootMp (%)
Plants
Q = diag (0.2 1 1 1) and ε = 0.5
1.15698
0.52966
0.86781
4.778
1.841
6.622
4.00
3.29
Plant.2
Q = diag (0.2 1 1 1) and ε = 2.5
1.45019
0.62500
1.30523
4.201
1.655
6.451
7.20
5.41
Plant.2
Q = diag (0.2 1 1 1) and ε = 3.0
1.79987
0.75000
1.78883
3.727
1.488
5.961
8.08
7.61
Plant.2
Q = diag (0.2 1 1 1) and ε = 4.0
2.31710
1.00000
2.00000
3.151
1.316
4.201
7.35
11.00
Plant.2
Q = diag (0.2 1 1 1) and ε = 5.0
2.46310
1.25000
1.00000
1.795
1.266
2.587
3.97
28.46
Plant.2Step response of the proposed method under variousε-stability (for plant.2).
(a)
(b)


### 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 Figures6, 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.
Methods
k
p
k
i
k
d
Rise timetr (sec)
Delay timetd (sec)
Peak timetp (sec)
Setting timets (sec)
Peak overshootMp (%)
Plants
Proposed method
0.13121
0.58138
0.58655
4.060
1.685
5.310
3.54
4.08
Plant.1
Z-N method
0.5971
1.0750
0.2687
1.371
0.877
1.643
6.41
11.08
Plant.1
IA-ISE method
0.4479
0.9994
0.3340
1.678
0.983
4.015
3.13
4.24
Plant.1
IA-ITSE method
0.7377
0.9886
0.3557
1.429
0.845
3.506
4.88
5.27
Plant.1
Proposed method
1.15698
0.52966
0.86781
4.778
1.841
6.622
4.00
3.29
Plant.2
Z-N method
2.8130
1.7194
1.1505
1.620
1.195
2.458
3.81
43.09
Plant.2
Gain and phase method
2.3870
0.8747
0.9476
1.938
1.302
2.478
3.34
15.02
Plant.2
A-H method
2.5778
1.4518
1.2941
1.763
1.286
2.537
4.43
33.65
Plant.2
Proposed method
0.32447
0.01537
1.45256
19.486
14.306
20.87
60.56
6.67
Plant.3
Z-N method
Z-N method (C-C)
0.3060
0.0205
1.1384
19.238
14.367
24.98
51.60
18.18
Plant.3
Z-N method (PRC)
0.3000
0.0150
0.0600
20.286
14.635
26.43
56.06
21.17
Plant.3
Cohen-Coon method
0.3958
0.0219
1.2178
17.252
13.400
21.43
67.05
38.06
Plant.3
A-H method (real PID controller form) with derivative filter constantTf=0.04947.Step response of plant.1 under different tuning methods.
(a)
(b)Step response of plant.2 under different tuning methods.
(a)
(b)Step response of plant.3 under different tuning methods.
(a)
(b)The plant.2 is a SOPTD plant with dead time ratioτ/T=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τ/T=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 Figures9, 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 Figures14, 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 h1(s)/Qic(s)=(0.4465s+0.3064)/(s2+1.340s+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 controllerkp*=10.000, ki*=0.32549, kd*=0.001, and ε=0 in set point 300 mm.Figure 15
Experimental response of the CWTS with the PID controllerskp*=10.000, ki*=0.10118, kd*=0.001, and ε=0 in set point 300 mm.Figure 16
Experimental response of the CWTS with the PID controllerskp*=10.000, ki*=0.22986, kd*=0.001, and ε=0 in set point 300 mm.Figure 17
Experimental response of the CWTS with the PID controllerskp*=6.000, ki*=0.13817, kd*=0.001, and ε=0 in set point 300 mm.Figure 18
Experimental response of the CWTS with the PID controllerskp*=10.000, ki*=0.32549, kd*=0.001, and ε=0 in set point 160 mm.Figure 19
Experimental response of the CWTS with the PID controllerskp*=10.000, ki*=0.10118, kd*=0.001, and ε=0 in set point 160 mm.Figure 20
Experimental response of the CWTS with the PID controllerskp*=10.000, ki*=0.22986, kd*=0.001, and ε=0 in set point 160 mm.Figure 21
Experimental response of the CWTS with the PID controllerskp*=6.000, ki*=0.13817, kd*=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 (kp=10.000,ki=0.32549,kd=0.001,ε=0;kp=10.000,ki=0.10118,kd=0.001,ε=0;kp=10.000,ki=0.22986,kd=0.001,ε=0;kp=6.000,ki=0.13817,kd=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 h1 is set at 300.00 mm and second group experiments’ setting point h1 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 h1, the solid magenta line is the dynamic response of the liquid level h2, 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 Figures18–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.


## 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 Figures6, 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.
Methods
k
p
k
i
k
d
Rise timetr (sec)
Delay timetd (sec)
Peak timetp (sec)
Setting timets (sec)
Peak overshootMp (%)
Plants
Proposed method
0.13121
0.58138
0.58655
4.060
1.685
5.310
3.54
4.08
Plant.1
Z-N method
0.5971
1.0750
0.2687
1.371
0.877
1.643
6.41
11.08
Plant.1
IA-ISE method
0.4479
0.9994
0.3340
1.678
0.983
4.015
3.13
4.24
Plant.1
IA-ITSE method
0.7377
0.9886
0.3557
1.429
0.845
3.506
4.88
5.27
Plant.1
Proposed method
1.15698
0.52966
0.86781
4.778
1.841
6.622
4.00
3.29
Plant.2
Z-N method
2.8130
1.7194
1.1505
1.620
1.195
2.458
3.81
43.09
Plant.2
Gain and phase method
2.3870
0.8747
0.9476
1.938
1.302
2.478
3.34
15.02
Plant.2
A-H method
2.5778
1.4518
1.2941
1.763
1.286
2.537
4.43
33.65
Plant.2
Proposed method
0.32447
0.01537
1.45256
19.486
14.306
20.87
60.56
6.67
Plant.3
Z-N method
Z-N method (C-C)
0.3060
0.0205
1.1384
19.238
14.367
24.98
51.60
18.18
Plant.3
Z-N method (PRC)
0.3000
0.0150
0.0600
20.286
14.635
26.43
56.06
21.17
Plant.3
Cohen-Coon method
0.3958
0.0219
1.2178
17.252
13.400
21.43
67.05
38.06
Plant.3
A-H method (real PID controller form) with derivative filter constantTf=0.04947.Step response of plant.1 under different tuning methods.
(a)
(b)Step response of plant.2 under different tuning methods.
(a)
(b)Step response of plant.3 under different tuning methods.
(a)
(b)The plant.2 is a SOPTD plant with dead time ratioτ/T=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τ/T=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 Figures9, 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 Figures14, 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 h1(s)/Qic(s)=(0.4465s+0.3064)/(s2+1.340s+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 controllerkp*=10.000, ki*=0.32549, kd*=0.001, and ε=0 in set point 300 mm.Figure 15
Experimental response of the CWTS with the PID controllerskp*=10.000, ki*=0.10118, kd*=0.001, and ε=0 in set point 300 mm.Figure 16
Experimental response of the CWTS with the PID controllerskp*=10.000, ki*=0.22986, kd*=0.001, and ε=0 in set point 300 mm.Figure 17
Experimental response of the CWTS with the PID controllerskp*=6.000, ki*=0.13817, kd*=0.001, and ε=0 in set point 300 mm.Figure 18
Experimental response of the CWTS with the PID controllerskp*=10.000, ki*=0.32549, kd*=0.001, and ε=0 in set point 160 mm.Figure 19
Experimental response of the CWTS with the PID controllerskp*=10.000, ki*=0.10118, kd*=0.001, and ε=0 in set point 160 mm.Figure 20
Experimental response of the CWTS with the PID controllerskp*=10.000, ki*=0.22986, kd*=0.001, and ε=0 in set point 160 mm.Figure 21
Experimental response of the CWTS with the PID controllerskp*=6.000, ki*=0.13817, kd*=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 (kp=10.000,ki=0.32549,kd=0.001,ε=0;kp=10.000,ki=0.10118,kd=0.001,ε=0;kp=10.000,ki=0.22986,kd=0.001,ε=0;kp=6.000,ki=0.13817,kd=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 h1 is set at 300.00 mm and second group experiments’ setting point h1 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 h1, the solid magenta line is the dynamic response of the liquid level h2, 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 Figures18–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 nth-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 q1, q2, q3, and q4. The peak overshoot and response velocity will be heightened by increasing weight factor q1 (with other factors fixed). Likewise, the overshoot and response velocity will be heightened by reducing weight factors q2, q3, and q4 (with q1 fixed). On the contrary, the overshoot and response velocity will be reduced by reducing weight factor q1 (with other factors fixed) or increasing weight factors q2, q3, and q4 (with q1 fixed). Through the study, for nth-order processes, the above conclusions still hold; it can be inferred that system performances are affected by the control weight factors q1,q2,…,qn+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.



---

*Source: 582879-2013-05-16.xml*