H∞ Optimal Filtering and Control of Wind Energy Conversion Systems Hoa M. Nguyen
D. Subbaram Naidu
S. Hossein Mousavinezhad
Department of Electrical Engineering School of Engineering Idaho State University Pocatello, Idaho 83209-8060 Email:
[email protected] Department of Electrical Engineering School of Engineering Idaho State University Pocatello, Idaho 83209-8060 Email:
[email protected] Department of Electrical Engineering School of Engineering Idaho State University Pocatello, Idaho 83209-8060 Email:
[email protected] Abstract—This paper presents a reduced-order H∞ optimal control for wind energy conversion systems. Two different timescale (slow and fast) dynamics of wind energy conversion systems are separated and processed independently using the singular perturbation theory. By using the decomposition technique, low-order, independent H∞ optimal filters and controllers are obtained, which provide computational advantages and enable implementations with different sampling rates. The control robustness and efficiency are shown by computer simulations.
I. I NTRODUCTION Wind energy has become the fastest growing energy among other renewable energy sources for the last few decades [1], [2], [3]. As wind energy conversion systems (WECSs) have been larger in sizes, the power conversion efficiency has become more important. Structural designs and operations have been optimized in order to increase the overall efficiency [4], [5]. One of the main control problems in WECSs is to find a good balance between two conflicting control goals which are the maximizing captured power from wind and minimizing fatigue loads caused by control effort [6], [7]. To solve the above control problem, the optimal control theory is a suitable method because it allows to find optimal solutions to control problems where control objectives and control effort are given in terms of cost functions [8]. Two common forms of the optimal control are Linear Quadratic Regulator (LQR) used for deterministic systems and Linear Quadratic Gaussian (LQG) used for stochastic systems. These optimal control versions have been presented to WECSs. For instance, in [9], the LQR was proposed to control the pitch system which reduces the power when wind speed exceeds the rated speed. Wu et al. in [10] also designed an LQR to regulate the output power of a wind turbine with doubly-fed induction generator based on a linear model derived by the exact linearization method which provides better performance than the approximate linearization method. Simulation results show better transient responses and fault ride-through ability of the LQR in comparison with PI control. The LQG, compared with the LQR, has more applications because it takes into account noises present in WECSs. For example, the LQG control was combined with a neural controller to regulate the output power and decrease cyclic loads on the
drive train of a wind turbine in [11]. Simulations show the maximum power tracking is obtained and mechanical loads are significantly reduced. Similarly, the control problems in [11] are also addressed in the report [12]. However, [12] is intended to study the effect of the LQG control technique in reducing fatigue loads for differently complex models of a wind turbine. Alternately, researchers in [13] investigated on the LQG control for multiple linear models obtained from different operating conditions of a wind turbine, in which each LQG controller was designed for each linear model and the control actions are switched based on an interpolation algorithm. However, none of the reports above and in [14] considers different time-scale dynamics in WECSs which contain both mechanical and electrical components, thereby exhibiting slow and fast phenomena. An interesting fact is such systems with slow and fast dynamics can be separated into slow and fast subsystems by singular perturbation methods [15], [16]. As a result, each subsystem can be analyzed and synthesized independently. To our best knowledge, there are not many applications of the singular perturbation theory on WECSs so far. Indeed, the technical note [17] is one of few studies that exploits the time-scale decomposition technique on WECSs. However, the paper is only focused on the LQR control under deterministic conditions which are not guaranteed in practice due to the presence of system noises and uncertainty. In this paper, H∞ optimal filters and controllers are designed for a singularly perturbed Wind Energy Conversion System (WECS) with a Permanent Magnet Synchronous Generator (PMSG). The advantage of the H∞ optimal control over the LQR and LQG control is that the H∞ controller is robust to internal and external noises. Additionally, the H∞ filtering does not require noises to be Gaussian stochastic processes as in the Kalman filtering. This paper is organized as follows. Section II describes a nonlinear of the PMSG-based WECS and its linear state space model. Section III presents the H∞ filtering and control methods for singularly perturbed systems. Simulations verifying the H∞ optimal filtering and control performance are given in Section IV.
978-1-4673-5208-6/13/$31.00 ©2013 IEEE
pitch angle β and the tip-speed ratio λ, where λ is defined as
II. M ODELING A. Nonlinear Model
ωr R , (2) V with ωr is the wind rotor speed. The torque coefficient CQ can be approximated by a polynomial of λ and β [6] λ=
In general, a WECS is a converter that transforms wind power into electrical power. Although a WECS involves many different devices, these devices can be grouped as functional blocks with respect to power flow as shown in Fig. 1. The aerodynamics block converts the power from wind into mechanical power (rotations) and these rotations are increased and transmitted to the generator block by the drive train block. Finally, the mechanical power is transformed into electrical power by the generator block.
CQ (λ, β)
=
= k0 + k1 λ + k2 λ2 + k3 λ3 + k4 λ4
WECS Block Diagram
5
1 ρπR3 CQ (λ, β)V 2 , 2
(4)
6
+k5 λ + k6 λ .
Though a WECS can operate over a wide range of wind speeds, basically it is only active in two regions or load conditions including partial load and full load as depicted in Fig. 2. In the partial load, the WECS is operating in such a way that maximum power is captured as the wind speed (V ) changes from the cut-in speed (Vcut−in) to the rated-speed (Vrated ). When V is greater than Vrated , the WECS switches to the full load operation where the output power is clipped at the rated power (Prated ). The WECS is inactive as V is less than Vcut−in or greater than the cut-out speed (Vcut−out ). For the control purpose, only aerodynamics, drive train dynamics, and generator dynamics are taken into account. The aerodynamics is the most complicated dynamics [18]. It takes the wind speed V and wind rotor speed ωr as inputs. The output is given in terms of aerodynamic torque Tr as follows Tr
(3)
It should be emphasized that for variable speed, variable pitch WECSs, there are two primary control strategies corresponding to two load conditions. In the partial load regime, the pitch angle is kept zero and the generator speed is adjusted to maintain the optimal tip-speed ratio as the wind speed changes. In the full load regime, the generator speed is kept constant and the pitch angle is controlled by the pitch system to reduce the aerodynamic power. In this paper, only the control in the partial load is considered, therefore, the torque coefficient in this region is only dependent on the tip-speed ratio and given as [7] CQ (λ)
Fig. 1.
c0 + c1 β + c2 β 2 + c3 β 3 + c4 λβ +c5 λβ 2 + c6 λβ 3 .
=
The drive train block consists of a low-speed shaft and a high-speed shaft connected to each other through a gearbox which increases the rotational speed. The drive train can be represented by a rigid or flexible model [7]. In this study, a flexible drive train is used as shown in Fig. 3 which has the model
(1) Fig. 3.
Flexible Drive Train Block
where ρ is the air density, R is the radius of the wind rotor plane, CQ is the torque coefficient given as a function of the ω˙ r
=
ω˙ g
=
T˙H
=
i 1 TH + Tr , ηJr Jr 1 1 TH − Tg , Jg Jg
−
iKg ωr − Kg ωg − Bg ( +
Fig. 2.
Wind Power versus Wind Speed
(5) (6) 1 i2 + )TH Jg ηJr
(7)
iBg Bg Tr + Tg , Jr Jg
where ωg is the generator speed, TH is the internal torque, Jr is the wind rotor inertia, Jg is the generator inertia, Kg is the stiffness coefficient of the high-speed shaft (the generator
shaft), Bg is the damping coefficient of the high-speed shaft, i is the gearbox ratio, and η is the gearbox efficiency. Essentially, asynchronous and synchronous generators are two primary types of generator which have been used in WECSs. Three popular generators are Squirrel Cage Induction Generator (SCIG), Doubly Fed Induction Generator (DFIG), and Permanent Magnet Synchronous Generator (PMSG) [19]. This study is focused on the PMSG which has the model in the direct and quadrature, (d, q), axes as follows [7] i˙ ds i˙ qs Tg
Rs ids + Ld Rs = − iqs − Lq = pφm iq ,
= −
pLq 1 iqs ωg − uds , Ld Ld p 1 (Ld ids − φm )ωg − uqs , Lq Lq
(8) (9) (10)
where ids and iqs are the d and q components of the stator current; Ld and Lq are the d and q components of the stator inductance; uds and uqs are the d and q components of the stator voltage; Rs is the stator resistance, p is the number of pole pairs, φm is the flux. The complete nonlinear model of a PMSG-based WECS is obtained by combining (5)-(9).
Bx =
where γ =
δ˙x = Ax δx + Bx δu ,
(11)
1 3 ¯ ¯2 2Jr ω ¯ r ρπR CQ (λ)γ V
Ax = iKg +
0 iBg 3 ¯ ¯2 2Jr ω ¯ r ρπR CQ (λ)γ V 0 0
− ηJi r −Bg
1 Jg 1 i2 Jg + ηJr
0 0
0 0 0 s −R Ld pLd ω ¯ − Lq g
− Lpq 0
0 0 −Kg pLq ¯ Ld iqs (Ld¯ids − φm )
,
(13)
FOR
A. H∞ Optimal Control Consider a singularly perturbed linear time-invariant system under disturbances: x˙ 1 x˙ 2
=
A1 x1 + A2 x2 + B1 u + D1 w,
(14)
=
A3 x1 + A4 x2 + B2 u + D2 w,
(15)
where xi ∈ Rni , i = 1, 2 are slow and fast state variables, u ∈ Rr is the control input, w ∈ Rr is the system disturbance, is a small positive singular perturbation parameter, A1 to A4 , B1 and B2 , D1 and D2 are constant matrices with appropriate dimensions. The cost function associated with the system (14)(15) is J=
Z
1 2
∞ 0
xT Qx + uT Rr u dt,
(16)
where Q ≥ 0, Rr > 0 are weighting matrices. The optimal control that assures a µ level of optimality is u∗
= =
T −R−1 r B Pr x T B1 x1 −1 −Rr Pr , 1 B x2 2
(17)
where µ is a real positive parameter that represents an optimal disturbance attenuation level, which is defined as
inf sup u
w
( √
J ||w| |
)
< µ,
(18)
and Pr is a positive semidefinite stabilizing solution of the following algebraic Riccati equation:
m − pφ Jg Bg pφm , Jg pLq ω ¯g Ld Rs − Ld
¯ 0 (λ) ¯ λC Q ¯ . CQ (λ)
This section briefly presents the H∞ optimal control and filtering for linear time-invariant singularly perturbed systems. The focus is on the decomposition of an H∞ optimal controller and filter into completely independent, slow and fast H∞ controllers and filters. More details are referred to [20] and references therein.
¯ and δu = u − u ¯ are variations of variables where δx = x − x in the neighborhood of the operating point. The system and control matrices are given as
0 0 0 0 − L1q
III. H∞ O PTIMAL C ONTROL AND F ILTERING S INGULARLY P ERTURBED S YSTEMS
B. Linear Model The function CQ (λ) in (4) reaches its maximum value at the so-called optimal tip-speed ratio λ∗ . WECSs are expected to operate optimally such that the λ∗ is maintained regardless the wind speed fluctuations. As observed from (2), to keep λ at the λ∗ as V changes, ωr must be adjusted by controlling Tg . Note that V is a time-varying parameter, therefore the nonlinear model (5)-(10) is linearized for different V . The nonlinear (5)-(10) has the state vector x = [ωr ωg TH ids iqs ]T and the input vector u = [uds uqs ]T . ¯ = [¯ Choosing an operating point with x ωr ω ¯ g T¯H ¯ids ¯iqs ]T ¯ = [¯ and u uds u¯qs ]T at a certain wind speed V¯ , the linearized state space model is obtained as follows
0 0 0 − L1d 0
(12)
1 AT Pr + Pr A + Q − Pr S − 2 Z Pr = 0, µ where
(19)
=
S1 Z1
= =
A1 1 A3
1 S2 1 2 S3
A2 S , S = 1 1T 1 A S 4 2 1 Z1 Z Q1 Q2 2 , Q = (20) 1 T 1 1 T 1 Z Z Q Q 2 2 3 2 3 T −1 T −1 T B1 R−1 r B1 , S2 = B1 Rr B2 , S3 = B2 Rr B2 T T T D1 D1 , Z2 = D1 D2 , Z3 = D2 D2 .
A = Z
with PF 1 to PF 3 are component matrices of the positive semidefinite stabilizing solution of the H∞ filter algebraic Riccati equation: 1 APF + PF AT − PF C T V −1 C − 2 GT RF G PF µ +DW DT = 0, (29) where
An algorithm to find the solution Pr of the H∞ algebraic regulator Riccati equation (19) in terms of slow and fast solutions is given in [20]. B. H∞ Optimal Filtering Consider a singularly perturbed linear time-invariant system under disturbances given in (14)-(15) with the system output: y = C1 x1 + C2 x2 + v,
(21)
where y ∈ Rp is the system output, v ∈ Rp is the measurement disturbance, C1 and C2 are constant matrices with appropriate dimensions. An H∞ filter is designed to estimate the state which is given in terms of a linear combination of x1 and x2 by z = G1 x1 + G2 x2 ,
(22)
where G1 and G2 are arbitrary matrices with appropriate dimensions. The H∞ filtering problem is to find an estimate ˆz of z such that the energy gain from the disturbances to the estimation errors, z − ˆ z, is less than a predefined level µ2 as [20] sup J < µ2 ,
(23)
w,v
where J is the performance measure for the infinite horizon estimation which is defined as
J = R∞ 0
R∞ 0
||z − ˆ z| |2RF dt
, ||w| |2W −1 + ||v| |2V −1 dt
(24)
with RF ≥ 0, W > 0, V > 0 are weighting matrices with appropriate dimensions. The H∞ optimal closed-loop filter driven by the control input and system output is given in [20] as ˆ˙ 1 x ˆ˙ 2 x
= =
ˆ 1 + (A2 − K1 C2 ) x ˆ2 (A1 − K1 C1 ) x +B1 u + K1 y, ˆ 1 + (A4 − K2 C2 ) x ˆ2 (A3 − K2 C1 ) x +B2 u + K2 y,
where filter gains K1 , K2 are given as K1 = PF 1 C1T + PF 2 C2T V −1 , K2 = PFT2 C1T + PF 3 C2T V −1 ,
(25) (26)
(27) (28)
PF
=
C
=
D1 PF 1 PF 2 , D = 1 PFT2 1 PF 3 D2 [C1 C2 ] , G = [G1 G2 ]
(30)
Note that the solution of the H∞ filter algebraic Riccati equation (29) is also obtained in terms of slow and fast solutions as in the H∞ regulator algebraic Riccati equation (19) using the dual property between the regulator and filter [20]. C. Decomposition of the H∞ Optimal Filter and Controller In order to decouple the H∞ optimal closed-loop filter given in (25)-(26) into two independent, reduced-order filters, a similarity transformation, TF , is proposed to transform the old state coordinate into a new state coordinate [20], which is x1 x1 xs xs = TF , or = T−1 , (31) F x2 x2 xf xf where TF =
I − HL −H L I
,
(32)
with I representing for identity matrices with appropriate dimensions, L and H are the solutions of the following Chang decoupling equations: (A4 − K2 C2 ) L − (A3 − K2 C1 )
−L [(A1 − K1 C1 ) − (A2 − K1 C2 ) L] = 0, (33) −H (A4 − K2 C2 ) + (A2 − K1 C2 ) − HL (A2 − K1 C2 ) + [(A1 − K1 C1 ) − (A2 − K1 C2 ) L] H = 0. (34) It should be noted that solutions L and H can be found by the numerical Newton algorithm given in [20]. Applying the state coordinate transformation (31) to the H∞ optimal closed-loop filter given in (25)-(26) yields independent, reduced-order, optimal slow and fast filters, which are ˆ˙ s ˆs x a1 a2 x −1 = T T F 1 F ˆf a3 1 a4 x ˆ˙ f x B1 K1 +TF 1 u + TF 1 y, B2 K2 ˆs As 0 x = ˆf 0 1 Af x Bs Ks + 1 u+ 1 y, (35) B K f f
where a1 = A1 − K1 C1 , a2 = A2 − K1 C2 a3 = A3 − K2 C1 , a4 = A4 − K2 C2 , As = a1 − a2 L, Af = a4 + La2 Bs = B1 − HLB1 − HB2 , Bf = B2 + LB1 Ks = K1 − HLK1 − HK2 , Bf = K2 + LK1 .
(36)
Similarly, applying the similarity transformation (31) to the H∞ optimal closed-loop control in (17) produces two independent, reduced-order slow and fast controls corresponding to two independent, reduced-order, slow and fast filters in (35), which are T ˆs B1 x −1 u∗ = − R−1 P T r F 1 r ˆf x B2 {z } | [Fs Ff ]
ˆ s + Ff x ˆf ) . = − (Fs x
(37)
The implementation diagram of the reduced-order H∞ optimal filtering and control is shown in Fig. 4 where two independent, reduced-order filters are used to estimate the slow ˆ s , and fast state vector, x ˆ f , based on the system state vector, x output, y, and control input, u. The final control is a composite of two reduced-order H∞ optimal controls, us and uf , as indicated in (37). IV. S IMULATION The proposed method is validated by numerical simulations on MATLAB R and SIMULINK R environments with system parameters given as: ρ = 1.25 kg/m3 , R = 2.5 m, i = 6, η = 1, Jr = 2.88 kgm2 , Jg = 0.22 kgm2 , Kg = 75 N m/rad, Bg = 0.3 kgm2 /s, Ld = Lq = 0.04156 H, p = 3, Rs = 3.3 Ω, φm = 0.4382 W b, k0 = 0.0061, k1 = −0.0013, k2 = 0.0081, k3 = −9.7477 × 10−4 , k4 = −6.5416 × 10−5 , k5 = 1.3027 × 10−5 , k6 = −4.54 × 10−7 , λ∗ = 7, = 0.04156. The nonlinear PMSG-based WECS is linearized at the wind speed V¯ = 10 m/s and the ¯ = [24.54 147.24 29.07 4.45 3.50]T where operating point x the maximum power conversion is achieved at the chosen wind speed. System matrices are obtained and chosen as follows −1.0133 0 −2.0833 A1 = 0 0 4.5455 , 448.176 −75 −5.1136
Fig. 4.
Reduced-order H∞ Optimal Filtering and Control
0 0 0 A2 = 0 −5.9755 , A3 = 0 0 1.7926 A4 =
−3.3 18.3583 −18.3583 −3.3
B2 = D2 =
C2 =
0 0
−1 0
, G1 =
0 −1
0.4370 0 0.7594 0
,
0 0 , B1 = D1 = 0 0 , 0 0
, C1 =
0 0.1 0
0
1 0
, G2 =
,
0 0
.
Weighting matrices are chosen by trial and error as Rr = 0.25I[2×2], RF = I[1×1] , W = 0.5I[2×2], V = I[1×1] , and the optimality level is specified as µ = 1. In the simulations, disturbances are assumed as random noises and their upper ¯ = 2, v ¯ = 1. Matrices of the H∞ optimal, bounds are w reduced-order slow and fast filters given in (35) and controllers given in (37) are obtained as
−1.0133 −0.0098 −2.0833 0 0.0345 4.5426 , as = 448.1760 −75.0219 −5.1128
0.0004 0.0003 0.0098 Bs = −0.3159 0.0560 , Ks = 0.0602 , 0.0771 −0.0670 −0.0066 af =
−3.3000 18.3693 −18.3583 −3.3040
Kf = 10−3 ×
Fs =
−0.9922 −0.1818
, Bf =
, Ff =
−1 0
0 −1
,
−0.5683 0.0003 −0.0001 −0.5689
0.6012 −1.5426 0.0187 −0.1048 0.3013 −0.0228
,
.
The comparisons between the full-order H∞ optimal controller (without using the time scale method) and reducedorder H∞ optimal controller are shown in Fig. 5 and Fig. 6, indicating that the reduced-order responses are really close to the full-order responses. Thus, this closeness verifies the efficiency of the proposed reduced-order control. The reducedorder filters’ performance is also shown in Fig. 7 which exhibits all states including slow and fast states are filtered out from the process and measurement disturbances. Additionally, the closeness between the full-order and reduced-order filters with respect to the generator speed is presented in Fig. 8.
V. C ONCLUSION
12 Full−order Reduced−order
Generator Speed (rad/s)
10
8
6
4
2
0
−2
0
2
4
6
8
10
In this paper, a reduced-order H∞ optimal control is proposed to regulate the maximum power conversion of a PMSGbased WECS. The key idea is based on the decoupling of slow and fast dynamics using the singular perturbation theory. As a result, the controllers are not only robust to disturbances, but also are of low order, which requires less computational effort compared to full-order controllers. The proposed control design also allows parallel processing at different time scales of different dynamics existing in WECSs.
Time (s)
R EFERENCES Fig. 5. System Output Responses of the Full-order and Reduced-order H∞ Optimal Control. 18 Full−order Reduced−order
16 14
d and q voltages (V)
12 10 8 6 4 2 0 −2 0
2
4
6
8
10
Time (s)
Fig. 6. Control Inputs of the Full-order and Reduced-order H∞ Optimal Control.
Wind rotor speed, wr Generator speed, wg Internal torque, TH d−component current, id q−component current, iq
40
Slow and Fast States
30 20 10 0 −10 −20 −30 0
2
4
6
8
10
Time (s)
Fig. 7.
Reduced-order Filtered States
11 Full−order Reduced−order
10
Generator speed, wg (rad/s)
9 8 7 6 5 4 3 2 1 0
0
2
4
6
8
Time (s)
Fig. 8.
Filtered Generator Speed
10
[1] M. Stiebler, Wind Energy Systems for Electric Power Generation. Springer-Verlag, 2008. [2] P. Musgrove, Wind Power. Cambridge University Press, 2010. [3] W. Shepherd and L. Zhang, Electricity Generation Using Wind Power. Singapore: World Scientific Publishing Co.Pte.Ltd, 2011. [4] H. Wagner and J. Mathur, Introduction to Wind Energy Systems: Basics, Technology and Operation. Springer-Verlag, Berlin Heidelberg, 2009. [5] D. Spera, Wind Turbine Technology: Fundamental Concepts of Wind Turbine Engineering, 2nd Edition. ASME Press, 2009. [6] F. Bianchi, H. Battista, and R. Mantz, Wind Turbine Control System: Principles, Modelling, and Gain Scheduling Design. Springer-Verlag, London, UK, 2007. [7] I. Munteanu, A. Bratcu, N. Cutuluslis, and E. Ceanga, Optimal Control of Wind Energy Systems: Toward a Global Approach. Springer, 2008. [8] D. Naidu, Optimal Control Systems. Chapman and CRC Press, 2003. [9] J. Li, H. Xu, L. Zhang, Zhuying, and S. Hu, “Disturbance accommodating LQR method based pitch control strategy for wind turbines,” in Intelligent Information Technology Application, 2008. IITA ’08. Second International Symposium on, vol. 1, Dec. 2008, pp. 766 –770. [10] F. Wu, X. Zhang, P.J., and M. Sterling, “Decentralized nonlinear control of wind turbine with doubly fed induction generator,” Power Systems, IEEE Transactions on, vol. 23, no. 2, pp. 613 –621, May 2008. [11] E. Muhando, T. Senjyu, N. Urasaki, H. Kinjo, and T. Funabashi, “Online WTG dynamic performance and transient stability enhancement by evolutionary LQG,” in Power Engineering Society General Meeting, 2007. IEEE, June 2007, pp. 1 –8. [12] S. Nourdine, H. Camblong, I. Vechiu, and G. Tapia, “Comparison of wind turbine LQG controllers designed to alleviate fatigue loads,” in Control and Automation (ICCA), 2010 8th IEEE International Conference on, June 2010, pp. 1502 –1507. [13] F. Lescher, J. Zhao, and A. Martinez, “LQG multiple model control of a variable speed pitch regulated wind turbine,” in 17th IMACS World Congress, Paris, France, 2005. [14] H. Nguyen and D. Naidu, “Advanced control strategies for wind energy systems: An overview,” in Power Systems Conference and Exposition (PSCE), 2011 IEEE/PES, March 2011, pp. 1 –8. [15] P. Kokotovic, H. K. Khali, and J. O’Reilly, Singular Perturbation Methods in Control: Analysis and Design, 1st ed. SIAM, 1986. [16] D. Naidu, Singular Perturbation Methodology in Control Systems, 1st ed. Peter Peregrinus Ltd., 1988. [17] H. Nguyen and D. Naidu, “Time scale analysis and control of wind energy conversion systems,” in Resilient Control Systems (ISRCS), 2012 5th International Symposium on, Aug. 2012, pp. 149 –154. [18] P. Moriarty and S. Butterfield, “Wind turbine modeling overview for control engineers,” in American Control Conference, 2009. ACC ’09., Jun. 2009, pp. 2090 –2095. [19] H. Li and Z. Chen, “Overview of different wind generator systems and their comparisons,” Renewable Power Generation, IET, vol. 2, no. 2, pp. 123 –138, Jun. 2008. [20] Z. Gajic and M. Lim, Optimal control of singularly perturbed linear systems and applications. 270 Madison Avenue, New York, NY 10016: Marcel Dekker, Inc., 2001.