自动化学报  2017, Vol. 43 Issue (8): 1434-1442   PDF    
基于二阶滑模控制的变速双馈风力发电系统最大功率点跟踪
刘向杰1, 王诚诚1, 韩耀振1,2     
1. 华北电力大学新能源电力系统国家重点实验室 北京 102206;
2. 山东交通学院信息科学与电气工程学院 济南 250357
摘要: 本文提出一种超螺旋二阶滑模控制方案同时实现双馈变速风力发电系统最大风能捕获和无功功率调节.通过设计两个二阶滑模控制器,实现控制目标,降低机械磨损,提高控制精度,通过调节发电机转子电压,跟踪风机最优转速和转子电流设定值,实现额定风速以下的最大风能捕获和无功功率调节.采用二次型李雅普诺夫函数确定控制参数范围、确保系统有限时间稳定性.1.5 MW风机系统仿真实验验证所提方案有效性.
关键词: 双馈异步发电机     最大功率点跟踪     二阶滑模控制     变速风电机组    
Second-order Sliding Mode Control of DFIG Based Variable Speed Wind Turbine for Maximum Power Point Tracking
Xiangjie Liu1, Chengcheng Wang1, Yaozhen Han1,2     
1. State Key Laboratory of Alternate Electrical Power System with Renewable Energy Sources, North China Electric Power University, Beijing 102206, China;
2. School of Information Science and Electrical Engineering, Shandong Jiaotong University, Jinan 250357, China
Abstract: This paper proposes a super-twisting second order sliding mode control scheme to maximize the wind energy capture of a doubly fed induction generator based variable speed wind turbine (VSWT) system, and minimize the reactive power simultaneously. Two second order sliding mode controllers are designed to achieve the control objectives, reduce mechanical stress and improve control accuracy. By regulating the generator rotor voltage, one controller makes the wind turbine rotor speed track the optimal speed, which can maximize power generation. The other maintains the rotor current at rated value to minimize the reactive power. A quadratic form Lyapunov function is adopted to determine the range of controller parameters and guarantee the finite time stability. Simulation results on a 1.5MW doubly fed induction generator (DFIG)-based variable speed wind turbine demonstrate the validity of the proposed control strategy.
Key words: Doubly fed induction generator (DFIG)     maximum power point tracking (MPPT)     second order sliding mode (SOSM) control     variable speed wind turbine (VSWT)    
1 Introduction

Recent years wind energy conversion systems have advanced rapidly, due to the increasing concern for the environment pollution caused by the traditional energy sources. Through tracking the changing wind speed by adjusting the shaft speed, the variable speed wind turbines (VSWTs) enable the turbine to operate at its maximum power coefficient over a wide range of wind speeds [1]. To achieve the major task of power efficiency maximization, the VSWT should track the maximum power point despite wind variations [2].

Originally, the VSWT control system is designed based on the linear system theory [3]-[5]. The linear control encountered difficulties since VSWT system is a complex and highly nonlinear system with strong coupling features and uncertainty in both the aerodynamics and the electrical parts. To account for the nonlinear behavior, and to deal with the recognized problem of parameter variations, various advanced control strategies have been proposed, such as neural networks [6], adaptive control [7], feedback linearization technique [8], predictive control [9] and sliding mode approach [10]-[11]. Among these control strategies, sliding mode control (SMC) has proven to be very robust with respect to system parameter variations and external disturbances, and thus suitable for realizing VSWT system control. However, the standard first-order sliding mode (FOSM) control generally shows significant drawback of chattering phenomenon. Some key components in the wind turbine (WT), such as gear box, would get damaged by the abrupt commutation in forces and torques [12].

Instead of influencing the first order time derivative of sliding variable, the second-order sliding mode (SOSM) control acts on the second order derivatives of the sliding surfaces, which can help reduce the chattering phenomenon and provide higher sliding precision. Particularly, the SOSM super-twisting algorithm has been an effective way of controlling the VSWT system, since it only requires measurement of the sliding variable without using information about derivatives of the sliding constraint. A number of contributions on this control strategy have recently appeared for the realization of maximum power point tracking (MPPT) of VSWT system. Valenciaga et al.[13] obtained the control objective of wind energy conversion optimization of brushless doubly fed reluctance machine (BDFRM) based VSWT. Beltran et al. [14] presented a SOSM controller that facilitates the generator torque tracking the optimal torque, and maximized the power extraction of doubly fed induction generator (DFIG) based VSWT. Evangelista et al. [15] synthesized a super-twisting sliding mode control with variable gain, and showed a better performance in terms of mechanical loads and power tracking. However, in the aforementioned literatures, the SOSM control mainly concentrated on the electrical side dynamic control. For achieving the control objective of maximum power extraction, the generator torque was set to track the optimal aerodynamic torque, and it means that the energy loss caused by friction and mechanical rotational inertia in the transmission system was neglected in such circumstance. In other words, the captured power from wind was assumed to be equal to generator output power. This assumption is not always consistent with the actual situation, especially in turbulent wind scenarios.

Therefore, this paper aims to investigate a new SOSM control method for the DFIG based VSWT system. A complete dynamic model of the VSWT system is established by integrating both aerodynamic and electrical parts, together with parameter uncertainties and perturbations. The control objective is to make the wind turbine rotor speed track the desired speed (the speed that is given by a MPPT) in spite of system uncertainties, and maintain the rotor current at rated value which can minimize the reactive power. This control strategy can provide a faster control action as the wind speed variations can be reflected instantaneously and significantly on the rotor speed reference signal, and yield more energy. Besides, by adopting quadratic Lyapunov function [16], the range of control parameters is obtained to guarantee the finite time convergence of the SOSM control system.

The article is organized as follows. The modeling of the VSWT system is presented in Section 2. Section 3 explains the problem formulation and the detailed control strategy. The comparative simulation results for the 1.5 MW three-blade DFIG-based WT system are investigated in Section 4. Finally, the conclusion is stated in Section 5.

2 VSWT System Modeling

The VSWT system is mainly composed of a turbine, a gearbox, and a generator, which are combined to convert wind energy into electric energy. According to different wind speeds, the VSWT system works at four different operating regions, as shown in Fig. 1. While the wind speed is below the cut-in speed $ v_{\rm cut\_in} $, the produced torque is insufficient to drive the turbine, thus it cannot realize the startup function. While the wind speed increases over the cut-in speed $ v_{\rm cut\_in} $ and below the rated speed $ v_{\rm rated} $, the VSWT system tries to capture the maximum wind energy by adjusting the rotor speed of WT. This region is called the partial load region. With the increasing wind speed satisfying $ v_{\rm rated} <v<v_{\rm cut\_off} $, the WT tries to maintain the rated output power by regulating the pitch angle, called the full load region. There is an additional region, called the shutdown region, where the turbine stops operation due to high wind speed.

Figure 1 WT operation regions.

Obviously, the control task in the partial load region is more challenging, since the dynamics in this region is quite complex. Thus the study in this paper concentrates on the modeling and control in this region. The mathematical model in this region consists of formulas accounting for the electrical dynamics of the rotor and the dynamics of the mechanical rotational speed, including also the damping and resistance uncertainties caused by long-term mechanical wear and ambient temperature changes [8], [14]. It is expressed as:

$ \begin{align} \label{eq1}\begin{cases} \dot {\omega }_w =\frac{T_w }{J}-\frac{K+\Delta K}{J}\omega _w +\frac{3L_m \phi _s n_p n_g }{2L_s J}I_{rq} \\[2mm] \dot {I}_{rd} =\frac{L_s (R_r +\Delta R_r )}{L_m ^2-L_r L_s }I_{rd} +\omega _1 I_{rq} -n_p n_g \omega _w I_{rq} \\[1mm] \qquad\ \ +~\frac{L_s }{L_r L_s -L_m ^2}U_{rd} \\[2mm] \dot {I}_{rq} =\frac{L_m \phi _s n_p n_g }{L_r L_s -L_m ^2}\omega _w -\omega _1 I_{rd} +\frac{L_s (R_r +\Delta R_r )}{L_m ^2-L_r L_s }I_{rq} \\[1mm] \qquad\ \ +~n_p n_g \omega _w I_{rd} +\frac{L_s }{L_r L_s -L_m ^2}U_{rq} -\frac{L_m \phi _s \omega _1 }{L_r L_s -L_m ^2} \end{cases} \end{align} $ (1)

where $ \Delta K $ and $ \Delta R_r $ represent the parameter uncertainties of the damping and resistance respectively, and the constraints $ \left| {\Delta K} \right|\le 0.5K $ and $ \left| {\Delta R_r } \right|\le 0.5R_r $ are satisfied. The other symbols in (1) are described as:

$ \omega _w $ wind turbine rotor speed;
$ I_{rd} $ the $d$-axis component of rotor current;
$ I_{rq} $ the $q$-axis component of rotor current;
$ U_{rd} $ the $d$-axis component of the rotor voltage;
$ U_{rq} $ the $q$-axis component of the rotor voltage;
$ J $ the inertia of the combined rotating parts;
$ K $ turbine total external damping;
$ n_g $ gearbox ratio;
$ \phi _s $ stator flux;
$ L_m $ mutual inductance;
$ L_s $ stator leakage inductance;
$ L_r $ rotor leakage inductance;
$ R_r $ rotor resistance;
$ \omega _1 $ synchronous speed;
$ n_p $pole pair number.

In order to concisely express (1), define

$ \begin{align*}&k_1 =\frac{1}{J}, \hspace{20mm} k_2 =\frac{K}{J} \\&\Delta k_2 ={\Delta K}{J}, \hspace{12.6mm}k_3 =\frac{3L_m \phi _s n_p n_g }{2L_s J}\\& k_4 =\frac{L_s R_r }{L_m ^2-L_r L_s }, \hspace{6.1mm}\Delta k_4 =\frac{\Delta R_r L_s }{L_m ^2-L_r L_s } \\ &k_5 =\omega _1, \hspace{20mm} k_6 =n_p n_g\\ & k_7 =\frac{L_s }{L_r L_s -L_m ^2}, \hspace{6.5mm} k_8 =\frac{L_m \phi _s n_p n_g }{L_r L_s -L_m ^2}\\& k_9 =\frac{L_m \phi _s \omega _1 }{L_r L_s -L_m ^2} \end{align*} $

this leads to

$ \begin{align} \label{eq2} \begin{cases} \dot {\omega }_w =k_1 T_w -(k_2 +\Delta k_2 )\omega _w +k_3 I_{rq} \\ \dot {I}_{rd} =(k_4 +\Delta k_4 )I_{rd} +k_5 I_{rq} -k_6 \omega _w I_{rq} +k_7 U_{rd} \\ \dot {I}_{rq} =(k_4 +\Delta k_4 )I_{rq} -k_5 I_{rd} +k_6 \omega _w I_{rd} +k_7 U_{rq}\\ \qquad\ \ +~k_8 \omega _w -k_9. \end{cases} \end{align} $ (2)

The VSWT system model (2) is composed of two parts, e.g., the turbine model and the DFIG model. The single mass model of the turbine is [8]:

$ \begin{align} \label{eq3} J\dot {\omega }_w =T_w -K\omega _w -n_g T_{em} \end{align} $ (3)

where $ T_w $ is the aerodynamic torque and $ T_{em} $ is the generator electromagnetic torque.

The nonlinear characteristics of $ T_w $ is established based on the aerodynamic principles [14], e.g.,

$ \begin{align} \label{eq4} T_w =\frac{P_w }{\omega _w }=\frac{1}{2}\rho \pi R^3\frac{C_p (\lambda, \beta )}{\lambda }v^2 \end{align} $ (4)

where $ P_w $ is the aerodynamic power and $ C_p (\lambda, \beta ) $ is the power coefficient relating to a nonlinear function affected by the blade pitch angle $ \beta $ and the tip-speed ratio $ \lambda $, defined as

$ \begin{align} \label{eq5} \lambda =\frac{\omega _w R}{v}. \end{align} $ (5)

The DFIG model in the synchronously rotating frame $dq$ is established as (6) based on the stator flux orientation, whose phase diagram is shown in Fig. 2.

Figure 2 Phase diagram of stator flux orientation.
$ \begin{align} \label{eq6} \begin{cases} U_{sd} =0 \\ U_{sq} =\omega _1 \phi _{sd} =U_s \\ U_{rd} =R_r I_{rd} +\left( {L_r -\frac{L_m ^2}{L_s }} \right)\dot {I}_{rd} -\left( {L_r -\frac{L_m ^2}{L_s }} \right)I_{rq} \omega _s \\ U_{rq} =R_r I_{rq} +\left( {L_r -\frac{L_m ^2}{L_s }} \right)\dot {I}_{rq} +\left( {L_r -\frac{L_m ^2}{L_s }} \right)I_{rd} \omega _s \\\qquad\ \ +~\frac{L_m \phi _S }{L_s }\omega _s \\ T_{em} =-\frac{3}{2}n_p \frac{L_m \phi _S }{L_s }I_{rq}. \end{cases} \end{align} $ (6)
3 The VSWT System Control Strategy 3.1 Problem Formulation

Considering (4) and (5) in Section 2, the aerodynamic power that can be captured by a wind turbine is:

$ \begin{align} \label{eq7} P_w =\frac{1}{2}\rho \pi R^2C_p (\lambda, \beta )v^3. \end{align} $ (7)

As the main control objective of the above VSWT system in the partial load region is to track the maximum power point and harvest more wind energy, the power conversion coefficient $ C_p (\lambda, \beta ) $ must reach the maximum value $ C_{p\max } $ in various wind speeds.

$ C_p (\lambda, \beta ) $ is generally expressed as [17]:

$ \begin{align} \label{eq8} C_p \left(\lambda, \beta \right)=c_1 \left(\frac{c_2 }{\Lambda }-c_3 \beta -c_4 \right)\times e^{\frac{-c_5 }{\Lambda }}+c_6 \lambda \end{align} $ (8)

with ${1}/{\Lambda }={1}/{(\lambda +0.08\beta )}-{0.035}/{(\beta ^3+1)} $. The coefficients $ c_1 -c_6 $ depend on the WT design characteristics. And the following coefficient values are adopted in this work: $ c_1$ $=$ $0.5176 $, $ c_2$ $=$ $116 $, $ c_3 =0.4 $, $ c_4 =5 $, $ c_5 =21 $, $ c_6=0.0068 $.

Actually, in the partial load region, the pitch angle $ \beta $ is fixed at zero, to achieve the maximum power conversion coefficient $ C_{p\max } $ [2]. Then the $ C_p -\lambda $ characteristic with a constant pitch angle $ \beta =0^\circ $ is shown in Fig. 3. The tip-speed ratio $ \lambda $ must be maintained at an optimal point for achieving the maximum power coefficient $ C_{p\max } $.

Figure 3 $ C_p-\lambda $ characteristic of wind turbine.

Considering Fig. 3 and the definition of $ \lambda $ in (5), it is thus necessary to adjust the rotor speed $ \omega _w $ according to the wind speed variation in the partial load region, to maintain the tip-speed ratio $ \lambda $ at an optimal point $ \lambda _{\rm opt} $. This could guarantee the maximum power coefficient $ C_{p\max } $, and from (7), when $ C_p $ is controlled at the maximum value, maximum output power $ P_{a\max } $ is extracted from wind. Thus, the optimal WT rotor speed $ \omega _{\rm opt} $ is given by:

$ \begin{align} \label{eq9} \omega _{\rm opt} =\frac{\lambda _{\rm opt} v}{R}. \end{align} $ (9)

In the DFIG based VSWT system (2), this maximum wind energy capture objective can be achieved by means of the rotor voltage regulation in the generator. The rotor voltage can also control the rotor current $ I_{rd} $, and thus regulate the stator-side reactive power $ Q_s $, which can be expressed as:

$ \begin{align} \label{eq10} Q_s =\frac{3}{2}(U_{sd} I_{sq} +U_{sq} I_{sd} )=\frac{3U_s (\phi _s -L_m I_{rd} )}{2L_s }. \end{align} $ (10)

Thus the designed SOSM control system should accomplish two major objectives. One is to maximize power extraction, by controlling the rotor speed to track the optimal rotor speed $ \omega _{\rm opt} $, which is regulated by $q$-axis component of generator rotor voltage. The other is to minimize the reactive power by controlling $ I_{rd} $ to follow an external reference, which is regulated by $d$-axis component of generator rotor voltage.

Here, $ Q_s $ is set to zero to ensure a unity power factor operation of the studied VSWT system which renders the rotor reference current $ I_{rdr} $ as

$ \begin{align} \label{eq11} I_{rdr} =\frac{U_s }{L_m \omega _1 }. \end{align} $ (11)
3.2 SOSM Controller Design

In order to achieve the control objectives in the partial load region, the nonlinear SMC scheme is presented in this section. The schematic diagram of the DFIG-based WT system is shown in Fig. 4.

Figure 4 The schematic diagram of the DFIG-based WT system.

In designing a general SMC for the VSWT system, the tracking error between the actual rotor speed and the desired value $ \omega _{\rm opt} $, and also the error between the actual and the desired $d$-axis component of rotor current, are defined as:

$ \begin{align} \label{eq12} \begin{cases} e_1 =\omega _w -\omega _{\rm opt} \\ e_2 =I_{rd} -I_{rdr}. \end{cases} \end{align} $ (12)

Then the sliding variables are defined as follows:

$ \begin{align} \label{eq13} \sigma _1 &=ce_1 +\dfrac{{d}}{{d}t}e_1\notag \\ &=\dot {\omega }_w +c\omega _w -c\omega _{\rm opt} -\dot {\omega }_{\rm opt} \end{align} $ (13)

where $ c>0 $.

$ \begin{align} \label{eq14} \sigma _2 =e_2 =I_{rd} -I_{rdr}. \end{align} $ (14)

Respecting the VSWT system model (2), the first-order derivatives of sliding variables are

$ \begin{align} \dot {\sigma }_1 =&\ \ddot {e}_1 +c\dot {e}_1 \notag\\ =&\ k_1 \dot {T}_w +(c-k_2 -\Delta k_2 )\dot {\omega }_w +k_3 k_7 U_{rq}\notag\\ &-\ddot {\omega }_{\rm opt} -c\dot {\omega }_{\rm opt} +k_3 ( (k_4 +\Delta k_4 )I_{rq}\notag \\ & -k_5 I_{rd} +k_6 \omega _w I_{rd} +k_8 \omega _w -k_9 ) \notag\\ =&\ G_1 +k_3 k_7 U_{rq} \end{align} $ (15)
$ \begin{align} \dot {\sigma }_2 =&\ \dot {e}_2 =\dot {I}_{rd} -\dot {I}_{rdr} \notag\\ =&\ (k_4 +\Delta k_4 )I_{rd} -\dot {I}_{rdr} +k_5 I_{rq} -k_6 \omega _w I_{rq} +k_7 U_{rd}\notag \\ =&\ G_2 +k_7 U_{rd} \end{align} $ (16)

where $ G_1 $ and $ G_2 $ contains the state variables, VSWT parameters and parameter perturbations:

$ \begin{align} G_1 =&\ k_1 \dot {T}_w +(c-k_2 -\Delta k_2 )\dot {\omega }_w -\ddot {\omega }_{\rm opt} -c\dot {\omega }_{\rm opt} \notag \\ & +k_3 \left( {(k_4 +\Delta k_4 )I_{rq} -k_5 I_{rd} +k_6 \omega _w I_{rd} +k_8 \omega _w -k_9 } \right) \end{align} $ (17)
$ \begin{align} \label{eq18} G_2 =&\ (k_4 +\Delta k_4 )I_{rd} -\dot {I}_{rdr} +k_5 I_{rq} -k_6 \omega _w I_{rq}. \end{align} $ (18)

In (13) and (15), since $ \frac{\partial \sigma _1 }{\partial U_{rq} }=0 $ and $ \frac{\partial \dot {\sigma }_1 }{\partial U_{rq} }\ne 0 $, the relative degree with respect to $ \sigma _1 $ equals to $ 1 $. Similarly, from (14) and (16), the relative degree with respect to $ \sigma _2 $ is $ 1 $ as well.

Using the standard SMC with the approaching law method [18], the so-called exponential approaching law is selected as

$ \begin{align} \label{eq19} \dot {\sigma }=-\varepsilon {\rm sgn}(\sigma )-\delta \sigma \end{align} $ (19)

where $ \sigma $ is the sliding variable, $ \varepsilon $ and $ \delta $ are the designed parameters.

Then the $dq$-axis components of rotor voltages are deduced as

$ \begin{align} \label{eq20} &U_{rq} =\frac{1}{k_3 k_7 }(-\varepsilon _1 {\rm sgn}(\sigma _1 )-\delta _1 \cdot \sigma _1 -G_1 ) \end{align} $ (20)
$ \begin{align} \label{eq21} &U_{rd} =\frac{1}{k_7 }(-\varepsilon _2 {\rm sgn}(\sigma _2 )-\delta _2 \cdot \sigma _2 -G_2 ). \end{align} $ (21)

Since the switching term $ \varepsilon {\rm sgn}(\sigma ) $ acts on the controller directly, it makes the control inputs $ U_{rd} $ and $ U_{rq} $ discontinuous and leads to intensive chattering effect.

Using the super-twisting algorithm [19], the control inputs $ U_{rd} $ and $ U_{rq} $ are reconstructed as

$ \begin{align} \label{eq22} \begin{cases} U_{rq} =u_1 +u_2 \\ u_1 =-\gamma _1 \vert \sigma _1 \vert ^\frac{1}{2}{\rm sgn}(\sigma _1 ) \\ \dot {u}_2 =-\phi _1 {\rm sgn}(\sigma _1 ) \end{cases} \end{align} $ (22)
$ \begin{align} \begin{cases} U_{rd} =u_3 +u_4 \\ u_3 =-\gamma _2 \vert \sigma _2 \vert ^\frac{1}{2}{\rm sgn}(\sigma _2 ) \\ \dot {u}_4 =-\phi _2 {\rm sgn}(\sigma _2 ) \end{cases} \end{align} $ (23)

where $ \gamma _1 $, $\phi _1 $, $\gamma _2 $, $\phi _2 $ are undetermined control parameters.

Notice that the discontinuity of super-twisting SOSM control appears only in its derivative term, such that the control inputs $ U_{rd} $ and $ U_{rq} $ in system (2) are actually continuous. Namely, the SOSM reduces the chattering effect while preserving all the advantages of standard SMC.

By using this SOSM scheme, it is very important to choose suitable control parameters for stable operation of WT system. Here a Lyapunov method is employed to determine the range of control parameters, and then guarantee stable operation of the WT system.

In order to ensure $ \sigma _1 \to 0 $, $ \dot {\sigma }_1 \to 0 $ in finite time, which means achieving a precise tracking of optimal rotor speed, the ranges of $ \gamma _1 $ and $ \phi _1 $ are needed to be addressed. Firstly, a state transformation is introduced as

$ \begin{align} \label{eq24} y=H_1 -k_3 k_7 \phi _1 \int_0^t {{\rm sgn}(\sigma _1 )} d\tau \end{align} $ (24)

and taking into consideration (22), then, (15) is deduced as

$ \begin{align} \label{eq25} \begin{cases} \dot {\sigma }_1 =-k_3 k_7 \gamma _1 \vert \sigma _1 \vert ^\frac{1}{2}{\rm sgn}(\sigma _1 )+y \\[1mm] \dot {y}=-k_3 k_7 \phi _1 {\rm sgn}(\sigma _1 )+\dot {G}_1. \end{cases} \end{align} $ (25)

As shown in (17), the perturbation $ G_1 $ is differentiable, and its first-order derivative is uniformly bounded, such that it satisfies

$ \begin{align} \label{eq26} \left| {\dot {G}_1 } \right|\le \Pi _1 {\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}\forall t>0 \end{align} $ (26)

where $ \Pi _1 $ is a positive constant.

Now, choose the quadratic Lyapunov function

$ \begin{align} \label{eq27} V(\sigma, y)=\zeta ^TP\zeta \end{align} $ (27)

where $ \zeta = [{\zeta _1 {\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}\zeta _2 }]^T=[{\vert \sigma _1 \vert ^{1/2}{\rm sgn}(\sigma _1 ){\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}y}]^T $, and $ P $ is a particular symmetric and positive definite matrix chosen as

$ P=\left[ {\begin{array}{l} 2k_3 k_7 \phi _1 +\frac{1}{2}k_3 ^2k_7 ^2\gamma _1 ^2{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}-\frac{1}{2}k_3 k_7 \gamma _1 \\[1mm] {\kern 1pt}-\frac{1}{2}k_3 k_7 \gamma {\kern 1pt}_1 {\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}1 \\ \end{array}} \right].{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt} $

From $ \frac{{ d}\left| x \right|}{{d}t}=\dot {x}{\rm sgn}(x) $ and (25), $ \dot {\zeta } $ is deduced as

$ \begin{align} \label{eq28} \dot {\zeta }=\frac{1}{\left| {\zeta _1 } \right|}\left(A\zeta +B\hat {\dot {G}}_1 \right) \end{align} $ (28)

where the matrix $ A=\left[{\begin{array}{l} -\frac{1}{2}k_3 k_7 \gamma _1 {\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}\frac{1}{2} \\ -k_3 k_7 \phi _1 {\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}0 \\ \end{array}} \right] $, $ B=\left[{0{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}1} \right]^T $, $ \hat {\dot {G}}_1$ $=$ $\left| {\zeta _1 } \right|\dot {G}_1 $.

Setting the matrix $ C=\left[{1{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}0} \right] $, the derivative of the Lyapunov function along the trajectories of the perturbed system (25) is deduced as

$ \begin{align} \dot {V}(\sigma, y)=&\ 2\dot {\zeta }^TP\zeta \notag\\ =&\ \frac{1}{\left| {\zeta _1 } \right|}(2\zeta ^TA^T+2\hat {\dot {G}}_1 B^T)P\zeta\notag \\ \le&\ \frac{1}{\left| {\zeta _1 } \right|}(2\zeta ^TA^TP\zeta +2\hat {\dot {G}}_1 B^TP\zeta +\Pi _1 ^2\left| {\zeta _1 } \right|-\hat {\dot {G}}_1 ^2)\notag \\ =&\ \frac{1}{\left| {\zeta _1 } \right|}(2\zeta ^TA^TP\zeta +2\hat {\dot {G}}_1 B^TP\zeta +\Pi _1 ^2\zeta ^TC^TC\zeta\notag\\ & -\hat {\dot {G}}_1 ^2)\notag \\ \le&\ \frac{1}{\left| {\zeta _1 } \right|}(\zeta ^TA^TP\zeta +\zeta ^TPA\zeta +\Pi _1 ^2\zeta ^TC^TC\zeta\notag \\ &+\zeta ^TPBB^TP\zeta )\notag \\ =&\ \frac{1}{\left| {\zeta _1 } \right|}\zeta ^T(A^TP+PA+\Pi _1 ^2C^TC+PBB^TP)\zeta. \end{align} $ (29)

Choose $ Q=-(A^TP+PA+\Pi _1 ^2C^TC+PBB^TP) $, then

$ \begin{align} \label{eq30} \dot {V}\le -\frac{1}{\left| {\zeta _1 } \right|}\zeta ^TQ\zeta. \end{align} $ (30)
$ \begin{align} Q&=-(A^TP+PA+\Pi _1 ^2C^TC+PBB^TP)\notag \\ &=\left[{\begin{array}{l} k_3 ^2k_7 ^2\Big(\dfrac{1}{2}k_3 k_7 \gamma _1 ^3+\gamma _1 \phi _1 -\dfrac{1}{4}\gamma _1 ^2\Big)-\Pi _1 ^2{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}\dfrac{1}{2}k_3 k_7 \Big(\gamma _1-k_3 k_7 \gamma _1 ^2\Big) \\[2mm] {\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}\dfrac{1}{2}k_3 k_7 \Big(\gamma _1 -k_3 k_7 \gamma _1 ^2\Big){\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}\dfrac{1}{2}{\kern 1pt}k_3 k_7 \gamma _1 -{\kern 1pt}1{\kern 1pt}{\kern 1pt}{\kern 1pt} \\ \end{array}} \right] \end{align} $ (31)

When $ Q $ is a symmetric and positive definite matrix, there exists $ \dot {V}\le -\frac{1}{\left| {\zeta _1 } \right|}\zeta ^TQ\zeta <0 $, which satisfies the Lyapunov theorem. Then the control system will converge to origin in finite time [16].

The matrix $ Q $ is given as (31), see the bottom of this page.

Set $ Q $ to be a positive definite matrix, and for finite time convergence characteristics the ranges of $ \gamma _1 $ and $ \phi _1 $ must be satisfied with

$ \begin{align} & {\begin{cases} \gamma _1 >\dfrac{2}{k_3 k_7 }\\ \phi _1 >\dfrac{k_3 k_7 \gamma _1 ^2}{4(k_3 k_7 \gamma _1 -2)}+\dfrac{\Pi _1 ^2}{k_3 k_7 \gamma _1 } \end{cases}} {\kern 1pt} \notag\\[2mm] &\quad\ \ \Rightarrow \begin{cases} \gamma _1 >\dfrac{4J(L_r L_s -L_m ^2)}{3n_g n_p L_m \phi _s } \\ \phi _1 >\dfrac{3L_m \phi _s n_p n_g \gamma _1 ^2}{12L_m \phi _s n_p n_g \gamma _1 -16J(L_r L_s -L_m ^2)}\\[4mm] \qquad\ \ +\dfrac{2J\Pi _1 ^2(L_r L_s -L_m ^2)}{3L_m \phi _s n_p n_g \gamma _1 }. \end{cases} \end{align} $ (32)

Similarly, the perturbation $ G_2 $ is differentiable as shown in (18), and its first-order derivative is uniformly bounded, so that it satisfies

$ \begin{align} \label{eq31} \left| {\dot {G}_2 } \right|\le \Pi _2 {\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}\forall t>0. \end{align} $ (33)

The ranges of $ \gamma _2 $ and $ \phi _2 $, which ensure $ \sigma _2 \to 0 $, $ \dot {\sigma }_2 \to 0 $ in finite time, are deduced as

$ \begin{align} \label{eq32}&\begin{cases} \gamma _2 >\dfrac{2}{k_7 } \\ \phi _2 >\dfrac{k_7 \gamma _2 ^2}{4(k_7 \gamma _2 -2)}+\dfrac{\Pi _2 ^2}{k_7 \gamma _2 } \end{cases} {\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt} \notag\\[2mm] &\quad\ \ \Rightarrow \begin{cases} \gamma _2 >\dfrac{2(L_r L_s -L_m ^2)}{L_s } \\ \phi _2 >\dfrac{L_s \gamma _2 ^2}{4L_s \gamma _2 -8(L_r L_s -L_m ^2)}+\dfrac{\Pi _2 ^2(L_r L_s -L_m ^2)}{L_s \gamma _2 }. \end{cases} \end{align} $ (34)

As long as the control parameters $ \gamma _1 $, $ \phi _1 $, $ \gamma _2 $ and $ \phi _2 $ fall into the range of (32) and (34), the whole DFIG-based WT system may operate steadily.

4 Simulation Results

The proposed SOSM was simulated on a 1.5 MW WT system. The system parameters based on a 1.5 MW wind turbine are shown in Table Ⅰ [14]. The parameter perturbations of resistances $ \Delta K $ and the turbine damping $ \Delta R_r $ in $ \Delta k_2 $ and $ \Delta k_4 $, are taken into account as $ \pm 20 \% $ of their nominal values and shown in Table Ⅰ. The simulation step size is set as 0.001 s. In designing the SOSM controller for WT system, the parameters are deliberately chosen according to the conditions (32) and (34), together with a thorough analysis and comprehensive computer simulation, to guarantee its finite time stability. Finally, these parameters are chosen as $ \phi _1 =20 $, $ \gamma _1 =1\times 10^5 $, $ \phi _2 =0.03 $, $ \gamma _2 =1000 $.

Table Ⅰ Characteristic of the Simulated Wind Turbine System

In order to verify the control effect under different operating environment, two case studies are considered in this section. Case 1 uses a stepwise-varying wind speed to test the validity of the proposed SOSM control strategy under a sudden wind change. A PID control and a standard first order sliding mode (FOSM) control are also constituted in this case, to show the robustness and chattering-free behavior of the SOSM controller. The effectiveness has been further demonstrated using the turbulent wind speed in Case 2, to test validity of the proposed SOSM control strategy under realistic VSWT conditions.

Generally, the cut-in wind speed of 1.5 MW WT is 3 m/s, the cut-out wind speed is 25 m/s, and the rated wind speed is in the range of 11 m/s to 13 m/s. Therefore, in these simulations, the wind speed varies within the range between 3 m/s and 12 m/s, making the WT system operates in the partial load region. Two case studies are described as follows.

4.1 Case 1: Stepwise Wind Speed

In order to show the controller performances under a sudden wind change, the fast step variations of the wind speed are used in this simulation, as shown in Fig. 5. For comparing purpose, a PID control and a FOSM control are also constituted. The PID control parameters are chosen as $ k_{p1} =1.2 $, $ k_{i1} =0.08 $, $ k_{d1} =0.07 $; $ k_{p2} =6.9 $, $ k_{i2} =20.6 $, $ k_{d2}=7.6 $, and the FOSM control parameters are $ \varepsilon _1 =5 $, $ \delta _1 $ $=$ $100 $; $ \varepsilon _2 =3 $, $ \delta _2 =10 $.

Figure 5 Stepwise wind speed profile.

Fig. 6 shows the regulation performances of the wind turbine rotor speed $ \omega _w $, the generator rotor current $ I_{rd} $ and the power coefficient $ C_p $. It can be observed from Fig. 6 (a) that all the three controllers can make the rotor speed track the optimal speed under the same wind speed variation. How-ever, the tracking accuracy of the rotor speed with the proposed SOSM control strategy is higher, and the response time is shorter. The resulting generator d-axis rotor current $ I_{rd} $ shown in Fig. 6 (b) demonstrates that the PID control scheme is not robust for parameter uncertainties. Since the wind turbine rotor speed can track its optimal value tightly, the power coefficient $ C_p $ will reach its maximum value $ C_{p\max } $. As shown in Fig. 6 (c), the SOSM control strategy gives a faster response during sudden change of the wind speed, and makes the $ C_p $ recover to $ C_{p\max } $ in 0.2 s after the sudden drop, yielding the maximum aerodynamic power.

Figure 6 Regulation performances of $ \omega _w $, $ I_{rd} $ and $ C_p . $

The actual controller output $ U_{rd} $ and $ U_{rq} $ are shown in Fig. 7. The switching term $ \varepsilon \cdot {\rm sgn}(\sigma ) $ appears directly in control input (20) and (21) for the FOSM controller, while the discontinuity in the SOSM controller (22) and (23) is hidden in their derivative term. Therefore, as is shown in Fig. 7, the chattering phenomenon, which is rather serious under FOSM method, is almost eliminated under the proposed SOSM control strategy. Actually, smooth $ U_{rd} $ and $ U_{rq} $ can reduce the mechanical stress on the VSWT system.

Figure 7 Comparison of the FOSM and the SOSM controller outputs.
4.2 Case 2: Randomly Varying Wind Speed

In order to evaluate the tracking performance, robustness and chattering-free behavior of the proposed SOSM control strategy under realistic VSWT conditions, a 10-min randomly varying wind speed is adopted in the simulation. As shown in Fig. 8, the wind speed is ranging between 3 m/s and 10 m/s. The evolution of the rotor speed and the optimal speed is depicted in Fig. 9 (a). Compared with the PID controller, the SOSM controller features more accurate and faster response in rotor speed tracking. $ I_{rd} $ is shown in Fig. 9 (b). The result shows that, the SOSM controller is robust for the internal disturbances and external disturbances, which are caused by parameter uncertainties and wind speed fluctuations respectively. And this yields a very good tracking performance of rotor current. Fig. 9 (c) shows the power conversion coefficient, since the SOSM control strategy gives an accurate tracking of the optimal rotor speed, the $ C_p $ successfully maintains at $ C_{p\max } $, which yields the maximum aerodynamic power. The three simulation results in Fig. 10 indicate that the control objectives of maximum power point tracking and minimized stator reactive power are achieved effectively.

Figure 8 Randomly varying wind speed profile.
Figure 9 Regulation performances of $ \omega _w $, $ I_{rd} $ and $ C_p . $
Figure 10 SOSM controller outputs.

The evolution of the SOSM controller output $ U_{rd} $ and $ U_{rq} $ are shown in Fig. 10. As can be seen from the figures, the SOSM controller does not cause any intense chattering, and the smooth features determined by the use of SOSM control strategy will reduce the mechanical stress on the system.

5 Conclusion

In this paper, a new SOSM control approach for DFIG-based VSWT system is proposed to achieve the objectives of maximum power point tracking and minimum stator reactive power. The two SOSM controllers are designed based on the super-twisting algorithm, which only require measurement of the sliding variables without using information about the time derivatives of the sliding constraint. Quadratic form Lyapunov function method is employed to choose controller parameters, and guarantee the finite time stabilization of closed-loop system. The controllers are simulated based on a complete model of the DFIG-based VSWT, which includes both the mechanical and the electric dynamics, together with parameter uncertainties and perturbations. A PID control method and a standard SMC are also carried out for comparison. The simulation results indicate that the proposed control strategy is rather suitable for controlling the DFIG-based VSWT system, and the control objectives are successfully achieved in both step wise-varying wind speed and a randomly varying turbulent wind speed. Compared with the PID controller and conventional first-order SMC, the proposed control strategy shows a higher robustness, and the chattering phenomenon of rotor control voltage is almost eliminated.

References
1
C. Ming and Z. Ying, "The state of the art of wind energy conversion systems and technologies: A review, " Energy Convers. Manage. , vol. 88, pp. 332-347, Dec. 2014. http://www.sciencedirect.com/science/article/pii/S0196890414007614
2
M. L. Corradini, G. Ippoliti, and G. Orlando, "Fully sensorless robust control of variable-speed wind turbines for efficiency maximization, " Automatica, vol. 49, no. 10, pp. 3023 -3031, Oct. 2013. http://www.sciencedirect.com/science/article/pii/S0005109813003762
3
F. Poitiers, T. Bouaouiche, and M. Machmoum, "Advanced control of a doubly-fed induction generator for wind energy conversion, " Electric Power Syst. Res. , vol. 79, no. 7, pp. 1085-1096, Jul. 2009. http://www.sciencedirect.com/science/article/pii/S0378779609000352
4
J. S. Wang, N. Tse, and Z. W. Gao, "Synthesis on PI-based pitch controller of large wind turbines generator, " Energy Convers. Manage. , vol. 52, no. 2, pp. 1288-1294, Feb. 2011. http://www.sciencedirect.com/science/article/pii/S0196890410004280
5
I. Munteanu, N. A. Cutululis, A. I. Bratcu, and E. Ceangǎ, "Optimization of variable speed wind power systems based on a LQG approach, " Control Eng. Pract. , vol. 13, no. 7, pp. 903-912, Jul. 2005. http://www.sciencedirect.com/science/article/pii/S0967066104002254
6
H. Chitsaz, N. Amjady, and H. Zareipour, "Wind power forecast using wavelet neural network trained by improved Clonal selection algorithm, " Energy Convers. Manage. , vol. 89, pp. 588-598, Jan. 2015. http://www.sciencedirect.com/science/article/pii/S0196890414008814
7
S. Abdeddaim, A. Betka, S. Drid, and M. Becherif, "Implementation of MRAC controller of a DFIG based variable speed grid connected wind turbine, " Energy Convers. Manage. , vol. 79, pp. 281-288, Mar. 2014. http://www.sciencedirect.com/science/article/pii/S0196890413007796
8
B. Boukhezzar and H. Siguerdidjane, "Nonlinear control with wind estimation of a DFIG variable speed wind turbine for power capture optimization, " Energy Convers. Manage. , vol. 50, no. 4, pp. 885-892, Apr. 2009. http://www.sciencedirect.com/science/article/pii/S0196890409000065
9
X. B. Kong and X. J. Liu, "Nonlinear model predictive control for DFIG-based wind power generation, " Acta Automat. Sin. , vol. 39, no. 5, pp. 636-643, May 2013. http://en.cnki.com.cn/Article_en/CJFDTOTAL-MOTO201305023.htm
10
Z. Y. Chen, M. H. Yin, C. X. Cai, B. Y. Zhang, and Y. Zou, "An accelerated optimal torque control of wind turbines for maximum power point tracking, " Acta Automat. Sin. , vol. 41, no. 12, pp. 2047-2057, Dec. 2015. http://www.aas.net.cn/EN/Y2015/V41/I12/2047
11
B. Beltran, T. Ahmed-Ali, and M. El Hachemi Benbouzid, "Sliding mode power control of variable-speed wind energy conversion systems, " IEEE Trans. Energy Convers. , vol. 23, no. 2, pp. 551-558, Jun. 2008. http://ieeexplore.ieee.org/document/4270775/
12
G. Bartolini, A. Ferrara, and E. Usai, "Chattering avoidance by second-order sliding mode control, " IEEE Trans. Automat. Control, vol. 43, no. 2, pp. 241-246, Feb. 1998. http://ieeexplore.ieee.org/xpls/abs_all.jsp?arnumber=661074
13
F. Valenciaga and C. A. Evangelista, "2-Sliding active and reactive power control of a wind energy conversion system, " IET Control Theory Appl. , vol. 4, no. 11, pp. 2479-2490, Nov. 2010. http://ieeexplore.ieee.org/xpls/abs_all.jsp?arnumber=5645800
14
B Beltran, M El Hachemi Benbouzid, and T Ahmed-Ali, "Second-order sliding mode control of a doubly fed induction generator driven wind turbine, " IEEE Trans. Energy Convers. , vol. 27, no. 2, pp. 261-269, Jun. 2012. http://ieeexplore.ieee.org/document/6130597/
15
C. Evangelista, F. Valenciaga, and P. Puleston, "Active and reactive power control for wind turbine based on a MIMO 2-sliding mode algorithm with variable gains, " IEEE Trans. Energy Convers. , vol. 28, no. 3, 682-689, Sep. 2013. http://ieeexplore.ieee.org/document/6582571/
16
J. A. Moreno and M. Osorio, "Strict Lyapunov functions for the super-twisting algorithm, " IEEE Trans. Autom. Control, vol. 57, no. 4, pp. 1035-1040, Apr. 2012. http://ieeexplore.ieee.org/document/6144710
17
J. Zaragoza, J. Pou, A. Arias, C. Spiteri, E. Robles, and S. Ceballos, "Study and experimental verification of control tuning strategies in a variable speed wind energy conversion system, " Renew. Energy, vol. 36, no. 5, pp. 1421-1430, May 2011. http://www.sciencedirect.com/science/article/pii/S0960148110005045
18
Gao W. B.. Variable Structure Control Theory. Beijing, China: China Science and Technology Press, 1990.
19
A. Levant, "Sliding order and sliding accuracy in sliding mode control, " Int. J. Control, vol. 58, no. 6, pp. 1247-1263, Dec. 1993.