The Credible Regions on the Approximate Stability Boundaries of ...
1486
IEEE TRANSACTIONS ON AUTOMATIC CONTROL, VOL. 52, NO. 8, AUGUST 2007
The Credible Regions on the Approximate Stability Boundaries of Nonlinear Dynamic Systems Yong Min, Lei Chen, Kaiyuan Hou, and Yonghua Song
Abstract—The stability boundary of a nonlinear dynamic system may often be constructed from stable manifolds of type-1 equilibrium points. These stable manifolds are typically approximated numerically, with only local validity. A numerical method is proposed, based on 1- and 2-D characteristic invariant manifolds (CIMs) which lie on the stable manifold and can be computed with greater accuracy than approximations to the entire manifold because of their lower dimension, to estimate Credible Regions in which the mismatch between the approximate and actual boundaries is below a threshold. Index Terms—Credible region, invariant manifold, stability boundary, stability region, stable manifold.
I. INTRODUCTION It has been proved in [1]–[3] that under certain assumptions the stability boundary of a nonlinear autonomous dynamic system equals to the union of the stable manifolds of the unstable equilibrium points (UEPs) lying on the boundary. A fairly large class of systems satisfies these assumptions, so their stability boundaries can be produced by the approximation of stable manifolds of type-1 UEPs.1 The problem of estimating the stable manifold of a hyperbolic equilibrium has been extensively explored during the past thirty or so years and is still under investigation. Most of the approaches developed so far are concentrated on constructing low-order approximations [4]–[13]. Unfortunately, the approximations obtained by these studies are all only valid locally in a neighborhood of the corresponding UEP, and no method has yet been suggested to determine the size of the neighborhood or to measure the approximate errors. This constitutes the main limitation of the stable manifold approximation type approaches because their results have to be confirmed by other type approaches. A recent work [14] proposes a new theoretical approach to the problem of finding invariant manifolds for nonlinear dynamic systems. It states that the local invariant manifolds passing through an equilibrium point can be determined via an associated system of singular first-order quasi-linear PDEs. Unlike previous studies, such as [4] and [12], which deal only with the stable or unstable manifolds of an equilibrium point, the main results of [14] is applicable to not only stable or unstable manifolds but also other invariant manifolds passing through an equilibrium point. This note proposes a numerical approach to determine the credible region around a type-1 UEP on the approximate stability boundary for the given threshold. The credible region is defined as the maximum neighborhood of the type-1 UEP on the stability boundary in which Manuscript received February 15, 2006; revised March 6, 2007 and March 18, 2007. Recommended by Associate Editor J. Berg. This work was supported by the Special Fund of the National Priority Basic Research of China (2004CB217904) and the National Science Foundation of China (50323002). Y. Min and L. Chen are with State Key Lab of Power Systems, Department of Electrical Engineering, Tsinghua University, Beijing 100084, China (e-mail: [email protected]). K. Hou was with the Department of Electrical Engineering, Tsinghua University, Beijing 100084, China. He is now with the Northeast China Grid Company Limited, Shenyang 110006, Liaoning Province, China. Y. Song is with the University of Liverpool, Liverpool L69 3GJ, U.K. Digital Object Identifier 10.1109/TAC.2007.902759 1An equilibrium point is called type-1 UEP if the Jacobian matrix of the system at this point has only one positive real part eigenvalue.
the approximate precision is under the given threshold. The spirit of the work is to find some 1- or 2-D characteristic invariant manifolds (CIMs), which are exactly located on the stability boundary. Based on the theory developed in [14], the Taylor series solution of relevant PDEs is adopted to approximate both the stable manifold and the CIMs passing through a type-1 UEP. Because the CIMs are only 1- or 2-D, it is possible to calculate them with a high precision by higher-order Taylor series expansion. The calculated results of CIMs are used as benchmarks to evaluate the mismatches between the approximate stability boundary and the true boundary. This technique gives a way to measure the approximate errors and determine the credible region. This note proceeds as follows. Section II introduces the fundamental theory of stability regions as well as the main results of [14]. Section III derives an associated system of PDEs to find the stable manifold passing through a type-1 UEP of the original system. The Taylor series solution of PDEs is provided to approximate the local stability boundary. In Section IV, the CIMs passing through equilibrium points are defined. Then the corresponding systems of PDEs are derived and the Taylor series solutions of PDEs are provided similarly. In Section V, the credible region around a type-1 UEP is defined and a procedure to estimate the credible region for a given approximate stability boundary is introduced. An example is given in Section VI. II. PRELIMINARY A. Stability Region and Stability Boundary Consider the following autonomous dynamic system: x _
!
=
f (x);
x
2
R
n
(1)
n R is assumed to be smooth. Let x0 be a hyperbolic where f : Rn equilibrium point of (1). Its stable and unstable manifolds are denoted as W s (x0 ) and W u (x0 ) respectively [1]. If xs is a Stable Equilibrium Point (SEP), A(xs ) = W s (xs ) is defined as the stability region surrounding xs and its boundary is denoted as @ A(xs ). It has been proved in [1], [2] that for systems satisfying the following conditions: 1) all the equilibrium points on @ A(xs ) are hyperbolic; 2) the stable and unstable manifolds of equilibrium points on @ A(xs ) satisfy the transversality condition; 3) every trajectory on @ A(xs ) converges to an equilibrium point; the stability boundary @ A(xs ) is the union of the stable manifolds @ A (xs ); i = of the UEPs on the boundary. Furthermore, let xi 1; 2; . . . ; be all the type-1 UEPs on the boundary, then @ A(xs ) = s s W (xi ), where A(xs ) and W (xi ) denote the closures of A(xs ) and W s (xi ) respectively. This topological property makes it possible to estimate the boundary of the stability region of a SEP by approximating the stable manifolds of all type-1 UEPs on the boundary.
2
B. Finding Invariant Manifolds Through Singular PDES This section introduces a new approach developed in [14] to find invariant manifolds for nonlinear real analytic dynamic systems. First, we rewrite system (1) as follows: y _
=
G(y; z );
2
z _
2
=
H (y; z )
(2)
m p U ;z U and U m ; U p being open where x = (y; z ), with y m p subsets of R and R , respectively, (m + p = n). Without loss of generality, we suppose that the origin (0,0) is an equilibrium point, which means G(0; 0) = 0, H (0; 0) = 0.
0018-9286/$25.00 © 2007 IEEE
IEEE TRANSACTIONS ON AUTOMATIC CONTROL, VOL. 52, NO. 8, AUGUST 2007
If (@G)=(@z )j(0;0) = 0, system (2) can be rewritten as follows:
y_ = Ay + g(y; z );
z_ = Bz + Cy + h(y; z )
1487
J 01 J1 , with J0 being the element at J2 J0 the bottom right corner of J (x ) such that J1 2 R( 01)21 ; J2 2 R12( 01) , and J 01 2 R( 01)2( 01) . Then taking notice of J (x )q = q , (7) is rewritten as follows: Rewrite
(3)
where g (y; z ) and h(y; z ) are analytic functions whose Taylor series expansions begin with terms of order higher than one. Then we consider the following first-order PDEs:
J (x
n
u)
n
as
n
u
(0) = 0
(4) where : Rm ! Rp is the unknown vector field of (4). The above system of PDEs is singular at the origin (0,0) because the principal part G(y; z ) = Ay + g (y; z ) vanishes at the origin. Only under some specific conditions, the existence and uniqueness of the solution can be ensured. So the following theorem is given in [14]. Theorem 1([14]): Suppose the origin (0,0) is an equilibrium point of system (2), let k1 ; k2 ; . . . ; km be the eigenvalues of the matrix A = (@G)=(@y )j(0;0) and 1 ; 2 ; . . . ; p be the eigenvalues of the matrix B = (@H )=(@z )j(0;0) . Assume the following: 1) (@G)=(@z )j(0;0) = 0; 2) The real parts of all the kj are different from zero and have the same sign; 3) kj and i are not related by any equation of the form m
mk j
j
=
;
i = 1; . . . ; p
i
(5)
j =1
where all the mj are nonnegative integers satisfying
0.
m j =1
m > j
Then there exists a neighborhood V Rm of y0 = 0 and a unique and locally analytic mapping : V ! Rp such that
S = f(y; z ) 2 V
2 R : z = (y); (0) = 0g p
(6)
is an analytic local invariant manifold of system (2) that passes through the origin, with z = (y ) being the solution of (4). The proof and discussions of Theorem 1 can be found in [14]. III. STABLE MANIFOLD OF TYPE-1 UEP Suppose that system (1) satisfies the conditions given in Section II.A and xu is a type-1 UEP of the system. Let J (xu ) be the Jacobian matrix at xu and u be its only eigenvalue with positive real part. The eigenvector of J (xu ) corresponding to u is denoted as2
q = [q1 ; . . . ; q 01 ; q n
Construct the matrix Q =
x
u)
I 01 n
0
n]
T
= [q^
q^ q
n
T
;q
n]
T
y_ = (J 01 0 q^J2 =q )y + g(y; z ) z_ = J2 y=q + z + h(y; z ) n
n
Q01 (x 0
, system (1) is converted to the following form:
x^_ = Q01 J (x )Qx^ + R(^x) u
(7)
x) denotes terms of order higher than one. where R(^ 2We assume q 6= 0 in the following deduction. In case where q = 0, we could find a nonzero element in q , say q . Then if we exchange x = f (x) with x = f (x) in system (1), q and q in q would be switched.
(8)
u
x1 ; x^2 ; . . . ; x^n01 ]T ; z = x^n ; g(y; z ) and h(y; z ) denote where y = [^ corresponding terms of order higher than one. Because u is the only eigenvalue with positive real part, according to Theorem 1, the unique and locally analytic solution z = (y ) of the following PDEs: @ ((J 0 q^J =q )y + g(y; z )) = J y=q + + h(y; z ); 01 2 2 @y (0) = 0 (9) n
n
n
u
coincides with the stable manifold in the neighborhood of xu . Furthermore, if xu 2 @A(xs ), the solution z = (y ) will coincide with @A(xs ) in the neighborhood of xu . Generally it is not possible to solve for an analytic expression of z = (y) from the PDEs (9). We use the following multivariate series solution to approximate it:
z=
1
A (y) = Ey + y F y + 1 1 1 i
i=1
T
(10)
where Ai (y ) is the Taylor term of ith order and is an ith order homogeneous polynomial. Theoretically, substituting (10) into (9) and equating the same order Taylor coefficients at both sides, we can solve for all Taylor coefficients by a recursion method [14]. However, if the system’s order is higher than 3, the calculation is very difficult and only the second-order Taylor series expansion can usually be available. IV. CIMS PASSING THROUGH EQUILIBRIUM POINTS A. Definitions of CIMs Let x0 be a hyperbolic equilibrium point of system (1) and J (x0 ) be the Jacobian matrix at x0 . We denote one of the real eigenvalues of J (x0 ) by and a pair of complex conjugate eigenvalues of J (x0 ) by 6 j! . The corresponding eigenvectors are denoted by and 6 j respectively. Consider the linear part of system (1) at x0
x_ = J (x0 )(x 0 x0 );
:
and define x ^ =
n
u
n
@ (Ay + g(y; (y))) = B(y) + Cy + h(y; (y)); @y
n
u n
x2R : n
(11)
We know from the theory of linear systems that in state space passing through x0 there exists a pair of linear trajectories in the direction of , and a set of trajectories forming a hyperplane which is spanned by and . Because the line or hyperplane formed by the trajectories is actually an invariant manifold and is corresponding to the eigenvector(s) of the system, we call it CIMs. The definition is as follows. Definition 1: For system (11), a pair of linear trajectories in the direction of is defined as the 1-D CIM of x0 relative to the eigenvalue , denoted by W1 (x0 ; ) or W1 (x0 ), and a hyperplane spanned by the vectors and is defined as the 2-D CIM of x0 relative to the eigenvalues 6 j! , denoted by W2 (x0 ; 6 j! ) or W2 (x0 ).
1488
IEEE TRANSACTIONS ON AUTOMATIC CONTROL, VOL. 52, NO. 8, AUGUST 2007
According to the Hartman-Grobman’s theorem, system (1) is topologically equivalent to its linear part (11) in the neighborhood of x0 . Therefore, W1 (x0 ) or W2 (x0 ) is uniquely homologized to a 1- or a 2-D manifold of system (1) respectively. Thus CIMs of the nonlinear system (1) are defined as follows. Definition 2: The 1-D manifold of system (1), homologizing to the CIM W1 (x0 ) of its linear part (11), is defined as the 1-D CIM of (1), still denoted by W1 (x0 ). The 2-D manifolds of system (1), homologizing to the CIM W2 (x0 ) of its linear part, is defined as the 2-D CIM of (1), still denoted by W2 (x0 ). According to the above definitions, it is obvious that if x0 is a type-1 UEP on @A(xs ), then W1 (x0 ; ) @A(xs ) if < 0 and W2 (x0 ; 6 j! ) @A(xs ) if < 0. B. CIMs in the Form of PDEs The CIMs of a linear system can be easily found from the eigenvectors of the Jacobian matrix. For a nonlinear system, the CIMs can be found via Theorem 1. For the 1-D CIM W1 (x0 ), suppose = [1 ; 2 ; . . . ; n ]T is the eigenvector of J T (x0 ) corresponding to (in other words, T is the left eigenvector of J (x0 )). Assume 1 6= 0. Let ^ = [2 ; . . . ; n ] and construct the ma1 ^ J1 Jw . Rewrite J (x0 ) as ; with trix T = Jv Jn01 0 In01 J1 being the element at the top left corner of J (x0 ) such that Jw 2 R12(n01) ; Jv 2 R(n01)21 , and Jn01 2 R(n01)2(n01) . Then taking notice of J T (x0 ) = , the original system (1) is converted to the following form:
v_ = v + gv (v; w) w_ = Jv v=1 + (Jn01 0 Jv ^=1 )w + hw (v; w)
(12)
where x ^ = T (x 0 x 0 ); v = x ^1 ; w = [^ x2 ; x^3 ; . . . ; x^n ]T ; gv (v; w) and hw (v; w) denote terms of order higher than one. The first two conditions of Theorem 1 are automatically satisfied. The third condition is generally true in the actual system. So a unique analytic local invariant manifold
S1
f
= (v; w )
2 V 2 Rn0
1
:
w = '(v ); '(0) = 0g
(13)
exists in a neighborhood V 2 R of v = 0, with w = '(v ) being the unique and locally analytic solution of the system of the following PDEs:
@' ^=1 )' + hw (v; w ); (v + gv (v; w )) = Jv v=1 + (Jn01 0 Jv @v (14) '(0) = 0: Theorem 2: The unique analytic local invariant manifold defined in (13) coincides with the 1-D CIM W1 (x0 ; ) in the neighborhood of x0 in (1). Proof: See Appendix. For the 2-D CIM W2 (x0 ; 6 j! ), suppose 6 j = T T are the eigenvectors [ 1 ; 2 ; . . . ; n ] 6 j : [ 1 ; 2 ; . . . ; n ] 0! ; = of J T (x0 ) corresponding to 6 j! . Let A0 = ! 1 2 3 1 1 1 n ; and ^ = . 1 2 3 1 1 1 n
From J T (x0 )( 6 j ) = ( 6 j! )( 6 j ) we have [ ; ^ ]J ( x 0 ) =
A0 [ ; ^]
(15)
^ Js J2 ], with J2 ]. Rewrite J (x0 ) as [ In02 Jr Jn02 being the 2 2 2 submatrix at the top left corner of J (x0 ) such that Js 2 R22(n02) ; Jr 2 R(n02)22 ; Jn02 2 R(n02)2(n02) . Then taking notice of (15), the original system (1) can be converted to the following form: Define P = [
0
r_ = A0 r + gr (r; s) s_ = Jr 01 r + (Jn02 0 Jr 01 ^)s + hs (r; s)
(16)
where x ~ = P (x 0 x0 ); r = [~ x1 ; x~2 ]T ; s = [~ x3 ; x~4 ; . . . ; x~n ]T ; gr (r; s) and hs (r; s) denote terms of order higher than one. So there exists a neighborhood V 2 R2 of r = 0 and a unique analytic local invariant manifold
S2
f
= (r; s)
2 V 2 Rn0
2
:
s=
(r );
(0) = 0
g
(17)
with s = (r) being the unique and locally analytic solution of the system of the following PDEs:
@ 01 01 ^) (A0 r + gr (r; s)) = Jr r + (Jn02 0 Jr @r (0) = 0:
+ hs (r; s);
(18)
Theorem 3: The unique analytic local invariant manifold defined in (17) coincides with the 2-D CIM W2 (x0 ; 6 j! ) in the neighborhood of x0 in system (1). Proof: See the Appendix. C. Series Solutions of CIMs The Taylor series of the unique and locally analytic solutions w = '(v ) of the system of (14) and s = (r) of the system of (18) are denoted as follows:
w= s=
1 j =1
aj v j
1
bjh r1j r2h :
(19) (20)
Substituting (19) and (20) into (14) and (18), respectively, we can solve for aj and bjh . Because there are only one variable v and two variables r1 ; r2 in (19) and (20), respectively, it is feasible to solve for the highorder Taylor coefficients of CIMs with the aid of a symbolic software package. These high-order Taylor series approximations of CIMs can be then used as benchmarks to test the accuracy of the approximate stability boundary. V. CREDIBLE REGION ON THE APPROXIMATE STABILITY BOUNDARY AROUND TYPE-1 UEP A. Credible Region and Its Definition The credible region on the approximate stability boundary around a type-1 UEP is a region in which the maximum mismatch between the actual boundary and the approximate boundary is less than a given tolerance. The definition is as follows.
IEEE TRANSACTIONS ON AUTOMATIC CONTROL, VOL. 52, NO. 8, AUGUST 2007
Definition 3: Suppose xu is a type-1 UEP on @A(xs ). Let 0 de^ denote the note the true stability boundary passing through xu and 0 approximate one. For any given positive number , the credible region ^ around xu with a credible level is defined as on 0
^ : jx 0 0 j g D ( ) = fx 2 0
(21)
where jx 0 0j denotes the shortest Euclidean distance from a point 2 0^ to the manifold 0. To make the above definition more practical, we consider a neighborhood B (") of xu . Without loss of generality, let B (") be a spherical neighborhood
x
B ( " ) = f x 2 R n : jx 0 x u j " g
as D2i (y; z ). D1i (y; z ) and D2i (y; z ) are then used as the indices to ^ and 0. estimate the approximation error between 0 For a given error threshold , define
X1 = f(y; z) 2 Rn : D1i (y; z) = ; i = 1; . . . ; kg; X2 = f(y; z) 2 Rn : D2i (y; z) = ; i = 1; . . . ; lg; D = fj(y; z) 0 xu j 2 R : (y; z) 2 X1 [ X2 g: The sets X1 and X2 contain all the points on the CIMs which have the ^ . The set D contains all the Euclidean distances given deviation from 0 from the points in X1 [ X2 to xu . Therefore, the minimum element in set D can be used as the maximum available radius "max of the credible region B (")
(22)
where " > 0 is the radius of B ("). ^ = 0^ \ B ("), then for a given > 0; ^ Define = 0 \ B (") and
is a credible region with credible level if there is a positive number " ^ , where @ ^ denotes the boundary such that jx 0 j for any x 2 @
^. of
Due to the difficulties in calculating the accurate 0 or , it is ob^ can not be exactly determined for a given tolerance . In vious that
the following part of this section, we will use the high-order series ap^. proximation of CIMs to estimate
B. Estimating the Credible Region via CIMS The algorithm of estimating the credible region via the concept of CIMs is based on the following two assumptions. 1) The approximation of CIMs is more accurate than that of the stability boundary. ^ appears only on one of 2) The maximum mismatch between 0 and 0 the CIMs. The first one is generally true but the validity of the second one is not easy to confirm. However, even if the second one is not true, the result obtained by the following method is still valuable for estimating ^. the mismatches between 0 and 0 Suppose that xs is a SEP in (1) and let J (xu ) be the Jacobian matrix at a type-1 UEP xu on @A(xs ). From Section IV we know there is an 1-D CIM W1i (xu ) = W1 (xu ; i ) for any negative real eigenvalue i of J (xu ), and a 2-D CIM W2i (xu ) = W2 (xu ; i 6 j!i ) for any pair of complex eigenvalues i 6j!i with i < 0. According to the definitions of CIMs, we know W1i (xu ) 0; i = 1; . . . ; k and W2i (xu ) 0; i = 1; . . . ; l, where k denotes the number of negative real eigenvalues and l denotes the number of pairs of complex conjugate eigenvalues with negative real parts. Therefore, W1i (xu ) and W2i (xu ) can be used to estimate the ap^ and 0 with the expansion of B ("). Beproximation error between 0 cause it is difficult to find actual W1i (xu ) and W2i (xu )), their high ^ 1i (xu ) and W ^ 2i (xu ) are used precision Taylor series approximation W instead. Based on the transformations defined in Section III and IV, a point ^ 1i (xu ) is converted to the state space of (8) as follows: (v; w) 2 W
y v = Q01 T 01 : z w
1489
"max = minfd 2 Dg:
(24)
Remark 1: The calculation of the coefficients of the terms with order higher than 7 in (19) and (20) may be involved (especially for bjh ) and, on the other hand, is usually not helpful to improve the precision of the approximation of CIMs because of the limitation of the region of convergence of the solutions of (14) and (18). Remark 2: The proposed algorithm of credible region can be applied to the high-order systems and can work with any other stability boundary approximation approaches as long as D1i (y; z ) and D2i (y; z ) can be calculated. VI. EXAMPLE In this section we give an example to describe the steps to determine the credible region. The approximate stability boundary in the example is calculated by the method developed in Section III. However, the results of other quadratic form approximations can also be used instead. Example 1: The state equations of a 2-machine power system are as follows:
!_ 1 = P1M + P1E sin( + 12 ) 0 !0 D1 !1 =TJ 1 !_ 2 = P2M 0 P2E sin( 0 12 ) 0 !0 D2 !2 =TJ 2 _ = !1 0 !2
(25)
where P1E = 0115:50; P2E = 0144:38; !0 = 100; P1M = 52:17; P2M = 10:94; 12 = 0:27; TJ 1 = 5; TJ 2 = 4; D1 = 0:04; D2 = 0:024. The SEP of the system is xs (0; 0; 0:196). There are two type-1 UEPs on the stability boundary of xs : xu1 (30:569; 30:569; 3:085); xu2 (30:569; 30:569; 03:199).
Now we take the following steps to calculate the credible region on the approximate stability boundary passing through xu1 . 1) Calculate the Jacobian matrix J (xu1 ) and its eigenvalues. It has two negative real eigenvalues: 1 02:228; 2 016:920, and one positive real eigenvalue: u 14:750. ^ passing through xu1 by a 2) Approximate the stability boundary 0 second-order Taylor series expansion
(23)
^ 1i (xu ) from 0^ is defined as The deviation of the point (v; w) 2 W D1i (y; z ) = jz 0 z^j, where z^ = m A i=1 i (y). The deviation of a point ^ 2i (xu ) from 0^ can be defined in the same way and denoted (r; s) 2 W
z = K11 y1 + K12 y2 + K21 y12 + K22 y1 y2 + K23 y22
(26)
where K11 = 03:091E002; K12 = 3:210E002; K21 = 06:939E006; K22 = 1:488E005; K23 = 08:019E006.
1490
IEEE TRANSACTIONS ON AUTOMATIC CONTROL, VOL. 52, NO. 8, AUGUST 2007
TABLE I ^ (x TAYLOR COEFFICIENTS OF W
) IN EXAMPLE 1
^ 11 (xu1 ) and W ^ 12 (xu1 ) passing 3) There are two 1-D CIMs W through xu1 corresponding to 1 and 2 respectively. We approximate each of them by a fifth-order Taylor expansion ^ 11 (xu1 ) : w1 = W w 2
=
^ 12 (xu1 ) : w1 = W w2 =
5
j =1 5
j =1 5
j =1 5
j =1
a1j vj b1j
vj
(27)
a2j vj b2j vj :
(28)
The Taylor coefficients are shown in Table I. ^ and CIMs be = 0:001. From 4) Let the error threshold between 0 ^ (26), (27), and (28) we can solve for four points on 0
Fig. 1. Credible region of x . x
constructed from the vertices x
;x
;x
and
way. Theoretically, other quadratic form approximations can also work along with CIMs. The concept of credible region is helpful when put the stability boundary approximation into practice. For example, in the calculation of the critical clearing time in power system stability analysis [12], the result of the stable manifold approximation type algorithm can be considered to be acceptable if the exit point of the fault-on trajectory is inside credible region of the specified type-1 UEP. APPENDIX
x11 (066:905; 067:500; 2:833) x21 (27:848; 33:876; 3:439) x12 (130:505; 131:094; 3:336) x22 (33:417; 27:417; 2:728): In this example, we use a quadrilateral neighborhood as the credible region of xu1 . It is constructed from the vertices x11 ; x21 ; x12 and x22 , ^ 11 (xu1 ) and W ^ 12 (xu1 ) as well as the as shown in Fig. 1. The CIMs W ^ are also shown in Fig. 1. Because j2 j approximate stability boundary 0 ^ 12 (xu1 ) is much larger than j1 j, the size of the credible region along W ^ 11 (xu1 ) and we can see the is significant smaller than that along W ^ from W ^ 12 (xu1 ) outside the credible region in obvious deviation of 0 Fig. 1. VII. CONCLUSION A numerical method to determine the credible region of the approximation of stability boundaries in nonlinear dynamic systems has been proposed. The basic idea is to define and calculate some 1- or 2-D CIMs which are exactly located on the stability boundary. Because the CIMs are only 1- or 2-D, it is possible to calculate them with a high precision by higher-order Taylor expansion. The calculated results of CIMs are then used as benchmarks to evaluate the mismatches between the approximate stability boundary and the true boundary. This technique gives a possible way to measure the approximate errors and determine a credible region in which the precision of the approximate stability boundary is under the given threshold. The concept and the algorithm of the credible region can also be applied to high-order systems. Furthermore, although the calculation of CIMs is based on the theory of singular PDEs, it does not necessarily require that the stable boundary to be tested is approximated in the same
Here we only give the Proof of Theorem 3. The Proof of Theorem 2 is similar. Proof of Theorem 3: Let 5 = @'=@r j0 denote the Jacobian matrix of the solution s = (r) of the system of (18). By differentiating both sides of (18) with respect to r , at the point r = 0 we get the following equation:
5A0 0 Jr 01 0 (Jn02 0 Jr 01 ^)5 = 0: Let 52
= [I2 ; 5T ]T , rewrite the equation in matrix form A0 0 5 = 5 2 A0 : Jr 01 (Jn02 0 Jr 01 ^) 2
The matrix on the left side is the Jacobian matrix of system (16). Denoting the first column of 52 by and the second column by , the above equation means that and are, respectively, the real part and imaginary part of a pair of eigenvectors of this Jacobian matrix corresponding to 6 j! . According to the definitions of CIMs in Section IV.A, because 52 also implies the tangent hyperspace of the manifold s = (r) at the origin, so in system (16) s = (r) is tangent to the hyperspace spanned by and at the origin. Transforming back to the original coordinates, we have proved that s = (r) coincides with the 2-D CIM W2 (x0 ; 6 j!) in the neighborhood of x0 in the original system (1).
REFERENCES [1] H. D. Chiang, M. W. Hirsch, and F. F. Wu, “Stability region of nonlinear autonomous dynamical system,” IEEE Trans. Autom. Control, vol. 33, no. 1, pp. 16–27, Jan. 1988.
IEEE TRANSACTIONS ON AUTOMATIC CONTROL, VOL. 52, NO. 8, AUGUST 2007
[2] J. Zaborzsky, G. Huang, B. Zheng, and T. C. Leung, “On the phaseportrait of a class of large nonlinear dynamic systems such as the power system,” IEEE Trans. Autom. Control, vol. 33, no. 1, pp. 4–15, Jan. 1988. [3] H. D. Chiang, C. C. Chu, and G. Cauley, “Direct stability analysis of electric power systems using energy functions: Theory, applications, and perspective,” in Proc. IEEE, Nov. 1995, vol. 83, no. 11, pp. 1497–1529. [4] S. Ushiki, “Analytic expressions of unstable manifolds,” Proc. Jpn Acad. Ser. A.-Math. Sci., vol. 56, no. 6, pp. 239–244, 1980. [5] F. M. A. Salam, A. Araposthathis, and P. P. Varaiya, “Analytic expressions for the unstable manifold at equilibrium points in dynamical systems of differential equations,” in Proc. 22nd IEEE Conf. Decision Control, Dec. 1983, vol. 3, pp. 1389–1392. [6] S. Saha, A. A. Fouad, W. H. Kliemann, and V. Vittal, “Stability boundary approximation of a power system using the real normal form of vector fields,” IEEE Trans. Power Syst., vol. 12, no. 2, pp. 797–802, 1997. [7] R. Qi, D. Cook, W. Kliemann, and V. Vittal, “Visualization of stable manifolds and multidimensional surfaces in the analysis of power system dynamics,” J. Nonlin. Sci., vol. 10, pp. 175–195, 2000. [8] Z. Jing, Z. Jia, and Y. Gao, “Research of the stability region in a power system,” IEEE Trans. Circuits Syst. I, Fundam. Theory Appl., vol. 50, no. 2, pp. 298–304, Feb. 2003. [9] P. A. Cook and A. M. Eskicioglu, “Transient stability analysis of electric power systems by the method of tangent hypersurfaces,” in Inst. Elect. Eng. Proc.-Gener. Transm. Distrib., Jul. 1983, vol. 130, no. 4, pp. 183–193. [10] J. Zaborzsky, G. Huang, B. Zheng, and T. Leung, “New results on stability monitoring on the large electric power systems,” in Proc. 26th IEEE Conf. Decision Control, 1987, pp. 30–40. [11] M. Djukanovic, D. Sobajic, and Y.-H. Pao, “Neural-net based tangent hypersurfaces for transient security assessment of electric power systems,” Int. J. Electr. Power Energy Syst., vol. 16, no. 6, pp. 399–408, 1994. [12] V. Venkatasubramanian and W. Ji, “Numerical approximation of (n 1)-dimensional stable manifolds in large systems such as the power system,” Automatica, vol. 33, no. 10, pp. 1877–1883, Oct. 1997. [13] D. Cheng, J. Ma, Q. Lu, and S. Mei, “Quadratic form of stable submanifold for power systems,” Int. J. Robust Nonlin. Contr., vol. 14, pp. 773–788, 2004. [14] N. Kazantzis, “Singular PDEs and the problem of finding invariant manifold for nonlinear dynamical system,” Phys. Lett. A., vol. 272, no. 4, pp. 257–263, Jul. 2000.
0
IPA Derivatives for Make-to-Stock Production-Inventory Systems With Lost Sales Yao Zhao and Benjamin Melamed
Abstract—This note applies the stochastic fluid model (SFM) paradigm to a class of single-stage, single-product make-to-stock (MTS) production-inventory systems with stochastic demand and random production capacity, where the finished-goods inventory is controlled by a continuous-time base-stock policy and unsatisfied demand is lost. This note derives formulas for infinitesimal perturbation analysis (IPA) derivatives of the sample-path time averages of the inventory level and lost sales with respect to the base-stock level and a parameter of the production rate process. These formulas are comprehensive in that they are exhibited for any initial inventory state, and include right and left derivatives (when they differ). The formulas are obtained via sample path analysis under Manuscript received March 1, 2006; revised November 16, 2006. Recommended by Associate Editor I. Paschalidis. The authors are with the Department of Management Science and Information Systems, Rutgers Business School, Newark and New Brunswick, NJ 08854 USA (e-mail: [email protected]; [email protected]).
1491
very mild regularity assumptions, and are inherently nonparametric in the sense that no specific probability law need be postulated. It is further shown that all IPA derivatives under study are unbiased and fast to compute, thereby providing the theoretical basis for online adaptive control of MTS production-inventory systems. Index Terms—Infinitesimal perturbation analysis (IPA), lost sales, make-to-stock (MTS), production-inventory systems, stochastic fluid models (SFMs).
I. INTRODUCTION This note derives infinitesimal perturbation analysis (IPA) derivatives of selected random variables for a class of make-to-stock (MTS) systems in stochastic fluid model (SFM) setting, where the traditional discrete arrival, service, and departure stochastic processes are replaced by corresponding stochastic fluid-flow rate processes. We henceforth refer to this approach as IPA-over-SFM. The IPA derivatives provide sensitivity information on system metrics with respect to control parameters of interest, and as such can serve as the theoretical underpinnings for online control algorithms. Comprehensive discussions of IPA derivatives and their applications can be found in Fu and Hu [3] and Cassandras and Lafortune [1]. The IPA-over-SFM approach has been successfully applied to theoretical studies of various queuing and production-inventory systems; see, e.g., [2], [6], [4], and [7]. These studies constrain the system to start from a prescribed initial inventory state, and only consider cases where the left and right IPA derivatives coincide. In contrast, Zhao and Melamed [8] considered any initial inventory state for MTS systems with backorders and derived sided IPA derivative formulas as needed. The goal of this note is to derive IPA derivatives for MTS systems with lost sales, and to show them to be unbiased. First, we derive IPA derivatives for the time averages of inventory level and lost sales with respect to the base-stock level for all initial inventory states, including sided derivatives when they differ. We are only aware of one note [6] addressing IPA-over-SFM queues with finite buffers, which can be used to model MTS systems with lost sales, though it constrains the initial condition to an empty buffer. Second, we derive IPA derivatives for the aforementioned metrics with respect to a production rate parameter, including sided derivatives when they differ. We point out that the assumptions in [6] preclude differing left and right IPA derivatives. As will become evident in the sequel, MTS systems with lost sales are also analytically more challenging than MTS systems with backorders, because the inventory state of the former has an extra boundary, a fact that results in more elaborate formulas. The computation of IPA derivatives for all initial conditions is motivated by two reasons. The first reason is to enable potential applications of IPA derivatives to online control of MTS systems driven by nonstationary processes among others. The intent is to adjust the system parameters over time according to the changing statistics of the underlying processes, but not necessarily to optimize it. Clearly, IPA-based online control applications mandate the computation of IPA derivatives for all initial states, as well as all sided derivatives when they differ, since a control action can change the system parameters at any state (which is then considered as the new initial state). It makes little sense to wait for the system to return to certain selected inventory states as this could suspend control actions over extended periods of time. The second reason is that the transient IPA derivatives computed here are exact and unbiased, whereas their asymptotic counterparts may not provide adequate approximations. Furthermore, in order to compute asymptotic IPA derivatives, we still need to obtain their transient counterparts before sending time to infinity. Digital Object Identifier 10.1109/TAC.2007.902760
0018-9286/$25.00 © 2007 IEEE
IEEE TRANSACTIONS ON AUTOMATIC CONTROL, VOL. 52, NO. 8, AUGUST 2007
The Credible Regions on the Approximate Stability Boundaries of Nonlinear Dynamic Systems Yong Min, Lei Chen, Kaiyuan Hou, and Yonghua Song
Abstract—The stability boundary of a nonlinear dynamic system may often be constructed from stable manifolds of type-1 equilibrium points. These stable manifolds are typically approximated numerically, with only local validity. A numerical method is proposed, based on 1- and 2-D characteristic invariant manifolds (CIMs) which lie on the stable manifold and can be computed with greater accuracy than approximations to the entire manifold because of their lower dimension, to estimate Credible Regions in which the mismatch between the approximate and actual boundaries is below a threshold. Index Terms—Credible region, invariant manifold, stability boundary, stability region, stable manifold.
I. INTRODUCTION It has been proved in [1]–[3] that under certain assumptions the stability boundary of a nonlinear autonomous dynamic system equals to the union of the stable manifolds of the unstable equilibrium points (UEPs) lying on the boundary. A fairly large class of systems satisfies these assumptions, so their stability boundaries can be produced by the approximation of stable manifolds of type-1 UEPs.1 The problem of estimating the stable manifold of a hyperbolic equilibrium has been extensively explored during the past thirty or so years and is still under investigation. Most of the approaches developed so far are concentrated on constructing low-order approximations [4]–[13]. Unfortunately, the approximations obtained by these studies are all only valid locally in a neighborhood of the corresponding UEP, and no method has yet been suggested to determine the size of the neighborhood or to measure the approximate errors. This constitutes the main limitation of the stable manifold approximation type approaches because their results have to be confirmed by other type approaches. A recent work [14] proposes a new theoretical approach to the problem of finding invariant manifolds for nonlinear dynamic systems. It states that the local invariant manifolds passing through an equilibrium point can be determined via an associated system of singular first-order quasi-linear PDEs. Unlike previous studies, such as [4] and [12], which deal only with the stable or unstable manifolds of an equilibrium point, the main results of [14] is applicable to not only stable or unstable manifolds but also other invariant manifolds passing through an equilibrium point. This note proposes a numerical approach to determine the credible region around a type-1 UEP on the approximate stability boundary for the given threshold. The credible region is defined as the maximum neighborhood of the type-1 UEP on the stability boundary in which Manuscript received February 15, 2006; revised March 6, 2007 and March 18, 2007. Recommended by Associate Editor J. Berg. This work was supported by the Special Fund of the National Priority Basic Research of China (2004CB217904) and the National Science Foundation of China (50323002). Y. Min and L. Chen are with State Key Lab of Power Systems, Department of Electrical Engineering, Tsinghua University, Beijing 100084, China (e-mail: [email protected]). K. Hou was with the Department of Electrical Engineering, Tsinghua University, Beijing 100084, China. He is now with the Northeast China Grid Company Limited, Shenyang 110006, Liaoning Province, China. Y. Song is with the University of Liverpool, Liverpool L69 3GJ, U.K. Digital Object Identifier 10.1109/TAC.2007.902759 1An equilibrium point is called type-1 UEP if the Jacobian matrix of the system at this point has only one positive real part eigenvalue.
the approximate precision is under the given threshold. The spirit of the work is to find some 1- or 2-D characteristic invariant manifolds (CIMs), which are exactly located on the stability boundary. Based on the theory developed in [14], the Taylor series solution of relevant PDEs is adopted to approximate both the stable manifold and the CIMs passing through a type-1 UEP. Because the CIMs are only 1- or 2-D, it is possible to calculate them with a high precision by higher-order Taylor series expansion. The calculated results of CIMs are used as benchmarks to evaluate the mismatches between the approximate stability boundary and the true boundary. This technique gives a way to measure the approximate errors and determine the credible region. This note proceeds as follows. Section II introduces the fundamental theory of stability regions as well as the main results of [14]. Section III derives an associated system of PDEs to find the stable manifold passing through a type-1 UEP of the original system. The Taylor series solution of PDEs is provided to approximate the local stability boundary. In Section IV, the CIMs passing through equilibrium points are defined. Then the corresponding systems of PDEs are derived and the Taylor series solutions of PDEs are provided similarly. In Section V, the credible region around a type-1 UEP is defined and a procedure to estimate the credible region for a given approximate stability boundary is introduced. An example is given in Section VI. II. PRELIMINARY A. Stability Region and Stability Boundary Consider the following autonomous dynamic system: x _
!
=
f (x);
x
2
R
n
(1)
n R is assumed to be smooth. Let x0 be a hyperbolic where f : Rn equilibrium point of (1). Its stable and unstable manifolds are denoted as W s (x0 ) and W u (x0 ) respectively [1]. If xs is a Stable Equilibrium Point (SEP), A(xs ) = W s (xs ) is defined as the stability region surrounding xs and its boundary is denoted as @ A(xs ). It has been proved in [1], [2] that for systems satisfying the following conditions: 1) all the equilibrium points on @ A(xs ) are hyperbolic; 2) the stable and unstable manifolds of equilibrium points on @ A(xs ) satisfy the transversality condition; 3) every trajectory on @ A(xs ) converges to an equilibrium point; the stability boundary @ A(xs ) is the union of the stable manifolds @ A (xs ); i = of the UEPs on the boundary. Furthermore, let xi 1; 2; . . . ; be all the type-1 UEPs on the boundary, then @ A(xs ) = s s W (xi ), where A(xs ) and W (xi ) denote the closures of A(xs ) and W s (xi ) respectively. This topological property makes it possible to estimate the boundary of the stability region of a SEP by approximating the stable manifolds of all type-1 UEPs on the boundary.
2
B. Finding Invariant Manifolds Through Singular PDES This section introduces a new approach developed in [14] to find invariant manifolds for nonlinear real analytic dynamic systems. First, we rewrite system (1) as follows: y _
=
G(y; z );
2
z _
2
=
H (y; z )
(2)
m p U ;z U and U m ; U p being open where x = (y; z ), with y m p subsets of R and R , respectively, (m + p = n). Without loss of generality, we suppose that the origin (0,0) is an equilibrium point, which means G(0; 0) = 0, H (0; 0) = 0.
0018-9286/$25.00 © 2007 IEEE
IEEE TRANSACTIONS ON AUTOMATIC CONTROL, VOL. 52, NO. 8, AUGUST 2007
If (@G)=(@z )j(0;0) = 0, system (2) can be rewritten as follows:
y_ = Ay + g(y; z );
z_ = Bz + Cy + h(y; z )
1487
J 01 J1 , with J0 being the element at J2 J0 the bottom right corner of J (x ) such that J1 2 R( 01)21 ; J2 2 R12( 01) , and J 01 2 R( 01)2( 01) . Then taking notice of J (x )q = q , (7) is rewritten as follows: Rewrite
(3)
where g (y; z ) and h(y; z ) are analytic functions whose Taylor series expansions begin with terms of order higher than one. Then we consider the following first-order PDEs:
J (x
n
u)
n
as
n
u
(0) = 0
(4) where : Rm ! Rp is the unknown vector field of (4). The above system of PDEs is singular at the origin (0,0) because the principal part G(y; z ) = Ay + g (y; z ) vanishes at the origin. Only under some specific conditions, the existence and uniqueness of the solution can be ensured. So the following theorem is given in [14]. Theorem 1([14]): Suppose the origin (0,0) is an equilibrium point of system (2), let k1 ; k2 ; . . . ; km be the eigenvalues of the matrix A = (@G)=(@y )j(0;0) and 1 ; 2 ; . . . ; p be the eigenvalues of the matrix B = (@H )=(@z )j(0;0) . Assume the following: 1) (@G)=(@z )j(0;0) = 0; 2) The real parts of all the kj are different from zero and have the same sign; 3) kj and i are not related by any equation of the form m
mk j
j
=
;
i = 1; . . . ; p
i
(5)
j =1
where all the mj are nonnegative integers satisfying
0.
m j =1
m > j
Then there exists a neighborhood V Rm of y0 = 0 and a unique and locally analytic mapping : V ! Rp such that
S = f(y; z ) 2 V
2 R : z = (y); (0) = 0g p
(6)
is an analytic local invariant manifold of system (2) that passes through the origin, with z = (y ) being the solution of (4). The proof and discussions of Theorem 1 can be found in [14]. III. STABLE MANIFOLD OF TYPE-1 UEP Suppose that system (1) satisfies the conditions given in Section II.A and xu is a type-1 UEP of the system. Let J (xu ) be the Jacobian matrix at xu and u be its only eigenvalue with positive real part. The eigenvector of J (xu ) corresponding to u is denoted as2
q = [q1 ; . . . ; q 01 ; q n
Construct the matrix Q =
x
u)
I 01 n
0
n]
T
= [q^
q^ q
n
T
;q
n]
T
y_ = (J 01 0 q^J2 =q )y + g(y; z ) z_ = J2 y=q + z + h(y; z ) n
n
Q01 (x 0
, system (1) is converted to the following form:
x^_ = Q01 J (x )Qx^ + R(^x) u
(7)
x) denotes terms of order higher than one. where R(^ 2We assume q 6= 0 in the following deduction. In case where q = 0, we could find a nonzero element in q , say q . Then if we exchange x = f (x) with x = f (x) in system (1), q and q in q would be switched.
(8)
u
x1 ; x^2 ; . . . ; x^n01 ]T ; z = x^n ; g(y; z ) and h(y; z ) denote where y = [^ corresponding terms of order higher than one. Because u is the only eigenvalue with positive real part, according to Theorem 1, the unique and locally analytic solution z = (y ) of the following PDEs: @ ((J 0 q^J =q )y + g(y; z )) = J y=q + + h(y; z ); 01 2 2 @y (0) = 0 (9) n
n
n
u
coincides with the stable manifold in the neighborhood of xu . Furthermore, if xu 2 @A(xs ), the solution z = (y ) will coincide with @A(xs ) in the neighborhood of xu . Generally it is not possible to solve for an analytic expression of z = (y) from the PDEs (9). We use the following multivariate series solution to approximate it:
z=
1
A (y) = Ey + y F y + 1 1 1 i
i=1
T
(10)
where Ai (y ) is the Taylor term of ith order and is an ith order homogeneous polynomial. Theoretically, substituting (10) into (9) and equating the same order Taylor coefficients at both sides, we can solve for all Taylor coefficients by a recursion method [14]. However, if the system’s order is higher than 3, the calculation is very difficult and only the second-order Taylor series expansion can usually be available. IV. CIMS PASSING THROUGH EQUILIBRIUM POINTS A. Definitions of CIMs Let x0 be a hyperbolic equilibrium point of system (1) and J (x0 ) be the Jacobian matrix at x0 . We denote one of the real eigenvalues of J (x0 ) by and a pair of complex conjugate eigenvalues of J (x0 ) by 6 j! . The corresponding eigenvectors are denoted by and 6 j respectively. Consider the linear part of system (1) at x0
x_ = J (x0 )(x 0 x0 );
:
and define x ^ =
n
u
n
@ (Ay + g(y; (y))) = B(y) + Cy + h(y; (y)); @y
n
u n
x2R : n
(11)
We know from the theory of linear systems that in state space passing through x0 there exists a pair of linear trajectories in the direction of , and a set of trajectories forming a hyperplane which is spanned by and . Because the line or hyperplane formed by the trajectories is actually an invariant manifold and is corresponding to the eigenvector(s) of the system, we call it CIMs. The definition is as follows. Definition 1: For system (11), a pair of linear trajectories in the direction of is defined as the 1-D CIM of x0 relative to the eigenvalue , denoted by W1 (x0 ; ) or W1 (x0 ), and a hyperplane spanned by the vectors and is defined as the 2-D CIM of x0 relative to the eigenvalues 6 j! , denoted by W2 (x0 ; 6 j! ) or W2 (x0 ).
1488
IEEE TRANSACTIONS ON AUTOMATIC CONTROL, VOL. 52, NO. 8, AUGUST 2007
According to the Hartman-Grobman’s theorem, system (1) is topologically equivalent to its linear part (11) in the neighborhood of x0 . Therefore, W1 (x0 ) or W2 (x0 ) is uniquely homologized to a 1- or a 2-D manifold of system (1) respectively. Thus CIMs of the nonlinear system (1) are defined as follows. Definition 2: The 1-D manifold of system (1), homologizing to the CIM W1 (x0 ) of its linear part (11), is defined as the 1-D CIM of (1), still denoted by W1 (x0 ). The 2-D manifolds of system (1), homologizing to the CIM W2 (x0 ) of its linear part, is defined as the 2-D CIM of (1), still denoted by W2 (x0 ). According to the above definitions, it is obvious that if x0 is a type-1 UEP on @A(xs ), then W1 (x0 ; ) @A(xs ) if < 0 and W2 (x0 ; 6 j! ) @A(xs ) if < 0. B. CIMs in the Form of PDEs The CIMs of a linear system can be easily found from the eigenvectors of the Jacobian matrix. For a nonlinear system, the CIMs can be found via Theorem 1. For the 1-D CIM W1 (x0 ), suppose = [1 ; 2 ; . . . ; n ]T is the eigenvector of J T (x0 ) corresponding to (in other words, T is the left eigenvector of J (x0 )). Assume 1 6= 0. Let ^ = [2 ; . . . ; n ] and construct the ma1 ^ J1 Jw . Rewrite J (x0 ) as ; with trix T = Jv Jn01 0 In01 J1 being the element at the top left corner of J (x0 ) such that Jw 2 R12(n01) ; Jv 2 R(n01)21 , and Jn01 2 R(n01)2(n01) . Then taking notice of J T (x0 ) = , the original system (1) is converted to the following form:
v_ = v + gv (v; w) w_ = Jv v=1 + (Jn01 0 Jv ^=1 )w + hw (v; w)
(12)
where x ^ = T (x 0 x 0 ); v = x ^1 ; w = [^ x2 ; x^3 ; . . . ; x^n ]T ; gv (v; w) and hw (v; w) denote terms of order higher than one. The first two conditions of Theorem 1 are automatically satisfied. The third condition is generally true in the actual system. So a unique analytic local invariant manifold
S1
f
= (v; w )
2 V 2 Rn0
1
:
w = '(v ); '(0) = 0g
(13)
exists in a neighborhood V 2 R of v = 0, with w = '(v ) being the unique and locally analytic solution of the system of the following PDEs:
@' ^=1 )' + hw (v; w ); (v + gv (v; w )) = Jv v=1 + (Jn01 0 Jv @v (14) '(0) = 0: Theorem 2: The unique analytic local invariant manifold defined in (13) coincides with the 1-D CIM W1 (x0 ; ) in the neighborhood of x0 in (1). Proof: See Appendix. For the 2-D CIM W2 (x0 ; 6 j! ), suppose 6 j = T T are the eigenvectors [ 1 ; 2 ; . . . ; n ] 6 j : [ 1 ; 2 ; . . . ; n ] 0! ; = of J T (x0 ) corresponding to 6 j! . Let A0 = ! 1 2 3 1 1 1 n ; and ^ = . 1 2 3 1 1 1 n
From J T (x0 )( 6 j ) = ( 6 j! )( 6 j ) we have [ ; ^ ]J ( x 0 ) =
A0 [ ; ^]
(15)
^ Js J2 ], with J2 ]. Rewrite J (x0 ) as [ In02 Jr Jn02 being the 2 2 2 submatrix at the top left corner of J (x0 ) such that Js 2 R22(n02) ; Jr 2 R(n02)22 ; Jn02 2 R(n02)2(n02) . Then taking notice of (15), the original system (1) can be converted to the following form: Define P = [
0
r_ = A0 r + gr (r; s) s_ = Jr 01 r + (Jn02 0 Jr 01 ^)s + hs (r; s)
(16)
where x ~ = P (x 0 x0 ); r = [~ x1 ; x~2 ]T ; s = [~ x3 ; x~4 ; . . . ; x~n ]T ; gr (r; s) and hs (r; s) denote terms of order higher than one. So there exists a neighborhood V 2 R2 of r = 0 and a unique analytic local invariant manifold
S2
f
= (r; s)
2 V 2 Rn0
2
:
s=
(r );
(0) = 0
g
(17)
with s = (r) being the unique and locally analytic solution of the system of the following PDEs:
@ 01 01 ^) (A0 r + gr (r; s)) = Jr r + (Jn02 0 Jr @r (0) = 0:
+ hs (r; s);
(18)
Theorem 3: The unique analytic local invariant manifold defined in (17) coincides with the 2-D CIM W2 (x0 ; 6 j! ) in the neighborhood of x0 in system (1). Proof: See the Appendix. C. Series Solutions of CIMs The Taylor series of the unique and locally analytic solutions w = '(v ) of the system of (14) and s = (r) of the system of (18) are denoted as follows:
w= s=
1 j =1
aj v j
1
bjh r1j r2h :
(19) (20)
Substituting (19) and (20) into (14) and (18), respectively, we can solve for aj and bjh . Because there are only one variable v and two variables r1 ; r2 in (19) and (20), respectively, it is feasible to solve for the highorder Taylor coefficients of CIMs with the aid of a symbolic software package. These high-order Taylor series approximations of CIMs can be then used as benchmarks to test the accuracy of the approximate stability boundary. V. CREDIBLE REGION ON THE APPROXIMATE STABILITY BOUNDARY AROUND TYPE-1 UEP A. Credible Region and Its Definition The credible region on the approximate stability boundary around a type-1 UEP is a region in which the maximum mismatch between the actual boundary and the approximate boundary is less than a given tolerance. The definition is as follows.
IEEE TRANSACTIONS ON AUTOMATIC CONTROL, VOL. 52, NO. 8, AUGUST 2007
Definition 3: Suppose xu is a type-1 UEP on @A(xs ). Let 0 de^ denote the note the true stability boundary passing through xu and 0 approximate one. For any given positive number , the credible region ^ around xu with a credible level is defined as on 0
^ : jx 0 0 j g D ( ) = fx 2 0
(21)
where jx 0 0j denotes the shortest Euclidean distance from a point 2 0^ to the manifold 0. To make the above definition more practical, we consider a neighborhood B (") of xu . Without loss of generality, let B (") be a spherical neighborhood
x
B ( " ) = f x 2 R n : jx 0 x u j " g
as D2i (y; z ). D1i (y; z ) and D2i (y; z ) are then used as the indices to ^ and 0. estimate the approximation error between 0 For a given error threshold , define
X1 = f(y; z) 2 Rn : D1i (y; z) = ; i = 1; . . . ; kg; X2 = f(y; z) 2 Rn : D2i (y; z) = ; i = 1; . . . ; lg; D = fj(y; z) 0 xu j 2 R : (y; z) 2 X1 [ X2 g: The sets X1 and X2 contain all the points on the CIMs which have the ^ . The set D contains all the Euclidean distances given deviation from 0 from the points in X1 [ X2 to xu . Therefore, the minimum element in set D can be used as the maximum available radius "max of the credible region B (")
(22)
where " > 0 is the radius of B ("). ^ = 0^ \ B ("), then for a given > 0; ^ Define = 0 \ B (") and
is a credible region with credible level if there is a positive number " ^ , where @ ^ denotes the boundary such that jx 0 j for any x 2 @
^. of
Due to the difficulties in calculating the accurate 0 or , it is ob^ can not be exactly determined for a given tolerance . In vious that
the following part of this section, we will use the high-order series ap^. proximation of CIMs to estimate
B. Estimating the Credible Region via CIMS The algorithm of estimating the credible region via the concept of CIMs is based on the following two assumptions. 1) The approximation of CIMs is more accurate than that of the stability boundary. ^ appears only on one of 2) The maximum mismatch between 0 and 0 the CIMs. The first one is generally true but the validity of the second one is not easy to confirm. However, even if the second one is not true, the result obtained by the following method is still valuable for estimating ^. the mismatches between 0 and 0 Suppose that xs is a SEP in (1) and let J (xu ) be the Jacobian matrix at a type-1 UEP xu on @A(xs ). From Section IV we know there is an 1-D CIM W1i (xu ) = W1 (xu ; i ) for any negative real eigenvalue i of J (xu ), and a 2-D CIM W2i (xu ) = W2 (xu ; i 6 j!i ) for any pair of complex eigenvalues i 6j!i with i < 0. According to the definitions of CIMs, we know W1i (xu ) 0; i = 1; . . . ; k and W2i (xu ) 0; i = 1; . . . ; l, where k denotes the number of negative real eigenvalues and l denotes the number of pairs of complex conjugate eigenvalues with negative real parts. Therefore, W1i (xu ) and W2i (xu ) can be used to estimate the ap^ and 0 with the expansion of B ("). Beproximation error between 0 cause it is difficult to find actual W1i (xu ) and W2i (xu )), their high ^ 1i (xu ) and W ^ 2i (xu ) are used precision Taylor series approximation W instead. Based on the transformations defined in Section III and IV, a point ^ 1i (xu ) is converted to the state space of (8) as follows: (v; w) 2 W
y v = Q01 T 01 : z w
1489
"max = minfd 2 Dg:
(24)
Remark 1: The calculation of the coefficients of the terms with order higher than 7 in (19) and (20) may be involved (especially for bjh ) and, on the other hand, is usually not helpful to improve the precision of the approximation of CIMs because of the limitation of the region of convergence of the solutions of (14) and (18). Remark 2: The proposed algorithm of credible region can be applied to the high-order systems and can work with any other stability boundary approximation approaches as long as D1i (y; z ) and D2i (y; z ) can be calculated. VI. EXAMPLE In this section we give an example to describe the steps to determine the credible region. The approximate stability boundary in the example is calculated by the method developed in Section III. However, the results of other quadratic form approximations can also be used instead. Example 1: The state equations of a 2-machine power system are as follows:
!_ 1 = P1M + P1E sin( + 12 ) 0 !0 D1 !1 =TJ 1 !_ 2 = P2M 0 P2E sin( 0 12 ) 0 !0 D2 !2 =TJ 2 _ = !1 0 !2
(25)
where P1E = 0115:50; P2E = 0144:38; !0 = 100; P1M = 52:17; P2M = 10:94; 12 = 0:27; TJ 1 = 5; TJ 2 = 4; D1 = 0:04; D2 = 0:024. The SEP of the system is xs (0; 0; 0:196). There are two type-1 UEPs on the stability boundary of xs : xu1 (30:569; 30:569; 3:085); xu2 (30:569; 30:569; 03:199).
Now we take the following steps to calculate the credible region on the approximate stability boundary passing through xu1 . 1) Calculate the Jacobian matrix J (xu1 ) and its eigenvalues. It has two negative real eigenvalues: 1 02:228; 2 016:920, and one positive real eigenvalue: u 14:750. ^ passing through xu1 by a 2) Approximate the stability boundary 0 second-order Taylor series expansion
(23)
^ 1i (xu ) from 0^ is defined as The deviation of the point (v; w) 2 W D1i (y; z ) = jz 0 z^j, where z^ = m A i=1 i (y). The deviation of a point ^ 2i (xu ) from 0^ can be defined in the same way and denoted (r; s) 2 W
z = K11 y1 + K12 y2 + K21 y12 + K22 y1 y2 + K23 y22
(26)
where K11 = 03:091E002; K12 = 3:210E002; K21 = 06:939E006; K22 = 1:488E005; K23 = 08:019E006.
1490
IEEE TRANSACTIONS ON AUTOMATIC CONTROL, VOL. 52, NO. 8, AUGUST 2007
TABLE I ^ (x TAYLOR COEFFICIENTS OF W
) IN EXAMPLE 1
^ 11 (xu1 ) and W ^ 12 (xu1 ) passing 3) There are two 1-D CIMs W through xu1 corresponding to 1 and 2 respectively. We approximate each of them by a fifth-order Taylor expansion ^ 11 (xu1 ) : w1 = W w 2
=
^ 12 (xu1 ) : w1 = W w2 =
5
j =1 5
j =1 5
j =1 5
j =1
a1j vj b1j
vj
(27)
a2j vj b2j vj :
(28)
The Taylor coefficients are shown in Table I. ^ and CIMs be = 0:001. From 4) Let the error threshold between 0 ^ (26), (27), and (28) we can solve for four points on 0
Fig. 1. Credible region of x . x
constructed from the vertices x
;x
;x
and
way. Theoretically, other quadratic form approximations can also work along with CIMs. The concept of credible region is helpful when put the stability boundary approximation into practice. For example, in the calculation of the critical clearing time in power system stability analysis [12], the result of the stable manifold approximation type algorithm can be considered to be acceptable if the exit point of the fault-on trajectory is inside credible region of the specified type-1 UEP. APPENDIX
x11 (066:905; 067:500; 2:833) x21 (27:848; 33:876; 3:439) x12 (130:505; 131:094; 3:336) x22 (33:417; 27:417; 2:728): In this example, we use a quadrilateral neighborhood as the credible region of xu1 . It is constructed from the vertices x11 ; x21 ; x12 and x22 , ^ 11 (xu1 ) and W ^ 12 (xu1 ) as well as the as shown in Fig. 1. The CIMs W ^ are also shown in Fig. 1. Because j2 j approximate stability boundary 0 ^ 12 (xu1 ) is much larger than j1 j, the size of the credible region along W ^ 11 (xu1 ) and we can see the is significant smaller than that along W ^ from W ^ 12 (xu1 ) outside the credible region in obvious deviation of 0 Fig. 1. VII. CONCLUSION A numerical method to determine the credible region of the approximation of stability boundaries in nonlinear dynamic systems has been proposed. The basic idea is to define and calculate some 1- or 2-D CIMs which are exactly located on the stability boundary. Because the CIMs are only 1- or 2-D, it is possible to calculate them with a high precision by higher-order Taylor expansion. The calculated results of CIMs are then used as benchmarks to evaluate the mismatches between the approximate stability boundary and the true boundary. This technique gives a possible way to measure the approximate errors and determine a credible region in which the precision of the approximate stability boundary is under the given threshold. The concept and the algorithm of the credible region can also be applied to high-order systems. Furthermore, although the calculation of CIMs is based on the theory of singular PDEs, it does not necessarily require that the stable boundary to be tested is approximated in the same
Here we only give the Proof of Theorem 3. The Proof of Theorem 2 is similar. Proof of Theorem 3: Let 5 = @'=@r j0 denote the Jacobian matrix of the solution s = (r) of the system of (18). By differentiating both sides of (18) with respect to r , at the point r = 0 we get the following equation:
5A0 0 Jr 01 0 (Jn02 0 Jr 01 ^)5 = 0: Let 52
= [I2 ; 5T ]T , rewrite the equation in matrix form A0 0 5 = 5 2 A0 : Jr 01 (Jn02 0 Jr 01 ^) 2
The matrix on the left side is the Jacobian matrix of system (16). Denoting the first column of 52 by and the second column by , the above equation means that and are, respectively, the real part and imaginary part of a pair of eigenvectors of this Jacobian matrix corresponding to 6 j! . According to the definitions of CIMs in Section IV.A, because 52 also implies the tangent hyperspace of the manifold s = (r) at the origin, so in system (16) s = (r) is tangent to the hyperspace spanned by and at the origin. Transforming back to the original coordinates, we have proved that s = (r) coincides with the 2-D CIM W2 (x0 ; 6 j!) in the neighborhood of x0 in the original system (1).
REFERENCES [1] H. D. Chiang, M. W. Hirsch, and F. F. Wu, “Stability region of nonlinear autonomous dynamical system,” IEEE Trans. Autom. Control, vol. 33, no. 1, pp. 16–27, Jan. 1988.
IEEE TRANSACTIONS ON AUTOMATIC CONTROL, VOL. 52, NO. 8, AUGUST 2007
[2] J. Zaborzsky, G. Huang, B. Zheng, and T. C. Leung, “On the phaseportrait of a class of large nonlinear dynamic systems such as the power system,” IEEE Trans. Autom. Control, vol. 33, no. 1, pp. 4–15, Jan. 1988. [3] H. D. Chiang, C. C. Chu, and G. Cauley, “Direct stability analysis of electric power systems using energy functions: Theory, applications, and perspective,” in Proc. IEEE, Nov. 1995, vol. 83, no. 11, pp. 1497–1529. [4] S. Ushiki, “Analytic expressions of unstable manifolds,” Proc. Jpn Acad. Ser. A.-Math. Sci., vol. 56, no. 6, pp. 239–244, 1980. [5] F. M. A. Salam, A. Araposthathis, and P. P. Varaiya, “Analytic expressions for the unstable manifold at equilibrium points in dynamical systems of differential equations,” in Proc. 22nd IEEE Conf. Decision Control, Dec. 1983, vol. 3, pp. 1389–1392. [6] S. Saha, A. A. Fouad, W. H. Kliemann, and V. Vittal, “Stability boundary approximation of a power system using the real normal form of vector fields,” IEEE Trans. Power Syst., vol. 12, no. 2, pp. 797–802, 1997. [7] R. Qi, D. Cook, W. Kliemann, and V. Vittal, “Visualization of stable manifolds and multidimensional surfaces in the analysis of power system dynamics,” J. Nonlin. Sci., vol. 10, pp. 175–195, 2000. [8] Z. Jing, Z. Jia, and Y. Gao, “Research of the stability region in a power system,” IEEE Trans. Circuits Syst. I, Fundam. Theory Appl., vol. 50, no. 2, pp. 298–304, Feb. 2003. [9] P. A. Cook and A. M. Eskicioglu, “Transient stability analysis of electric power systems by the method of tangent hypersurfaces,” in Inst. Elect. Eng. Proc.-Gener. Transm. Distrib., Jul. 1983, vol. 130, no. 4, pp. 183–193. [10] J. Zaborzsky, G. Huang, B. Zheng, and T. Leung, “New results on stability monitoring on the large electric power systems,” in Proc. 26th IEEE Conf. Decision Control, 1987, pp. 30–40. [11] M. Djukanovic, D. Sobajic, and Y.-H. Pao, “Neural-net based tangent hypersurfaces for transient security assessment of electric power systems,” Int. J. Electr. Power Energy Syst., vol. 16, no. 6, pp. 399–408, 1994. [12] V. Venkatasubramanian and W. Ji, “Numerical approximation of (n 1)-dimensional stable manifolds in large systems such as the power system,” Automatica, vol. 33, no. 10, pp. 1877–1883, Oct. 1997. [13] D. Cheng, J. Ma, Q. Lu, and S. Mei, “Quadratic form of stable submanifold for power systems,” Int. J. Robust Nonlin. Contr., vol. 14, pp. 773–788, 2004. [14] N. Kazantzis, “Singular PDEs and the problem of finding invariant manifold for nonlinear dynamical system,” Phys. Lett. A., vol. 272, no. 4, pp. 257–263, Jul. 2000.
0
IPA Derivatives for Make-to-Stock Production-Inventory Systems With Lost Sales Yao Zhao and Benjamin Melamed
Abstract—This note applies the stochastic fluid model (SFM) paradigm to a class of single-stage, single-product make-to-stock (MTS) production-inventory systems with stochastic demand and random production capacity, where the finished-goods inventory is controlled by a continuous-time base-stock policy and unsatisfied demand is lost. This note derives formulas for infinitesimal perturbation analysis (IPA) derivatives of the sample-path time averages of the inventory level and lost sales with respect to the base-stock level and a parameter of the production rate process. These formulas are comprehensive in that they are exhibited for any initial inventory state, and include right and left derivatives (when they differ). The formulas are obtained via sample path analysis under Manuscript received March 1, 2006; revised November 16, 2006. Recommended by Associate Editor I. Paschalidis. The authors are with the Department of Management Science and Information Systems, Rutgers Business School, Newark and New Brunswick, NJ 08854 USA (e-mail: [email protected]; [email protected]).
1491
very mild regularity assumptions, and are inherently nonparametric in the sense that no specific probability law need be postulated. It is further shown that all IPA derivatives under study are unbiased and fast to compute, thereby providing the theoretical basis for online adaptive control of MTS production-inventory systems. Index Terms—Infinitesimal perturbation analysis (IPA), lost sales, make-to-stock (MTS), production-inventory systems, stochastic fluid models (SFMs).
I. INTRODUCTION This note derives infinitesimal perturbation analysis (IPA) derivatives of selected random variables for a class of make-to-stock (MTS) systems in stochastic fluid model (SFM) setting, where the traditional discrete arrival, service, and departure stochastic processes are replaced by corresponding stochastic fluid-flow rate processes. We henceforth refer to this approach as IPA-over-SFM. The IPA derivatives provide sensitivity information on system metrics with respect to control parameters of interest, and as such can serve as the theoretical underpinnings for online control algorithms. Comprehensive discussions of IPA derivatives and their applications can be found in Fu and Hu [3] and Cassandras and Lafortune [1]. The IPA-over-SFM approach has been successfully applied to theoretical studies of various queuing and production-inventory systems; see, e.g., [2], [6], [4], and [7]. These studies constrain the system to start from a prescribed initial inventory state, and only consider cases where the left and right IPA derivatives coincide. In contrast, Zhao and Melamed [8] considered any initial inventory state for MTS systems with backorders and derived sided IPA derivative formulas as needed. The goal of this note is to derive IPA derivatives for MTS systems with lost sales, and to show them to be unbiased. First, we derive IPA derivatives for the time averages of inventory level and lost sales with respect to the base-stock level for all initial inventory states, including sided derivatives when they differ. We are only aware of one note [6] addressing IPA-over-SFM queues with finite buffers, which can be used to model MTS systems with lost sales, though it constrains the initial condition to an empty buffer. Second, we derive IPA derivatives for the aforementioned metrics with respect to a production rate parameter, including sided derivatives when they differ. We point out that the assumptions in [6] preclude differing left and right IPA derivatives. As will become evident in the sequel, MTS systems with lost sales are also analytically more challenging than MTS systems with backorders, because the inventory state of the former has an extra boundary, a fact that results in more elaborate formulas. The computation of IPA derivatives for all initial conditions is motivated by two reasons. The first reason is to enable potential applications of IPA derivatives to online control of MTS systems driven by nonstationary processes among others. The intent is to adjust the system parameters over time according to the changing statistics of the underlying processes, but not necessarily to optimize it. Clearly, IPA-based online control applications mandate the computation of IPA derivatives for all initial states, as well as all sided derivatives when they differ, since a control action can change the system parameters at any state (which is then considered as the new initial state). It makes little sense to wait for the system to return to certain selected inventory states as this could suspend control actions over extended periods of time. The second reason is that the transient IPA derivatives computed here are exact and unbiased, whereas their asymptotic counterparts may not provide adequate approximations. Furthermore, in order to compute asymptotic IPA derivatives, we still need to obtain their transient counterparts before sending time to infinity. Digital Object Identifier 10.1109/TAC.2007.902760
0018-9286/$25.00 © 2007 IEEE
