Direct Voltage Control of a Doubly Fed Induction Generator by Means of Optimal Strategy

Abstract

The major objective of the investigation reported in this article is to demonstrate the feasibility of controlling a Doubly Fed Induction Generator actuating directly on the rotor voltage produced by the Rotor Side Converter, as its reference value may be determined analytically, after definition of the control objective. Two usual objectives are here considered: maximum power extraction from wind (MPPT) and stator reactive power equal to zero. This last objective defines the reference slip to be considered in the formulation of developed power that, jointly with the reactive power equation, forms the system to calculate the rotor reference voltages. The process is completed by specifying the desired dynamical response. Thus, the angular velocity of the rotor should quickly reach its reference value, which requires maximal power acceleration at the beginning, but respects the restriction that no overshoot should be allowed. This is achieved by means of a constrained optimization process solved in real time. Following recent trends, only measurements obtained from stator (voltages and currents) sensors are used. This way, angular velocity and rotor currents are estimated in real time. An algorithm for inductance estimation is also included, which prevents deviations of nominal values that could lead to false reference voltages.

1. Introduction

Although many wind turbine technologies are in use [1], and some others are in development [2,3], nowadays, the Doubly Fed Induction Generator (DFIG) is still used in about 50% of all Wind Energy Conversion Systems (WECSs) in operation. The predominance of this technology, in comparison to others, is due to its capacity of maximizing energy production in a wide range of wind conditions, allied to a good performance in dynamical response to input variations, including undesirable disturbance, in the Power System [1]. Thus, the focus of the present work is how to control a DFIG efficiently, without relying on parameters that can be corrupted over time. Therefore, it is noteworthy that the proposed control strategy is specific to DFIG, and there is no guarantee of the efficiency of its application in other technologies. The WECS based on DFIG possesses the stator connected directly to the grid and is projected considering that about 30% of active nominal power may be transferred through the rotor circuit, which is connected to the grid via back-to-back converters. The Grid-Side Converter (GSC) is responsible for maintaining constant voltage at the DC-link and controlling the reactive power at its output. Thus, active power may flow from the rotor to stator via GSC. The Rotor-Side Converter (RSC) controls active and reactive powers at the stator [4] by means of velocity control. Both stator and rotor circuits deliver power to the grid, but the direction of power flow via the rotor depends on the wind velocity and, consequently, on the angular velocity of the rotor.

Several strategies for controlling RSC are found in technical literature. The most usual method adopts a Field-Oriented Control (FOC), using linear PI controllers [5,6,7,8,9,10,11]. By this proposition, a direct and quadrature (d-q) decomposition for the rotor currents is accomplished, so that active and reactive powers flowing through the gap to the stator are controlled by d-q components. The FOC usually adopts two PI control loops: the inner one for controlling the rotor current and the outer for controlling the voltage applied to the rotor circuit. Owing to the d-q decomposition, the control of DFIG may be accomplished in decoupled loops, that is, the stator active power is controlled by variables of the q-axis, and the stator terminal voltage or reactive power is controlled by the control loop of the d-axis [7].

Another widely studied control approach is the Direct Power Control (DPC) [12,13,14]. This approach is based on Direct Torque Control (DTC) [12], which proposes firstly the elimination of the current control loop, and secondly a control of the stator active and reactive powers, accomplished by selecting directly the voltage applied to the rotor via RSC [13]. For this purpose, the states of stator active and reactive powers, as well the stator flux position, are measured or estimated. Using a hysteresis control allied to a switching table, the states of the RSC switches are defined for further application to the rotor circuit [13,14].

Other authors adopt different approaches like Predictive Control [15], Robust Control [16], Sliding-mode Control (SMC) [17], or even combined techniques as in [18], aiming to overcome specific deficiencies or to increase the efficiency inherent in the traditional algorithms of FOC or DPC. Propositions of optimal controllers supported by metaheuristics may be also found in the literature, as in [19,20]. A formulation of an optimization nonlinear process for determining the optimal values of rotor voltages is presented in [21]. The constraints are defined in order to respect regions where steady-state stable operation of DFIG occurs, limited by optimal values of rotor voltages. The authors of [22] present a formulation to determine analytically the values of voltages to be applied at the rotor circuit in a steady state. A nonlinear system composed of three equations is built, consisting of two equations corresponding to the real part and imaginary part of the rotor current, and the third one defines the electromagnetic torque. However, the system is solved symbolically using proprietary software. Moreover, the proposals of [21,22] do not incorporate dynamical simulations.

This article proposes a new control strategy for RSC. GSC is considered to work optimally and efficiently, that is, its power factor is maintained unitary and the dc battery voltage is controlled, so that active power can flow from or to the rotor. An equation system derived from the equivalent circuit is adopted to represent the desired operating point of DFIG, and the solution of this system defines the value of the steady-state rotor supply voltages. In this system, the developed power ($P_D$), which is the mechanical power converted to electrical by the generator, must be equal to the maximum power provided by the turbine from the wind, ignoring mechanical losses, and the stator reactive power ($Q_S)$ must equal zero. At steady state, the wind turbine is considered to work at the point of maximum wind power extraction, so the steady-state reference speed is known. It is also considered that the wind turbine is connected to the 13.8 kV bus of a relatively weak 69/13.8 kV grid represented by its equivalent of Thévenin; therefore, the stator voltage may vary during transient states. For this reason, the equations system must be solved in real time. The unknown variables are the real and imaginary parts of the complex rotor voltage ($V_{Rr}$ and $V_{Ri}$, respectively). After determining the rotor reference voltage, it is possible to outline a strategy for achieving it. Direct application of the steady-state reference voltage to the machine rotor can cause undesirable transients, such as large variations in electromagnetic torque and, consequently, very high transient currents. For this reason, a new strategy is proposed here to achieve rotor voltage reference values, without violating existing limits in other variables. This strategy is designed so that the speed of evolution of the control variables ($V_{Rr}$ and $V_{Ri}$) is governed by the dynamic evolution of the rotor speed (${\omega }_m$). Thus, the rate of evolution of the control variables in the direction of the reference must be controlled. An optimization process was formulated to minimize the convergence time of the rotor speed, restricting transient variations in net generated power, in order to accelerate the machine without exceeding the DFIG current nominal limits.

Magnetic saturation affects the inductance values of the DFIG, and consequently the control performance [23]. However, the effect of this saturation mainly influences the magnetizing inductance, since the path of the magnetic leakage flux is mainly air; thus, the variations of the leakage inductances are irrelevant to the control [24]. Changes in inductance values may lead to a different operating point than expected, and for this reason a method for estimating the magnetizing inductance under steady-state conditions is proposed here. In addition, in order to eliminate sensor dependence for evaluating rotor quantities, mainly because of the disadvantages of using sensors (robustness, maintenance, costs and cabling [25,26]), a solution is proposed here to estimate currents and the rotor speed by means of stator quantities.

The main advantages of the proposed strategy are the accuracy of the steady-state response and that a single adjustment of a single control parameter is valid for the entire wind speed range, as reference values are known in advance, and the real-time optimization process adapts to actual conditions. Through this new control strategy, which avoids the knowledge of controller data (in this case, unnecessary), it would be easier to simulate the dynamic response. Sensitive data from manufacturers would no longer be needed. This would constitute a new dimension in the relationship between manufacturers, consultants, entrepreneurs and agencies for the regulation and operation of the electrical system. The motivation of our strategy was inspired to solve such problems, which are reported in textbooks, for example, in reference [27] (p. 132). If manufacturers make it possible to define a single control parameter at the level of supervisory control system [27], the wind farm operator may extend it, easily, to each wind profile along the day or for different seasons. The proposed method to estimate the slip and the magnetizing inductance only by means of stator variables makes it possible to minimize the sensors needed to control the DFIG. Another advantage is the online correction of the magnetizing inductance, which is shown to be a critical parameter to reach the unitary power factor operation.

2. Materials and Methods

The Wind Energy Conversion System used to specify, analyze and test the control strategy proposed in this work is sketched in Figure 1. It consists of a wind turbine equipped with a generator of type DFIG, connected to an infinite bus of 13.8 kV, by a line whose impedance is 0.3943 + 1.6564j p.u. in a 100 MVA basis. The GSC is assumed to work ideally, assuring constant voltage at the DC link, and unity power factor at its output (grid side). The RSC is a pulse width modulation (PWM) device of the voltage source type. Its dynamical response is not considered, once the objective is to investigate the electromechanical transients of the machine, whereby harmonics caused by converter switching are not considered. The machine adopted for simulations was a 2 MW induction generator. The system was simulated in the Scilab environment using the 4th order Runge–Kutta numerical method to solve differential equations.

2.1. Wind Turbine

The mechanical power extracted from the wind by the turbine, $P_t$, is given in watts by

$$P_t=0.5\rho \pi R^2{C}_p\left(\beta ,\lambda \right)V_w^3$$

(1)

where $\rho$ is the air density; $R$ is the radius, in meters, of the circle described by the blades; $V_w$ is the wind velocity, in m/s; and $C_p\left(\lambda ,\beta \right)$ is the power coefficient, which determines the capacity of the turbine by extracting power from the wind. It is a function of the pitch angle $\beta$ and of the tip-speed ratio $\lambda$ [28]. The tip-speed ratio $\lambda$ is calculated by

$$\lambda =\frac{{\omega }_{t}R}{V_w}$$

(2)

where ${\omega }_t$ is the rotational velocity of the turbine. WECS using DFIG requires a mechanical velocity multiplier (gearbox-$gb$) in order to convert the low rotational velocity of the turbine and the high mechanical torque into quantities compatible to the DFIG operation. Thus, the relationship between the turbine velocity and the DFIG rotor velocity is given by

$${\omega }_t=\frac{{\omega }_m}{gb}$$

(3)

where $w_m$ is the mechanical speed of DFIG.

The mechanical torque, $T_m$, imposed by the turbine to the generator may be obtained from

$$T_m=\frac{P_t}{{\omega }_m}=\frac{0.5\rho \pi R^3{C}_p\left(\beta ,\lambda \right)V_w^2}{\lambda gb}$$

(4)

The curve of the power coefficient $C_p\left(\lambda ,\beta \right)$ may assume different forms for different turbines. In this work, the adopted curve $C_p\left(\lambda ,\beta \right)$ is given in [29].

The maximum power extraction from wind occurs when the power coefficient assumes its maximal value, corresponding to $\beta =0$ and $\lambda ={\lambda }_{otm}=7.206$. To attend this condition, the reference velocity of the DFIG rotor, for different wind velocities, is given by

$${\omega }_m^{*}=\frac{{\lambda }_{otm}{V}_{w}gb}{R}$$

(5)

In order to maintain the turbine in Maximum Power Point Tracking (MPPT), the angular velocity must be controlled so that the power coefficient is continuously maximal, corresponding to an optimal tip-speed ratio. This must occur from the moment the generator begins to produce power, at a cut-in wind velocity, until the rotor reaches the maximal admitted velocity for the turbine. When this occurs, the control objective is no longer to maintain the optimal tip-speed ratio but to assure the turbine velocity is no higher than the maximal allowed value. This causes a reduction of the power coefficient, reducing the amount of power extracted from the wind. An increase in the wind velocity would lead the turbine to produce more power than the nominal power of DFIG. From this moment on, the pitch control acts on the pitch angle of the blades, so that the absorbed power by the turbine equals the nominal value of DFIG [22,29].

The turbine used in the present work was taken from [29]. Its characteristics are shown in

Table 1. Parameters of the wind turbine.

Parameter	Value	Unit
Diameter	75	m
Gearbox ratio	1:100	-
Inertia constant	2.5	s
Nominal rotor speed	18	RPM
Nominal wind speed	12	m/s

.

2.2. DFIG Model

The DFIG was modeled by both its steady-state equations and its representation in the direct and quadrature d-q coordinates. Equivalent circuit equations were used to find the rotor steady-state values and to formulate the constrained nonlinear optimization problem, which defines the strategy of evolution of the control variables to their reference value. The DFIG dynamic equations resulting from the d-q reference frame transformation were used to calculate estimated mechanical velocity, rotor currents and magnetizing inductance values. The parameters of the generator used in the simulations can be seen in

Table 2. Parameters of the Doubly Fed Induction Generator (DFIG).

Parameter	Value	Unity
Number of poles	4	-
Frequency	50	Hz
Generator speed	1000–1900	RPM
Inertia constant	0.5	s
Mutual inductance	3.0	p.u.
Stator leakage inductance	0.1	p.u.
Rotor leakage inductance	0.08	p.u.
Stator resistance	0.01	p.u.
Rotor resistance	0.01	p.u.

[29].

2.2.1. Equivalent Circuit Equations

In [21,22] the authors demonstrate the steady-state equations in their trigonometric form. In this work, to simplify the algebraic solution of the equation system as well as the restricted nonlinear optimization, the equations are presented in their algebraic form.

The equivalent circuit of DFIG can be seen in Figure 2.

Considering ${\underset{\_}{V}}_S=V_S+j0$, the complex currents of stator (${\underset{\_}{I}}_S$) and rotor (${\underset{\_}{I}}_R^{\prime })$, obtained from the equivalent circuit, are

$${\underset{\_}{I}}_S=\frac{-V_s\left(R_R^{\prime }+js{X}_R^{\prime }\right)+j{\underset{\_}{V}}_R^{\prime }{X}_m}{\left(R_s{R}_R^{\prime }-s{X}_s{X}_R^{\prime }+s{X}_m^2\right)+j\left(R_R^{\prime }{X}_s+R_{s}s{X}_R^{\prime }\right)},$$

(6)

$${\underset{\_}{I}}_R^{\prime }=\frac{{\underset{\_}{V}}_R^{\prime }\left(R_S+j{X}_S\right)-j{V}_{s}s{X}_m}{\left(R_s{R}_R^{\prime }-s{X}_s{X}_R^{\prime }+s{X}_m^2\right)+j\left(R_R^{\prime }{X}_s+R_{s}s{X}_R^{\prime }\right)},$$

(7)

where, $R_S$ and $X_s$ are the stator resistance and reactance, respectively; $s$ is the slip; $R_R^{\prime }$ and $X_R^{\prime }$ are the rotor resistance and reactance referred to the stator, respectively; and $X_m$ is the magnetizing reactance. The parameters $X_s$ and $X_R^{\prime }$ are defined as

$$X_S=X_{lS}+X_m,$$

(8)

$$X_R^{\prime }=X_{lR}^{\prime }+X_m,$$

(9)

where, $X_{lS}$ and $X_{lR}^{\prime }$ are the stator and rotor leakage reactances, respectively.

Using Equations (5) and (6), it is possible to calculate the stator and rotor powers, as well the developed power $P_D$. In order to simplify the equations, the variables below are introduced:

$$A_D=\left(R_s{R}_R^{\prime }-s{X}_s{X}_R^{\prime }+s{X}_m^2\right),$$

(10)

$$B_D=\left(R_R^{\prime }{X}_s+R_{s}s{X}_R^{\prime }\right).$$

(11)

The complex powers of stator ($S_S$) and rotor ($S_R$) are calculated by

$$S_S=3V_S{\underset{\_}{I_S}}^{*}=3\left(P_S+j{Q}_S\right),$$

(12)

$$S_R=3{\underset{\_}{V}}_R^{\prime }{\underset{\_}{I}}_R^{\prime }{}^{*}=3\left(P_R+j{Q}_R\right).$$

(13)

Defining ${\underset{\_}{V}}_R^{\prime }=V_{Rr}+j{V}_{Ri}$ and substituting Equations (5) and (6) in (11) and (12) obtains each power equation (active and reactive) as a function of three variables: slip, real and imaginary parts of the rotor voltage, as presented in Equations (13)–(16).

$$P_S=\frac{3}{A_D^2+B_D^2}\left[V_S{X}_m\left(B_D{V}_{Rr}-A_D{V}_{Ri}\right)-V_S^2\left(R_R^{\prime }{A}_D+s{X}_R{B}_D\right)\right],$$

(14)

$$Q_S=\frac{3}{A_D^2+B_D^2}\left[V_S{X}_m\left(-A_D{V}_{Rr}-B_D{V}_{Ri}\right)+V_S^2\left(-R_R^{\prime }{B}_D+s{X}_R{A}_D\right)\right],$$

(15)

$$P_R=\frac{3}{A_D^2+B_D^2}\left[V_{Rr}^2\left(R_S{A}_D+X_S{B}_D\right)+V_{Ri}^2\left(R_S{A}_D+X_S{B}_D\right)-V_{S}s{X}_m\left(B_D{V}_{Rr}+A_D{V}_{Ri}\right)\right],$$

(16)

$$Q_R=\frac{3}{A_D^2+B_D^2}\left[V_{Rr}^2\left(R_S{B}_D-X_S{A}_D\right)+V_{Ri}^2\left(R_S{B}_D-X_S{A}_D\right)+V_{S}s{X}_m\left(A_D{V}_{Rr}-B_D{V}_{Ri}\right)\right].$$

(17)

The developed power in the machine is given by [30]

$$P_D=-\frac{3\left(1-s\right)}{s}{R}_R^{\prime }{\underset{\_}{I}}_R^{\prime }{}^2+\frac{3\left(1-s\right)}{s}{P}_R,$$

(18)

Considering yet

$$I_R^{\prime }{}^2={\left|{\underset{\_}{I}}_R^{\prime }\right|}^2={\underset{\_}{I}}_R^{\prime }{\underset{\_}{I}}_R^{\prime }{}^{*},$$

(19)

and, substituting Equation (6) in (18) results in

$${\left|{\underset{\_}{I}}_R^{\prime }\right|}^2=\frac{\left(R_S^2+X_S^2\right)(V_{Rr}^2+V_{Ri}^2)-2V_{S}s{X}_S{X}_m{V}_{Rr}-2V_{S}s{X}_m{R}_S{V}_{Ri}+V_S^2{s}^2{X}_m^2}{A_D^2+B_D^2}.$$

(20)

Finally, the equation for the developed power of a DFIG, in steady state, may be obtained by substitution of Equations (15) and (19) into (17):

$$\begin{gathered} P_D=\frac{3\left(1-s\right)}{A_D+j{B}_D}\left[R_S{X}_m^2(V_{Rr}^2+V_{Ri}^2)+V_S{X}_m\left(2R_R^{\prime }{X}_S-B_D\right)V_{Rr}+V_S{X}_m\left(2R_S{R}_R^{\prime }-A_D\right)V_{Ri} \\ -V_S^{2}s{R}_R^{\prime }{X}_m^2\right] \end{gathered}$$

(21)

The stator power may be so calculated, using the energy conservation principle, as

$$P_S=P_D+P_R-P_{R,loss}-P_{S,loss}$$

(22)

where:

$$P_{S,loss}=R_S{\left|\underset{\_}{I_S}\right|}^2$$

(23)

$$P_{R,loss}=R_R^{\prime }{\left|{\underset{\_}{I}}_R^{\prime }\right|}^2$$

(24)

2.2.2. Dynamic Model of DFIG

The dynamical modeling of DFIG is adopted in the synchronously rotating d-q reference frame, as presented in [28]:

$$V_{Sd}=R_S{I}_{Sd}-{\omega }_e{\lambda }_{Sq}+\frac{d{\lambda }_{Sd}}{dt}$$

(25)

$$V_{Sq}=R_S{I}_{Sq}+{\omega }_e{\lambda }_{Sd}+\frac{d{\lambda }_{Sq}}{dt}$$

(26)

$$V_{Rd}=R_R{I}_{Rd}-\left({\omega }_e-\frac{P}{2}{\omega }_m\right){\lambda }_{Rq}+\frac{d{\lambda }_{Rd}}{dt}$$

(27)

$$V_{Rq}=R_R{I}_{Rq}+\left({\omega }_e-\frac{P}{2}{\omega }_m\right){\lambda }_{Rd}+\frac{d{\lambda }_{Rq}}{dt}$$

(28)

where ${\omega }_e$ is the synchronous electric frequency of the grid; $P$ is the number of poles; $R_S$ and $R_R$ are the stator and rotor resistances; $I_{Sd}$, $I_{Sq}$, $I_{Rd}$ and $I_{Rq}$ are the stator and rotor currents; ${\lambda }_{Sd}$, ${\lambda }_{Sq}$, ${\lambda }_{Rd}$ and ${\lambda }_{Rq}$ are the stator and rotor linkage fluxes; $V_{Sd}$, $V_{Sq}$, $V_{Rd}$ and $V_{Rq}$ are the stator and rotor voltages, respectively. The equations relating currents and fluxes are

$${\lambda }_{Sd}=L_S{I}_{Sd}+L_m{I}_{Rd},$$

(29)

$${\lambda }_{Sq}=L_S{I}_{Sq}+L_m{I}_{Rq},$$

(30)

$${\lambda }_{Rd}=L_R{I}_{Rd}+L_m{I}_{Sd},$$

(31)

$${\lambda }_{Rq}=L_R{I}_{Rq}+L_m{I}_{Sq},$$

(32)

where, $L_m$ is the mutual inductance, $L_S$ is the stator proper inductance and $L_R$ is the rotor proper inductance. Parameters $L_s$ and $L_R$ are defined as

$$L_S=L_{lS}+L_m,$$

(33)

$$L_R=L_{lR}+L_m,$$

(34)

where, $L_{lS}$ is the stator leakage inductance and $L_{lR}$ is the rotor leakage inductance.

The electromagnetic torque may be given by

$$T_e=\frac{3}{2}\frac{P}{2}\left({\lambda }_{Rq}{I}_{Rd}-{\lambda }_{Rd}{I}_{Rq}\right),$$

(35)

and the equation to describe the dynamical behavior of the machine may be expressed using the moment of inertia as

$$J\frac{d{\omega }_m}{dt}=T_m-T_e,$$

(36)

where $J$ is the moment of inertia, whose relationship with the inertia constant $H$, in seconds, is [31]

$$J=\frac{2H}{{\omega }_0^2}{S}_b,$$

(37)

where ${\omega }_0$ is the synchronous velocity in mechanical radians per second and $S_b$ is the apparent power in volt-amperes.

2.3. Calculation of Reference Values for Control Variables

The operating point adopted as the basis for the control efficiency study proposed in this article corresponds to the maximum utilization of wind power and unit power factor in the stator. As explained in Section 2.1, for the turbine to extract maximum wind power, the power coefficient must assume the value $\lambda ={\lambda }_{otm}=7.206$. Therefore, Equation (4) calculates the reference speed and the corresponding reference slip in steady state.

DFIG’s power generation and speed control are done by RSC, which is a voltage source PWM converter. Thus, it is necessary to define the voltage that the converter must supply to the DFIG rotor circuit to control it at the desired operating point. In most studies found in the literature, the rotor reference voltage is given in polar form; in the present work, the reference voltage will be calculated in its Cartesian form. Thus, the rotor reference voltage will have the following form:

$${\underset{\_}{V}}_R^{\prime }{}^{*}=V_{Rr}^{*}+j{V}_{Ri}^{*}$$

(38)

This way, to keep DFIG operating at the desired operating point, the values of $V_{Rr}^{*}$ and $V_{Ri}^{*}$ must be known.

To extract maximum wind power, it is necessary that all mechanical energy absorbed by the turbine be transformed into electrical energy in the generator, disregarding the rotational losses, at the optimal power coefficient value. The second condition imposed on the control is that the stator reactive power is zero. Therefore, the equations system that defines the desired steady-state operating point for the wind turbine is

$$\{\begin{array}{c}{P}_D=P_t\\ Q_S=0\end{array}$$

(39)

Equations (14) and (20) will form the system proposed in Equation (37) to find the values of $V_{Rr}^{*}$ and $V_{Ri}^{*}$ that ensures the DFIG will work at the desired operating point in steady state. Thus, the proposed equation system is

$$\{\begin{array}{c}\frac{3\left(1-s^{*}\right)}{A_D^2+B_D^2}[R_S{X}_m^2(V_{Rr}^{*}{}^2+V_{Ri}^{*}{}^2)+V_S{X}_m\left(2R_R{X}_S-B_D\right)V_{Rr}^{*}+V_S{X}_m\left(2R_S{R}_R-A_D\right)V_{Ri}^{*}-V_S^2{s}^{*}{R}_R{X}_m^2]-P_t=0\\ V_S{X}_m\left(-A_D{V}_{Rr}^{*}-B_D{V}_{Ri}^{*}\right)+V_S^2\left(-R_R{B}_D+s^{*}{X}_R{A}_D\right)=0\end{array},$$

(40)

$$A_D=\left(R_s{R}_R-s^{*}{X}_s{X}_R+s^{*}{X}_m^2\right),$$

(41)

$$B_D=\left(R_R{X}_s+R_s{s}^{*}{X}_R\right).$$

(42)

In addition to the DFIG equivalent circuit parameters, the variables that form the equation system above are the reference slip, $s^{*}$, which is known for the desired operating point; the stator voltage magnitude ($V_S$); and the control variables of RSC, $V_{Rr}^{*}$ and $V_{Ri}^{*}$. As can be seen from Figure 1, the stator voltage will vary transiently according to the net current generated by the DFIG. Thus, the use of voltage sensors to measure the amplitude and phase of the stator voltage is necessary for the solution of the equation system.

Having all the data, the system of equations has a simple solution obtained analytically by solving the system defined in Equation (38) (e.g., by variable elimination), thus finding the desired values. Observe that the $P_D$ equation has a degree of 2, which results in two solutions. The choice of the appropriate solution is done by analyzing the magnitude of the current corresponding to each of these solutions. The value of the current is calculated by Equation (19), and the solution chosen will correspond to a current smaller than the magnitude of the DFIG rotor-rated current. The solution [$V_{Rr}^{*}$,$V_{Ri}^{*}$] is found considering ${\underset{\_}{V}}_S=V_S+j0$. As shown in the system proposed by Figure 1, the voltage ${\underset{\_}{V}}_S$ has a nonzero phase angle, so it is necessary that the actual value of its angle be compensated at the end of the process of calculating the rotor voltage reference values. In the simulations, the stator voltage angle was obtained by the solution of the Thevenin equivalent circuit, but in practice, as it is not possible to know the equivalent of Thevenin at all times, it is necessary to use a Phase-Locked Loop (PLL) algorithm.

2.4. Velocity of the Control Variables Trajectory

The desired operating point is defined by the measured wind velocity, which determines the reference velocity ${\omega }_m^{*}$, according to Equation (4). A change in wind velocity requires new settings of the control variables. The new reference values for these variables, $V_{Rr}^{*}$ and $V_{Ri}^{*}$, may be calculated directly by solving Equation (38). Thereafter, a strategy is required for defining the transition from initial steady-state values of these variables ($V_{Rr}^i$ and $V_{Ri}^i$) until the corresponding final values ($V_{Rr}^f$ and $V_{Ri}^f$), where $V_{Rr}^f=V_{Rr}^{*}$ and $V_{Ri}^f=V_{Ri}^{*}$. A sudden change of the control variables to their new reference values may cause undesirable or even damaging transients to the network state. For exemplary purposes, a wind velocity step from 7 m/s to 7.5 m/s was simulated, with application of the corresponding rotor voltage step. The resulting rotor current is shown in Figure 3, presenting a transient value high above the rated value.

The current transient caused by application of the rotor voltage step reaches, approximately, two times the rotor nominal current. A current of such magnitude could damage the circuits of the RSC converter. For this reason, the trajectory of control variables towards their reference value must be carefully designed, respecting physical limitations of generator components.

The basic idea of the strategy proposed here is to link the time variations of the control variables, ($\frac{d{V}_{Rr}}{dt}$ and $\frac{d{V}_{Rr}}{dt}$), to the time variations of the machine speed ($\frac{d{\omega }_m}{dt}$) in order to define the evolution of the control variables, assuring their convergence compatible to the convergence of the rotor velocity, when steady state is reached. This way, major disturbances are avoided.

Firstly, it will be stated that no overshot in the time response of the rotor velocity is desirable. Thus, this velocity is assumed to change from a previous reference to the next one according to a time-varying exponential function:

$${\omega }_m=\left({\omega }_m^i-{\omega }_m^f\right)e^{-f_{\omega }t}+{\omega }_m^f,$$

(43)

where, ${\omega }_m^i$ is the initial value of the velocity, ${\omega }_m^f$ is its final value (${\omega }_m^f={\omega }_m^{*}$) and $f_{\omega }$ is a frequency that defines the convergence of the rotor velocity. As the equations from the equivalent circuit presented in Section 2.2.1 contain the variable slip, instead of velocity, a simple relationship between $\frac{d{\omega }_m}{dt}$ and $\frac{ds}{dt}$ can be obtained:

$$\frac{ds}{dt}=\frac{d}{dt}\left(\frac{{\omega }_s-{\omega }_m}{{\omega }_s}\right),$$

(44)

$$\frac{ds}{dt}=-\frac{1}{{\omega }_s}\frac{d{\omega }_m}{dt},$$

(45)

Combination of Equations (41) with (43) leads to

$$\frac{ds}{dt}=\frac{1}{{\omega }_s}{f}_{\omega }\left({\omega }_m^i-{\omega }_m^f\right)e^{-f_{\omega }t},$$

(46)

Evaluating Equation (44) in $t=0$ gives

$${\left[\frac{ds}{dt}\right]}_{t=0}=\frac{1}{{\omega }_s}{f}_{\omega }\left({\omega }_m^i-{\omega }_m^f\right),$$

(47)

$${\left[\frac{ds}{dt}\right]}_{t=0}=\frac{1}{{\omega }_s}{f}_{\omega }\left[\left(1-s^i\right){\omega }_s-{\omega }_m^f\right],$$

(48)

Equation (46) may be rearranged as

$${\left[\frac{ds}{dt}\right]}_{t=0}=f_{\omega }\left(s^f-s^i\right).$$

(49)

where $s^f$ is the final slip, which must equal the new reference slip $s^{*}$; and $s^i$ is the initial slip, which could equal the previous reference value, if the machine was previously in a steady state. In order to establish a relationship between the time variations of the control variables with the corresponding slip time variations, Equation (48) is proposed. The basic idea of these equations is to search the aims of maintaining both power conversion and reactive control dynamically.

$$\{\begin{array}{c}\frac{d\left(P_D-P_t\right)}{dt}=0\\ \frac{d{Q}_S}{dt}=0\end{array}.$$

(50)

Unlike the system of equations seen in Equation (38), where the reference slip is defined by wind speed, the system of equations proposed here takes into account slip variations that are proportional to the rotor speed variations. Adopting linearization at a time $t=t_0$:

$$\frac{d\left(P_D-P_t\right)}{dt}=\frac{\partial \left(P_D-P_t\right)}{\partial s}\frac{ds}{dt}+\frac{\partial \left(P_D-P_t\right)}{\partial V_{Rr}}\frac{d{V}_{Rr}}{dt}+\frac{\partial \left(P_D-P_t\right)}{\partial V_{Ri}}\frac{d{V}_{Ri}}{dt}=0,$$

(51)

$$\frac{d{Q}_S}{dt}=\frac{\partial Q_S}{\partial s}\frac{ds}{dt}+\frac{\partial Q_S}{\partial V_{Rr}}\frac{d{V}_{Rr}}{dt}+\frac{\partial Q_S}{\partial V_{Ri}}\frac{d{V}_{Ri}}{dt}=0,$$

(52)

thus, the equations system is defined by

$$\{\begin{array}{c}\frac{\partial \left(P_D-P_t\right)}{\partial s}\frac{ds}{dt}+\frac{\partial \left(P_D-P_t\right)}{\partial V_{Rr}}\frac{d{V}_{Rr}}{dt}+\frac{\partial \left(P_D-P_t\right)}{\partial V_{Ri}}\frac{d{V}_{Ri}}{dt}=0\\ \frac{\partial Q_S}{\partial s}\frac{ds}{dt}+\frac{\partial Q_S}{\partial V_{Rr}}\frac{d{V}_{Rr}}{dt}+\frac{\partial Q_S}{\partial V_{Ri}}\frac{d{V}_{Ri}}{dt}=0\end{array}.$$

(53)

where, the numerical values of the partial derivatives $\frac{\partial \left(P_D-P_t\right)}{\partial V_{Rr}}$, $\frac{\partial \left(P_D-P_t\right)}{\partial V_{Ri}}$, $\frac{\partial Q_S}{\partial s}$ and $\frac{\partial Q_S}{\partial V_{Rr}}$ as well $\frac{\partial \left(P_D-P_t\right)}{\partial s}$ and $\frac{\partial Q_S}{\partial s}$ may be calculated at each time $t=t_0$ using the expressions

$$\frac{\partial \left(P_D-P_t\right)}{\partial V_{Rr}}=\frac{3\left(1-s\right)}{A_D^2+B_D^2}\left[V_S{X}_m\left(2R_R{X}_S-B_D\right)+2V_{Rr}{R}_S{X}_m^2\right],$$

(54)

$$\frac{\partial \left(P_D-P_t\right)}{\partial V_{Ri}}=\frac{3\left(1-s\right)}{A_D^2+B_D^2}\left[V_S{X}_m\left(2R_S{R}_R-A_D\right)+2V_{Ri}{R}_S{X}_m^2\right],$$

(55)

$$\frac{\partial \left(P_D-P_t\right)}{\partial s}=\frac{D{P}_D\frac{\partial N{P}_D}{\partial s}-N{P}_D\frac{\partial D{P}_D}{\partial s}}{D{P}_D^2},$$

(56)

where, $N{P}_D$ and $D{P}_D$ are, respectively, the numerator and denominator of the developed power from Equation (20). From Equation (54), the partial derivative of the developed power $P_D$ is obtained as a function of slip, s. In this calculation, $A_D$ and $B_D$ are considered as defined in Equations (9) and (10). Thus, the solution of Equation (54) is

$$\begin{array}{ll}N{P}_D& =3\left(1-s\right)[R_S{X}_m^2(V_{Rr}^2+V_{Ri}^2)+V_S{X}_m\left(R_R{X}_S-R_{S}s{X}_R\right)V_{Rr}\\ & +V_S{X}_m\left(R_S{R}_R+s\left(-X_m^2+X_S{X}_R\right)\right)V_{Ri}-V_S^{2}s{R}_R{X}_m^2]\end{array},$$

(57)

$$\begin{array}{ll}\frac{\partial N{P}_D}{\partial s}& =-3[R_S{X}_m^2(V_{Rr}^2+V_{Ri}^2)+V_S{X}_m\left(R_R{X}_S-2R_{S}s{X}_R\right)V_{Rr}\\ & +V_S{X}_m\left(R_S{R}_R+2s\left(-X_m^2+X_S{X}_R\right)\right)V_{Ri}-2V_S^{2}s{R}_R{X}_m^2]\\ & +3\left[-V_S{X}_m{R}_S{X}_R{V}_{Rr}+V_S{X}_m\left(-X_m^2+X_S{X}_R\right)V_{Ri}-V_S^2{R}_R{X}_m^2\right]\end{array}$$

(58)

$$D{P}_D=A_D^2+B_D^2$$

(59)

$$\frac{\partial D{P}_D}{\partial s}=2A_D\left(-X_s{X}_R+X_m^2\right)+2B_D{R}_s{X}_R.$$

(60)

The partial derivative expressions of the stator reactive power are as follows:

$$\frac{\partial Q_S}{\partial V_{Rr}}=\frac{1}{A_D^2+B_D^2}\left[-V_S{X}_m{A}_D\right],$$

(61)

$$\frac{\partial Q_S}{\partial V_{Ri}}=\frac{1}{A_D^2+B_D^2}\left[-V_S{X}_m{B}_D\right],$$

(62)

$$\frac{\partial Q_S}{\partial s}=\frac{D{Q}_S\frac{\partial N{Q}_S}{\partial s}-N{Q}_S\frac{\partial D{Q}_S}{\partial s}}{D{Q}_S^2},$$

(63)

where, $N{Q}_S$ and $D{Q}_S$ are, respectively, the numerator and denominator of the stator reactive power Equation (14). $A_D$ and $B_D$ are considered as shown in Equations (9) and (10). Thus, the solution of Equation (61) is

$$\begin{gathered} N{Q}_S=V_S{X}_m\left[\left(R_s{R}_R-s{X}_s{X}_R+s{X}_m^2\right)V_{Rr}-\left(R_R{X}_s+R_{s}s{X}_R\right)V_{Ri}\right]+V_S^2(-R_R^2{X}_s-s^2{X}_s{X}_R^2 \\ +s^2{X}_m^2{X}_R) \end{gathered}$$

(64)

$$\frac{\partial N{Q}_S}{\partial s}=V_S{X}_m\left[\left(-X_s{X}_R+X_m^2\right)V_{Rr}-R_s{X}_R{V}_{Ri}\right]+V_S^2\left(-2s{X}_s{X}_R^2+2s{X}_m^2{X}_R\right)$$

(65)

$$D{Q}_S=A_D^2+B_D^2.$$

(66)

$$\frac{\partial D{Q}_S}{\partial s}=2A_D\left(-X_s{X}_R+X_m^2\right)+2B_D{R}_s{X}_R.$$

(67)

After evaluating all partial derivatives and knowing the variables $V_{Rr}$, $V_{Ri}$ and $s$ at each instant, the solution of the system of equations proposed in Equation (51) is simply given by

$$\left[\begin{array}{c}\frac{d{V}_{Rr}}{dt}\\ \frac{d{V}_{Ri}}{dt}\end{array}\right]={\left[\begin{array}{cc}\frac{\partial \left(P_D-P_t\right)}{\partial V_{Rr}}& \frac{\partial \left(P_D-P_t\right)}{\partial V_{Ri}}\\ \frac{\partial Q_S}{\partial V_{Rr}}& \frac{\partial Q_S}{\partial V_{Ri}}\end{array}\right]}^{-1}\left[\begin{array}{c}\frac{\partial \left(P_D-P_t\right)}{\partial s}\\ \frac{\partial Q_S}{\partial s}\end{array}\right]\frac{ds}{dt}.$$

(68)

$$\left[\begin{array}{c}\frac{d{V}_{Rr}}{dt}\\ \frac{d{V}_{Ri}}{dt}\end{array}\right]=\left[\begin{array}{c}{\xi }_{Vr}\\ {\xi }_{Vi}\end{array}\right]\frac{ds}{dt}$$

(69)

where $\frac{ds}{dt}$ is calculated by Equation (47).

2.5. Optimizing Rotor Speed Convergence

Calculation of the convergence speeds of the control variables shown in Equation (66), Section 2.4, depends on the calculation of $\frac{ds}{dt}$. In Equation (47) the variable $f_{\omega }$ is responsible for defining the convergence speed of the rotor (mechanical) speed to its reference value, according to the exponential proposed in Equation (41). In order to maximize $f_{\omega }$ to lead to a fast convergence of the rotor speed, as well as not violating DFIG limits, a nonlinear constrained optimization technique is proposed in this work. Nonlinear constrained programming problems propose the maximization or minimization and has the following generic form [32]:

$$\{\begin{array}{c}\min f\left(x\right),x\in {\Re }^n\\ s.t.\{\begin{array}{c}{h}_i\left(x\right)=0,i=1,\dots ,m<n\\ g_j\left(x\right)\ge 0,j=1,\dots ,p\end{array}\end{array}$$

(70)

The problem consists of an objective function $f\left(x\right)$, respecting limits imposed by constraint functions, which may be equality $h_i\left(x\right)$ or inequality $g_j\left(x\right)$. Depending on the proposed problem and the nature of its component functions, different solution techniques may be applied. For the solution of the problem proposed in this work, the Lagrange multiplier method can be applied because it is a relatively simple problem, where the solution of the problem can be found analytically without the need for concern with numerical convergence, as in iterative methods. The first-order condition for the constrained problem using the Lagrange multiplier method is then given by [32].

$$\nabla ℒ\left(x\right)=\nabla f\left(x\right)+\lambda \nabla \mathrm{h}\left(x\right)+\mu \nabla g\left(x\right)=0.$$

(71)

where $ℒ\left(x\right)$ is the Lagrange Augmented Function (Lagrangian) and $\mu$ is the Lagrange multiplier.

The optimization problem of variable $f_{\omega }$ will be defined as follows for rising wind speed:

$$\{\begin{array}{c}\max f_{\omega }\left(V_{Rr},V_{Ri}\right)\\ s.t.P_{liq}\left(V_{Rr},V_{Ri}\right)\ge P_{lim}\end{array},$$

(72)

where $f_{\omega }$ can be defined by Equations (34) and (41) for $t=0$ as follows:

$$f_{\omega }=\frac{T_m-T_e}{J\left({\omega }_m^f-{\omega }_m^i\right)},$$

(73)

Torque $T_m$ is calculated by Equation (3), and $T_e$ is obtained from the developed power equation $P_D$, as shown below:

$$T_e=\frac{P_D}{{\omega }_m}=\frac{P_D}{\left(1-s\right){\omega }_s},$$

(74)

$$T_e=\frac{3\left[R_S{X}_m^2(V_{Rr}^2+V_{Ri}^2\right)+V_S{X}_m\left(2R_R{X}_S-B_D\right)V_{Rr}+V_S{X}_m\left(2R_S{R}_R-A_D\right)V_{Ri}-V_S^{2}s{R}_R{X}_m^2]}{{\omega }_s\left(A_D^2+B_D^2\right)}$$

(75)

The net power is the difference between the stator and the rotor active power, $P_{liq}=P_S-P_R$. Thus, from Equations (13) and (15), $P_{liq}$ is defined as

$$\begin{array}{ll}{P}_{liq}& =\frac{3}{A_D^2+B_D^2}[-\left(R_S{A}_D+X_S{B}_D\right)\left(V_{Rr}^2+V_{Ri}^2\right)+\left(1+s\right)V_S{X}_m{B}_D{V}_{Rr}\\ & +\left(-1+s\right)V_S{X}_m{A}_D{V}_{Ri}+V_S^2\left(R_R{A}_D+s{X}_R{B}_D\right)]\end{array}.$$

(76)

The limitation of the net power generated is justified by the fact that, to produce rotor acceleration, it is necessary that there is a difference between the mechanical torque and the electromagnetic torque, as can be observed by Equation (34). Thus, as the electromagnetic torque is directly proportional to the developed power, as shown in Equation (72), it is necessary to limit the distance of the electromagnetic torque, in order to respect an allowable limit of net power that already discounted the losses in the DFIG. This power limit will be proportional to the available wind power

$$P_{lim}=\alpha P_t,$$

(77)

the constant $\alpha$ is chosen by experimental simulations, and it is different for rising or falling wind speed, reflecting operational conditions of the site where the turbine will be installed.

The only constraint proposed in the problem is an inequality constraint, and for the solution of the first-order condition it must be considered active, that is

$$g\left(V_{Rr},V_{Ri}\right)=P_{liq}\left(V_{Rr},V_{Ri}\right)-P_{lim}=0.$$

(78)

The first-order condition for the optimization problem proposed in Equation (70) is:

$$\nabla ℒ\left(V_{Rr},V_{Ri},\mu \right)=\nabla f_{\omega }\left(V_{Rr},V_{Ri}\right)+\mu \nabla g\left(V_{Rr},V_{Ri}\right)=0$$

(79)

To find the stationary points of the function presented in Equation (77), it is necessary to assemble the system of equations that define the first-order condition for the Lagrange function and the problem constraint [32], as shown below:

$$\nabla f_{\omega }\left(V_{Rr},V_{Ri}\right)=\left[\begin{array}{c}\frac{\partial f_{\omega }}{\partial V_{Rr}}\\ \frac{\partial f_{\omega }}{\partial V_{Ri}}\end{array}\right]$$

(80)

$$\nabla g\left(V_{Rr},V_{Ri}\right)=\left[\begin{array}{c}\frac{\partial g}{\partial V_{Rr}}\\ \frac{\partial g}{\partial V_{Ri}}\end{array}\right]$$

(81)

$$\{\begin{array}{c}\frac{\partial f_{\omega }}{\partial V_{Rr}}+\mu \frac{\partial g}{\partial V_{Rr}}=0\\ \frac{\partial f_{\omega }}{\partial V_{Ri}}+\mu \frac{\partial g}{\partial V_{Ri}}=0\\ g\left(V_{Rr},V_{Ri}\right)=0\end{array}$$

(82)

where:

$$\frac{\partial f_{\omega }}{\partial V_{Rr}}=\frac{-1}{J\left({\omega }_m^f-{\omega }_m^i\right)}\frac{3}{{\omega }_s\left(A_D^2+B_D^2\right)}\left[V_S{X}_m\left(2R_R{X}_S-B_D\right)+2R_S{X}_m^2{V}_{Rr}\right]$$

(83)

$$\frac{\partial f_{\omega }}{\partial V_{Ri}}=\frac{-1}{J\left({\omega }_m^f-{\omega }_m^i\right)}\frac{3}{{\omega }_s\left(A_D^2+B_D^2\right)}\left[V_S{X}_m\left(2R_S{R}_R-A_D\right)+2R_S{X}_m^2{V}_{Ri}\right]$$

(84)

$$\frac{\partial g}{\partial V_{Rr}}=\frac{3}{\left(A_D^2+B_D^2\right)}\left[\left(1+s\right)V_S{X}_m{B}_D-2\left(R_S{A}_D+X_S{B}_D\right)V_{Rr}\right]$$

(85)

$$\frac{\partial g}{\partial V_{Ri}}=\frac{3}{\left(A_D^2+B_D^2\right)}\left[\left(-1+s\right)V_S{X}_m{A}_D-2\left(R_S{A}_D+X_S{B}_D\right)V_{Ri}\right]$$

(86)

Therefore, System (80) can be rewritten as

$$\{\begin{array}{ll}\frac{-1}{{\omega }_{s}J\left({\omega }_m^f-{\omega }_m^i\right)}[V_S{X}_m\left(2R_R{X}_S-B_D\right)& +2R_S{X}_m^2{V}_{Rr}\left]+\mu \right[\left(1+s\right)V_S{X}_m{B}_D-2\left(R_S{A}_D+X_S{B}_D\right)V_{Rr}]=0\\ \frac{-1}{{\omega }_{s}J\left({\omega }_m^f-{\omega }_m^i\right)}[V_S{X}_m\left(2R_S{R}_R-A_D\right)& +2R_S{X}_m^2{V}_{Ri}\left]+\mu \right[\left(-1+s\right)V_S{X}_m{A}_D-2\left(R_S{A}_D+X_S{B}_D\right)V_{Ri}]=0\\ & P_{liq}\left(V_{Rr},V_{Ri}\right)-P_{lim}=0\end{array}$$

(87)

The solution of the system is obtained analytically by manipulating these equations. Dividing the first one by the second leads to

$$\begin{array}{ll}2V_{Rr}[V_S{R}_S{X}_m^3{A}_D& \left(-1+s\right)+V_S{X}_m\left(2R_S{R}_R-A_D\right)\left(R_S{A}_D+X_S{B}_D\right)]\\ & -2V_{Ri}\left[V_S{R}_S{X}_m^3{B}_D\left(1+s\right)+V_S{X}_m\left(2R_R{X}_S-B_D\right)\left(R_S{A}_D+X_S{B}_D\right)\right]\\ & =V_S^2{X}_m^2{B}_D\left(2R_S{R}_R-A_D\right)\left(1+s\right)+V_S^2{X}_m^2{A}_D\left(2R_R{X}_S-B_D\right)\left(1-s\right)\end{array}$$

(88)

Thus, the Lagrange multiplier is eliminated, and an equation is obtained from which it is possible to explicit $V_{Rr}$ as a function of $V_{Ri}$. Entering the following constants:

$$K_{Vr}=2\left[V_S{R}_S{X}_m^3{A}_D\left(-1+s\right)+V_S{X}_m\left(2R_S{R}_R-A_D\right)\left(R_S{A}_D+X_S{B}_D\right)\right],$$

(89)

$$K_{Vi}=2\left[V_S{R}_S{X}_m^3{B}_D\left(1+s\right)+V_S{X}_m\left(2R_R{X}_S-B_D\right)\left(R_S{A}_D+X_S{B}_D\right)\right],$$

(90)

$$K_i=V_S^2{X}_m^2{B}_D\left(2R_S{R}_R-A_D\right)\left(1+s\right)+V_S^2{X}_m^2{A}_D\left(2R_R{X}_S-B_D\right)\left(1-s\right),$$

(91)

one obtains

$$V_{Rr}=\frac{K_{Vi}{V}_{Ri}+K_i}{K_{Vr}}.$$

(92)

Substitution of Equation (90) in constraint $g\left(V_{Rr},V_{Ri}\right)$ leads to

$$\begin{array}{ll}[-& (R_S{A}_D+X_S{B}_D)\left({\left(\frac{K_{Vi}{V}_{Ri}+K_i}{K_{Vr}}\right)}^2+V_{Ri}^2\right)+\left(1+s\right)V_S{X}_m{B}_D\left(\frac{K_{Vi}{V}_{Ri}+K_i}{K_{Vr}}\right)\\ & +\left(-1+s\right)V_S{X}_m{A}_D{V}_{Ri}+V_S^2\left(R_R{A}_D+s{X}_R{B}_D\right)]-P_{lim}\frac{A_D^2+B_D^2}{3}=0\end{array}$$

(93)

This polynomial has two roots, which represent two possible solutions for $V_{Ri}^{otm}$. Replacing these solutions into Equation (90), yields two solution pairs [$V_{Rr}^{otm},V_{Ri}^{otm}$], known as stationary (or critical) points of the problem.

The first analysis of the critical points is done by observing the signal of the Lagrange multiplier $\mu$, calculated by the following equation, taken from the system presented in Equation (85):

$$\mu =\frac{1}{{\omega }_{s}J\left({\omega }_m^f-{\omega }_m^i\right)}\left[\frac{V_S{X}_m\left(2R_R{X}_S-B_D\right)+2R_S{X}_m^2{V}_{Rr}^{otm}}{\left(1+s\right)V_S{X}_m{B}_D-2\left(R_S{A}_D+X_S{B}_D\right)V_{Rr}^{otm}}\right],$$

(94)

When the constraint is of inequality type, the sign of the Lagrange multiplier expresses the need to use the imposed constraint to limit the objective function. If the problem is maximization and the inequality constraint is of type “greater or equal”, a positive sign of the multiplier indicates that the constraint is naturally limiting the objective function because its gradients are opposite to ensure the first-order condition seen in Equation (77). If the problem is minimization, with the inequality constraint of type “less than or equal”, a negative sign of the Lagrange multiplier leads to the same conclusion. When, in the cases presented, the signs are opposite, it means that the objective function is being forced to go to the constraint limit, in which case it would be unwanted.

After applying the first-order condition and calculating the critical points, it is necessary to analyze the second-order condition and the Hessian matrix of the problem. The second-order condition is applied as shown below:

$${\nabla }^2ℒ\left(V_{Rr},V_{Ri},\mu \right)={\nabla }^2{f}_{\omega }\left(V_{Rr},V_{Ri}\right)+\mu {\nabla }^{2}g\left(V_{Rr},V_{Ri}\right),$$

(95)

and the Hessian matrix is presented below:

$$H=\left[\begin{array}{cc}\omega I-{\nabla }^2ℒ\left(V_{Rr}^{otm},V_{Ri}^{otm}\right)& \nabla g\left(V_{Rr}^{otm},V_{Ri}^{otm}\right)\\ {\nabla }^{\mathrm{T}}g\left(V_{Rr}^{otm},V_{Ri}^{otm}\right)& 0\end{array}\right].$$

(96)

The analysis of the eigenvalues of the Hessian matrix is necessary to define the convexity of the region around each critical point and to verify if they are minimum or maximum points [31]. The eigenvalues are obtained from the Hessian determinant as shown below:

$$det\left[\begin{array}{lll}{\omega }^2-\frac{\partial ℒ\left(V_{Rr}^{otm},V_{Ri}^{otm}\right)}{\partial V_{Rr}^2}& -\frac{{\partial }^2ℒ\left(V_{Rr}^{otm},V_{Ri}^{otm}\right)}{\partial V_{Rr}\partial V_{Ri}}& \frac{\partial gℒ\left(V_{Rr}^{otm},V_{Ri}^{otm}\right)}{\partial V_{Rr}}\\ -\frac{{\partial }^2ℒ\left(V_{Rr}^{otm},V_{Ri}^{otm}\right)}{\partial V_{Ri}\partial V_{Rr}}& \omega -\frac{{\partial }^2ℒ\left(V_{Rr}^{otm},V_{Ri}^{otm}\right)}{\partial V_{Ri}^2}& \frac{\partial g\left(V_{Rr}^{otm},V_{Ri}^{otm}\right)}{\partial V_{Ri}}\\ \frac{\partial g\left(V_{Rr}^{otm},V_{Ri}^{otm}\right)}{\partial V_{Rr}}& \frac{\partial g\left(V_{Rr}^{otm},V_{Ri}^{otm}\right)}{\partial V_{Ri}}& 0\end{array}\right]=0$$

(97)

If all the eigenvalues of the Hessian matrix, analyzed for a certain critical point, have a positive sign, then the region around the critical point is convex; therefore, this is a minimum. If negative, the critical point is a maximum of the objective function [32].

The objective function and the constraint function of the problem proposed by Equation (70) are both second-order polynomials. Thus, it is expected that the critical points found are one maximum and one minimum. Therefore, the critical point that solves this problem will be one that maximizes the objective function. Similarly, when the wind is decreasing, the problem becomes minimization, as shown below:

$$\{\begin{array}{c}\min f_{\omega }\left(V_{Rr},V_{Ri}\right)\\ s.t.P_{liq}\left(V_{Rr},V_{Ri}\right)\le P_{lim}\end{array},$$

(98)

All mathematical inferences presented throughout this section are valid for this formulation as well.

In the second-order condition, for the problem presented in Equation (70), we search for the critical point that generates the eigenvalue of the Hessian matrix, with negative sign, which maximizes the objective function. For the problem presented in Equation (96) the eigenvalue of the Hessian matrix should be positive.

At the end of the optimization process of $f_{\omega }$, we arrive at a pair $V_{Rr}^{otm}$ and $V_{Ri}^{otm}$. Therefore,

$$f_{\omega }^{otm}=\frac{T_m-T_e\left(V_{Ri}^{otm},V_{Ri}^{otm}\right)}{J\left({\omega }_m^f-{\omega }_m^i\right)}.$$

(99)

and

$${\frac{ds}{dt}}^{otm}=f_{\omega }^{otm}\left(s^f-s^i\right).$$

(100)

The optimal values of the derivatives of the control variables $V_{Rr}$ and $V_{Ri}$ at $t=0$ are determined from Equations (67) and (98).

$$\left[\begin{array}{c}{\frac{d{V}_{Rr}}{dt}}^{otm}\\ {\frac{d{V}_{Ri}}{dt}}^{otm}\end{array}\right]=\left[\begin{array}{c}{\xi }_{Vr}\\ {\xi }_{Vi}\end{array}\right]{\frac{ds}{dt}}^{otm}$$

(101)

As for the rotor speed, the exponential behavior for the convergence of control variables was also adopted. Therefore, generally,

$$V_{Rr}=\left(V_{Rr}^i-V_{Rr}^f\right)e^{-f_{vr}t}+V_{Rr}^f,$$

(102)

$$V_{Ri}=\left(V_{Ri}^i-V_{Ri}^f\right)e^{-f_{vi}t}+V_{Ri}^f,$$

(103)

where, $f_{vr}$ and $f_{vi}$ are the frequencies that define the convergence velocity, and $V_{Rr}^f$ and $V_{Ri}^f$ are the final values of the control variables (i.e., their reference values).

To find the optimal convergence frequencies $f_{vr}$ and $f_{vi}$, a relationship to the optimal velocity convergence frequency is established. Deriving $V_{Rr}$ and $V_{Ri}$:

$$\frac{d{V}_{Rr}}{dt}=-f_{vr}\left(V_{Rr}^i-V_{Rr}^f\right)e^{-f_{vr}t}$$

(104)

$$\frac{d{V}_{Ri}}{dt}=-f_{vi}\left(V_{Ri}^i-V_{Ri}^f\right)e^{-f_{vi}t}$$

(105)

Substituting Equations (102) and (103) into Equation (99) gives

$$-f_{vr}^{otm}\left(V_{Rr}^i-V_{Rr}^f\right)e^{-f_{vr}^{otm}t}={\xi }_{Vr}{\frac{ds}{dt}}^{otm}$$

(106)

$$-f_{vi}^{otm}\left(V_{Ri}^i-V_{Ri}^f\right)e^{-f_{vi}^{otm}t}={\xi }_{Vi}{\frac{ds}{dt}}^{otm}$$

(107)

Substituting Equation (98) into Equations (104) and (105) and analyzing them at $t=0$, one concludes

$$f_{vr}^{otm}=\frac{-{\xi }_{Vr}{f}_{\omega }^{otm}\left(s^f-s^i\right)}{\left(V_{Rr}^i-V_{Rr}^f\right)}$$

(108)

$$f_{vi}^{otm}=\frac{-{\xi }_{Vi}{f}_{\omega }^{otm}\left(s^f-s^i\right)}{\left(V_{Ri}^i-V_{Ri}^f\right)}$$

(109)

2.6. Speed Estimation and Magnetizing Inductance Correction

Magnetic saturation affects the inductance value of the DFIG, and consequently the control performance [23]. However, the effect of this saturation mainly influences the magnetizing inductance since the path of the magnetic leakage flux is mainly air; thus, the variations of the leakage inductances are irrelevant to the control [24]. Changes in inductance values may lead to a different operating point than expected. Still, considering that the magnetizing inductance is the parameter that has the highest value, compared to the other DFIG parameters, it is important to identify its value to ensure control efficiency. Because of these facts, a strategy is proposed here for estimating the value of the magnetizing inductance by means of measuring stator voltages and currents and estimating rotor slip and rotor currents. The proposed procedure considers that the identification of magnetizing inductance should not necessarily be performed at each step of the control. An error in its value would lead the DFIG to an erroneous steady-state operating point. This way, correction of the magnetizing inductance value should be performed when detecting the achievement of this condition, with differences in the expected values for variables such as, for example, reactive power. For this reason, the proposed strategy is based on the dynamic equations of DFIG, considering that, for steady state, one has $\frac{d\lambda }{dt}=0$. It is assumed that the d axis of the Park transform is aligned with the stator voltage, which means that $V_{Sd}=V_S$ and $V_{Sq}=0$, where $V_{Sd}$ and $V_{Sq}$ are the stator voltages at the d and q coordinates. From these considerations and from Equations (27) to (34), the problem can be solved as follows:

$$V_{Sd}=-R_S{I}_{Sd}-{\omega }_e\left(-L_S{I}_{Sq}+L_m{I}_{Rq}\right),$$

(110)

$$0=-R_S{I}_{Sq}+{\omega }_e\left(-L_S{I}_{Sd}+L_m{I}_{Rd}\right),$$

(111)

$$V_{Rd}=R_R{I}_{Rd}-s{\omega }_e\left(L_R{I}_{Rq}-L_m{I}_{Sq}\right),$$

(112)

$$V_{Rq}=R_R{I}_{Rq}+s{\omega }_e\left(L_R{I}_{Rd}-L_m{I}_{Sd}\right).$$

(113)

It should be noted that Equations (27) to (34) assume that the directions of stator and rotor currents are entering the DFIG. The stator current direction assumed for the mathematical equation of the strategy proposed here, as shown in Figure 2, comes out of the generator, so the stator current signals have been changed to suit the calculation.

From Equations (108) and (109), the value of the rotor currents can be estimated as follows:

$${\tilde{I}}_{Rd}=\frac{R_S{\widehat{I}}_{Sq}+{\omega }_e{\tilde{L}}_S{\widehat{I}}_{Sd}}{{\omega }_e{\tilde{L}}_m}$$

(114)

$${\tilde{I}}_{Rq}=\frac{-{\widehat{V}}_{Sd}-R_S{\widehat{I}}_{Sd}+{\omega }_e{\tilde{L}}_S{\widehat{I}}_{Sq}}{{\omega }_e{\tilde{L}}_m}$$

(115)

The values of the variables represented with “~” are estimated, whereas the variables represented with “^” have their values measured. Therefore, ${\tilde{I}}_{Rd}$ and ${\tilde{I}}_{Rq}$ are the estimated values of rotor currents, and ${\widehat{V}}_{Sd}$, ${\widehat{V}}_{Sq}$, ${\widehat{I}}_{Sd}$ and ${\widehat{I}}_{Sq}$ are the measured values of stator voltages and currents, respectively. ${\tilde{L}}_m$ and ${\tilde{L}}_S$ are the estimated values of the magnetizing inductance and stator self-inductance, respectively. Initially, their nominal values are considered for the calculation.

As previously mentioned, variations in the dispersion inductances are irrelevant in the control, so it is possible to state that variations in the DFIG magnetizing inductance are reflected in equal proportions in the rotor and stator proper inductances, as follows:

$$\frac{L_m}{{\tilde{L}}_m}=\frac{L_S}{{\tilde{L}}_S}$$

(116)

where $L_m$ and $L_S$ are the nominal values provided by the manufacturer. Substituting Equation (31) into (114) results in

$$\frac{L_{lS}+L_m}{L_m}=\frac{{\tilde{L}}_S}{{\tilde{L}}_m}$$

(117)

$$\left(1+k_{lS}\right){\tilde{L}}_m={\tilde{L}}_S$$

(118)

where,

$$k_{lS}=\frac{L_{lS}}{L_m}$$

(119)

Similarly to the rotor:

$$\left(1+k_{lR}\right){\tilde{L}}_m={\tilde{L}}_R$$

(120)

$$k_{lR}=\frac{L_{lR}}{L_m}$$

(121)

and, therefore, it is possible to modify Equations (112) and (113) to be exclusively dependent on ${\tilde{L}}_m$,

$${\tilde{I}}_{Rd}=\frac{R_S{\widehat{I}}_{Sq}+{\omega }_e\left(1+k_{lS}\right){\tilde{L}}_m{\widehat{I}}_{Sd}}{{\omega }_e{\tilde{L}}_m},$$

(122)

$${\tilde{I}}_{Rq}=\frac{-{\widehat{V}}_{Sd}-R_S{\widehat{I}}_{Sd}+{\omega }_e\left(1+k_{lS}\right){\tilde{L}}_m{\widehat{I}}_{Sq}}{{\omega }_e{\tilde{L}}_m}.$$

(123)

With the estimated values of rotor currents by Equations (120) and (121), it is possible to estimate the slip and magnetizing inductance by Equations (110) and (111), respectively, as shown below:

$$\tilde{s}=\frac{-V_{Rd}+R_R{\tilde{I}}_{Rd}}{{\omega }_e(\left(1+k_{lR}\right){\tilde{L}}_m{\tilde{I}}_{Rq}-{\tilde{L}}_m{\widehat{I}}_{Sq})},$$

(124)

$$\begin{array}{ll}{\tilde{L}}_m^2[\tilde{s}{\omega }_e^2{\widehat{I}}_{Sd}& -\tilde{s}{\omega }_e^2(1+k_{lr})(1+k_{ls}){\widehat{I}}_{Sd}]\\ & +{\tilde{L}}_m[{\omega }_e{V}_{Rq}-{\omega }_e{R}_R(1+k_{ls}){\widehat{I}}_{Sq}-\tilde{s}{\omega }_e(1+k_{lr})R_S{\widehat{I}}_{Sq}]+R_R({\widehat{V}}_{Sd}\\ & +R_S{\widehat{I}}_{Sd})=0,\end{array}$$

(125)

and the values of $V_{Rd}$ and $V_{Rq}$ are found by Equations (100) and (101).

The slip estimation process is performed in real time, and its estimated value is used in the optimization process presented in Section 2.4 and Section 2.5. The estimation of magnetizing inductance is not performed in real time, as the current transient would generate unrealistic values for ${\tilde{L}}_m$. The strategy proposes that the magnetizing inductance value should be calculated when the machine is in a steady state and the stator reactive power is not close to zero, considering a tolerance. Failure to achieve this objective of steady-state control indicates that the value of the magnetizing inductance used in calculating the reference values of the control variables is different from the actual value. It is expected that, under this condition, the estimated rotor and slip current values are wrong. Upon detecting this condition, the value of ${\tilde{L}}_m$ will be corrected using Equation (123), and the corrected value will be used for the calculation of the analytic solution and optimal strategy routines.

There are two values for ${\tilde{L}}_m$ possible for the solution of Equation (123). To achieve a less turbulent convergence toward the actual inductance value, the value closest to the value calculated in the previous step is always chosen.

3. Results

The DFIG RSC control methodology proposed in this article was analyzed by applying wind speed steps, as shown in Figure 4. This wind profile was adopted with the intention of demonstrating the effectiveness of the proposed method in low and high wind speeds, as well as near zero slip. The wind speed ranging from 5.5 to 10.0 m/s was adopted to respect the generator speed limits presented in Table 2, corresponding to a slip range of 30% to −30%, approximately. In order to ensure that the control assumptions are respected, the convergence of velocity to its reference value should be analyzed, as well as the behavior of the control variables.

For each wind speed, a reference rotor speed for the maximum wind power extraction was defined, as shown in Equation (4). The control variables $V_{Rr}$ and $V_{Ri}$ had their reference values calculated over time by Equation (38) and changed due to stator voltage variation, caused by transitory currents. The exponentials presented in Equations (100) and (101) defined the trajectory of $V_{Rr}$ and $V_{Ri}$ to their reference values. Figure 5 and Figure 6 show this behavior.

The optimization process was proposed to define an optimal rotor convergency speed and, by this, to calculate a trajectory for the control variables that respects some power limits. The strategy intended to accelerate the DFIG rotor by maintaining a difference between the mechanical torque and the electromagnetic torque, and this difference depended on the value that was chosen for the constant $\alpha$, as shown at Equation (75). The value of $\alpha$ has been set to 0.85 for rising speed wind and at 1.15 for falling speed wind, and the effect of this strategy can be observed by the torques in Figure 7.

For the effect of the distance between the torques, the DFIG rotor accelerated or deaccelerated to reach the speed corresponding to the maximum power extraction point operation. Figure 9 shows the rotor speed behavior and its convergence to the reference values. The value of the estimated rotor speed was accurate, and the difference between both curves can hardly be observed in Figure 8.

To verify the effectiveness of the proposed strategy at steady state, the developed, active and wind powers were analyzed. The behavior of the torques directly influenced the active power output, as can be analyzed from Figure 9.

The net active power moved away from its previous value, increasing or decreasing, because of the behavior of the electromagnetic torque, but it quickly converged to its steady-state value. The developed power had a similar behavior, and transiently, the value of $\lambda$ was not optimal. At convergence, the value of the developed power was equal to the value of the power delivered by the wind. Net active power was lower due to stator and rotor losses.

Similar to net active power, due to machine acceleration, the value of stator reactive power did not remain constant, as can be seen from Figure 10.

Since there was no voltage control, transient reactive power variations caused variations in stator voltage. Despite the reactive power peaks, the stator voltage remained within the limits imposed by the international grid codes, transiently [33], as shown in Figure 11.

The slip estimation strategy based on the rotor current estimation, using only the stator variables, proved to be efficient and accurate, as shown in Figure 8. In order to analyze the algorithm to correct the value of the magnetizing inductance online, the following situation was supposed. The value of this parameter, used to solve the Runge–Kutta method, was changed to 70% of its nominal value. Thus, control algorithms (composed of analytical solutions to the control variables), the optimal strategy for rotor speed convergency and estimation of the rotor currents and slip all became dependent on a value of the magnetizing inductance that was different from the real value at the DFIG. Since its value was corrected only in a steady-state situation, the control will perform transiently and give wrong information. To test the algorithm, a variation of wind speed from 7.5 to 8 m/s was used, and the result is shown in Figure 12.

As should be expected, there was an error between the actual rotor speed and its estimated value, caused by the wrong value of the magnetizing inductance. In addition, the rotor speed would not converge to its optimal reference value, since the reference values for the control variables were miscalculated. As can be observed in Figure 13, the first convergence value, approximately between 3.3 and 4.4 s, was very close to the optimal value. This means the tip-speed ratio was near its optimal value, and the effect of the wrong value of the magnetizing inductance had a small effect on tracking the maximum power point; however, the value of the reactive power at the stator was significantly affected, as shown in Figure 13. Between 3.3 and 4.4 s the reactive power consumed from the grid was about 280 kvar.

At 4.4 s, the proposed algorithm to correct the value of magnetizing inductance started. This algorithm rapidly converged to a value very close to the real one. The control variables were quickly corrected too, as shown in Figure 14, Figure 15 and Figure 16. Then, the rotor speed and the stator reactive power converged to its desired value, as shown in Figure 12 and Figure 13.

In order to test the proposed control in situations near to real conditions, a set of real measurement data of wind speed in a real farm was obtained. A selected range of 10 min, characterized by the mean value and standard deviation, was used to compose 600 points following a normal distribution. This variation is presented in Figure 17.

Observe that the variations incorporated ascendant and descendant ramps, differently from the previous simulations. The resulting speeds obtained by applying the new control strategy are sketched in Figure 18. Once more, the estimated rotor speed superposed the real speed, attesting the efficiency of the rotor speed estimation.

The efficiency of extracting power from wind may be evaluated by calculating the total optimal energy available and the total energy converted for the period considered in Figure 18, as suggested by [15]. Proceeding this way, the obtained value was about 98%, which was 5% higher than the best result presented in that reference. However, it is important to point out that different wind speed profiles were used in both simulations; the wind speed profile of Figure 17 better reflected the wind behavior from the specific site where the measurement data were obtained. The corresponding voltages were such that the difference between the maximal and minimal voltages by this simulation corresponded to a range of about 0.5%. No violations in variable limits were detected.

In order to give another idea about the efficiency of the proposed control strategy, the variation of the power coefficient with time is sketched in Figure 19. In this figure, a comparison between the optimal value, the mean value and the actual value of the power coefficient was registered for a time period of 300 s. The mean value corresponded to 99.77% of the optimal value, with a standard deviation of 0.19%.

4. Discussion

In previous works, objectives are set, and the control variables match their final values by approximation as objectives are achieved. In the present work, final values of the control variables (references) are calculated by means of an analytical procedure to solve a system of nonlinear equations. If the terminal voltage cannot be considered constant because of high equivalent Thévenin impedance, the reference rotor voltages must be recalculated using information acquired from a voltage sensor and a PLL.

A new DFIG control strategy was proposed and presented using a nonlinear, real-time optimization process, formulated from the equivalent circuit equations. In addition, a solution for the machine slip estimation and magnetizing inductance estimation for online correction of values was presented. The proposed control strategy performed the DFIG rotor acceleration optimally, respecting the limits of net power injected into the grid, by means of controlling the trajectory of the control variables that defined the RSC voltage. The electromagnetic and mechanical torques moved apart for a brief period to cause acceleration of the rotor, which caused the injected net power to rise or fall for a few seconds. Control has proven to be effective over a wide range of wind speeds, and the only parameter to be defined, in order to reach this purpose, was the net active power limit. Rotor currents were transiently controlled well, and reactive power peaks caused mild overvoltages or undervoltages well below the limits imposed by actual regulations. The strategy to estimate rotor currents, slip and the magnetizing inductance obtained from the dynamical formulation of DFIG has proven to be effective. It used only stator voltage and current sensors, being independent of rotor and speed sensors, thus reducing the number of sensors in the system.

Table 3. Comparison of proposed Direct Voltage Control with other methods.

Method	${\mathbf{L}}_{\mathbf{m}}\mathbf{Estimation}$	Rotor Sensors Dependency	Active Power Flow Inversion	Respect Voltage Limits	Convergency Time to 2 m/s Step	PI Gains Dependency for RSC Control
[8]	✕	✓	✕	✓	>10 s	✓
[9] (a *)	✕	✓	✓	Not considered	<1 s	✓
[9] (b *)	✕	✓	✕	Not considered	<5 s	✓
[10]	✕	✓	✓	Not considered	<1 s	✓
DVC	✓	✕	✕	✓	<5 s	✕

* (a): Rotor speed control; (b): Power control.

presents a synthetic comparation of the Direct Voltage Control (DVC), here proposed, with other methods found in the literature.

5. Conclusions

The major objective of the present work was to control a DFIG efficiently, without relying on parameters that can be corrupted over time like PI gains. An optimization process was proposed to define an optimal rotor convergency speed and thus define the trajectory of the control variables, respecting some power limits. The strategy intended to accelerate the DFIG rotor by maintaining a difference between the mechanical torque and the electromagnetic torque, and this difference depended on the allowed limit in net active power. This limit defines the momentaneous acceleration of the machine, in response to wind speed variations. When the wind speed changes, a new mechanical torque is imposed. The difference to the instantaneous electromagnetic torque, resulting in the acceleration torque, may not be enough to cause a fast convergence to the new reference speed. This way, the objective of optimization was to maximize acceleration, respecting limits on other variables.

Implementation of the optimization process is dependent on the instantaneous values of rotor speed. In order to avoid the necessity of using a speed sensor, a process of rotor speed estimation was proposed. As presented in Figure 8, the estimated speed resulted very close to the actual rotor speed, obtained by solving the differential equation set.

Simulation results demonstrated that too high acceleration torques (or powers) may cause violations both in currents or in voltage limits (see Figure 3). This way, optimization had to be constrained. This was implemented by limiting the net generated active power as a percentage of the previous reference power. For both wind speed profiles used (step change and near to real conditions) as example in this article, a variation of ±15% was adjusted so that no violations in currents or voltages magnitudes could have resulted. Considering all simulations, the highest resulting voltage amplitude was about 1.5% above the nominal value, as shown in Figure 11.

From the results of simulations performed for a real wind speed profile, it was observed (see Figure 18) that the proposed control was efficient in keeping the rotor speed close to the reference speed. Furthermore, the time variation of the power coefficient, as shown in Figure 19, presented an average very close to the optimal value.

Variations in the magnetizing inductance value have been shown to cause DFIG to operate at a different operating point than desired, implying small errors in rotor steady state velocity but a significant error in reactive power. The proposed algorithm quickly detected errors in the magnetizing inductance and corrected its value with good precision, causing the generator to work as expected, according to the right trajectory of control. The main advantages of the proposed method can be observed by comparing it other methods with respect to some desired features, as compiled in Table 3. From this table, it can be observed (third and fifth columns) that requiring fast convergence leads to a compromise in active power flow. It must be pointed out that oscillations in power flow may cause voltage fluctuations and, therefore, should be avoided. In general, the DFIG Direct Voltage Control with an optimum rotor acceleration strategy and parameter estimation proved to be effective in reaching the proposed objectives. The proposed technique is planned to be checked in an experimental system in a further work.