Computation of stabilizing PI and PID controllers by using Kronecker summation method

Energy Conversion and Management 50 (2009) 1821–1827

Contents lists available at ScienceDirect

Energy Conversion and Management journal homepage: www.elsevier.com/locate/enconman

Computation of stabilizing PI and PID controllers by using Kronecker summation method Jian’an Fang, Da Zheng *, Zhengyun Ren College of Information Science and Technology, Donghua University, Shanghai 201620, China

a r t i c l e

i n f o

Article history: Received 15 December 2007 Received in revised form 3 August 2008 Accepted 8 March 2009 Available online 11 April 2009 Keywords: Stabilization PI control PID control Relative stabilization Kronecker summation Coefficient uncertainty

a b s t r a c t In this paper, a new method for computation of all stabilizing PI controllers is given. The method is based on the plant model in time domain, and by using the extraordinary feature results from Kronecker sum operation, an explicit equation of control parameters defining the stability boundary in parametric space is obtained. Beyond stabilization, the method is used to shift all poles to a shifted half plane that guarantees a specified settling time of response. The stability regions of PID controllers are given in (kp, ki), (kp, kd) and (ki, kd) plane, respectively. The proposed method is also used to compute all the values of a PI controller stabilizing a control system with uncertain parameters. The proposed method is further extended to determine stability regions of uncertain coefficients of the system. Examples are given to show the benefits of the proposed method. Ó 2009 Published by Elsevier Ltd.

1. Introduction Despite continual advances in control theory and development of advanced control strategies, the proportional, integral, and derivative (PID) control algorithm still finds wide applications in industrial process control systems. It has been reported [1] that 98% in of the control loops in the pulp and paper industries are controlled by proportional integral (PI) controllers. Moreover, more than 95% of the controllers used in process control applications are of the PID type [2]. The popularity among industrial practitioners stems from the facts that the PID control structure is simple and its principle is easy to understand and that the PID controllers are deemed to be satisfactory and robust for a vast majority of processes. The primary problem associated with the use of PID controllers is tuning, that is, the determination of PID controller parameters for satisfactory control performance. Since the primary requirement of the candidate PID controller parameters is that of making the closed-loop system stable, it is often desired to construct the complete sets of stabilizing PID parameters. For instance, as fuzzy PID controllers have been used more and more widely in various control systems [3–5], it is necessary to known the complete sets of stabilizing PID parameters. And as the PID optimization technique develops [6–9], with the complete sets of stabilizing PID controller parameters being available for a given process, it can avoid the time consuming stability check for * Corresponding author. E-mail address: [email protected] (D. Zheng). 0196-8904/$ - see front matter Ó 2009 Published by Elsevier Ltd. doi:10.1016/j.enconman.2009.03.007

each set of PID controller parameters in the search process and thereby to save the controller tuning time. Up to now, there has been a great amount of research work on the determination of stabilizing sets of PI (proportional integral) and PID (proportional integral derivative) controllers [10–12]. A complete analytical solution based on the generalized version of the Hermite–Biehler theorem has been proposed [10] for computation of all stabilizing constant gain controllers for a given plant. In Ref. [11], the computation of all stabilizing PI and PID controllers for a given plant by linear programming has been proposed. This approach, besides its numerical efficiency, has also revealed important structural properties of PI and PID controllers. It shows that for a fixed proportional gain, the set of stabilizing integral and derivative gains lie in a convex set. Such an approach can deal with systems that are open loop stable or unstable, minimum or non-minimum phase. However, the computation time for this approach increases in an exponential manner as the order of the system increases, which is a disadvantage of the method. An alternative approach for fast computation of stabilizing PI and PID controllers based on the use of Nyquist plot has been proposed in Refs. [12,13]. Some fast approaches based on the gridding of frequency have been given in Refs. [14,15]. A stability boundary locus approach for the design of PI and PID controllers has been proposed in Ref. [16]. A parameter space approach using the singular frequency concept has been given in Ref. [17] for design of robust PID controllers. Usually, controller design in classical control engineering is based on a plant with fixed parameters. However, in the real world most practical system models are not known exactly, meaning that

1822

J. Fang et al. / Energy Conversion and Management 50 (2009) 1821–1827

the plant contains uncertainties. Thus, sometimes it is desirable to know the permissible varying range of plant parameters once the controller parameters are fixed. That is, to determine the stability regions of plant parameters. Many approaches to determine the stabilizing sets of PI/PID parameters employ the specific form of PI/PID controllers or the characteristic equation of PI/PID controlled system, thus can not be used to determine the stability regions of plant uncertain parameters. In this paper, a new method is given for computation of stabilizing PI and PID controllers in the parameter plane. The novel approach makes use of the extraordinary feature of the Kronecker summation operation and we obtain the explicit equation that PI and PID controllers parameters corresponding to the boundary of stability region must satisfy. The proposed method is also used for computation of stabilizing PI and PID controllers for relative stabilization. Actually, the proposed method has a wider application. It is further extended to determine stability regions of uncertain parameters in coefficient space. The paper is organized as follows: the proposed method is presented in Section 2. In Section 3, the computation of PI controllers for relative stabilization is given. In Section 4, the proposed method is used to determine the stability region of PID controllers. The computation of PI controllers for interval plant stabilization is given in Section 5. In section 6, the proposed method is extended to the determination of stability regions of uncertain coefficients. Concluding remarks are given in Section 7.

The characteristic equation of the closed-loop system is

CEðsÞ ¼ sDðsÞ þ ðkp s þ ki ÞNðsÞ ¼ fn ðkp ; ki Þsn þ    þ f1 ðkp ; ki Þs þ f0 ðkp ; ki Þ

ð3Þ

Transform Eq. (3) to differential equation matrix and define

x_ 1 ¼ x2 x_ 2 ¼ x3 .. . f ðk ;k Þ f ðk ;k Þ f ðkp ;ki Þ x_ n ¼  f0 ðkp ;ki Þ x1  f1 ðkp ;ki Þ x2      n1 f ðk ;k Þ n

p

n

i

p

n

i

p

i

yields

X_ ¼ AX

ð4Þ T

where X_ ¼ ½x_ 1 ; . . . ; x_ n  ; X_ ¼ ½x1 ; . . . ; xn 

2 6 6 6 6 6 A¼6 6 6 6 6 4

T

3

0

1

0

0



0

0 .. . .. .

0

1

0



0 .. .

0 .. .

1 .. .

 .. .

0 .. .

0

0

0

 0 1 f ðkp ;ki Þ        n1 fn ðkp ;k Þ

f ðk ;k Þ

f ðk ;k Þ

f ðk ;k Þ

 f0n ðkpp ;ki Þ  f1n ðkpp ;ki Þ  f2n ðkpp ;ki Þ i

i

i

0

i

7 7 7 7 7 7 7 7 7 7 5 nn

The relation between Eqs. (3) and (4) is given by the following equation: 2. Stabilization using a PI controller

CEðsÞ ¼ fn ðkp ; ki Þ detðsI  AÞ ¼ 0

In early work, many methods proposed to compute the stabilizing sets of PI/PID controllers are based on determining the stability boundary of parameters, the essence of which is to find all the values of parameters which will render pure imaginary roots. Here we study an alternative procedure: Kronecker sum method. Kronecker summation of two matrices and its properties: In matrix algebra [18,19] the Kronecker sum of square matrices M1(n1  n1) and M2(n2  n2) is defined as M 1  M 2 ¼ M1  In2 þ In1  M 2 , where M 1 2 Rn1 n1 ; andM 2 2 Rn2 n2 . Here  denotes the Kronecker summation and  the Kronecker product operations. The most critical feature of the Kronecker summation of M1 and M2 is that this new square matrix M1  M 2 2 Rðn1 ;n2 Þðn1 ;n2 Þ has n1  n2 eigenvalues which are indeed pair-wise combinatoric summations of the n1 eigenvalues of M1 and n2 eigenvalues of M2 [18]. That is, the Kronecker sum operation, in fact, induces the ‘‘eigenvalue addition” character to the matrices. We take advantage of this feature to obtain the equation that all the values of (kp, ki) that render pure imaginary roots must satisfy. Consider the single input, single output (SISO) control system of Fig. 1 where

GðsÞ ¼

NðsÞ DðsÞ

ð1Þ

is the plant to be controlled and C(s) is a PI controller of the form

CðsÞ ¼ kp þ

ki kp s þ ki ¼ s s

ð2Þ

ð5Þ

It can be seen from Eqs. (3) and (5) that s is the root of Eq. (3) as well as the eigenvalues of matrix A. Due to the fact that A is a constant matrix, the complex conjugates of s also satisfy Eq. (5).

detðs I  AÞ ¼ 0

ð6Þ

Therefore, if s = jx is the root of Eq. (3), it must be the eigenvalue of matrix A. And at the same time, s = jx is also the root of Eq. (3) and the eigenvalue of matrix A. Since the sum of two eigenvalues s = jx and s = jx is zero then the Kronecker sum of two matrices must be singular when such hkp ; ki ; xi correspondence occurs. That is,

ACE ¼ det½A  A ¼ 0

ð7Þ

Note that Eq. (7) does not contain s and is the function of (kp, ki). For every value (kp, ki) that satisfy Eq. (7), Eq. (3) will have a pair of conjugate imaginary roots or roots at origin. As we know, the imaginary axis and the origin are the only places that the stability shift of the system will occur. Thus, Eq. (7) defines the boundary of the stability region in kp–ki plane. The stability boundary may divide the parametric plane into several separate regions. To determine the actual stability regions, choose one point in one separate region, which will result in a polynomial, determine the stability posture of the polynomial by calculating the roots of polynomial. The corresponding system is stable if the polynomial has no right half plane (RHP) root. Then by executing the D-subdivision method [20], the prospective region is a stable one. Repeating the same procedure and test other separated regions and all stability regions can be determined. 2.1. Example 1 Consider the control system of Fig. 1 with transfer function

GðsÞ ¼ Fig. 1. A SISO control system.

NðsÞ s3 þ 4s2  s þ 1 ¼ DðsÞ s5 þ 2s4 þ 32s3 þ 14s2  4s þ 50

ð8Þ

J. Fang et al. / Energy Conversion and Management 50 (2009) 1821–1827

The closed-loop characteristic equation is 6

5

4

3.1. Example 2 3

CEðsÞ ¼ s þ 2s þ ðkp þ 32Þs þ ð4kp þ ki þ 14Þs þ ð4ki  kp  4Þs2 þ ðkp  ki þ 50Þs þ ki

Consider the function of Manutec robot [21]

ð9Þ GðsÞ ¼

Refer to Eq. (4), transform Eq. (9) into the form of differential equation matrix yields

2 6 6 6 6 A¼6 6 6 6 4

0

0

0

0

0

0

1

0

0

0

0

0

1

0

0

0

0

0

1

0 7 7 7 0 7 7 0 7 7 7 1 5

0 0 0 0 ðkp  ki þ 50Þ ð4ki  kp  4Þ ð4kp þ ki þ 14Þ ðkp þ 32Þ 2

When q = 0, the closed-loop characteristic equation is

Refer to Eq. (7), execute Kronecker sum operation yields

n 4 3 2 2 ACE ¼ ki 10kp þ ð465  5ki Þkp þ ð35ki  1125ki  4386Þkp 3

2

CEðsÞ ¼ s5 þ 87s4 þ 1970s3 þ 3600s2 þ 10kp s þ 10ki

ð14Þ

Replace s with s + q in Eq. (14) yields

4

þð30ki  2125ki þ 43480ki  217964Þkp þ 5ki o 3 2 235ki  6626ki þ 42974ki  4308200

f qÞ ¼ s5 þ ð5q þ 87Þs4 þ ð10q2 þ 348q þ 1970Þs3 CEðs;

The explicit equation Eq. (10) defines the stability boundary in kp–ki plane. After generating the stability boundary defined by Eq. (10), we choose one point in one separated region, calculate the roots of the characteristic equation. If the characteristic equation has no right half plane (RHP) root, the system is stable. According to D-subdivision method, the prospective region is a stable one. Repeating

6 6 A¼6 6 4

ð13Þ

3

1

ki

2

10 sðs þ 2Þðs þ 40Þðs þ 45Þ

0

0

1823

þ ð10q3 þ 522q2 þ 5910q þ 3600Þs2 þ ð10kp þ 5q4 þ 348q3 þ 5190q2 þ 7200qÞs þ 10qkp þ 10ki þ q5 þ 87q4 þ 1970q3 þ 3600q2

ð15Þ

Refer to Eq. (4), transform Eq. (15) into the form of differential equation matrix yields

3 0 1 0 0 0 7 0 0 1 0 0 7 7 0 0 0 1 0 7 5 0 0 0 0 1 5 4 3 2 4 3 2 3 2 2 10qkp þ 10ki þ q þ 87q þ 1970q þ 3600q ð10kp þ 5q þ 348q þ 5190q þ 7200qÞ ð10q þ 522q þ 5910q þ 3600Þ ð10q þ 348q þ 1970Þ ð5q þ 87Þ

the same procedure and tested other separated regions and all stability regions can be determined. The stability boundary and actual stability region (shaded) of the example system is shown in Fig. 2. 3. Relative stabilization In this section, the approach is further developed for relative stabilization. In the analysis and design of control systems, it is important to shift all poles of the characteristic equation of a control system to a desired region in the complex plane, for example, to a shifted half plane that guarantees a specified settling time of the response. The aim of this section is to find all values of kp and ki that put all the closed-loop poles to the left of s = q(q = constant). Using s + q instead of s in Eqs. (1) and (2)

~ Nðs þ qÞ NðsÞ ~ GðsÞ ¼ Gðs þ qÞ ¼ ¼ ~ Dðs þ qÞ DðsÞ

ACE ¼ ð10qkp þ 10ki þ q5 þ 87q4 þ 1970q3 þ 3600q2 Þ n 2  25ð4q þ 87Þ2 kp þ ½ð200q þ 4350Þki þ 20ðq þ 1Þðq þ 20Þ   2q þ 45Þð12q3 þ 348q2  1970q  83859Þkp 16qðq þ 21Þð2q þ 470ð2q þ 85Þðq þ 1Þ2 ðq þ 20Þ2  o2  2q þ 45Þ2 ¼ 0

ð16Þ

The stability boundary in kp–ki plane is defined by Eq. (16). Plot the stability boundary corresponding to different q in kp–ki plane. The stability region for q = 0 and q = 0.2 is shown is Fig. 3. 4. Design of PID controllers

ð11Þ

and

kp s þ kp q þ ki ~ CðsÞ ¼ Cðs þ qÞ ¼ sþq

Refer to Eq. (7), execute Kronecker sum operation yields

ð12Þ

In this case, the characteristic equation of closed-loop system still have the form of Eq. (3). And the Kronecker sum operation can be used to obtain the equation defining the boundary of stability region.

Consider the hydro-turbine governing system [8] shown in Fig. 4. The transfer function of the hydraulic pressure servo system is

G1 ðsÞ ¼

1 ð1 þ T S sÞ

ð17Þ

The transfer function of the hydro-turbine system is

G2 ðsÞ ¼

ey  ðeqy eh  eqh ey ÞT w s 1 þ eqh T w s

ð18Þ

1824

J. Fang et al. / Energy Conversion and Management 50 (2009) 1821–1827

Fig. 2. Stability boundary and actual stability region (shaded).

Fig. 5. Stability region in kp–ki plane for fixed kd.

Fig. 3. Stability region in kp–ki plane.

Fig. 6. Stability region in kp–kd plane for fixed ki.

optimizing procedure, the lower and upper bounds of three PID parameters should be set beforehand. With larger searching ranges of three PID parameters, it is possible to find more optimal control parameters. Thus, it is necessary to know the permissible ranges of three PID parameters. System parameters are TS = 5.1, Ta = 3.9, Tw = 0.365, ey = 0.313, eh = 1.3028, eqy = 1.0045, eqh = 0.3843 and en = 0.5. The characteristic equation of the closed-loop system is

CEðsÞ ¼ 2:79s4 þ ð20:7947  0:2935kd Þs3 þ ð1:313kd þ 6:5201  0:2935kp Þs2 þ ð0:5 þ 1:313kp  0:2935ki Þs þ 1:313ki

Refer to Eq. (4), transform Eq. (21) into the form of differential equation matrix yields

Fig. 4. The system structure of the hydro-turbine governing system.

2 6 6 A¼6 4

0

1

0

0

0

0

1

0

0

0

i  1:313k 2:79



ð0:5þ1:313kp 0:2935ki Þ 2:79

0 

ð1:313kd þ6:52010:2935kp Þ 2:79

1

ki þ kd s s

7 7 7 5

Refer to Eq. (7), execute Kronecker sum operation yields

ð19Þ

The transfer function of the PID governing system is

GC ðsÞ ¼ kp þ

3

dÞ  ð20:79470:2935k 2:79

The transfer function of the generator and load is

1 G3 ðsÞ ¼ T a s þ en

ð21Þ

ð20Þ

The PID parameters are optimized using deterministic-chaoticmutation evolutionary programming (DCMEP). When executing the

n 2 ACE ¼ ki ð1:6202  0:0143kd Þkp  ð21:6447 þ 0:498ki

þ4:2176kd Þkp þ ð0:0032ki þ 0:0639kd Þkp kd þ ð0:0243kd o2 1:0835ki  1:604Þkd þ ð0:0304ki þ 76:6622Þki ð22Þ After obtaining Eq. (22) defining stability boundary in kp–ki–kd plane, we can also obtain the equation defining stability boundary

1825

J. Fang et al. / Energy Conversion and Management 50 (2009) 1821–1827

Fig. 8. Position control system.

N1 ðsÞ ¼ z0 þ z1 s þ z2 s2 þ z3 s3 þ     2 s2 þ p  3 s3 þ    D1 ðsÞ ¼ p0 þ p1 s þ p N2 ðsÞ ¼ z0 þ z1 s þ z2 s2 þ z3 s3 þ    1 s þ p 2 s2 þ p3 s3 þ    D2 ðsÞ ¼ p0 þ p and

ð24Þ 2

3

N3 ðsÞ ¼ z0 þ z1 s þ z2 s þ z3 s þ    0 þ p1 s þ p2 s2 þ p  3 s3 þ    D3 ðsÞ ¼ p N4 ðsÞ ¼ z0 þ z1 s þ z2 s2 þ z3 s3 þ    0 þ p 1 s þ p2 s2 þ p3 s3 þ    D4 ðsÞ ¼ p Taking all combinations of the Ni(s) and Dj(s) for i, j = 1, 2, 3, 4, yields the following 16 Kharitonov plants family

GK ðsÞ ¼ Gij ðsÞ ¼

Fig. 7. (a) Stability region (shaded) in ki–kd plane for fixed kp. (b) Stability region (shaded) in ki–kd plane for fixed kp.

in kp–ki plane, kp–kd plane and ki–kd plane for fixed kd, ki and kp, respectively. The stability region in kp–ki plane for fixed kd is shown in Fig. 5. The stability region in kp–kd plane for fixed ki is shown in Fig. 6. The stability region in ki–kd plane for fixed kp is shown in Fig. 7a and b. Fig. 7a and b shows that for fixed kp, the stability region in ki–kd plane is a convex polygon, which accord with the result proved in Ref. [11].

Ni ðsÞ Dj ðsÞ

ð25Þ

where i, j = 1, 2, 3, 4. Define the set S(C(s)G(s)) that contains all the values of the parameters of the controller C(s) that stabilize G(s), and then, the set of all the stabilizing values of the parameters of a PI controller that stabilize the interval plant of Eq. (23) can be written as:

SðCðsÞGðsÞÞ ¼ SðCðsÞGK ðsÞÞ ¼ SðCðsÞG11 ðsÞÞ \ SðCðsÞG12 ðsÞÞ    SðCðsÞG44 ðsÞÞ where GK(s) represents the sixteen Kharitonov plants family given in Eq. (24).

5. Interval plant stabilization

5.1. Example 3

There are some important results in the literature about stabilization of interval systems. For example, it was shown that a constant gain controller stabilizes an interval plant family if and only if it stabilizes a set of eight of the extreme plants that can be obtained from the upper and lower values of uncertain parameters [22]. And a first order controller stabilizes an interval plant if it stabilizes the set of extreme plants [23]. A first order controller stabilizes an interval plant if and only if it simultaneously stabilizes the sixteen Kharitonov plants family [24]. Someone used the generalized version of the Hermite–Biehler theorem has been used for the stabilization of interval systems [25]. In this section, to characterize all the parameters of a first order controller that stabilize an interval plant, the Kronecker summation method is used to find all the values of the parameters of a PI controller for which the given interval plant is Hurwitz stable. Consider a unity feedback system with a PI controller of Eq. (2) and an interval plant

Consider the position control system with the block diagram of Fig. 8, where the motor is assumed to have the transfer function [26]

GðsÞ ¼

NðsÞ zm sm þ zm1 sm1 þ    þ z0 ¼ DðsÞ pn sn þ pn1 sn1 þ    þ p0

ð23Þ

j ; j ¼ 0; 1; . . . ; n. Let where zi 2 ½zi ; zi ; i ¼ 0; 1; . . . ; m and pj 2 ½pj ; p the Kharitonov polynomials associated with N(s) and D(s) be, respectively:

GðsÞ ¼

Km Km ¼ sðJs þ bÞðLf s þ Rf Þ Lf Js3 þ ðbLf þ JRf Þs2 þ bRf s

ð26Þ

Rf and Lf are the resistance and inductance of the field, Km is the motor constant, J is the inertia, b is the viscous friction and the controller is a PI controller of the form Eq. (2). To eliminate the steady state error of a ramp input or a steady state error to a torque disturbance, an integrator is needed here. The characteristic equation of the closed-loop system is 2

CEðsÞ ¼ Lf Js4 þ ðLf b þ Rf JÞs3 þ Rf bs þ K m kp s þ K m ki

ð27Þ

Refer to Eq. (4), transform Eq. (27) into the form of differential equation matrix yields

2 6 6 A¼6 6 4

0

1

0

0

0

0

1

0

0

0

 KLm Jki f

K k  Lm J p f

3

0 1 Rf b ðLf bþRf JÞ LJ  LJ f

f

7 7 7 7 5

1826

J. Fang et al. / Energy Conversion and Management 50 (2009) 1821–1827

Refer to Eq. (7), execute Kronecker sum operation yields 2

2

2

2

ACE ¼ 16K 3m ki ½K m Lf Jkp  ðLf Rf b þ JbRf Þkp þ ðb L2f þ 2bLf Rf J þ R2f J 2 Þki =L7f J 7

ð28Þ

The major uncertainty is assumed to be in the parameters Km, b, Lf and J. In nominal case, Km = 60  103, Rf = 1.2, b = 2.5  103, Lf = 1.3  102 and J = 2  103. Assume now that the parameters Km, b, Lf and J may vary by 10%, 15%, 20% and 40% around their nominal values, respectively, then

 m  ¼ ½54  103 ; 66  103 ; K m 2 ½K m ; K  ¼ ½2:1  103 ; 2:9  103 ; b 2 ½b; b

ð29Þ

Lf 2 ½Lf ; Lf  ¼ ½1:04  102 ; 1:56  102 ; J 2 ½J; J ¼ ½1:2  103 ; 2:8  103 

Fig. 11. Stability region in TS–Tw plane for the example system.

Therefore, by over-bounding the coefficients of Eq. (27), since the uncertain parameters of Eq. (27) are multi-linearly dependent on the above parameters, the transfer function of the motor can be written in the following form:

GðsÞ ¼

z0 p3 s3 þ p2 s2 þ p1 s þ p0

ð30Þ

 ¼  m  ¼ ½54  103 ; 66  103 , p ¼ 0; p 2 ½bRf ; bR where z0 2 ½K m ; K f 0 1 3 3    ½2:52  10 ; 3:48  10 , p 2 ½bLf þ JRf ; bLf þ JRf  ¼ ½1:4618  103 ; 2

3:4052  103  and p3 2 ½Lf J; Lf J ¼ ½1:25  105 ; 4:37  105  And in nominal case, z0 = 60  103, p0 = 0, p1 = 3  103, p2 = 2.4325  103, p3 = 2.6  105. The stabilizing parameters can be calculated for all eight Kharitonov plants (i = 1–2 and j = 1–4). Now the characteristic equation of the closed-loop system can be rewritten as

CEðsÞ ¼ p3 s4 þ p2 s3 þ p1 s2 þ ðp0 þ z0 kp Þs þ z0 ki

Refer to Eq. (4), transform Eq. (31) into the form of differential equation matrix yields

2

4.5

G24

4 3.5

0

6 6 0 A¼6 6 0 4  zp03ki

G12

G14

G22

3

1

0

0

1

0 

0

ðp0 þz0 kp Þ p3

 pp13

0

3

7 0 7 7 1 7 5

 pp23

Refer to Eq. (7), execute Kronecker sum operation yields

ki

2.5 2

2

Gno

ACE ¼

1.5

G11

1

16z0 ki ½p3 z20 kp þ ð2z0 p0 p3  z0 p1 p2 Þkp þ z0 p22 ki þ p20 p3  p0 p1 p2  p73

¼0

0.5

ð32Þ

G21 0

ð31Þ

0

G23

2

G13

4

6

8

10

12

14

16

18

kp

Fig. 9. The stability region in nominal case and of eight Kharitonov plant.

The stability boundary in nominal case and in all eight Kharitonov plants are shown in Fig. 9. And the stability region (shaded) that guarantees stability under parameter varying is shown in Fig. 10.

6. Determination of stability region of uncertain coefficients In the above sections, the proposed method is used to determine stabilizing sets of PID controller parameters. Actually, the proposed method has a wider application. It can be used to determine the stability region of other uncertain coefficients. Consider the example case studied in Section 4. Suppose system parameters are Ta = 3.9, ey = 1.313, eh = 1.3028, eqy = 1.0045, eqh = 0.3843, en = 0.5. And the PID parameters are chosen as kp = 40, ki = 10, kd = 30. Now determine the stability region in TS–Tw plane. The characteristic equation of the closed-loop system is

CEðsÞ ¼ 1:4988T S T w s4 þ ð0:1921T S T w þ 3:9T S  22:6235T w Þs3 þ ð0:5T S  31:9709T w þ 43:29Þs2 þ ð53:02  8:0408T w Þs þ 13:13 Fig. 10. The stability region (shaded) that guarantee stability under parameter varying case.

ð33Þ

Refer to Eq. (4), transform Eq. (33) into the form of differential equation matrix yields

J. Fang et al. / Energy Conversion and Management 50 (2009) 1821–1827

2 6 6 A¼6 6 4

0

1

0

0

0

0

1

0

0 13:13  1:4988T STw

0 w  53:028:0408T 1:4988T S T w

3

0 1 0:1921T S T w þ3:9T S 22:6235T w w þ43:29  0:5T S 31:9709T  1:4988T S T w 1:4988T S T w

From the above matrix A, it is obvious that TSTw – 0. Thus, there is two singular boundaries in TS–Tw plane: TS = 0, Tw = 0. Refer to Eq. (7), execute Kronecker sum operation yields

n ACE ¼ ð167:0557T S þ 20452ÞT 3w þ ð4:4214T 2S  7360:2T S

138920ÞT 2w þ ð106:4272T 2S þ 35248T S þ 182600ÞT w o þ338:7105T 2S  31478T S =ðT 7S T 7w Þ

ð34Þ

The stability boundaries and the actual stability region are shown in Fig. 11. 7. Conclusions In this paper, a new approach has been proposed for computation of stability boundary of PI and PID controller parameters that guarantee stability. The approach is based on the ‘‘eigenvalue addition” feature result from Kronecker sum operation by which the explicit equation that defines the stability boundaries in parametric space can easily be obtained. Thus, the proposed method is not confined to the determination of stabilizing sets of PID controller parameters. It can be used to determine the stability regions of uncertain parameters in coefficient space. The examples given clearly show the value of the method presented. Acknowledgement This work was supported by the National Natural Science Foundation of China (Grant Nos. 60674088, 60874113), Ph.D. Programs Foundation of Ministry of Education of China (Grant No. 200802550007) and Innovation Program of Shanghai Municipal Education Commission (Grant No. 09ZZ66). References [1] Bialkowski WL. The control handbook, In: Levine WS, editor; 1996. p. 1219–42. [2] Astrom KJ, Hagglund T. PID controllers: theory, design, and tuning. 2nd ed. NC, USA: Research Triangle Park; 1995. [3] El-Sherbiny MK, El-Saady G, Yousef AM. Efficient fuzzy logic load–frequency controller. Energy Convers Manage 2002;43:1853–63. [4] Yesil E, Guzelkaya M, Eksin I. Self tuning fuzzy PID type load and frequency controller. Energy Convers Manage 2004;45:377–90.

1827

7 7 7 7 5

[5] Guesmi K et al. Systematic design approach of fuzzy PID stabilizer for DC–DC converters. Energy Convers Manage 2008. doi:10.1016/j.enconman.2008. 03.01. [6] Ngamroo I. An optimization technique of robust load–frequency stabilizer for superconducting magnetic energy storage. Energy Convers Manage 2005;46(18–19):3060–90. [7] Pothiya S, Ngamroo I. Optimal fuzzy logic-based PID controller for load– frequency control including superconducting magnetic energy storage units. Energy Convers Manage 2008. doi:10.1016/j.enconman.2008.03.01. [8] Jiang Chuanwen, Ma Yuchao, Wang Chengmin. PID controller parameters optimization of hydro-turbine governing systems using deterministic-chaoticmutation evolutionary programming (DCMEP). Energy Convers Manage 2006;47:1222–30. [9] Zhuang M, Atherton DP. Automatic tuning of optimum PID controllers. IEE Proc Part D 1993;140:216–24. [10] Ho MT, Datta A, Bhattacharyya SP. A new approach to feedback stabilization. In: Proceedings of the 35th CDC; 1996. p. 4643–8. [11] Ho MT, Datta A, Bhattacharyya SP. A linear programming characterization of all stabilizing PID controllers. In: Proceedings of American control conference; 1997. [12] Munro N, Söylemez MT. Fast calculation of stabilizing PID controllers for uncertain parameter systems. In: Proc symposium on robust control, Prague; 2000. [13] Soylemez MT, Munro N, Baki H. Fast calculation of stabilizing PID controllers. Automatica 2003;39:121–6. [14] Shafiei Z, Shenton AT. Frequency domain design of PID controllers for stable and unstable systems with time delay. Automatica 1997;33:2223–332. [15] Huang YJ, Wang YJ. Robust PID tuning strategy for uncertain plants based on the Kharitonov theorem. ISA Trans 2000;39:419–31. [16] Tan N, Kaya I, Atherton DP. Computation of stabilizing PI and PID controllers. In: Proceedings of the 2003 IEEE international conference on the control application (CCA2003); 2003. p. 876–81. [17] Ackermann J, Kaesbauer D. Design of robust PID controllers. In: European control conference; 2001. p. 522–7. [18] Bernstein DS. Matrix mathematics. Princeton University Press; 2005. p. 7. [19] Brewer JW. Kronecker products and matrix calculus in system theory. IEEE Trans Circ Syst 1978;CAS-25:772–81. [20] El’sgol’ts LE, Norkin SB. Introduction to the theory and application of differential equations with deviating arguments. New York: Academic Press; 1973. [21] Dorf RC, Bishop RH. Modern control systems. New Jersey: Prentice Hall; 2001. [22] Ghosh BK. Some new results on the simultaneous stabilization of a family of single input, single output systems. Syst Control Lett 1985;6:39–45. [23] Hollot CV, Yang F. Robust stabilization of interval plants using lead or lag compensators. Syst Control Lett 1990;14:9–12. [24] Barmish BR, Holot CV, Kraus FJ, Tempo R. Extreme points results for robust stabilization of interval plants with first order compensators. IEEE Trans Automat Control 1993;38:1734–5. [25] Ho MT, Datta A, Bhattacharyya SP. Design of P, PI and PID controllers for interval plants. In: Proceedings of American control conference, Philadelphia; June 1998. [26] Tan N, Kaya I, Atherton DP, et al. Computation of stabilizing PI and PID controllers using the stability boundary locus. Energy Convers Manage 2006;47:3045–58.