Miguel Jiménez
email: miguel.jimenez@lss.supelec.fr

Carrizosa " Gilney

Abdelkrim Damm
email: abdelkrim.benchaib@alstom.com

Pedro Benchaib

Mariana Alou

Françoise Netto

Lamnabhi-Lagarrigue

Bilinear and nonlinear control algorithms for a DC/DC converter for Multi-Terminal HVDC networks

Keywords: Control of bilinear systems, nonlinear control, DC/DC converters, Multi-Terminal HVDC

come   L'archive ouverte pluridisciplinaire

Bilinear and nonlinear control algorithms for a DC/DC converter for Multi-Terminal HVDC networks M. Jiménez Carrizosa * G. Damm * * A. Benchaib * * * P. Alou * * * * M. Netto † and F. Lamnabhi-Lagarrigue *

INTRODUCTION

The increasing development of electricity needs has brought current AC grids close to their saturation. In the same time, the increase of renewable sources, with their intermittent and variable nature, and their ever increasing penetration in the classic AC power networks makes this problem even more delicate. It is in this context where Multi-Terminal High Voltage Direct Current (MT-HVDC) grids have gained importance with respect to AC grids. Some of their advantages is that they do not have reactive power, so their transmission capacities are increased, they also produce fewer losses than AC grids, the weight and cost of components are smaller (fewer cables), and the efficiency is higher (see [START_REF] Arrillaga | Flexible Power Transmission. The HVDC Options[END_REF]; [START_REF] Asplund | Ultra high voltaje trasmission[END_REF]). It is clear than in MT-HVDC the DC/DC converters play a decisive role, since energy is generated at low voltages and highly efficient devices are needed to increase these voltages to higher values or to interconnect DC Grids at different voltage levels.

The control of DC/DC converters has been highly studied over the last decades, but especially in low voltage and low power applications. However the same philosophy of these techniques may be adapted for high or medium

This work is supported by WINPOWER project (ANR-10-SEGI-016).

voltage and power applications. These techniques can be addressed from different viewpoints. One of them is to use the average model in each switching cycle of the control variables (see [START_REF] Sira-Ramirez | Dynamical adaptive pulse-width-modulation control of dc-to-dc power converters: a backstepping approach[END_REF]). To achieve this average model, pulse width modulation (PWM) is commonly used. Whatever the control objective, the average model of DC/DC converters via PWM have a structure of bilinear system.

A system is called bilinear if it is described by linear differential equations in which the control inputs appear as coefficients. The study of bilinear systems began in the 60's as a gateway between the linear and nonlinear systems. A large amount of definitions and properties have been formulated until our days with respect to bilinear systems and their control, see: [START_REF] Mohler | Bilinear Control Processes[END_REF][START_REF] Mohler | Nonlinear systems. Applications to Bilinear Control[END_REF]; [START_REF] Landau | On the optimal regulator problem and stabilization of bilinear systems[END_REF]; [START_REF] Elliot | Bilinear Control Systems[END_REF]. However there are not yet general results to study their stability and controllability as in the case of linear systems. There exists only global proofs in the case of 2 nd order bilinear systems (see [START_REF] Koditschek | The controllability of planar bilinear systems[END_REF]; [START_REF] Barras | Controllability of two-dimensional bilinear systems[END_REF]; [START_REF] Ayala | Controllability of bilinear systems on the projective space[END_REF]; [START_REF] Koditschek | Stabilizity of 2 n d order bilinear systems[END_REF]). Due to this fact, it is very usual to treat bilinear systems as nonlinear systems and apply all the results of nonlinear control, since the theory of nonlinear control has global results on stability and controllability (see [START_REF] Isidori | Nonlinear control systems[END_REF]; [START_REF] Khalil | Nonlinear systems[END_REF]). In this paper, both a bilinear controller and a nonlinear controller based on Lyapunov theory are designed and compared.

The paper is outlined as follows: Section II presents the modeling of the dual active bridge which is the considered converter. Section III presents the control of a single DAB converter cell, by either a bilinear and a nonlinear controls. Section IV presents the simulations of these controllers. In Section V a comparison of both controllers is shown. Finally, in Section VI the conclusions are explained.

DAB

The dual active bridge, DAB (see figure 1), is a widely used DC/DC converter for high voltage (more than 2 kV) or medium (around 1-2 kV) applications. It is composed of two full bridge converters, each one connected to one side of the transformer, which provides galvanic insulation necessary to protect against unforeseen circumstances (see G. Ortiz et al. (2010)). An inductance L, which is used as an energy transfer device, is placed in series with the transformer. Usually in each bridge it can be used either two phases (high voltage applications) as three phases (medium voltage). In this work we have chosen the last, in order to apply the dq transformation. DAB is shown in figure 1, the objective of this paper is to implement a control algorithm taking into account its three state variables. Thereby improving the performance of the control compared to previous results in the literature as [START_REF] Jiménez | Dc transformer for dc/dc connection in hvdc network[END_REF].

Averaged dynamic model

Thanks to the well known technique Pulse Width Modulation (PWM), it is possible to obtain an averaged dynamic model from the system shown in figure 1 which is shown in figure 2, notice that all devices are in the primary side of the transformer, where n is the transformer ratio. This way of proceeding, in order to achieve an average model, is standard, as for example in [START_REF] Barragán-Villarejo | Dynamic modelling and control of a shunt-series power flow controller based on ac-link[END_REF]. It must be emphasized that using this technique, it appears higher switching losses, so its use for high voltage levels is not recommended. The converter can then be modeled as in Figure 2, where u 1a = u C1 m 1 cos(ωt + δ) and u 2a = u C2 m 2 cos(ωt), with 0 < m 1 < 1, 0 < m 2 < 1 and δ is the phase for AC left side. So using the dq transformation (see [START_REF] Park | Two reaction theory of synchronous machines[END_REF]), we obtain that:

                   u 1d = 3 2 u C1 m 1d -1 < m 1d = m 1 cos(δ) < 1 u 1q = - 3 2 u C1 m 1q -1 < m 1q = m 1 sin(δ) < 1 u 2d = 3 2 u C1 m 2d 0 < m 2d = m 2 < 1 u 2q = 0 (1)
where u 1d and u 1q , are the dq voltage components of the three phases voltage system u 1a , u 1b , and u 1c (AC left side), and in analogous form for u 2d and u 2q (AC right side).

It should be noted that u 2q = 0, and this is due to our choice of the voltage phase in the right side of the converter was zero. This is in concordance with the real behavior of the system, because when we work in dq frame, a phase lock loop (PLL) is necessary in order to generate a signal whose phase is related to the phase of an input signal (see [START_REF] Mendel | IEEE Transactions on Circuits and Systems II[END_REF]). Usually in electronic devices this technique is used to detect the frequency of the voltage and current and to detect the phase between them. So a voltage reference is needed, and if we select voltage u 2 as reference, its phase is zero, and consequently m 2q = 0, and therefore u 2q = 0 .

On the other hand, and taking into account the power equalities

u C1 I 1 = u 1a i La + u 1b i Lb + u 1c i Lc and u C2 I 2 = u 2a i La + u 2b i Lb + u 2c i 2c
, and carry out the Park's transformation we obtain that:

       I 1 = 3 2 m 1d i Ld + m 1q i Lq I 2 = 3 2 m 2d i Ld (2)
The system shown in figure 2 can then be described by:

                       d dt i Ld = 1 L -R • i Ld + ω • L • i Lq + 3 2 (m 1d • uc1 -m 2d • uc2) d dt i Lq = 1 L -ω • L • i Ld -R • i Lq + 3 2 m1q • uc1 d dt u C1 = - 1 C1 3 2 (m 1d • i Ld + m1q • i Lq ) - V1 -u C1 R1 d dt u C2 = n 2 C2 3 2 m 2d i Ld - u C2 -V2 • n R2 • n 2 (3) where i Ld , i Lq , u C1 , u C2 are state variables, m 1d , m 2d , m 1q are control variables and V 1 , V 2 are external variables.
The systems equations shown in (3) may be represented by:

d dt x = Ax + 5 i=1 u i B i x + Cu (4)
where

A =       - R L ω 0 0 -ω - R L 0 0 0 0 - 1 R1C1 0 0 0 0 - 1 R2C2       ; B 1 = 3 2     0 0 1 L 0 0 0 0 0 - 1 C1 0 0 0 0 0 0 0    
19th IFAC World Congress Cape Town, South Africa. August 24-29, 2014

B 2 = 3 2     0 0 0 - 1 L 0 0 0 0 0 0 0 0 n 2 C2 0 0 0     ; B 3 = 3 2     0 0 0 0 0 0 1 L 0 0 - 1 C1 0 0 0 0 0 0     B 4 = B 5 = 0 C =     0 0 0 0 0 0 0 0 0 0 0 0 0 1 R1C1 0 0 0 0 0 n R2C2     x =   i Ld i Lq u C1 u C2   u =    m 1d m 2d m1q V1 V2   

CONTROL OBJECTIVE

As shown in system equations (3), we have three variables on which we may act (m 1d , m 1q , m 2d ) and four state variables, in principle we may only control three of these state variables, leaving a free behavior in the other variable. The control objective is then to assure that each state variable which could be controlled tracks its reference.

Under-actuation

Since we may only arbitrarily control at most three state variables, we will now study whether these input variables we may control these three states. Controllability of general bilinear systems is still an open problem, there exist sufficient conditions for planar bilinear systems (scalar control) (see [START_REF] Koditschek | The controllability of planar bilinear systems[END_REF] [START_REF] Ayala | Controllability of bilinear systems on the projective space[END_REF]), but unfortunately there exits only partial criteria to assure controllability of the systems for the case n ≥ 3 (see [START_REF] Barras | Controllability of two-dimensional bilinear systems[END_REF]). Necessary conditions for bilinear controllability systems are shown in [START_REF] Mohler | Nonlinear systems. Applications to Bilinear Control[END_REF].

In our case the system is not fully actuated, so we can separate it into two parts, one actuated and one non-actuated. This separation will help us to study the controllability of the actuated part.

ẋ = ẋa ẋna = A a A 12 A 21 A na • x a x na + + 5 i=1 u i B ia B i12 B i21 B ina x a x na + C a C na u (5)
where x a denotes the actuated states, x na denotes the nonactuated state, A a denotes the actuated part of matrix A, A na denotes the non-actuated part of A and the same reasoning for B and C .

Three variables (m 1d , m 1q , m 2d ) are used for controlling three state variables, leaving the last free. The reasoning of the choice of the control variables is the following: since i Lq is closely linked with the reactive power consumed by the inductance in the AC step of the converter, and since this power decreases the transmitted active power, it is clearly a variable whose control is important and whose reference in most cases will be to make it as small as possible or even zero. Between the three state variables that remains (i Ld , u C1 and u C2 ), we will always chose a voltage as variable to control in order to guarantee the voltage stability of the converter. As we have chosen u C2 as phase reference for PLL, we will chose this variable as control objective, remaining the other two variables to pick one in order to complete the set. We will present results with both control objectives, that means, either we control the trio (i Ld , i Lq , u C2 .) or the trio (i Lq , u C1 , u C2 ). In both cases we can easily verify the necessary conditions explained above (see [START_REF] Mohler | Nonlinear systems. Applications to Bilinear Control[END_REF]) for bilinear controllability systems.

Bilinear system stability and control law

Bilinear systems' stability has been studied extensively over the past decades (see [START_REF] Mohler | Nonlinear systems. Applications to Bilinear Control[END_REF]). In [START_REF] Koditschek | Stabilizity of 2 n d order bilinear systems[END_REF] the stability of second order bilinear systems under certain conditions is proven. In [START_REF] Landau | On the optimal regulator problem and stabilization of bilinear systems[END_REF] a quadratic feedback control u(x):

u i = -α[B i x + c i ] T Px (6) with i = 1, • • • ,
m stabilizes the system presented in (3) in the origin, by means of the quadratic Lyapunov function, V = x T Px. Where P is positive definite and α > 0.

Indeed, if we derive the Lyapunov function, we obtain

V = x T [PA + A T P]x -2α m i=1 x T P(B i x + c i )x 2 (7)
where it is assumed that

   (B 1 x + c 1 ) T Px . . . (B m x + c m ) T Px    = 0 (8)
∀x such that x = 0 and x(PA + A T P)x ≤ 0 in order to guarantee V ≤ 0.

Since the system is not fully actuated, and we have separated our system into two parts, we may apply the control strategy presented in 3.2 to the controllable part.

Applying the quadratic feedback control explained above and taking into account section 3.1, we will only act over i Lq , u C1 and u C2 . So x a = [i Lq , u C1 ,u C2 ]. If xa is an equilibrium state of the controlled variables corresponding to an input ū, making the change of variables x a = xa +x a , and u = ū + ũ, in order to achieve that the state variables track their references, and taking into consideration that

B(u, x) = ( m i=1 u i B i ) x = ( n i=1 x i H i ) u, with a suited H i ∈ R n×m , we obtain that: B(u, x a ) = B(ū + ũ, xa + xa ) = (9) = B(ū, xa ) + B(ū, xa ) + B(ũ, xa ) + B(ũ, xa ) and consequently from (5): ẋa = Ãa • xa + B(ũ, xa ) + B na (ũ, x na ) + Ca • ū (10)
where B(ũ, xa ) and B na (ũ, x na ) ∈ R 5 → R 3 and:

Ãa = A a + m i=1 ūi B ai Ca = C a + n i=1
x ai H ai • (11) Proceeding in analogous form showed above, we are going to prove that the following control (12) stabilizes the system ( 10 and Ãa in our case has the following form if we control the trio i Ld , i Lq , u C2

Ãa =       - R L ω - ū2 √ 3 L √ 2 -ω - 1 R 1 C 1 0 n 2 ū2 √ 3 C 2 √ 2 0 - 1 R 2 C 2       (15)
which is Hurwitz ∀ ū2 ∈ R. Or if we control the trio i Lq , u C1 , u C2 then Ãa has the form:

Ãa =       - R L ū3 √ 3 L √ 2 0 - ū3 √ 3 C 1 √ 2 - 1 R 1 C 1 0 0 0 - 1 R 2 C 2       (16)
which is Hurwitz ∀ ū3 ∈ R. So the Lyapunov function V < λ Ṽ , with λ a positive constant, and consequently the control law shown in (12) stabilizes exponentially the system. (Note that for the other combinations of actuated variables its corresponding Ãa are always Hurwitz). Now remains to study what happens with the nonactuated variable. But in our case, due to it being a first order differential equation with a constant negative value in the coefficient of the state variable, and as we have controlled the remaining states, this free variable is also exponentially stable.2 Theorem 1. The converter described by (3), with arbitrary initial conditions, under control law shown in (12), will exponentially tracks its references.

Nonlinear system stability and control law

As discussed above, bilinear systems are a class of nonlinear systems, but bilinear systems paradoxically have not global stability specified results. Therefore, they are commonly dealt with non-linear control techniques. In this section we will use Lyapunov theory to develop another controller to our system.

If we take a glance on system (3), we observe that the variable i Lq may be controlled through the control variable m 1q . Something similar happens with the variable u C2 , which may be controlled by m 2d . So it is clear that both variables can be controlled independently. With respect to the other variables, i Ld and u C1 , the same reasoning explained in 3.1 is valid, that means, we may only control one of them, and the other has a free behavior. In the section above we have shown a control for the trio i Lq , u C1 and u C2 . The present section will develop a controller for the trio i Ld , i Lq and u C2 , although one could have chosen the other trio. Let's start with the following equation:

d dt i Lq = -ω • i Ld - R L • i Lq + m 1q √ 3 L √ 2 • u c1 (17)
Considering the following tracking error ĩLq = i Lq -īLq (18) with the associated Lyapunov function

Q i Lq = 1 2 (i Lq -īLq ) 2 (19)
Its time derivate is:

Qi Lq = ĩLq • -ω • i Ld - R L • i Lq + m 1q √ 3 L √ 2 • u c1 -iLq
(20) in order to achieve negative derivative of the Lyapunov function:

-ω • i Ld - R L • i Lq + m 1q √ 3 L √ 2 • u c1 -iLq = -α 1 • ĩLq (21)
with α 1 positive constant. So, the time derivative becomes:

Qi Lq = -α 1 • ĩ2 Lq (22)
which applying the Lyapunov theory one can guarantee its exponential stability.

The control law for this variable becomes:

m 1q = L √ 2 u C1 √ 3 ω • i Ld + R L • i Lq + iLq -α 1 • ĩLq (23)
If we apply the same reasoning to the equation 1 :

d dt u c2 = n 2 C 2 3 2 m 2d • i Ld - u C2 -V 2 • n n 2 R 2 (24)
and considering the following tracking error ũC2 = u C2 -ūC2 (25) with the associated Lyapunov function

Q u C2 = 1 2 (u C2 -ūC2 ) 2 (26)
we obtain the following control law which guaranteed the stability of the variable u C2 , where α 2 is a positive constant:

m 2d = C 2 √ 2 n 2 • i Ld √ 3 u C2 -V 2 • n R 2 C 2 -uC2 -α 2 ũC2 (27)
At this point, considering that m 1d is the only variable that remains to be defined, we will use it to stabilize i Ld proceeding in analogous form to the following equation:

d dt i Ld = - R L •i Ld +ω•i Lq + m 1d √ 3 L √ 2 •u c1 - m 2d √ 3 L √ 2 •u c2 (28)
and considering the following tracking error ĩLd = i Ld -īLd (29) with the associated Lyapunov function

Q i Ld = 1 2 (i Ld -īLd ) 2 (30)
1 In the following it is important to remark that in practice the three states can never remain equal to zero.

19th IFAC World Congress Cape Town, South Africa. August 24-29, 2014

we obtain the following control law which guaranteed the stability of the variable i Ld , where α 3 is a positive constant:

m 1d = L √ 2 u C1 √ 3 R L i Ld -ω • i Lq + m 2d √ 3 L √ 2 u C2 + iLd -α 3 ĩLd (31)
Now we will study the behavior of the free state in order to verify if and where it converges. To carry out this task, we will apply the well known theory of zero dynamics (see [START_REF] Isidori | Nonlinear control systems[END_REF] and [START_REF] Chen | Analysis of a control strategy for a multiterminal vsc-hvdc system[END_REF]). To analyze it, we first divide the state variables (considering the error variables) in two parts:

η = u C1 ξ = [ ĩLd , ĩLq , ũC2 ] t (32)
So the whole system may be considered as in the normal form:

η = f 1 (η, ξ, u) ξ = f 2 (η, ξ, u) (33)
with u = f 3 (η, ξ) (34) given the equations ( 23), ( 27) and ( 31).

So given u such that ξ → ξ ( ξ could be zero or not) then the behavior of the system (33) is governed by the differential equation:

η = f 1 (η, ξ) (35)
If we apply this result to system (3), the behavior of the free variable u C1 = η could be written as:

uC1 = k 1 u C1 + k 2 • u C1 + k 3 ( 36
)
where

k 1 = - 1 C 1 R ī2 Ld + ī2 Lq + 1 n 2 ūC2 -V 2 • n R 2 ūC2
(37)

k 2 = - 1 C 1 R 1 , k 3 = V 1 C 1 R 1 (38)
If we study the behavior of equation ( 36), calling f (u C1 ) = k1 u C1 + k 2 • u C1 + k 3 , and taking into account that in steady state by energy balance it is always true that:

ūC1 -V 1 R 1 ūC1 = R ī2 Ld + ī2 Lq + 1 n 2 ūC2 -V 2 • n R 2 ūC2
(39) then the derivate of f is always negative and consequently equation ( 36) is globally exponentially stable, so the whole system is stabilized by the control law u = f 3 (η, ξ) given by equations ( 23), ( 27) and (31). We could obtain an analogous result, if we had chosen u C1 instead of i Ld proceeding in a similar form. Theorem 2. The converter described by (3) under control laws ( 23), ( 27) and (31) (with α 1 , α 2 and α 3 positive constants) is globally exponentially stabilized to an equilibrium point given by its references.

The proof of theorem 2 is composed by the steps of the current subsection.

SIMULATIONS

The controllers studied in the previous sections are tested in MATLAB-Simulink.

Bilinear control

In this section it is shown the simulations of control for the trio i Lq , u C1 and u C2 . The parameter values of the simulated example are presented in table 1. In figure 3 the behavior of i Lq is shown. We observe as the variable tracks its references (100 A until 0.4 s and 85 A for the rest) satisfactorily. In figure 4 the behavior of u C1 and u C2 are shown. We observe as the u C1 tracks its reference (990 V until 0.5 s and 980 V for the rest) in an appropriate way, where the voltage has a small oscillation band, caused by the size of the capacitors. For u C2 , we observe as the variable tracks its references (10015 V until 0.8 s and 10005 V for the remaining) satisfactorily. Observing figure 4 it is clear that the power direction flows from left side to right side of the converter, so i Ld , the free variable, must be positive. This is in according with the simulations shows in figure 5

Nonlinear control

In this section it is shown the simulations for the control of the trio i Ld , i Lq and u C2 , leaving u C1 as free variable. The parameter values of the simulated example are the same as above. In figure 6 the behavior of i Ld and i Lq are shown. We observe as in figure 6a the variable tracks its references satisfactorily(-200 A until 0.8 s and -250 A for the remaining). In figure 6b the variable tracks its references satisfactorily(-15 A until 1.2 s and -45 A for the rest). In figure 7b the behavior of u C2 is shown. We observe as the variable tracks its references (9990 V until 1.6 s and 9995 V for the remaining) satisfactorily. We observe that due to the current has negative value and the voltage in capacitor C 2 is less than V 2 = 10000 V, the power direction flows from right side to left side of the converter. So u C1 (figure 7a), the free variable, must be greater than V 1 =1000 V. We can check it showing figure 7a.

COMPARISON

In view of the results, we can say that both controllers provide interesting results, and they are able to correctly stabilize the system. Due to its simplicity, bilinear control is easier to implement, and may be preferable for real applications. On the other hand, fully nonlinear approaches have a deeper established theory, with broader tools. For this reason it can be more interesting for future works.

CONCLUSIONS

The main objective of this paper is to provide a control for a DC/DC converter used in a multi-terminal DC grid in order to integrate renewable energy. In order to solve the problem, two viewpoints have been addressed and compared. Due to the intrinsic nature of this type of converters, the problem leads to a bilinear approach in a natural way, so a solution supported by bilinear theory has been presented. Due to the system being not completely actuated, a solution for the actuated subsystem has been developed. In a second step a nonlinear control has been developed. A detailed stability analysis of the zero dynamics shows that our system is exponentially stable. Both control algorithms provide suitable results in simulations, so both will be considered for subsequent test bed experimentations.

Fig. 1 .

 1 Fig. 1. Dual Active Bridge.

Fig. 2 .

 2 Fig. 2. DAB average model.

  ) in the desired point given by the references, by means of the quadratic Lyapunov function, Ṽ = xT Px. ũi = -α[B ai xa + cai + b nai x na ] T Px a (12) with i = 1, • • • , 3 where cai are the columns of the matrix Ca , b nai are the columns of matrix B na , and P is positive 19th IFAC World Congress Cape Town, South Africa. August 24-29, 2014definite symmetric matrix of appropriate dimension (in this case the identity matrix of order 3).Proof.-System (10) can be rewritten asẋa = Ãa • xa + m i=1 (ũ i [B ai xa + cai + b nai x na ]) (13)and taking into account that given column vectors a and b then a T P b = b T P a because P is symmetric. Applying this property:

Fig. 3 .

 3 Fig. 3. i Lq behavior (bilinear control).

Fig. 4 .

 4 Fig. 4. a) u C1 b) u C2 behavior (bilinear control).

Fig. 5 .

 5 Fig. 5. i Ld behavior (bilinear control).

Fig. 6 .

 6 Fig. 6. a) i Ld b) i Lq behavior (nonlinear control).

Fig. 7 .

 7 Fig. 7. a) u C1 b) u C2 behavior (nonlinear control).

Table 1 .

 1 Simulation values.

	n	Transformer ratio	11/100
	R total-i	Total resistance branch i	0.022 Ω
	R 1	Resistance source 1	0.001 Ω
	R 2	Resistance source 2	0.1 Ω
	V 1	Input voltage	1000 V
	V 2	Output voltage	10000 V
	L i	Inductance branch i	0.01 H
	C 1	Input capacitor	1 mF
	C 2	Output capacitor	20 µ F
	f	Switching frequency	1000 Hz

19th IFAC World Congress Cape Town, South Africa. August 24-29, 2014


Fernando Jaramillo-Lopez

Gilney Damm
email: gilney.damm@ibisc.fr

Godpromesse Kenne
email: gokenne@yahoo.com

Francoise Lamnabhi-Lagarrigue

Laia

Adaptive control scheme for maximum power point tracking of a photovoltaic system connected to the grid

An adaptive control scheme for maximum power point tracking of a single-phase single-stage photovoltaic system connected to the grid is presented.

The maximum power point depends on temperature and solar irradiance, ambient conditions that are time-varying and difficult to measure.

Two solutions are presented. Each solution derive an estimator that approximate three different parameters. These parameters are functions of solar irradiance and temperature. In this manner, we eliminate the necessity of climatic sensors. The first solution, uses an adaptive estimator that is able to estimate constant parameters, and the second one uses a sliding mode estimator that is capable of estimate time-varying parameters.

A complete analysis was done taking into account the nonlinearities showed by the closed-loop system. The Lyapunov redesign technique was used to derive a controller that gives globally asymptotically stable trajectories of the closed-loop system. Computer simulations are presented to compare the performance of both estimators and also to show the good performance of the controller.

I. INTRODUCTION.

New policies and regulations have been developed to face the growing energy needs and climate change. These facts have stimulated the interest on renewable energy sources. Solar photovoltaic is one of them. Photovoltaic systems (PVS) converts sunlight directly into electricity by means of a semiconductive process.

Grid-connected PVS usually consists of a photovoltaic panel or array and a power conditioning system (PCS).

In [START_REF] Soeren | A Review of Single-Phase Grid-Connected Inverters for Photovoltaic Modules[END_REF] and [START_REF] Xue | Topologies of Single-Phase Inverters for Small Distributed Power Generators: An Overview[END_REF], the authors analyze advantages and disadvantages of various single-phase grid-connected inverters topologies. The single-phase single-stage central inverter DC/AC topology was chosen in this work, because of its simplicity and low price.

The output power of a photovoltaic array is function of irradiance and temperature. To increase the efficiency of the overall system, PVS always needs to work in its maximum power point, to deliver the maximum amount of energy (Fig. 1). Hence, an algorithm that can follow these power changes is needed. This is the maximum power point tracking (MPPT) algorithm. There are several MPPT algorithms like perturb and observe, incremental conductance, extremum-seeking among others. They compute the value of the voltage corresponding to the maximum power point. This voltage is then used as a reference value in the controller. In [START_REF] Femia | Optimization of Perturb and Observe Maximum Power Point Tracking Method[END_REF] and [START_REF] Esram | Comparison of photovoltaic array maximum power point tracking techniques[END_REF] the authors present some improvements and comparisons between these algorithms.

Previous works on MPPT control of PVS were done by splitting the problem in two parts: a control capable of tracking the MPP and a control capable of deliver a sinusoidal current in phase with the grid voltage.

Solutions like feedback linearization and sliding mode techniques were applied to the first part.

Several solutions have been developed for the second part, like the P+ resonant controller [START_REF] Daniel | Stationary Frame Current Regulation of PWM Inverters With Zero Steady-State Error[END_REF] and the adaptive resonant controller [START_REF] Timbus | Adaptive Resonant Controller for Grid-Connected Converters in Distributed Power Generation Systems[END_REF]. However these solutions are linear based techniques. Approximations that doesn't take into account the nonlinear behavior of the PV array.

An interesting solution that uses the nonlinear model of PVS was presented in [START_REF] Meza | Passive P-Control of a Grid-Connected Photovoltaic Inverter[END_REF]. In this work, the MPPT and the unity power factor tasks are satisfied with a single passive control, based on the work of [START_REF] Sanders | Lyapunov-based control for switched power converters[END_REF]. MPP tracking is accomplished for time-varying environmental conditions, with the help of an estimator that compute a reference for the controller. However, global stability with the use of this estimator is not proved.

In [START_REF] Meza | Lyapunov-Based Control Scheme for Single-Phase Grid-Connected PV Central Inverters[END_REF] an attempt to prove global stability for the whole system was made. In this case, an estimation of the parameter that depends on irradiance and temperature is used together with the Lyapunov-based controller to prove global stability. However, in that work it is assumed that another two parameters that depends on temperature are known.

In this paper, we are presenting two solutions that uses an adaptive scheme control that is capable of achieving MPP tracking for changing environmental conditions and deliver unity power factor current to the grid. In the first solution, an adaptive estimator is designed, and in the second one, a sliding mode estimator capable of estimate time-varying parameters. Both, estimate parameters that depends on irradiance and temperature. No assumption about known parameters was made and global stability is proved in both cases, taking into account the nonlinear model and the estimates of the parameters. The main contribution consists in the full analysis -nonlinear model, controller and estimators-of the system, increasing of the robustness and elimination of the necessity of irradiance and temperature sensors with these estimators.

II. PROBLEM FORMULATION AND MATHEMATICAL

MODEL OF THE PHOTOVOLTAIC SYSTEM.

The system analyzed in this paper is shown in figure 2 and consists in a PV array of solar cells and a DC/AC PCS. The states of the system are given by the capacitor voltage (x 1 ) and the inductor current (x 2 ). The PCS is a full-bridge inverter driven by a bipolar PWM scheme. The PWM gives two discrete complementary signals s and s which turn on and off the four switches in the PCS. The PWM block is fed by the control signal u ∈ [-1, 1]. A. Mathematical model.

The PV array is composed by PV cells arranged in series and parallel. The PV cell model used in this work is the single-diode model with no resistors. In [START_REF] Gow | Development of a photovoltaic array model for use in power-electronics simulation studies[END_REF] and [START_REF] Saloux | Explicit model of photovoltaic panels to determine voltages and currents at the maximum power point[END_REF], the authors present some models and give explicit relations to get the electric characteristics which are functions of irradiance and temperature. A comparative analysis was made in [START_REF] Bennett | Photovoltaic model and converter topology considerations for MPPT purposes[END_REF], and it is shown that all models have no significant differences for MPPT purposes.

The I -V characteristic curve of a PV array with identical cells is given by:

I = I ph (G)N p -I o (T )N p exp qV nN s k B T -1 (1)
where I ph is the photocurrent, that depends on irradiance G, I o is the saturation current, that depends on temperature T , q is the absolute value of electron's charge, n is the quality factor of the diode, N p is the number of cells connected in parallel, N s is the number of cells connected in series, k B is Boltzmann's constant and T is the temperature of the P-N junction. V is the capacitor voltage x 1 . This equation can be written in a simpler way:

I(x 1 ) = c 1 -c 2 e c3x1 (2)
where c 1 is function of irradiance and temperature and c 2 and c 3 are functions of temperature.

The model of the whole system is given by:

ẋ = f (t, x) + G(t, x)u (3)
where

f (t, x)   I(x1) C -vg(t) L   , G(t, x)   -x2 C x1 L   (4)
C and L are known values of the capacitor and the inductor respectively, and v g the grid voltage. The signals that are measured are the states x ∈ R 2 , the voltage grid v g ∈ R and the PV array current I ∈ R, which is a common practice in this type of circuits. u is the control signal composed by the terms:

u = u n + w + δ(t, x) (5)
δ(t, x) is an uncertain term that satisfies the matching condition (i.e. it enters at the same point that the control signal u), w is the term that will be derived to compensate it and u n is the control component that turn the nominal system globally asymptotically stable. In section IV, these components are explained in detail.

The nominal system is the system without the uncertain term δ(t, x), hence:

ẋ = f (t, x) + G(t, x)u n (6)
where f (t, x) and G(t, x) are defined in the same way that before.

The main tasks that the signal control u must fulfill are: 1) To track the maximum power point of the PV array, despite of changes in the environmental variables irradiance and temperature. 2) To deliver a current in phase with the grid voltage (i.e.

unity power factor). The controller needs two references to accomplish the main tasks of the system. The first reference is the voltage DC value in the capacitor x1 * ∈ R + and is given by the MPPT algorithm. The second reference is the current value in the inductor x 2 * , which is taken from [START_REF] Meza | Lyapunov-Based Control Scheme for Single-Phase Grid-Connected PV Central Inverters[END_REF]:

x 2 * 2 v g (t) x1 * I * (x 1 * , c i ) A 2 (7)
where I * is the PV current (2) evaluated at x1 * and v g is assumed to be sinusoidal with constant amplitude A and frequency ω, v g = Asin(ωt).

The parameters c i in I * are unknown. Therefore, two different estimators were designed and presented here. They allow to calculate the full I * expression.

Hence, the second reference becomes:

x2 * 2 v g (t) x1 * Î * (x 1 * , ĉi ) A 2 (8)
where (•) indicate the estimate of (•).

III. ESTIMATORS DESIGN.

A reparametrization was made in order to be able to design the next two estimators.

From ( 2) and ( 6) the derivative of the current I is:

İ = -c 2 c 3 ẋ1 e c3x1 = 1 C (ux 2 -I)(θ 1 -θ 2 I) = Φ T θ (9)
where:

θ c 1 c 3 c 3 , Φ T 1 C ux 2 -I, -ux 2 I + I 2 (10)
also consider the following estimator and observer errors respectively, θ = θθ, Ĩ = Î -I for the next subsections.

A. Adaptive estimator.

An adaptive estimator [START_REF] Jean- | Applied Nonlinear Control[END_REF] capable of estimate constant parameters is designed in this subsection.

The estimator is given by:

İ = -λ( Î -I) + Φ T θ (11) θ = -ΓΦ( Î -I) (12)
where λ > 0 and Γ is a square diagonal matrix 2 x 2 with positive values. Let us assume the following: Assumption 1-A. The parameters vector θ is constant, i.e. θ = 0 Lemma 1. Consider the observer [START_REF] Saloux | Explicit model of photovoltaic panels to determine voltages and currents at the maximum power point[END_REF] and the estimator [START_REF] Bennett | Photovoltaic model and converter topology considerations for MPPT purposes[END_REF], satisfying assumption 1-A, then the product Φ T θ converges to zero. Proof. Let's first rewrite the estimation errors as:

İ = -λ Ĩ + Φ T θ (13) θ = -ΓΦ Ĩ (14)
Consider the following Lyapunov candidate function

V = 1 2 | Ĩ| 2 + 1 2 θT Γ -1 θ β( Ĩ, θ) (15)
Its time derivate is

V = -λ| Ĩ| 2 + ĨT Φ T θ -θT Φ Ĩ = -λ| Ĩ| 2 (16)
By Barbalat's lemma, we can conclude that:

Ĩ → 0 as t → ∞
Therefore, it is clear by ( 13) that Φ T θ → 0 as t → ∞ Furthermore, if the regressor Φ(t) satisfies the persistent excitation condition [START_REF] Kumpati | Stable Adaptive Systems[END_REF]:

t+T0 t Φ(τ )Φ T (τ )dτ ≥ αI m ∀ t ≥ t 0 ( 17
)
where I m is the identity matrix and t 0 , T 0 and α are positive constants, then we can in addition conclude that θ → 0 as t → ∞ .

B. Sliding mode estimator.

In [START_REF] Kenne | Nonlinear systems time-varying parameter estimation: Application to induction motors[END_REF] a new framework for nonlinear systems timevarying parameter estimation using sliding mode techniques was proposed.

The following notation is used in this subsection:

|x T | G = (|x 1 |, |x 2 |, ... |x n |)
diag(A) is the column vector whose elements are the diagonal elements of a given square diagonal matrix A.

The function sign(•) : R n → R n is defined as

sign(x T ) = (sign(x 1 ), sign(x 2 ), ... sign(x n )) (18)
with

x T = (x 1 , x 2 , ... x n ) and sign(x i )      = -1 for x i < 0 = [-1, 1]
for

x i = 0 = 1 for x i > 0
Now, let us assume the following: Assumption 1-B. The derivatives of the parameters vector θ are bounded and these bounds are known: θi (t) ≤ µ θi where µ θi are known positive numbers and i = 1, 2. In order to design the time-varying parameter estimator consider the following adaptive observer for the current I:

İ = v I + Φ T θ (19)
where

v I = -K I sign( Î -I) (20)
thus İ = v I + Φ T θ. For K I sufficiently chosen large

K I > |Φ T θ| (21)
and assuming that θ is bounded (the proof will be shown later) then a sliding mode regime occurs on the manifold Ĩ = 0 and 0 = v Ieq + Φ T θ then

v Ieq = -Φ T θ (22)
The following approximation is used (see the work of Utkin [START_REF] Utkin | Sliding Modes in Optimization and Control[END_REF]):

v Ieq = 1 1 + τ s v I ( 23
)
where s is the Laplace operator and τ → 0 is a positive constant. Hence θ = -(ΦΦ T ) -1 Φv Ieq (24)

and finally θ = -K θ sign( θ)

Let us choose the gain matrix for the estimator such that:

diag(K θ ) i > µ θi (26)
Lemma 2. Consider the observer given by ( 19) and (20); the estimator given by ( 23)-( 25); satisfying assumption 1-B, ( 21) and (26); then the estimated parameters θ will converge to their real values θ(t) in finite time.

Proof. Consider the following Lyapunov candidate function

V = 1 2 θT θ β( θ) (27)
Its time derivative is:

V = θT θ = -| θT | G diag(K θ ) -θT θ ≤ -| θT | G (diag(K θ ) -µ θ ) (28)
and due to the inequality (26), V is negative-definite for θ = 0. Therefore the estimation error θ will converge to 0 in finite time.

C. Persistent excitation condition.

It is not necessary that the regressor Φ(t) satisfies the persistent excitation condition [START_REF] Khalil | Nonlinear Systems 2nd ed[END_REF], because the estimates are just used in the calculation of I * (x 1 * , c i ). The only necessary condition is that the estimators approximate the current I by means of its dynamics: İ( Ĩ, Φ T θ) → İ(x 1 , c i ) and it is satisfied if:

Φ T θ → Φ T θ IV. CONTROLLER DESIGN.
In order to accomplish the control objectives, the problem has been divided in two parts:

1) To find the control portion that renders the nominal system (6) globally asymptotically stable (GAS). 2) To find the control portion that compensates the uncertain term δ(t, x) and gives the whole system (3) GAS.

A. Controlling the nominal system.

In this part of the analysis it is considered that the x 2 * reference is well known, i.e. Ĩ(t) = 0.

According to the idea given in [START_REF] Sanders | Lyapunov-based control for switched power converters[END_REF], the control is decomposed in two parts: one for the steady-state stage u * n and another for the dynamic stage e un .

The reference system is:

ẋ * = f (t, x * ) + G(t, x * )u * n (29)
Hence, the control component for the steady state is:

u * n = L ẋ2 * + v g (t) x 1 * (30)
Consider the tracking errors and the variable control:

e = x -x * , e un = u n -u * n ė = ẋ -ẋ * = I(x1)-I * C 0 + -e2u * n C -(e2+x2 * )eun C e1u * n L + (e1+x1 * )eun

L

Let us propose the Lyapunov candidate function

V n = C 2 e 2 1 + L 2 e 2 2
Its time derivative is:

Vn = Ce 1 ė1 + Le 2 ė2 = e 1 (I -I * ) -e 1 [
From the fact that the first term of the right side of (32) is < 0 ∀ e 1 = 0, the nominal system is GAS.

B. Controlling the whole system.

The control part of the nominal system was established in (30) and (31). These equations are functions of ẋ2 * and x 2 * , and these values are not exactly known. Therefore, the control for the whole system is

u = u * n (x 2 * ) + e un (x 2 * ) + ∆u * + ∆e u = u * n + e un + δ(t, x) (33)
where

δ(t, x) ∆u * + ∆e u = L∆ ẋ2 * x 1 * + e 1 ∆x 2 * (34)
is the uncertain term, and ∆ ẋ2 * , ∆x 2 * are the perturbations due to the transients in the observer, i.e. Ĩ = 0.

In order to derive the controller for the whole system the next assumption is needed. Assumption 2. The uncertain term is bounded and known, i.e.

||δ(t, x)|| 2 < ξ
where ξ is a known positive number. Since this uncertain term satisfies the matching condition, the Lyapunov redesign technique [START_REF] Khalil | Nonlinear Systems 2nd ed[END_REF] is used to derive a term w that will compensate it.

u = u n + δ(t, x) + w (35)
where

u n = u * n + e un
The closed-loop system becomes

ẋ = f (t, x) + G(t, x)u n + G(t, x)[w + δ(t, x)] (36)
and

ė = f (t, x) + G(t, x)u n -ẋ * + G(t, x)[w + δ(t, x)]
Consider the Lyapunov candidate function for the whole system:

V = C 2 e 2 1 + L 2 e 2 2 + β (37)
where β was defined for each estimator in ( 15) and ( 27).

Let us omit the arguments of the functions. The time derivative of V is:

V = ∂V ∂e {f + Gu n -ẋ * } + ∂V ∂e {G[w + δ]} + β = Vn + β + ∂V ∂e Gw + ∂V ∂e Gδ (38)
We can derive the w term as follows:

∂V ∂e Gw + ∂V ∂e Gδ ≤ ∂V ∂e Gw + ξ ∂V ∂e G
2 and now we can choose:

w -ξ ∂V ∂e G || ∂V ∂e G|| 2 = -ξ -e 1 x 2 + e 2 x 1 | -e 1 x 2 + e 2 x 1 | (39)
Theorem 1. Consider the system given by ( 3)-( 4); the estimator given by ( 11)-( 12)|( 19)-( 25); the controller given by ( 30), ( 31) and (39); satisfying assumptions 1-A|1-B and 2, then, the whole system is globally asymptotically stable.

Proof. Using (39) one can write:

∂V ∂e Gw + ∂V ∂e Gδ ≤ -ξ ∂V ∂e G 2 + ξ ∂V ∂e G 2 = 0 (40)
then, from ( 16), ( 28), (32), ( 38) and (40), V is negative definite and the whole system is GAS with the controller given by u n + w.

V. SIMULATION RESULTS.

Numerical simulations were made in the Simulink/Matlab platform to verify the performance of the estimators and the controller. Two scenarios were simulated. In the first one, realistic variations were applied to the irradiance signal and a ramp to the temperature signal, for both estimators. Graphs a) and b) in figure 3 show these variations.

The adaptive estimator is not able to track step changes since it was designed for constant parameters, and the second scenario was simulated just for the sliding mode estimator. It consists in step changes in both variables, irradiance and temperature. Its values were changed by 50%. This is one of the worst conditions for MPPT in PVS, and because of this, a standard simulation scenario.

For both simulations, the initial conditions in the plant were zero and in the estimators 1x10 -5 .

Figures 3 and4 show the results for the first scenario. The graphs c) and d) in figure 3 show the first reference x 1 * , given by the MPPT algorithm, and the capacitor voltage Graphs c) and d) in figure 4 show the output power in both estimators. The results are similar. This is the maximum power that the PVS can deliver with the irradiance and temperature given.

Figure 5 shows the results for the second scenario. The step changes in irradiance and temperature for the sliding mode estimator are shown in graphs a) and b) in figure 5. The values were changed by 50%.

The graphs c), d) and e) in the same figure show that the controller and the estimator present a good performance, even under these demanding conditions. Some noticeable ripple appears after the step variations, due not only to the variation speed but also to the variation level (50%).

The output power is shown in graph e) of the same figure. It clearly follows the environmental variations.

VI. CONCLUDING REMARKS AND FUTURE RESEARCH.

Two solutions for the MPPT of PVS were presented. For the first adaptive scheme control (ASC), an adaptive estimator was designed. It is capable of estimate constant parameters. For the second ASC, a sliding mode estimator was designed. It is capable of estimate time-varying parameters.

A Lyapunov function that proves GAS of the system was derived. The analysis includes the dynamics of the estimators and the uncertainty in the second reference x 2 * . Even when GAS was demonstrated for the system including estimators and perturbations, there is one subject that can enrich the analysis, it is the inclusion of the MPPT algorithm in the closed-loop system analysis. We are currently investigating how to achieve this. Numerical simulations were made to verify the performance of the solutions. Both ASC have good performance under realistic conditions. The sliding mode ASC has the better performance. It can works well even under very demanding conditions.

2 Fig. 1 .

 21 Fig. 1. MPP of a typical PV array for different environmental conditions. a) Voltage vs. current curve at T = 25 o C and different irradiance values. b) Voltage vs. power curve at G = 1000W/m 2 and different temp. values.

Fig. 2 .

 2 Fig. 2. Photovoltaic system analyzed in this paper.

e 2 u

 2 * n + (e 2 + x 2 * )e un ] +e 2 [e 1 u * n + (e 1 + x 1 * )e un ] = -c 2 e 1 (e c3x1 -e c3x1 * ) -e un (e 1 x 2 * -e 2 x 1 * ) Consider the controller for the error dynamics of the nominal system: e un e 1 x 2 * -e 2 x 1 * (31) thus Vn = -c 2 e 1 (e c3x1 -e c3x1 * ) -(e 1 x 2 * -e 2 x 1 * ) 2

x 1 Fig. 3 .

 13 Fig. 3. a) Irradiance signal b) Temperature signal c) x 1 * and x 1 signals in the adaptive estimator d) x 1 * and x 1 signals in the sliding mode estimator

Fig. 4 .

 4 Fig. 4. a) I * and Î * signals in the adaptive estimator b) I * and Î * signals in the sliding mode estimator c) Theoretical maximum power and output power in the adaptive estimator case d) Theoretical maximum power and output power in the sliding mode estimator case

Fig. 5 .

 5 Fig. 5. Step changes and results got in the numerical simulation for the sliding mode estimator. a) Irradiance signal b) Temperature signal c) x 1 * and x 1 signals d) I * and Î * signals e) Theoretical maximum power and output power

Acknowledgements

The work of F. Jaramillo has been supported by CONACyT, Mexico and the French Embassy in Mexico. Also, the authors give thanks to Dr. R. Ortega for his contribution with the adaptive estimator.

* This research has received partially support from European Union Seventh Framework Programme [FP7/2007-2013] under grant agreement no.257462 HYCON2 Network of excellence. * F. Lamnhabi and F.


Fernando Jaramillo-Lopez

Godpromesse Kenne
email: gokenne@yahoo.com

Gilney Damm
email: gilney.damm@ibisc.fr

Francoise Lamnabhi-Lagarrigue

Adaptive control scheme for grid-connected photovoltaic systems with unknown bounds

Keywords: Adaptive control, Parameter estimation, Global stability, Solar cells, Chattering

An adaptive control scheme for maximum power point tracking of a single-phase grid-connected photovoltaic system is presented. The difficulty on design a controller that may operate a photovoltaic system on its maximum power point (MPP) is that, this MPP depends on temperature and solar irradiance, ambient conditions that are time-varying and difficult to measure. A solution using an on-line sliding mode estimator is presented. It estimates three different parameters that depends on solar irradiance and temperature, eliminating the necessity of having any sensor for these environmental variables. It is capable of estimate time-varying parameters. A complete analysis was done taking into account the non-linearities showed by the closed-loop system. A variation of the Lyapunov redesign technique was used to find a controller that give globally asymptotically stable trajectories of the closed-loop system. An adaptive law was found to substitute a perturbation bound and also to eliminate possible chattering due to the discontinuous controller term. Computer simulations are presented to show the good performance of the controller.

INTRODUCTION.

The compromising situation of the environment due to the pollution, and the high costs of fossil fuels have originated new policies and regulations that have stimulated the interest on alternative energy sources. Developing countries have increased in an important way the penetration of these energy sources (IEA [2012]).

One of the most widely used renewable energy systems are the photovoltaic systems (PVS), which converts sunlight directly into electricity by means of a semiconductive process. PVS usually consists of a photovoltaic panel or array, a power conditioning system (PCS), and occasionally, an energy storage element. Most of the grid-connected PVS (GPVS) consist of the first two elements.

Two of the topologies used for low and medium power levels are the two-stage DC/DC + DC/AC, and the single-stage DC/AC converter. In [START_REF] Jain | A Single-Stage Grid Connected Inverter Topology for Solar PV Systems With Maximum Power Point Tracking[END_REF], the authors show the main advantages of the single-stage against the two-stage topology. In Kjaer et al. [2005] and Xue et al. [2004], the authors analyze advantages and disadvantages of various single-phase grid-connected This research has received partially support from European Union Seventh Framework Programme [FP7/2007[FP7/ -2013] ] under grant agreement no.257462 HYCON2 Network of excellence.

inverters topologies. The single-phase single-stage central inverter DC/AC topology was chosen in this work, because of its simplicity and low price.

The output power of a PV array is function of irradiance and temperature. So, changes in these environmental conditions will bring different power values. To increase the efficiency of the overall system, PVS always needs to work in its maximum power point (MPP), to deliver the maximum amount of energy (Fig. 1). Hence, an algorithm that can follow these power changes is needed. This is the maximum power point tracking (MPPT) algorithm. There are several MPPT algorithms like perturb and observe, incremental conductance, extremum-seeking among others. They compute the value of the voltage corresponding to the MPP. This voltage is then used as a reference value in the controller. In Femia et al. [2005] andEsram et al. [2007], the authors present some improvements and comparisons between these algorithms.

Previous works on MPPT control of PVS were done by splitting the problem in two parts: a control capable of tracking the MPP and a control capable of deliver a sinusoidal current in phase with the grid voltage -unity power factor-.

Solutions like feedback linearization and sliding mode techniques were applied to the first part. Several solutions have been developed for the second part, like the P+ resonant controller (see Zmood et al. [2003]) where a better alternative to the PI controller is introduced, and the adaptive resonant controller (Timbus et al. [2006]) which is an improve of the proportional resonant controller. However, these solutions are linear based techniques. Approximations that doesn't take into account the nonlinear behavior of the PV array.

An interesting solution that uses the nonlinear model of PVS was presented in [START_REF] Meza | Passive P-Control of a Grid-Connected Photovoltaic Inverter[END_REF]. In this work, the MPPT and the unity power factor tasks are satisfied with a single passive control, based on the work of [START_REF] Sanders | Explicit model of photovoltaic panels to determine voltages and currents at the maximum power point[END_REF]. MPPT is accomplished for time-varying environmental conditions, with the help of an estimator that compute a reference for the controller. However, global stability with the use of this estimator is not proved.

In Meza et al. [2012] an attempt to prove global stability for the whole system was made. In this case, an estimation of the parameter that depends on irradiance and temperature is used together with the Lyapunov-based controller to prove global stability. However, in that work it is assumed that another two parameters that depends on temperature are known, so, a temperature sensor and a very precise knowledge of the part of the model that depends on temperature are needed.

In El Fadil and Giri [2011] the authors solved the climatic sensorless problem. The main difference with the present work is the application: in that work they use the PVS to charge a battery. The objective of deliver a current in phase with the grid voltage is not present. Another important difference is that the analysis was made considering the time-varying parameters as constants.

In this paper, we are presenting a solution that uses an adaptive scheme control that is capable of achieving MPPT for changing environmental conditions and deliver unity power factor current to the grid. In this solution, a sliding mode estimator capable of estimate time-varying parameters was designed and it estimates parameters that depends on irradiance and temperature. No assumption about known parameters were made and global stability taking into account the nonlinear model and the estimates of the parameters was proved. The main contribution consists in the full analysis -nonlinear model, controller and estimator-of the system, increasing of the robustness and elimination of the necessity of irradiance and temperature sensors with this estimator.

PROBLEM FORMULATION AND

MATHEMATICAL MODEL.

The system analyzed in this paper is shown in figure 2 and consists in a photovoltaic array of solar cells and a DC/AC PCS. The states of the system are given by the capacitor voltage (x 1 ) and the inductor current (x 2 ). The PCS is a full-bridge inverter driven by a bipolar PWM scheme. The PWM gives two discrete complementary signals s and s which turn on and off the four switches in the PCS. The PWM block is fed by the control signal u ∈ [-1, 1].

Mathematical model.

The PV array is composed by PV cells arranged in series and parallel. The PV cell model used in this work is the single-diode model with no resistors. In Gow et al. [1999] and Saloux et al. [2011], the authors present some models and give explicit relations to get the electric characteristics which are functions of irradiance and temperature. In [START_REF] Aggarwal | Solar cell array parameters[END_REF], the authors studied the PV array modeling problem for nonidentical PV cells. A comparative analysis was made in Bennett et al. [2012], and it is shown that all models have no significant differences for MPPT purposes.

The I -V characteristic curve of a PV array with identical cells is given by:

I = I ph (G)N p -I o (T )N p exp qV nN s k B T -1 (1)
where I ph is the photocurrent, that depends on irradiance G, I o is the saturation current, that depends on temperature T , q is the absolute value of electron's charge, n is the quality factor of the diode, N p is the number of cells connected in parallel, N s is the number of cells connected in series, k B is Boltzmann's constant, T is the temperature of the P-N junction and V is the capacitor voltage x 1 .

This equation can be written in a simple way:

I(x 1 ) = c 1 -c 2 e c3x1 (2)
where c 1 is function of irradiance and temperature and c 2 and c 3 are functions of temperature.

The model of the whole system is given by: ẋ

= f (t, x) + G(t, x)u (3) where f (t, x)      I(x 1 ) C -v g (t) L      , G(t, x)     -x 2 C x 1 L     (4)
C and L are the known values of the capacitor and the inductor respectively, and v g the grid voltage. The signals that are measured are the states x ∈ R 2 , the voltage grid v g ∈ R and the PV array current I ∈ R, which is a common practice in this type of circuits.

u is the control signal composed by the terms:

u = u n + w + δ(t, x) (5) δ(t, x
) is an uncertain term that satisfies the matching condition (i.e. it enters at the same point that the control signal u), w is the term that will be derived to compensate the uncertain term and u n is the control component that turn the nominal system globally asymptotically stable. In section IV, these components are explained in detail.

The nominal system is the system without the uncertain term δ(t, x), hence: ẋ = f (t, x) + G(t, x)u n (6) where f (t, x) and G(t, x) are defined in the same way that before.

The main tasks that the signal control u must fulfill are:

(1) To track the MPP of the PV array, despite of the changes in irradiance and temperature. (2) To deliver a current in phase with the grid voltage.

Another important issue that must be taken into account, is the total harmonic distortion (THD) injected in the grid. One way to reduce this is computing the controller references at the ending of each voltage grid cycle.

The controller needs two references to accomplish the main tasks of the system.

The first reference is the voltage DC value in the capacitor x1 * ∈ R + (7) that provides the MPP and is given by a MPPT algorithm.

The second reference is the current value in the inductor x 2 * , which is taken from Meza et al. [2012]:

x 2 * 2 v g (t) x1 * I * (x 1 * , c i ) A 2 (8)
and is the inductor current value that achieves MPPT, and that has to be in phase with the grid voltage. I * is the PV current (2) evaluated at the voltage x1 * , v g is assumed to be sinusoidal with constant amplitude A and frequency ω, v g = Asin(ωt).

However, the parameters c i in I * are unknown. Therefore, an estimator was designed and it is presented here. It allows to calculate the complete I * expression.

Hence, the second reference becomes:

x2 * 2 v g (t) x1 * Î * (x 1 * , ĉi ) A 2 (9)
where (•) indicate the estimate of (•).

3. ESTIMATOR DESIGN.

Reparametrization.

A reparametrization was made in order to be able to design the next estimator.

From ( 2) and ( 6) the derivative of the current I is:

İ = -c 2 c 3 ẋ1 e c3x1 = 1 C (ux 2 -I)(θ 1 -θ 2 I) = Φ T θ (10)
where:

θ c 1 c 3 c 3 , Φ T 1 C ux 2 -I, -ux 2 I + I 2 (11)
also consider the following estimator and observer errors respectively θ = θ -θ, Ĩ = Î -I.

Sliding mode estimator.

In [START_REF] Kenne | Nonlinear systems time-varying parameter estimation: Application to induction motors[END_REF] a new framework for nonlinear systems time-varying parameter estimation using sliding mode techniques was proposed.

The following notation is used in this subsection:

|x T | G = (|x 1 |, |x 2 |, ... |x n |)
diag(A) is the column vector whose elements are the diagonal elements of a given square diagonal matrix A.

The function sign(•) : R n → R n is defined as sign(x T ) = (sign(x 1 ), sign(x 2 ), ... sign(x n )) (12) with

x T = (x 1 , x 2 , ... x n ) and sign(

x i ) =    -1 for x i < 0 [-1, 1] for x i = 0 1 for x i > 0
Now, let us assume that Assumption 1. The derivatives of the parameters vector θ are bounded and these bounds are known: θi (t) ≤ µ θi where µ θi are known positive numbers and i = 1, 2.

In order to design the time-varying parameter estimator consider the following adaptive observer for the current I:

İ = v I + Φ T θ (13) where v I = -K I sign( Î -I) thus İ = v I + Φ T θ (14)
For K I sufficiently chosen large

K I > |Φ T θ|
and assuming that θ is bounded (the proof will be shown later) then a sliding mode regime occurs on the manifold Ĩ = 0 and 0 = v Ieq + Φ T θ then

v Ieq = -Φ T θ (15)
The following approximation is used (see the work of [START_REF] Utkin | Sliding Modes in Optimization and Control[END_REF])

v Ieq = 1 1 + τ s v I ( 16
)
where s is the Laplace operator and τ → 0 is a positive constant.

Hence

θ = -Φ(Φ T Φ) -1 v Ieq (17) and finally θ = -K θ sign( θ) (18)
Choosing the gain matrix for the estimator such that:

diag(K θ ) i > µ θi (19)
then, the estimated parameters θ will converge to their real values θ(t) in finite time.

Proof. Let us consider the following Lyapunov candidate function

V = 1 2 θT θ β( θ) (20) V = θT θ = θT (-K θ sign( θ) -θ) ≤ -| θT | G (diag(K θ ) -µ θ ) (21)
From ( 19), one can see that V is negative-definite for θ = 0. Therefore, the estimation error θ will converge to 0 in finite time. 222 4. CONTROLLER DESIGN.

In order to accomplish the control objectives, the problem has been divided in two parts:

(1) To find the control portion that renders the nominal system (6) globally asymptotically stable.

(2) To find the control portion that compensates the uncertain term δ(t, x) and gives the whole system (3) globally asymptotically stable.

Controlling the nominal system.

In this part of the analysis it is considered that the x 2 * reference is well known, i.e. the observer error is zero Ĩ(t) = 0. The term δ(t, x) contains the perturbations due to Ĩ(t) = 0 and the second portion of the control handle it.

According to the idea given in [START_REF] Sanders | Explicit model of photovoltaic panels to determine voltages and currents at the maximum power point[END_REF], the control is decomposed in two parts: one for the steadystate stage u * n and another for the dynamic stage e un .

The reference system is: ẋ * = f (t, x * ) + G(t, x * )u * n (22) Hence, the control component for the steady state is:

u * n = L ẋ2 * + v g (t) x 1 * (23)
Consider the following tracking errors and control variable respectively:

e = x -x * , e un = u n -u * n ė = ẋ-ẋ * = I(x 1 ) -I * C 0 +    - e 2 u * n C - (e 2 + x 2 * )e un C e 1 u * n L + (e 1 + x 1 * )e un L   
Let us propose the Lyapunov candidate function (25) From the fact that the first term of the right side of ( 25) is < 0 ∀ e 1 = 0, the nominal system is globally asymptotically stable (GAS). Thus, the control component e un will drive the error dynamics e asymptotically to 0.

V n = C 2 e 2 1 + L

Controlling the whole system.

The control part of the nominal system was established in ( 23) and ( 24). These equations are functions of ẋ2 * and x 2 * , and these values are not exactly known. Therefore, the control for the whole system is

u = u * n (x 2 * ) + e un (x 2 * ) + ∆u * + ∆e u = u * n (x 2 * ) + e un (x 2 * ) + δ(t, x) (26)
where

δ(t, x) ∆u * + ∆e u = L∆ ẋ2 * x 1 * + e 1 ∆x 2 * (27)
is the uncertain term, and ∆ ẋ2 * , ∆x 2 * are the perturbations due to the transients in the observer, i.e. Ĩ = 0.

In order to find the controller for the whole system the next assumption is needed.

Assumption 2. The uncertain term is bounded, i.e. ||δ(t, x)|| 2 < ξ where ξ is a unknown positive number.

Since this uncertain term satisfies the matching condition, a variation of the Lyapunov redesign technique [START_REF] Khalil | Nonlinear Systems 2nd ed[END_REF]) is used here to compensate it. Now, a term w will be added to the controller to compensate the uncertain term:

u = u n + δ(t, x) + w (28) where u n = u * n + e un The closed loop system becomes ẋ = f (t, x) + G(t, x)[u * n + e un ] + G(t, x)[w + δ(t, x)] (29) and ė = f (t, x) + G(t, x)[u * n + e un ] -ẋ * + G(t, x)[w + δ(t, x)]
Consider the following Lyapunov candidate function for the whole system:

V = C 2 e 2 1 + L 2 e 2 2 + β( θ) (30)
where β( θ) was defined in (20).

Let us omit the arguments of the functions for the sake of brevity. The time derivative is:

V = ∂V ∂e {f + G(u * n + e un ) -ẋ * } + ∂V ∂e {G(w + δ)} + β = Vn + β + ∂V ∂e Gw + ∂V ∂e Gδ (31)
To compensate the perturbation, the following term is considered:

w - ξ ∂V ∂e G || ∂V ∂e G|| 2 (32)
with the adaptive law:

ξ = -k ξ ξ for | Ĩ| ≤ α for | Ĩ| > ( 33
)
where k ξ , α > 0 and ≥ 0.

First, let us analyze the system setting = 0. Then, when Ĩ > 0, ξ(t) = αt + ξ(0), the perturbation term δ will be bounded and > 0 and the next inequality will hold in finite time:

ξ ∂V ∂e G 2 > ∂V ∂e Gδ (34)
Thus, from ( 21), ( 25), ( 31), ( 32) and ( 34), V will be negative definite and the whole system will be globally asymptotically stable (GAS) with the controller given by u n + w. Now, when Ĩ = 0 the perturbation term δ will be zero also, and the adaptive term ξ will converge to zero, "turning off" in this way, the discontinuous term w, and the whole system will be GAS.

In practice, the observer error Ĩ can't be exactly zero for a large time, even when the system variations are small. This is why a small interval needs to be setting. Then, when | Ĩ| ≤ , the adaptive term ξ will converge to zero and the whole system will continue being GAS due that the perturbation term δ will be small enough that will be dominated by the negative terms in ( 21) and ( 25). In this way, we have GAS for the whole system and chattering due to the discontinuous term w can be avoided.

SIMULATION RESULTS.

Numerical simulations were made in the Simulink/Matlab platform to verify the performance of the adaptive controller. The simulation consists in step changes in the irradiance and temperature variables. Its values were changed by 50%. This is one of the worst conditions for MPPT in PVS. The initial conditions in the plant and in the adaptive term ξ were zero and in the estimator was 1x10 -5 . The MPPT algorithm used in these simulations was the extremum seeking algorithm (ESA). The initial condition in this algorithm has to be set > 0 to avoid singularities in (23). Figures 3 and4 show the results.

The step changes in irradiance and temperature are shown in graphs a) and b) in figure 3. Graph d) in figure 4 shows the output power. This is the maximum power that the PVS can deliver with the irradiance and temperature given. It clearly follows the environmental variations. At time 4 and 7 sec. we can see small variations due to the stepping changes in the temperature.

In graphs c) and d) in figure 3, the adaptive bound ξ and the observer error Ĩ are shown. The adaptive bound starts growing when the observer error is bigger than the interval. These perturbations occur at time 5 and 11 sec. due to the irradiance variations. The perturbations due to the temperature variations are so small that they fall into the interval. Once the observer error fall into the interval, the adaptive bound ξ decrease exponentially to zero. The discontinuous controller term w is different from zero only when medium and big perturbations occur, and it will be zero again as soon as the perturbations cease. In this way, the chattering problem is avoided.

CONCLUSION.

An adaptive control scheme for maximum power point tracking of a single-phase grid-connected PVS was presented. The PVS driven by this control law is capable of deliver the maximum power available under variable environmental conditions.

An estimator was designed using sliding mode techniques. It is capable of estimate time-varying parameters that depends on irradiance and temperature.

A Lyapunov function that proves globally asymptotically stability of the system was found. The analysis includes the non-linear model of the system, the dynamics of the estimator and the uncertainty in the second reference x 2 * .

An adaptive law was proposed to substitute the perturbation bound. This adaptive law helps the system to be GAS under perturbations and it is zero when the perturbations cease. The discontinuous controller term w is driven by this adaptive law and the chattering problem is avoided in this manner.

Numerical simulations were made to verify the performance of the solution. It was showed the good performance of the system even under very demanding conditions.

Fig. 1 .

 1 Fig. 1. MPP of a typical PV array for different environmental conditions a) T = 25 o C and different irradiance values. b) G = 1000W/m 2 and different temp. values.

Fig. 2 .

 2 Fig. 2. Photovoltaic system analyzed in this paper.

  Ce 1 ė1 + Le 2 ė2 = e 1 (I -I * ) -e 1 [e 2 u * n + (e 2 + x 2 * )e un ] +e 2 [e 1 u * n + (e 1 + x 1 * )e un ] = -c 2 e 1 (e c3x1 -e c3x1 * ) -e un (e 1 x 2 * -e 2 x 1 * ) and choose the controller for the error dynamics of the nominal system: e un e 1 x 2 * -e 2 x 1 * (24) then Vn = -c 2 e 1 (e c3x1 -e c3x1 * ) -(e 1 x 2 * -e 2 x 1 * ) 2

  Fig. 3. a) Irradiance signal b) Temperature signal c) Adaptive bound ξ d) Observer error

Acknowledgements

The work of F. Jaramillo has been supported by CONA-CyT, Mexico and the French Embassy in Mexico.


Micha L Pióro
email: michal.pioro@tele.pw.edu.pl

Ilya Kalesnikau

Michael Poss
email: michael.poss@lirmm.fr

Path Generation for Affine Flow Thinning

Keywords: Survivable networks, variable link capacity, linear and mixed-integer programming, multicommodity flows, affine routing, path generation

Flow thinning (FT) is a traffic routing and protection strategy for communication networks whose links experience fluctuations in available capacity (as, e.g., in wireless networks). To cope with this phenomenon, end-to-end traffic demands are assigned dedicated tunnels (e.g., MPLS tunnels) whose nominal capacity is subject to thinning in order to account for variable capacity of the links, fluctuating below their nominal values. Consequently, the instantaneous traffic sent between the demand's end nodes must accommodate to the current total capacity available on its dedicated tunnels. In the paper we present a path generation-based approach to network dimensioning for a practical version of FT applying affine flow thinninga fairly complicated issue not yet considered. We derive a relevant pricing problem and present a numerical study illustrating efficiency of the optimization algorithm.

Introduction

Flow thinning (FT), a concept introduced in [START_REF] Pióro | Optimizing flow thinning protection in multicommodity networks with variable links capacity[END_REF], is a traffic routing and protection strategy applicable to communication networks with variable capacity of links. In such networks the links do not attain their nominal (maximum) capacity simultaneously, so that in a typical network state only some links are fully available while on each of the remaining links only a fraction of its maximum capacity is usable. Each end-to-end traffic demand is assigned a set of logical tunnels whose total capacity is dedicated to carry the demand's traffic. The nominal (i.e., maximum) capacity of the tunnels, supported by the nominal (maximum) link capacity, is subject to state-dependent thinning to account for variable capacity of the links, fluctuating below the maximum. Accordingly, the capacity available on the tunnels is also fluctuating below their nominal levels and hence the instantaneous traffic sent between the demand's end nodes must accommodate to the current total capacity available on its dedicated tunnels. FT (a generalization of the demand-wise shared protection strategy [START_REF] Koster | Demand-wise shared protection for meshed optical networks[END_REF][START_REF] Wessäly | Demand-wise shared protection revisited: A new model for survivable network design[END_REF] devised for total link failures) is relevant, among other applications, for wireless mesh networks utilizing MPLS tunnels.

The multicommodity flow optimization problem for FT is N P-hard and its non-compact linear programming formulation requires path generation (PG) (see [START_REF] Ahuja | Network Flows: Theory, Algorithms, and Applications[END_REF][START_REF] Orlowski | Complexity of column generation in network design with path-based survivability mechanisms[END_REF]) using a pricing problem to find the paths to be added to the problem in the consecutive steps of the PG algorithm. In this paper we consider an affine version of FT, referred to as AFT, where the capacity of each tunnel is thinned according to a tunnel-dependent linear function whose arguments are link failure coefficients. In particular, we study the pricing problem (PP) for AFT -the issue not investigated in [START_REF] Pióro | Optimizing flow thinning protection in multicommodity networks with variable links capacity[END_REF]. Our motivation is that PP for the affine case turns out to be more complicated than its counterpart presented in [START_REF] Pióro | Optimizing flow thinning protection in multicommodity networks with variable links capacity[END_REF] for the general form of FT where thinning is not restricted. In the next sections we formulate the affine version of the FT optimization problem, derive the relevant PP, and present a numerical study illustrating efficiency of the PG algorithm, and effectiveness of the considered version of affine FT as compared to the effectiveness of FT.

Notation

The (undirected) network graph G = (V, E) is composed of the set of nodes V and the set of links E, where each link represents an undirected pair {v, w} of nodes for some v, w ∈ V, v = w. We mention here that considerations of the paper are applicable to directed networks as well.

The nominal capacity of link e ∈ E is denoted by y 0 e and its unit capacity cost by ξ(e) (cost of the network is given by C = e∈E ξ(e)y 0 e ). Yet, the capacity of links is variable, and typically not all the links reach their maximum (nominal) capacity simultaneously. Thus, at any time instant, nominal link capacities y 0 e , e ∈ E, are typically available only for a subset of links, and the remaining links have the available capacity reduced. To consider this, we specify a set of link availability states S. Each state s ∈ S is characterized by the link availability coefficients α(e, s) (0 ≤ α(e, s) ≤ 1), e ∈ E. By definition, the capacity of link e available in state s is equal to y s e := α(e, s)y 0 e . In the sequel we use link failure coefficients β(e, s) := 1 -α(e, s) instead of alphas. Traffic demands are represented by the set D. Each demand d ∈ D is characterized by an undirected pair {o(d), t(d)} composed of the demand's end nodes (origin and termination), and given traffic volumes h(d, s) to be realized in the considered link availability states s ∈ S. Traffic volumes and link capacities are expressed in the same units. Each demand d is assigned a set of admissible paths P(d) (called the path-list) composed of selected elementary paths between o(d) and t(d). Paths in P(d), used to realize the demand (traffic) volumes, are assigned nominal flows x 0 dp , p ∈ P(d), which are optimization variables. Each x 0 dp specifies the nominal capacity (expressed in the same units as link capacity and demand volume) reserved on the tunnel realized by path p ∈ P(d). The set of all admissible paths is denoted by P := d∈D P(d). The maximum path-lists, i.e., path-lists P(d) containing all the elementary paths between o(d) and t(d), will be denoted by P(d), d ∈ D, with P := d∈D P(d).

Since we assume elementary paths, the paths can be identified with the sets of links they traverse: p ⊆ E, p ∈ P. The given sets of admissible paths define the link-path incidence coefficients δ(e, d, p), e ∈ E, d ∈ D, p ∈ P d , where δ(e, d, p) = 1 if e ∈ p, and δ(e, d, p) = 0 otherwise. The sets of admissible paths P(d) are parameters in the problem formulations considered in the sequel, although in general we assume that all possible elementary paths can potentially be used if this is required to achieve the optimum.

In general, not all nominal path-flows x 0 dp , d ∈ D, p ∈ P(d), can be realized when the network is in state s ∈ S, as the available link capacity y s := (y s e , e ∈ E) is decreased with respect to the nominal link capacity y 0 . To account for that, the nominal flows must be thinned (decreased) to fit the reduced link capacity. Still, the thinned flows must be sufficient to carry the (reduced) demands h(d, s), d ∈ D. The thinned nominal path-flows for state s ∈ S are denoted by x s dp , d ∈ D, p ∈ P(d). These flows are reserved on the admissible paths for the duration of the state. As only thinning of nominal flows is allowed, the following inequalities must hold: e (e ∈ E). In the formulation, E(d, p) ⊆ E is a (predefined) subset of links whose failure coefficients influence the linear (affine to be precise) thinning formula for path p ∈ P(d). Dual variables indicated in the square brackets will be used in [START_REF] Orlowski | Complexity of column generation in network design with path-based survivability mechanisms[END_REF].

x s dp ≤ x 0 dp , d ∈ D, p ∈ P(d), s ∈ S.

ATOP -Affine Thinning Optimization Problem

Problem P(P):

C(P) = min e∈E ξ(e)y 0 e (1a) [π 0 e ≥ 0] d∈D p∈P(d) δ(e, d, p)x 0 dp ≤ y 0 e , e ∈ E (1b) [λ s d ≥ 0] p∈P(d) x s dp ≥ h(d, s), d ∈ D, s ∈ S (1c) [π s e ≥ 0] d∈D p∈P(d) δ(e, d, p)x s dp ≤ α(e, s)y 0 e , e ∈ E, s ∈ S(e) (1d) [ϕ s dp ] x s dp = x 0 dp -e∈E(d,p) β(e, s)z e dp , d ∈ D, p ∈ P(d), s ∈ S ( 1e
)
y 0 ≥ 0, x 0 ≥ 0, x s ≥ 0, z ≥ 0 and continuous. ( 1f
)
Objective (1a) minimizes the cost of links. The first capacity constraint (1b) makes sure that the nominal link loads do not exceed the nominal link capacities. Next, the demand constraint (1c) assures that in each state s ∈ S, the thinned nominal flows are sufficient to realize the volume of each demand d ∈ D assumed for this state. Then, the second capacity constraint (1d) ensures that the capacity of each link e ∈ E available in each state s ∈ S(e) is not exceeded (where S(e) denotes the set of states in which α(e, s) < 1). Finally, equalities (1e) ensure thinning since, by (1f), all z e dp are nonnegative. Note that constraint (1d) is, for each link e ∈ E, written down only for s ∈ S(e), and not for all s ∈ S. This is because for any state s ∈ S \ S(e) (in which the capacity of the considered link e is fully available, α(e, s) = 1) constraint (1d) is implied by (1b) and (1e). Note also that in general feasibility of ATOP is assured only if α(e, s) > 0 for all e ∈ E and s ∈ S.

The thinning formula (1e) can be generalized in a natural way to:

x s dp = z 0 dp -e∈E(d,p) β(e, s)z e dp , d ∈ D, p ∈ P(d), s ∈ S, (2)
where all z e dp are continuous and unconstrained in sign (note that z 0 dp provides an additional "degree of freedom" in the thinning formula), and x s dp ≤ x 0 dp , d ∈ D, p ∈ P d , s ∈ S. As discussed in [START_REF] Pióro | Optimizing flow thinning protection in multicommodity networks with variable links capacity[END_REF], typical examples of the sets E(d, p) are E(d, p) = p, E(d, p) = v∈V(p) δ(v) (where V(p) is the set of nodes traversed by path p, and δ(v) is the set of links incident with v), and E(d, p) = E.

Finally, the most general thinning is achieved when only the restriction

x s dp ≤ x 0 dp , d ∈ D, p ∈ P(d), s ∈ S (3)
is assumed and no affine formula as ( 2) is applied. The so obtained generalization of ATOP was studied in [START_REF] Pióro | Optimizing flow thinning protection in multicommodity networks with variable links capacity[END_REF] under the name FTOP.

We end this section with giving the formulation dual to ATOP (1):

Problem D(P):

W (P) = max d∈D s∈S h(d, s)λ s d (4a) π 0 e + s∈S(

Path Generation and the Pricing Problem

The link-path LP formulation (1) of ATOP is noncompact because of exponentially many routing paths in the maximum path-list P. Although the complexity of ATOP is open as yet, we may suspect that the problem is difficult, taking into account that we are not aware of its compact LP formulations, and that FTOP -the generalization of ATOP assuming (3) instead of (1e)is N P-hard [START_REF] Pióro | Optimizing flow thinning protection in multicommodity networks with variable links capacity[END_REF]. Since in practice formulation (1) cannot be solved directly (as we are not able to include all elementary paths in the path-lists), to consider all paths in P we need to apply path generation (PG, see [START_REF] Ahuja | Network Flows: Theory, Algorithms, and Applications[END_REF][START_REF] Orlowski | Complexity of column generation in network design with path-based survivability mechanisms[END_REF][START_REF] Pióro | Routing, Flow, and Capacity Design in Communication and Computer Networks[END_REF]) -a classical technique in multicommodity flow networks related to column generation in linear programming (cf. [START_REF] Minoux | Mathematical Programming: Theory and Algorithms[END_REF]). With PG, starting from some initial path-lists P(d), d ∈ D, we iteratively generate new paths, one per demand, and add those paths that may improve the solution to the path-lists -in effect, we are adding the corresponding variables and constraints to the problem formulation (i.e., we generate new columns and rows in ( 1)). This is done by solving an appropriate pricing problem (PP) using, as parameters, optimal dual variables, i.e., an optimal solution (λ * , π * ) of the current dual D(P). PP for ATOP turns out to be more complicated than its FTOP counterpart derived in Section 4 of [START_REF] Pióro | Optimizing flow thinning protection in multicommodity networks with variable links capacity[END_REF]. PP for ATOP is derived below.

Consider a given demand d ∈ D and a given path q ∈ P(d) \ P(d). Introducing a new path q ∈ P(d) \ P(d) to the dual (4) introduces both new (we call them local in the sequel) dual variables ϕ s , s ∈ S, and new constraints (4c)-(4e). It can happen that some of the new constraints (4c)-(4e) will be violated by the optimal dual solution (λ * , π * ) of D(P) for any values of the local variables. In such a case some of the new constraints will separate (λ * , π * ) from the dual polyhedron of D(P ∪ {q}). The sum of such violations is expressed through the quantity P (q) defined by the following LP.

P (q) = min Z + s∈S Y s + e∈E(d,q) X e (5a) [g ≥ 0] Z ≥ s∈S ϕ s -|q| 0 (5b) [a s ≥ 0] Y s ≥ λ s d * -|q| s -ϕ s , s ∈ S (5c) [b e ≥ 0] X e ≥ -s∈S β(e, s)ϕ s , e ∈ E(d, q) (5d) Z ∈ R + ; Y s ∈ R + , s ∈ S; X e ∈ R + , e ∈ E(d, q); ϕ s ∈ R, s ∈ S, (5e)
where |q| 0 := e∈q π 0 e (nominal dual length of q), and |q| s := e∈q∩E(s) π s e , s ∈ S (state-dependent dual length of q).

Since we are interested in maximizing the value of P (q), we cannot use a minimization objective as in (5a). Therefore we consider the problem dual to (5), formulated using the dual variables specified in the square brackets in constraints (5b)-(5d) (in ( 6) dual variables a s , s ∈ S, are eliminated):

P (q) = max -g|q| 0 + s∈S (g -e∈E(d,q) β(e, s)b e )(λ s d * -|q| s )) (6a)
g ≤ 1; b e ≤ 1, e ∈ E(d, q); e∈E(d,q) β(e, s)b e ≤ g, s ∈ S (6b)

g ∈ R + ; b e ∈ R + , e ∈ E. (6c)
PP consists in finding, for each d ∈ D, a path q that maximizes (6a). PP is derived from (6) by introducing binary variables u e , e ∈ E, and constraints (7b) (so that the set {e ∈ E : u e = 1} specifies a path), and by eliminating bi-linearities b e u e , gu e , b e u e u e that will appear. For undirected networks and E(d, p) = p, the resulting PP is as follows (below the superscript * in π * , λ * is skipped, and δ(v) denotes the set of links incident to v ∈ V):

max -e∈E π 0 e G e +
e∈δ(v) u e = 1, v ∈ {o(d), t(d)}; e∈δ(v) u e = 2z v , v ∈ V \ {o(d), t(d)} (7b) g ≤ 1; b e ≤ 1, e ∈ E;
u e ∈ {0, 1}, e ∈ E; z v ∈ {0, 1}, v ∈ V \ {o(d), t(d)}. ( 7h
)
An optimal solution u * of ( 7) defines the path q := {e ∈ E : u * e = 1} that should be added to ATOP (provided P (q) > 0).

Note that any feasible set of links U := {e ∈ E : u e = 1} will necessarily contain an elementary path between o(d) and t(d) and, possibly, a set of disjoint isolated loops. The path is elementary because variables z v , v ∈ V, are binary (note that assuming these variables to be integer-valued would allow paths to contain loops as well). Such loops may appear in optimal PP solutions for several reasons. One of them is that there is a tradeoff related to the number of elements in the set U: the more elements the more flexible formulae (1e), and the less elements the lower link capacities.

To eliminate the (isolated) loops we introduce the bi-directed version G = (V, A) of the original undirected network graph G = (V, E), where the set of (directed) arcs is defined as A := {(v, w), (w, v) : {v, w} ∈ E}, i.e., each undirected link e ∈ E is substituted by two oppositely directed arcs a (e) and a (e). Below, δ -(v) denotes the set of all arcs incoming to node v, and δ + (v) -the set of all arcs outgoing from node v (v ∈ V). Adding aggregated arc-flow variables f a ≥ 0, a ∈ A, and the following constraints to [START_REF] Pióro | Routing, Flow, and Capacity Design in Communication and Computer Networks[END_REF] will eliminate the loops in question.

a∈δ -(t(d)) f a = v∈V\{o(d),t(d)} z v (8a) a∈δ + (o(d)) f a = a∈δ -(o(d)) f a (8b) a∈δ + (v) f a = a∈δ -(v) f a + z v , v ∈ V \ {o(d), t(d)} (8c) f a (e) + f a (e) ≤ (|V| -2)u e , e ∈ E. (8d)
In the next section we will also consider the case E(d, p) = E, d ∈ D, p ∈ P(d). Then formulation (7) becomes a bit simpler since variables B e become equal to b e , and variables T e e express the product u e •b e rather than u e •u e •b e .

Numerical Study

Below we present numerical results for a network with 12 nodes, 18 (undirected) links and 66 (undirected) demands. The network is fully described in SNDlib (sndlib.zib.de, see [START_REF] Orlowski | SNDlib 1.0 -survivable network design library[END_REF]) under the name polska. We use the data from the first instance (i.e., polska-D-B-M-N-C-A-N-N) and assume the link unit capacity cost ξ e equal to its <module cost> divided by its <module capacity> for the first pair of these values. Thus, for the first link (between Gdansk and Warsaw) ξ 1 = 156/155 = 1.0645, and the entire vector ξ of the unit capacity costs is as follows: 1.0645, 1.7548, 1.0064, 1.2000, 1.7548, 1.5290, 1.3419, 1.1677, 1.3419, 1.6129, 2.0903, 2.0903, 1.6129, 1.0645, 1.9677, 0.9161, 1.2581, 1.8968 (the order of links is the same as in the SNDlib file).

In Table 1 we present the results illustrating time efficiency of the PG algorithm for ATOP and FTOP, and cost efficiency of ATOP as compared with FTOP (referred to as F-T). Two cases of ATOP are considered: A-S with E(d, p) = p, and A-F with E(d, p) = E. All the reported calculations were executed on a PC-class computer (Windows 10 64-bit, 8 GB RAM, Processor Intel Core i5-3210M, 4 logical processors, 2.5GHz) using CPLEX 12.4.0.0.

Table 1 shows the results for three failure scenarios (in all of them, the failing links loose half of their capacity, i.e, if link e fails in state s, its link availability coefficient is α(e, s) = 0.5):

• Single link (SL) failure scenario S1 contains the state s(0) (all links fully available), and all the states with exactly one link failing (note that |S1| = |E| + 1). 100% traffic protection is assumed so that h(d, s) := h(d), d ∈ D, s ∈ S1, and h(d) are the demand volumes described in SNDlib.

• Double link (DL) failure scenario S2 contains S1 (with no reduction of demand volumes) and all the states with exactly two links failing (|S2| = |S1| + |E|(|E| -1)/2). For the states in S2 \ S1 the demand volumes are reduced to 95% of the volumes assumed for the states in S1.

• Triple link (TL) failure scenario S3 contains S2, and all the states with exactly three links failing (|S3| = |S2| + |E|(|E| -1)(|E| -2)/6). For S2 the demand reductions are as in DL, and for S3 \ S2 the demand volumes are reduced to 90% of the S1 volumes.

For PG, the initial path list contains one path for each demand -the shortest path with respect to the link unit costs (66 paths in total). Note that assigning the entire demand volume h(d, s(0)) for each demand d ∈ D to such a shortest path will result in the optimally dimensioned network when protection is not considered. For the examined network example this cost is equal to 30275.5.

In Table 1 the columns presenting the results describe: initial cost (C 0 ), final cost (C * ), cost increase of ATOP with respect to FTOP (∆C), final number of paths |P * |, number of iterations (iter), total computation time (total), time per iteration (t/iter), percentage of the total computation time spent in the master problem (t/MP), and percentage of the total computation time spent in the pricing problem (t/PP). It should be noted that the computational time for A-S is usually greater than for A-F due to the fact that the PP corresponding to A-S contains additional binary variables with respect to the PP corresponding to A-F because E(d, p) = p for A-S and E(d, p) = E for A-F. For the same reason, the network cost achieved for A-F is always less than or equal to that for A-S. (In the reported cases these costs are equal but this is not the case in general). Finally, let us note that, again for the same reason, A-S needs to be considered as its potential network implementation is much simpler than for A-F. A general conclusion from the presented study is that the PG process is pretty effective, which is especially important with respect to PP -a MIP formulation of a potentially hard optimization problem.

Comments and Concluding Remarks

Although affine flow thinning (AFT) requires more link capacity than the original version (FT) with no constraints on thinning (except (3)), the affine versions are important for three reasons. First, as FT, AFT must be aware of the current link states α(•, s); yet it will use the (optimized) affine functions for thinning the nominal flows. Such thinning, contrary to FT, is applicable to all possible states, also to the states not in the set S assumed for optimization (traffic efficiency of such extended thinning is addressed in Section 8.2 of [START_REF] Pióro | Optimizing flow thinning protection in multicommodity networks with variable links capacity[END_REF]).

Second, AFT in its version with the thinning dependent only on the failure coefficients β along the path realizing the given flow (and for that matter, in the version taking into account the failure coefficients of all the links e belonging to the sets δ(v) for all the nodes v traversed by the considered path) is much more practical than FT in the network implementation aspect. In such a case disseminating link capacity information is, as discussed in Section 9.2 of [START_REF] Pióro | Optimizing flow thinning protection in multicommodity networks with variable links capacity[END_REF], not a major issue, as the path-originating node thins the corresponding path-flow only on the basis of the current links' capacities along the path.

Third, from the computational efficiency viewpoint, the main interest in AFT stems from its capability of handling large sets S of link availability states, provided S is given as the set of vertices (extreme points) of a polytope (called uncertainty polytope) described by a polynomial number of variables and constraints. This is explained in Section 7.3 of [START_REF] Pióro | Optimizing flow thinning protection in multicommodity networks with variable links capacity[END_REF].

In the paper we have discussed an important issue in optimization of affine thinning -path generation. We have derived an appropriate formulation of the pricing problem and shown by numerical examples that the path-generation resolution process of the basic optimization problem is quite effective. An issue left for future studies is path generation for the general affine thinning case defined by (2) (which requires less link capacity than the two versions considered in Table 1), and also for the case when the link availability states are specified by means of the uncertainty polytope.

For

  given admissible path-lists P(d), d ∈ D, the considered problem (ATOP) is denoted by P(P) (where P = d∈D P(d)). ATOP can be formulated as the following link-path linear programme (LP) involving variables x 0 dp (d ∈ D, p ∈ P(d)), x s dp (d ∈ D, p ∈ P(d), s ∈ S), z e dp (d ∈ D, p ∈ P(d), e ∈ E) and y 0

  e∈E β(e, s)B e ≤ g, s ∈ S (7c) B e ≤ u e , B e ≤ b e , B e ≥ b e + u e -1, e ∈ E (7d) G e ≤ g, G e ≤ u e , G e ≥ g + u e -1, e ∈ E (7e) T e e ≤ u e , T e e ≤ u e , T e e ≤ b e , T e e ≥ b e + u e + u e -2, e , e ∈ E (7f) g ∈ R + ; b e , B e , G e ∈ R + , e ∈ E; T ee ∈ R + , e, e ∈ E (7g)

1 This work was supported by NCN (Poland) [grant 2015/17/B/ST7/03910].