Hassan El Fadil

Fouad Giri
email: fouad.giri@unicaen.fr

Fatima-Zahra Z Chaoui

Ouadia El Magueri

O El Magueri

Accounting for Input Limitation in the Control of Buck Power Converters

Keywords: Absolute stability, DC-DC power converters, input limitation, voltage regulation

We are considering the problem of controlling a DC-DC switched power converter of the Buck type. The converter involves an inherent control limitation; accordingly the control signal (duty ratio) can only take values in the interval (0 1). In the relevant literature, such a physical control limitation is generally not taken into account when designing the converter regulators. This is only dealt with in the control implementation stage, placing an isolated limiter between the (linear) controller and the controlled system. Furthermore, the presence of such a limiter is generally ignored when analyzing the closed-loop control system. In the present paper, the control signal limitation is dealt with using a nonlinear regulator that involves an internal limiter. The resulting closed-loop control system is shown to be equivalent to a nonlinear feedback involving a linear dynamics block in closed-loop with a nonlinear static element. Using absolute stability tools, sufficient conditions are established for the involved feedback to be 2 -stable. If these conditions are respected when choosing the control design parameters then the regulator meets its objectives (closed-loop stability and output reference tracking). It is worth noting that, though the focus is made on a specific power converter, the paper includes an important theoretical dimension that may be of general interest.

I. INTRODUCTION

T HERE are three main types of switched power con- verters, namely boost, buck, and buck-boost. These have recently received an increasing deal of interest both in power electronics works and in automatic control applications. This is due to their wide applicability domain, e.g., domestic equipments, communication systems, computers, industrial electronics, battery-operating, embedded equipments, uninterruptible power sources. From an automatic control viewpoint, a switched power converter constitutes a challenging case study as it is variable-structure and nonlinear. Its rapid structure variation is generally coped with using averaged models [START_REF] Sira-Ramirez | Control Design Techniques in Power Electronics Devices[END_REF], [START_REF] Sira-Ramirez | Passivity-based controllers for the stabilisation of DC-to-DC power converters[END_REF], [START_REF] Rashid | Power Electronics Handbook[END_REF]. Based on these average models, different nonlinear controllers have been developed using passivity techniques [START_REF] Sira-Ramirez | Passivity-based controllers for the stabilisation of DC-to-DC power converters[END_REF], feedback linearization, flatness methods [START_REF] Gensior | On differential flatness, trajectory planning, observers, and stabilization for DC-DC converters[END_REF], sliding mode control ( [START_REF] Tan | A pulse-widthmodulation based sliding mode controller for buck converters[END_REF], [START_REF] Nguyen | Indirect implementation of sliding mode control law in buck-type converters[END_REF], [START_REF] Gupta | Frequency-domain characterization of sliding mode control of an inverter used in DSTATCOM application[END_REF]) and backstepping control technique [START_REF] Fadil | Nonlinear and adaptive control of buck power converters[END_REF]- [START_REF] Alvarez-Ramirez | Current-mode control of DC-DC power converters: A backstepping approach[END_REF]. In all works the proposed controllers are designed to achieve closed-loop global stability and voltage reference tracking. However, in these studies it was never accounted for the limitation of the control signal (namely, the duty ratio ). As a matter of fact, the control is not allowed to be outside the interval (0, 1), due to the technological nature of the controlled circuits. Therefore, the usual practice consists, when it comes to implementing the controllers, in simply placing an isolated limiter between the designed controller and the controlled system (see Fig. 1). Unfortunately, the nonlinear effect of such a limiter is never taken into account when analyzing the closed-loop system ([2]- [START_REF] Alvarez-Ramirez | Current-mode control of DC-DC power converters: A backstepping approach[END_REF]). As a consequence, the aforementioned stability results may loose their global nature, the controller transient performances may deteriorate and the output-reference tracking objective may not be achieved.

On the other hand, the problem of controlling linear systems subject to input saturation constraint has received a great deal of interest over the last two decades [START_REF] Glattfelder | Control Systems With Input and Output Constraints[END_REF]. Most theoretical results have focused on (local/global) stabilization. This was carried out following two main approaches. The first one consists in enforcing the control signal to stay all time within the allowed limits so that the closed-loop system stays in a region of linear behavior, see e.g., [START_REF] Benzaouia | The regulator problem for linear systems with saturations on the control and its increments or rate: An LMI approach[END_REF], [START_REF] Liu | Contrained control of positive discrete-time systems with delays[END_REF] and reference list therein. It is proved that linear regulators can be designed so that to ensure the positive invariance of this linear region (or a subset of it including all admissible initial conditions). Such a positive invariant set is then a local region of stability. The second research direction to handle input constraints is one where the control signal is allowed to saturate. During the time intervals of control saturation, the closed-loop system is no longer linear and the controlled system is steered in open loop. Then, global asymptotic stabilization is only possible for stable systems [START_REF] Glattfelder | Control Systems With Input and Output Constraints[END_REF]. The problem of output-reference (in presence of control input limitation) has not been so deeply investigated. When the controlled system is of type-1 and the reference signal is step-like, then the tracking problem can be transformed into a disturbance-free regulation problem and existing solutions can be applied, see e.g., [START_REF] Chaoui | Adaptive control of input-constrained type-1 plants: Stabilization and tracking[END_REF]. For non-type-1 systems, an integrator should be incorporated in the regulator to make the tracking objective achievable. But, the presence of an integrator generally results in large control actions which, due to the control limitation, may lead to undesirable oscillatory behavior. To avoid such behavior, many authors have proposed linear regulators together with ad hoc "anti-windup" devices (examples of such devices can be found in [START_REF] Glattfelder | Control Systems With Input and Output Constraints[END_REF]). Unfortunately, there is generally no formal proof that the proposed anti-windup regulators actually ensure the tracking objective for a specific class of reference signals.

The present paper is precisely focusing on the tracking issue in presence of control limitation. The focus is made on dc-to-dc Buck power converters operating according to the PWM principle. Our approach consists in first designing a linear control law that achieves the tracking objective in the absence of control limitation. Then, an adequate anti-windup device is incorporated, in the above control law, leading to a nonlinear regulator. It is shown that the resulting closed-loop control system is equivalent to a nonlinear feedback consisting of a linear dynamics block in closed-loop with a nonlinear static element. Sufficient conditions for -stability of this feedback are then established using tools from the absolute stability theory (circle criterion, Barbalat's lemma ). The main condition is the realpositivity (RP) of a specific transfer function that involves, on one hand, the poles of the controlled system and, on the other hand, those of the underlying linear closed-loop system (i.e., the linear closed-loop system obtained when the controller never saturates). In fact, the RP condition defines a neighborhood of the controlled system poles in which should be assigned those of the closed loop. Such a requirement can always be satisfied through an adequate choice of the controller parameters. Finally, it is shown that if the reference signal is slowly varying (in a well defined sense) then the proposed nonlinear controller stops saturating after finite transient-periods (following the reference signal changes) and the system output tracks asymptotically its reference. Additional simulation results show that the developed regulator has quite interesting robustness capability e.g., it preserves a reasonable level of performances when facing uncertain load changes and transient discontinuous conduction mode operation.

The paper is organized as follows: in Section II, the Buck converter is described and the control objective is formulated; Sections III and IV are devoted to the regulator design and theoretical analysis; the closed-loop performances are illustrated in Section V. A conclusion and a reference list end the paper.

II. CONVERTER MODELING AND CONTROL OBJECTIVE

A. Converter Modeling

A Buck converter is constituted of power electronic components connected together, as shown in Fig. 2. It operates according to the so-called Pulse Width Modulation (PWM) principle and is described by the following instantaneous model:

(1a)
where

. Resorting to averaging techniques (see, e.g., [START_REF] Sira-Ramirez | Control Design Techniques in Power Electronics Devices[END_REF], [START_REF] Sira-Ramirez | Passivity-based controllers for the stabilisation of DC-to-DC power converters[END_REF]), it is shown that the circuit, supposed to be operating As averaging is performed over cutting periods, the state variables and denote the average input current and average output capacitor voltage , respectively. The control signal (said duty ratio) is in turn the average value of the actual binary-type control . As a consequence of the binary nature of , the duty ratio turns out to be a continuous function of time that only takes values in the real interval i.e., , all the time. By definition, the system output is:

(2b)
In the sequel, the symbol ' ' refers generally to the time-derivation operator, i.e.,

. Then, for a given signal , one has:

(3)

Using this notation, the state-space model (2a) may be given the following input-output representation1

where

(5a) (5b)
The parameters , , are supposed to be constant and known during all stages of the regulator design (Section III) and closed-loop theoretical analysis (Section IV). Nevertheless, the effect of uncertain load changes and transient discontinuous conduction mode operation will be investigated when evaluating the regulator performances in Section V.

B. Control Objective

The control problem at hand is to design, for the considered circuit, a controller that meets the following.

1) The control action, i.e., the generated duty ratio function , should satisfy the following inequalities [START_REF] Fadil | Nonlinear and adaptive control of buck power converters[END_REF] where are any real constants such that .

2) All the closed-loop signals should remain bounded.

3) The output tracking error should converge to zero, where is a given bounded reference signal. This represents the desired trajectory of the output . But, for this tracking objective to be achievable, the signal has to be compatible with the control constraint (6), i.e.,: [START_REF] Sira-Ramirez | Dynamical adaptive pulse-width modulation control of DC to DC power converters: A backstepping approach[END_REF] Note that and this quantity is nothing other than the static gain of the system (4). The reference signal is further supposed to be slowly varying in the sense that:

(8)
Remark 2.1: 1) Condition ( 7)-( 8) will prove to be crucial in achieving the tracking objective in presence of the control limitation [START_REF] Fadil | Nonlinear and adaptive control of buck power converters[END_REF].

In other words, a reference signal that does not satisfy the double condition ( 7)-( 8) is not admissible in presence of the constraint (6) (though it is admissible in the constraint-free case).

2) Theoretically, one may let and in ( 6). However, it is not recommended in practice to impose to the duty ratio the extreme values 0 and 1. Therefore and should be given values close to (but different from) 0 and 1, respectively.

III. REGULATOR DESIGN

The first step in the regulator design consists in transforming the initial tracking problem into a regulation problem. To this end, introduce the following tracking error: [START_REF] Khalil | Nonlinear Systems Analysis[END_REF] Operating on and using (4) gives the following fictive system: [START_REF] Anderson | Optimal Control: Linear Quadratic Methods[END_REF] The latter is viewed as a system of input and output . The quantity , which belongs to due to (8), stands as a disturbance in the new system [START_REF] Anderson | Optimal Control: Linear Quadratic Methods[END_REF]. Enforcing the output of the initial system (4) to track a reference , amounts to enforcing the output of the new system (10) to vanish asymptotically. If the control input were not constrained, it would have been possible to achieve the new regulation objective using a linear regulator of the form . The differential operators and can be determined by linear design methods (e.g., pole-placement, linear-quadratic, see e.g., [START_REF] Anderson | Optimal Control: Linear Quadratic Methods[END_REF]). Since the control signal is constrained to stay in the interval , we consider the following nonlinear controller:

(11a) (11b)
The real coefficients in (12a)-(12b) will be defined later. The closed-loop system thus obtained is represented by Fig. 3. Note the presence of the nonlinear element in the controller block.

Remarks 3.1: a) Suppose that the computed control action satisfy the constraint all the time. Then, it follows from (13) that: and, consequently, the control law (11a) reduces to a standard linear regulator of the form [START_REF] Gupta | Frequency-domain characterization of sliding mode control of an inverter used in DSTATCOM application[END_REF] The closed-loop control system then turns out to be a standard linear feedback as shown by (Fig. 4). Notice that ( 14) contains an integral action . Moreover, the above control law is nothing other than a PID controller. Indeed, using (12a)-(12b), it can be easily checked by direct substitution that the transfer function may be given the following form with b) Contrary to ( 14), the nonlinear controller (11a)-(11b) involves two control variables, namely and . The Fig. 5. Closed-loop system obtained with the standard controller ( 14) and an isolated control limiter.

first represents the control action that is suggested by the underlying linear controller (11a). The second is the control action that is actually applied to the controlled system. Equation (11a) shows that the suggested control action accounts for the actual control action . Consequently, the nonlinear controller involves a feedback around the saturation element (see Fig. 3). This contrasts with the usual practice (Fig. 1) that consists in just placing the saturation element outside the loop, as this is illustrated by Fig. 5. In the later, there is no compensation of the disturbing effect due to the control limitation. It is clear that if the controller (11a)-(11b) never saturates (i.e., and ) it will coincide with the linear controller [START_REF] Gupta | Frequency-domain characterization of sliding mode control of an inverter used in DSTATCOM application[END_REF]. When this happens, the closed loop systems (Figs. 3 and4) will be identical. c) As this was mentioned in Section I, controllers like (11a)-(11b) that attempt to compensate the effect of control limitation are usually referred to anti-windup controllers (see e.g., [START_REF] Glattfelder | Control Systems With Input and Output Constraints[END_REF]). In this respect, the controller of Fig. 5 involves no anti-windup device Many design techniques may be used to obtain the operators and in (11a). Specifically, all linear design methods can be used e.g., pole-placement and linear-quadratic [START_REF] Anderson | Optimal Control: Linear Quadratic Methods[END_REF]. For the sake of simplicity, let us consider the pole placement technique. Accordingly, the above operators are uniquely obtained solving the following Bezout equation: [START_REF] Chaoui | Adaptive control of input-constrained type-1 plants: Stabilization and tracking[END_REF] where (16a) (16b) are any Hurwitz polynomials whose choice will be subject to conditions précised later 2 . For now, let us just note that for these polynomials to be Hurwitz, their coefficients must be positive (because of the Routh criterion). That is, one has: (16c)

IV. ANALYSIS OF THE CLOSED-LOOP CONTROL SYSTEM

The closed-loop control system of Fig. 3, consisting of the controlled system (4) and the regulator (11a)-(11b), will now be analyzed. The outcome of such an analysis will help making 2 Recall that the operators C (s) ; 3(s) are said Hurwitz if their zeroes have negative real parts. Then, the transfer functions 1=C (s) ; 1=3(s) are asymptotically stable, (see e.g., [START_REF] Glattfelder | Control Systems With Input and Output Constraints[END_REF]). suitable choices of the design parameters, namely the operators and , so that the control objective of Section II can be achieved. To this end, the closed-loop system will first be given a feedback representation that allows application of absolute stability theorems [START_REF] Khalil | Nonlinear Systems Analysis[END_REF].

A. A Suitable Feedback Representation of the Closed-Loop System

Combining ( 10) with (11a) so that to eliminate yields:

(17) Equation ( 17), together with the fact that , can be given the feedback representation of Fig. 6, with

(18) (19)
and is a (nonlinear) operator that maps into .

B. Conicity of the Nonlinear Element

The following analysis aims at establishing a sufficient condition (on the polynomial ) that ensures the -stablity of the above feedback. The first step is to show that belongs to the sector . Proposition 1: The nonlinear element , in the feedback of Fig. 6, belongs to the sector in the sense that: , where and Proof: It follows from ( 13) that, for all :

(20)

On the other hand, one has from (11b): C.

-Stability of the Feedback of Fig. 6 Proposition 1 will now be used to establish a sufficient condition, on the operator , that ensures the -stability of the feedback system of Fig. 6.

Proposition 2: Let the polynomial , in ( 15) and (16a), be chosen so that:

(24)
which, in view of (18), amounts to (25) Then, the feedback of Fig. 6 is -stable and, consequently (26)

Proof: As belongs to the sector and satisfies the positive real property (24), it follows from the circle criterion (see e.g., [9, p. 265]) that the feedback of Fig. 6 is -stable. On the other hand, as is Hurwitz and is a proper transfer function 3 , one has from ( 8) and ( 19) that In the light of Proposition 2, it is seen that all signals of the control system (Fig. 3) vary slowly, in the sense that their derivatives belong to . This will be used to establish a more crucial result, namely the fact that the deviation , between the computed and applied control, belongs to . This is the subject of the following proposition.

Proposition 3: Consider the control system of Fig. 3 where: i) is any Hurwitz operator satisfying inequality (25) ii) and are solutions of (15) iii)

is any Hurwitz operator of the form (16b) 3 Recall that a transfer function is proper if the degree of its numerator is equal or less than the degree of the denominator. An improper transfer function is not physically realizable.

Then, one has (32)

Proof: In view of (12b), (11a) can be rewritten as follows: where is an exponentially vanishing term due to initial conditions. Using (11b), (61) can be rewritten as follows:

(62) Let denote any state-space realization of the (strictly proper) transfer function . Then (62) can be given the following state-space representation:

(63a) (63b) with (63c) Equations (63a)-(63c) imply (64)
The function is Lipschitz because the operator is so. Furthermore, is piecewise continuous as and are so. Then, it follows from the existence theorem (see, e.g., [9, p. 93]) that, for initial conditions, the differential equation (64) has a unique solution . Furthermore, and its derivative are both continuous. Then, it follows from (63b) that the signal is in turn continuous. The continuity of is an obvious consequence of the continuity of . Since the signal is continuous and belongs to it converges to zero (by Barbalat's Lemma, see, e.g [9, p.323]). Finally, it follows from (4) that is continuous and in view of (10) so is . But this also belongs to . Then, using again Barbalat's Lemma, the signal in turn vanishes asymptotically. Theorem 1 is established.

Remarks 4.1:

1) The fact that converges to zero means that the tracking objective is perfectly reached asymptotically.

2) As vanishes asymptotically, the nonlinear controller (11a)-(11b) stops saturating (for sufficiently large ). Consequently, the controller behaves asymptotically like the underlying linear controller (defined by ( 16)). In practice, the transient period (before the nonlinear controller stops saturating) will last a finite time after each change in the output reference signal.

3) Conditions ( 7), [START_REF] Alvarez-Ramirez | Current-mode control of DC-DC power converters: A backstepping approach[END_REF], and (25) turned out to be crucial in achieving the above two results. Condition [START_REF] Sira-Ramirez | Dynamical adaptive pulse-width modulation control of DC to DC power converters: A backstepping approach[END_REF] ensures that the tracking objective is achievable with a control signal that takes its steady-state values in the interval . Condition [START_REF] Alvarez-Ramirez | Current-mode control of DC-DC power converters: A backstepping approach[END_REF] stipulates that the output reference signal should be slowly time-varying. This is necessary to limit the control activity. Finally, the positive real condition (25) defines a neighborhood of the controlled system poles (i.e., the zeros of ) in which those of the closed-loop system (i.e., the zeros of ) should be assigned. This means that the closed-loop system is not allowed to be too much faster than the controlled system. Doing so, one limits the involved control activity and avoids long saturation periods. 4) In practice, an appropriate choice of the polynomial consists in using the following general rule:

with It is clear that if then and condition (25) is satisfied because then one has: But, letting means that the closed-loop system is not allowed to be faster than the controlled system. For the closed-loop to be faster, the parameter should be positive because then the zeros of are obtained left shifting those of with . More specifically, if has a zero at then will have a zero at . The larger is the faster the closed-loop system (provided that condition (25) holds).

5)

can also be chosen of the form with . This choice is compatible with the fact that corresponds (in the state-space interpretation of the controller (11a)-(11b)) to the poles of the observer (see e.g., [START_REF] Glattfelder | Control Systems With Input and Output Constraints[END_REF]). Indeed, it is well known that the observer poles (zeros of ) should preferably be faster that the regulation poles (zeros of ). 6) It is worth noticing that in the usual control practice the input constraints ( 6) are essentially carried out limiting the time-variation of the output reference signal and slowing down the closed-loop system rapidity. However, these requirements are generally met using heuristic rules (generally try-an-error). One major motivation of the present work is the development of a theoretical background to that usual practice and formulates rigorously accurate requirements on the allowed reference time-variation and the closed-loop rapidity.

V. R EGULATOR EVALUATION

A. Circuit Simulation

The controlled Buck circuit is simulated by its instantaneous model defined by (1a) with the following characteristics: , H, F, V. The switching frequency (used in the PWM process) is fixed bearing in mind several considerations. First, the averaged model (2a) is obtained supposing that is sufficiently large compared to the cut-off frequency of the underlying low-pass filter see e.g., [START_REF] Rashid | Power Electronics Handbook[END_REF]. In the present case, the underlying low-pass filter is a second order whose cut-off frequency is equal to . Therefore, the above averaging requirement amounts to let which numerically yields Hz. Furthermore, in practical applications where highly silent circuits are required, the switching frequency should be above the audible frequency i.e., kHz. On the other hand, to minimize the current ripple in the coil, the switching frequency must be chosen large enough. In the light of the above discussion, it is seen that a suitable choice of the switching frequency can be 200 kHz.

B. Regulator Design

The regulator design is based on the average model equations (2a)-(2b) and the associated transfer function which, in view of the above circuit characteristics, turns out to be the following:

(65)
The next step in the regulator design is to solve the Bezout [START_REF] Chaoui | Adaptive control of input-constrained type-1 plants: Stabilization and tracking[END_REF] to get the operator and . In view of (5a), [START_REF] Benzaouia | The regulator problem for linear systems with saturations on the control and its increments or rate: An LMI approach[END_REF], and (16a)-(16c), it can be easily checked that the solution of ( 15) is defined by

(66a) (66b) (66c) (66d)
In the above equations the parameters and are freely chosen. But, the parameters and should be chosen bearing in mind the positive real condition (25). To this end, one gets from (5a) and (15a) (67) For the positive real condition (25) to hold, it is sufficient to choose and such that: and either (68a) or and . But, the last two conditions are in turn equivalent to the following single inequality:

C. Experimental Setup

The experimental control setup is described by Fig. 7 and is implemented using the MATLAB/SIMULINK software.

D. Simulation Results

The behavior of the closed-loop control system of (Fig. 7) is illustrated by Figs. 8910111213141516, where the parameters that characterizes the saturation function [START_REF] Liu | Contrained control of positive discrete-time systems with delays[END_REF] are given the values 0.05 and 0.95, respectively. Fig. 8 shows the closed-loop behavior in ideal conditions (no model uncertainty) when the output reference trajectory is a signal switching between 9 V and 15 V. It is clearly seen that the output voltage tracks perfectly its reference after transient periods (following the reference signal changes).

Fig. 9 shows the computed and applied control signals, and . is seen that the former tracks the latter, after finite transient periods. That is, the controller stops saturating after these periods, confirming thus Theorem 1.

Remark 5.1: It is readily seen from Fig. 8 that the closed-loop settling time is nearly 0.5 ms. This is much smaller than the settling time of the (open loop) controlled converter. Indeed, step responses of the (average) transfer function (65) show that the converter settling time is nearly equal to 10 ms. That is, the closed loop control system is much more rapid than the (open loop) controlled circuit. It may be possible to further increase the closed-loop rapidity using different polynomials and . However, there is a limit imposed by conditions (68a)-(68b). Indeed, if these are violated the regulator will be too rapid and may generate an excessive control activity. The control limitation (13) will then cause the saturation of the regulator, if not all time, at least over long periods. Consequently, there will be no guarantee that the tracking performances of Theorem 1 will be preserved. Hence, a suitable choice of the design parameters is one which leads to the best compromise between the involved control activity and the achievement of perfect asymptotic tracking.

2) Regulator Sensitivity to Load Uncertainty (Figs. 10 and 11): In practical situations the load resistance is changing and sometimes its value may not be precisely known. The point is that such features are not accounted for in the control model (2a). Therefore, it is important to see whether the controller (11a)-(11b) is robust with respect to such load uncertainty. Fig. 10 illustrates the response of the closed-loop system to a reference step of 12 V, in presence of a load change of 50% of the nominal value. Specifically, the load switches from 1.5 to 0.75 and vice-versa. While the load resistance is changing, its value in the controller design is kept constant all the time (it is set to the nominal value 1.5 ). It is seen from Fig. 10 that the proposed regulator compensates very well the effect of load uncertainty and variations. Fig. 11 shows that the involved control activity remains acceptable and, again, the controller stops saturating after these transient periods.

3) Regulator Behavior in Presence of Discontinuous Conduction Mode Operation (Figs. 12 and13): In practice, dc-dc converters may enter into a discontinuous conduction mode operation. This means that, in each switching period, the current vanishes during a very short time interval. The point is that such phenomenon is not accounted for in the control model (2a), based on in the regulator design. Therefore, it is of interest to check whether the regulator (11a)-(11b) preserves its performances when it faces such converter behavior. To push the converter into discontinuous mode operation, a sudden and drastic change of the load is produced at time instant s (Fig. 12). This results in a drastic decrease of the current that makes the converter operate in discontinuous mode during an interval following the sudden load change. This is illustrated making a zoom on the inductor current during the interval Fig. 13. Fig. 12 shows that such disturbing phenomenon produces only a temporary and slight deviation of the output voltage from its reference value.

4) Regulator Tracking Behavior When Condition (25) is Not Satisfied (see Figs. [START_REF] Gupta | Frequency-domain characterization of sliding mode control of an inverter used in DSTATCOM application[END_REF][START_REF] Chaoui | Adaptive control of input-constrained type-1 plants: Stabilization and tracking[END_REF][START_REF] Rashid | Power Electronics Handbook[END_REF]: Here, the values of the coefficients and have been chosen in such a way that the real positive condition (25) is not satisfied. Then the regulator performances may deteriorate drastically. Figs. [START_REF] Chaoui | Adaptive control of input-constrained type-1 plants: Stabilization and tracking[END_REF] and 16 illustrate such deterioration when (72)

The simulation conditions are similar to those of Sections V-A-C. It is seen from Fig. 14 that is actually negative for some values of . Figs. 15 and 16 show that the regulator behaviour deteriorates significantly (compared to Figs. 8 and9). In particular, it is observed that: (i) the reference-output tracking performance is not satisfactory; (ii) the controller never stops saturating.

VI. CONCLUSION

The control limitation ( 6) is inherent to all power converters operating according to the PWM principle, whatever the control design approach. In the present paper, such issue has been dealt with, for converters of the Buck type, using the nonlinear controller (11a)-(11b). This involves an internal loop around the saturation function , see Fig. 3. It contrasts with previous solutions (Fig. 5) that consist in first designing a (linear) controller as if there was no control limitation and then incorporating an isolated limiter, placed between the controller and the controlled system. The point is that the effect of such a limiter is also ignored, in previous works, when analyzing the closed-loop control system.

The control solution developed in the present paper has been analyzed using tools from the absolute stability theory. The key step is to have shown that the resulting (nonlinear) closed loop system (Fig. 3) gives rise to the feedback system of Fig. 6. This is composed of a linear direct path and a nonlinear feedback path of the conic type. Then, -stability of this feedback system is guaranteed provided that satisfies the real positive condition (25). In fact, this defines a neighborhood of the controlled systems poles (zeros of ) in which the closed-loop poles (zeros of ) should be assigned. For the sake of simplicity, the pole placement design technique has been resorted to meet such requirement. But other linear design techniques could as well be used. The stability of the feedback of Fig. 6 makes it possible to establish the tracking objective for a class of slowly varying reference signals defined by ( 7)- [START_REF] Alvarez-Ramirez | Current-mode control of DC-DC power converters: A backstepping approach[END_REF]. The proposed regulator shows quite interesting robustness properties especially against uncertain load changes and discontinuous conduction mode operation.

Fig. 1 .

 1 Fig. 1. Basic control setup of Buck converters.

Fig. 2 .

 2 Fig. 2. Buck converter circuit.

Fig. 3 .

 3 Fig. 3. Closed-loop system obtained with the nonlinear controller (11a)-(11b) that incorporates a control limiter.

Fig. 4 .

 4 Fig.4. Closed-loop system obtained with the standard linear controller (14) that involves no control limitation.

Fig. 6 .

 6 Fig. 6. Nonlinear feedback relating the signals sv and s.

  which, together with (20), implies that (23) This proves Proposition 1.

  follows that all internal signals of the feedback belong to . This particularly means that (28a) (28b) Noting that is Hurwitz, it follows from (10), (27), and (28a) that (29) Also, from (8) one gets (30) which together with (8) and (29), gives (31) This ends the proof of Proposition 2

  ) and (29) and the fact that is Hurwitz, it follows from (34) that (35) On the other hand, (4) can be transformed as follows: 28a)) and the transfer function is Hurwitz and strictly proper, it follows from (38) that . Then, using (31), it follows from (37) that (39)On the other hand, multiplying (36and the applied controls, both vanish asymptotically Proof: Let us first show that vanishes asymptotically. As , it is sufficient to show that and are continuous. From (17) one gets (61)

Fig. 7 .

 7 Fig. 7. Experimental setup used for regulator evaluation.

1 )Fig. 8 .

 18 Fig.8. Tracking behavior of the controller in presence of a time-varying output reference signal y (t) switching between 9 and 15 V.

Fig. 9 .

 9 Fig. 9. Comparison between the computed duty ratio (t) and the applied value (t) in presence of a time-varying input reference.

Fig. 10 .

 10 Fig.10. Robustness of the controller with respect to load resistance changes.

Fig. 11 .

 11 Fig. 11. Comparison between the computed duty ratio (t) and the applied value (t) in presence of load resistance changes.

Fig. 12 .

 12 Fig. 12. Regulator behaviour in presence of discontinuous conduction mode.

Fig. 13 .

 13 Fig.[START_REF] Liu | Contrained control of positive discrete-time systems with delays[END_REF]. Inductor current i during the discontinuous conduction period.

Fig. 14 .

 14 Fig. 14. Real part of C (j!)=A(j!) versus !.

Fig. 15 .

 15 Fig. 15. Deterioration of tracking performances when condition (25) is not satisfied.

Fig. 16 .

 16 Fig.[START_REF] Rashid | Power Electronics Handbook[END_REF]. Control activity when condition (25) is not satisfied: (t) computed value of the duty ratio, (t) the applied value.

When clear from the context, the symbol 's' will also used to designate the Laplace Transform variable. For instance, the transfer function of the system (4) is simply denoted B(s)=A(s).

With the above notation, (41) can be rewritten The rest of the proof consists in showing that . To this end let us consider the three cases suggested by [START_REF] Liu | Contrained control of positive discrete-time systems with delays[END_REF]. . Equations (11b) and ( 13) yield and, consequently (53) Case 3

. Equation (11b) and ( 13) give that and . Also, (43a) implies that Propositions 3 and 4 are key ingredients to get the main results of this paper i.e., the tracking objective and the fact that the control signal stops saturating asymptotically. This is properly formulated in the following theorem.

Theorem 1: Consider the closed-loop system sketched by Fig. 3 andcomposed of: i) the Buck power converter of Fig. 1, represented by ( 4) and subject to the control saturation constraint (6); ii) the nonlinear controller (11a)-(11b), where satisfies conditions ( 7)-( 8) and the operator satisfies the real positive condition (25).


Carlos L Pieck

Estanislao L Jablonski

Jose M Parera

Roger Frety

Frédéric Lefebvre

Frederic Lefebvret

Characterization of residual coke during burning

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

The carbonaceous deposit, or coke, formed on the naphtha re-forming catalyst during operation has a large molecular weight, containing polyaromatic hydrocarbons with paraffinic substitutions, and is poorly organized into pseudographitic structure [START_REF] Biswas | The Reformer Lineout Phenome non and its Fundamental Importance to Catalyst Deactivation[END_REF][START_REF] Parera | Coking on Bifunctional Catalysts[END_REF][START_REF] Caruso | Crystallinity of Coke on Platinum-Rhenium/ Alumina Reforming Catalyst during the Commercial Cycle[END_REF]. The degree of condensation of coke and its nature depend on the catalyst, type of feed, and operation parameters in the re-forming process [START_REF] Beltramini | In fluence of Pt Concentration on Activity and Combustion of Coke on Pt/ Al203[END_REF]Barbier et al., 1985a,b;[START_REF] Franck | Deactivation of Reforming Catalysts[END_REF][START_REF] Parmaliana | Platinum Reforming Catalysts: Effect of Chlorine Content on Coking and Catalyst Self-Regeneration[END_REF][START_REF] Zhorov | Kinetics of Platinum Catalyst Deactivation by Coke Deposits[END_REF].

It has been shown that coke consists basically of two types of carbonaceous deposits: one of a low polymeriza tion degree (H/C atomic ratio about 1) deposited on the metal, and the other of a higher degree of polymerization (H/C ratio about 0.5) deposited on the support [START_REF] Parera | The Influence of Coke Deposition on the Func tions of a Pt/ Al203-Cl Bifunctional Catalyst[END_REF][START_REF] Barbier | Selective Poisoning by Coke Formation on Pt/ Al203[END_REF]Barbier et al., , 1985a,b),b). Coke is elimi nated from the catalyst by burning with a diluted oxygen gas. The mechanism governing the burning process of the coke deposited in pelleted catalysts changes with tem perature [START_REF] Pieck | The Burning of Coke on Pt-Re/ Al203[END_REF]. At high temperatures, the burning follows the shrinking core model; while at a low temperature, burning occurs according to the different reactivity of coke components and oxygen pressure.

In this paper, the coke remaining after partial burning at different temperatures was studied. The coke of com mercially used catalyst pellets was partially burnt at various temperatures with analysis of the gaseous effluent during burning. After this partial burning, the catalyst pellets were finely ground and submitted to a tempera ture-programmed oxidation (TPO). After the coke was partiall y burned, the remaining coke was concentrated by acid elimination of the support and studied by X-ray diffraction (XRD) electron diffraction (ED), IR, 13C CP MAS NMR, electron spectroscopy for chemical analysis (ESCA), electron paramagnetic resonance (EPR), and elemental C and H analysis.

Experimental Section

Catalyst. A Pt-Re-S/ Al20 3 catalyst (1.5-mm-diameter and 5-mm-length extruded pellets) coked during a com mercial re-forming operation was used. When the opera tion was started, the catalyst had a specific surface area t lnstituto de Investigaciones en Catalisis y Petroquimica (IN CAPE).

i Conventionne a l'Universite Claude Bernard Lyon I.

of 187 m 2 g-1, a pore volume of 0.51 cm3 g-1, and the fol lowing composition by weight: 0.3% Pt, 0.3% Re, 0.96% Cl, and 0.024 % S. The catalyst was in operation for 6 months with a paraffinic naphtha at a pressure of 1.35 MPa and during this period accumulated 12.2% C in weight.

Coke Partial Burning. Coke partial burning was performed in equipment already described [START_REF] Sad | Evaluation of Activity, Selectivity and Stability of Catalysts for Naphtha Reforming[END_REF] that was modified to allow the entrance of the ox idant mixture. The procedure was as following: 4 g of the coked catalyst pellets were heated up to the burning tem perature in a nitrogen flow, and then nitrogen was replaced by the oxidant mixture (1.9 vol % 02 in N2) during 2 h.

The burning temperature varied between 623 and 823 K, the gas pressure was 0.1 MPa, and the gas flow rate was 125 mL(STP) min-1• Temperature-Programmed Oxidation. A Shimadzu Thermoanalyzer DT 30 was used, charging 6 mg of fine ly-ground catalyst. The heating rate was 20 K min-1, and the flow rate of the previously dried oxidant mixture (1.9 vol % 02 in N2) was 50 mL min-1• The same catalyst completely free of coke was used as reference.

Coke Concentration. For attainment of adequate signal levels, the concentration of coke was increased by dissolving alumina in an acid mixture, following a proce dure proposed by other authors [START_REF] Magnoux | New Technique for the Characterization of Carbonaceous Compounds Responsible for Zeolite Deactivation[END_REF]. A 1-g amount of the partially burnt catalyst ground to 25-80 mesh was attacked under vigorous stirring with a solution of 2 mL of 37 wt % hydrochloric acid and 8 mL of 40% hydrofluoric acid at room temperature for 6 h. After 1 night at room temperature, the suspension was filtered and the solid washed with distilled water, until the acids were eliminated, and dried at room temperature. The acid attack was repeated until no alumina was detected in the coke by XRD. Probably, this acid treatment causes no modification of the coke structure since coke in zeolites have been shown to be unaffected by zeolite dissolution in HF [START_REF] Magnoux | New Technique for the Characterization of Carbonaceous Compounds Responsible for Zeolite Deactivation[END_REF].

X-ray Diffraction. XRD spectra were obtained using a Rich Seifert ISO-Debyeflex 2002 powder diffractometer at 40 kV and 20 mA with nickel-filtered Cu Ka radiation (wavelength = 1.54 A). With the 28 values of the 002 plane peaks and the peak width measured at half the maximum height taken from the X-ray spectra, characteristic pa rameters of coke can be calculated.

Mean interlayer spacing of the aromatic system can be calculated with the Bragg equation:

d = _x_ (1) 2 sin 8
where X (wavelength) = 1.54 A.

Mean thickness of the aromatic layers can be calculated with the Scherrer equation (Anderson, 1985):

L = kX
(2) Bj cos 8where k = Scherrer constant for which generally a value of 0.90 is taken for crystalline materials (Anderson, 1985) and Bj = peak width measured at half the maximum, in radians.

Infrared Spectroscopy. A Perkin-Elmer (Model 580B) double-beam spectrophotometer was used. The instrument range is 180-4000 cm-1• Coke samples were pressed into disks with KBr.

Nuclear Magnetic Resonance. 13C CP-MAS NMR spectra were recorded on a Brucker MSL-300 spectrometer equipped with a double bearing probe head and operated at 75.47 MHz. The spinning frequency was typically 3 kHz, and in order to prevent a superposition of spinning side bands and aliphatic peaks, the TOSS sequence allowing a suppression of the former-was used. The contact time for cross polarization was 1 ms (the 90° value for 1H is 6.2 µs), and the delay between each scan was 10 s. Depending on the amount of the product, 360-10000 scans were accumulated.

Electron Diffraction. Electron diffraction patterns were recorded with a JEOL 100 CS instrument. Samples were suspended in water and mounted on a 400-mesh copper grid with a nitrocellulose film.

Chemical Analysis by ESCA. XPS spectra were re corded with a Shimadzu ESCA-750 electron spectrometer using magnesium Ka X-rays. The instrument working range is 0-1150-eV binding energy. Deconvolution of data was performed through a Shimadzu ESCAP AC 760 data system. For this analysis, samples were pressed at 0.35 MPa before the measurement.

Hydrogen-to-Carbon Ratio. Coke concentrated by the dissolution of alumina, using the acid treatment quoted above, was analyzed microgravimetrically for carbon and hydrogen.

In some cases, the H/C ratio of the burnt coke was obtained by chromatographic analysis of the combustion gases during the burning, using the method described elsewhere (Barbier et al., 1985a,b). From the coke com bustion reaction, the H/C ratio can be deduced:

CH,+ ( 1 + �)02--+ C02 + �H20 (3) H ( 02 consumed ) C = x = 4 C02 produced -1 (4)
In most cases, quantities of 02 consumed and C02 pro duced are similar; then their ratio is near 1. Therefore, when eq 4 is used, great errors are possible. Another problem of this method is the consumption of oxygen without combustion when burning starts at a low tern• perature: when oxygen is fed, combustion starts slowly and a certain amount of oxygen is consumed by adsorption on the equipment and on coke. As little C02 is evolved, the calculation gives a higher H/C ratio than the correct value.

At higher temperatures, the oxygen adsorbed or reacted with coke produces C02 without gaseous 02 consumption, giving the calculation a lower H/C ratio than the correct value. For this reason, only the average value including all the combustion gases can be taken as a correct value. Electron Paramagnetic Resonance (EPR). A Bruker Model ER-200 spectrometer operating in band X (9.7 GHz) at room temperature was used.

Results and Discussion

Figure 1 shows the TPO of the original coked catalyst, A, and of samples B and C obtained when previous TPO runs were stopped at 698 and 773 K, respectively. The original coked catalyst presents two burning zones. The first zone-at 573-723 K-corresponds to the burning of the coke over the metallic function (Pt-Re). The second zone corresponds to the burning of the coke over the acid function (Al203-Cl) that is more polymerized and poorer in hydrogen [START_REF] Parera | The Influence of Coke Deposition on the Func tions of a Pt/ Al203-Cl Bifunctional Catalyst[END_REF][START_REF] Barbier | Selective Poisoning by Coke Formation on Pt/ Al203[END_REF]Barbier et al., , 1985a,b),b). In sample B, all the coke on the metal and a very small fraction on the support were eliminated. In sample C, the burning of coke over the support continued. This shows that the TPO is selective; the coke on the metal is burnt first, and the burning continues with coke increas ingly difficult to burn each time (more polymerized) and located on the support. TPO is performed with a fine powder and there is no mass-transfer limitation, at least when the temperature is not very high.

In the commercial burning of the coked catalyst pellets, the process starts at a low temperature and with a very low oxygen concentration; consequently more hydrogen ated coke (located on the metal) is burned. When this burning is already finished, temperature and oxygen con centration increase slowly, thus allowing the burning of more coke. The burning starts selectively, controlled by the chemical step, but when the temperature is high, the controlling step can change. The oxygen diffusion can become the limiting step, and the burning will be nonse lective, with all the coke of external shells of the catalyst being eliminated first. When the experiment is performed isothermally, as in our case, a selective burning can be obtained at a low temperature and a nonselective one at a high temperature.

Table I shows the carbon remaining on the catalyst after burning during 2 hat various temperatures, the H/C at omic ratio of the residual coke, the interlayer spacing of the aromatic layers, and the mean thickness of these layers obtained by XRD in the coke concentrated after elimi nation of alumina by acid treatments.

During the 2-h burning, the catalyst bed temperature remained constant because of the great thermal capacity of the system: the catalyst bed has a small ring-shaped cross section (between thermocouple well and reactor wall), the reactor wall is 2 cm thick, and the gas has 98.2% N2•

T he measured burning temperature is the average parti cle-gas temperature, but the actual temperature inside the catalyst particle could be higher. The rate of burning is rather small because in all cases a great amount of carbon remains after the 2-h burning.

The original coke has an H/C ratio of 0.56, and the residual coke, when the partial burning is performed at 673 K, has a ratio of 0.4. At higher temperatures, the residual coke has the same H/C ratio as that of the original coke.

This can be explained considering that coke burning is selective at low temperatures (up to 723 K), burning first species richer in hydrogen, leaving on the catalyst a heavier coke Oower H/C ratio). [START_REF] Sotirchos | Combustion of Coke Deposits in a Catalyst Pellet[END_REF] found that hydrocarbonaceous species are five times more reactive to oxygen than the graphitic ones. At high burning tem peratures (773-823 K), all coke components are simulta neously burnt, and the residual coke has a H/C ratio similar to values of the original coke. This agrees with the results of [START_REF] Pieck | The Burning of Coke on Pt-Re/ Al203[END_REF]. who stated that, at low tem peratures, the burning is nearly uniform in all the catalyst particles (controlled by the chemical step) and, at high temperatures, the burning follows the shrinking core model (controlled by the oxygen diffusion step).

It can be see n (Table I) that the mean interlayer spacing is practicall y independent of the burning temperature and larger than the one of graphite (d002 = 3.354 A). The interlayer spacing is larger than that in graphite because coke is a disordered carbonaceous material formed by polyring aggregates with alkyl chains, as quoted before [START_REF] Bakulin | Transformation of Cyclic Hydrocarbons on an Aluminoplatinum Catalyst[END_REF][START_REF] Parera | Mechanism of Coke Formation During Naphtha Reforming[END_REF]. The mean thickness of the aromatic layers as a function of burning temperature has a maximum.

Figures 2 and3 show the electron diffraction patterns of the sample without burning treat ment and of the sample partially burnt at 673 K, respectively. In a comparison of all diagrams, the sample without treatment shows a n amorphous structure, while the residual coke of the burning at low temperatures presents ordered structures more evident in the sample burnt at 673 K. This means that, at a low temperatures, the smaller and less ordered coke particles burn first. At a high burning temperature and because of the shrinking core mode of burning, the residual coke has a structure similar to that of the original material before treatment.

Figure 4 shows m (1000-1800-cm-1) spectra of the ori ginal coke and of that from the partially burnt samples. Using the equipment computer, it was given the same intensity to the peak at 1680-1480 cm-1 in all the spectra in order to take it as a reference to compensate differences in coke concentration and thickness of the disks. The absorption band appearing between 1480 and 1680 cm-1 in all samples can be ascribed to the stretching of the C=C bond belonging to olefins, aromatic rings, and poly aromatics. The peak at about 1400 cm-1 can be ascribed to the CH bond stretching in C-(CHs)2 or C-(CHs)3 groups.

The absorption peak at about 1200 cm-1 is due to the stretching of the CH group in the aromatic plane [START_REF] Rouxhet | Characterization of Kero gens and of Their Evolution by Infrared Spectroscopy[END_REF].

In a comparison of the absorption spectra of all samples, it can be seen that there is no absorption due to C-OH (3400-3600 cm-1, not shown in Figure 4, and 1100 cm-1) and to C=O (1710 cm-1) in the sample without burning treatment. Meanwhile, samples partially burnt present an absorption peak approximately at 1715 cm-1 due to the group C=O and another peak at 1100 cm-1• This last one can be due to the presence of OH, because a broad peak in the range 3400-3600 cm-1 also appears. It can be seen that the intensity of the peak produced by C=O has a maximum value at 673 K, which corresponds to an in termediate coke burning. Then, coke burning follows the following consecutive reaction coke + 02 -+ oxidized coke -C02 + H20

(5)

which is similar to the model proposed by Furimsky et al.

( 1988) for coke deposited on nickel molybdate hydro treating catalyst.

The peak in the range 1480-1680 cm-1 moves to a higher frequency from the original sample to the one burnt at 673 K, moving then to lower frequencies at high burning tem peratures (723-823 K). This could indicate that the structure of residual coke should be more organized in the sample burnt at 673 K because of its highest vibration energy, being in agreement with the results of transmission electron microscopy (TEM) analysis above.

Solid-state NMR spectra show only three broad ab sorptions around 125, 73, and 25 ppm attributed to olefinic and aromatic carbons (125 ppm), carbon singly bonded to oxygen (75 ppm), and -CH2CH3 peaks (25 ppm). The main difference between spectra of coke samples is the relative amount of oxygen bonded to carbon: this peak can be very easily observed in samples burned at 623 and 673 K. In the same range of temperatures, IR data show the greatest absorption peak due to the C=O group, the position of the maximum not being the same. In NMR, the maximum is in the sample burnt at 623 K, while in IR spectra, the C=O group at 1715 cm-1 shows a maximum in the sample burnt at 673 K. This may be due to a different structure of the oxidized species observed in each case. If there is less hydrogen in the neighborhood of the C-0 group, the intensity measured by CP-MAS will be lower. No carbon in the C=O group was detected in NMR, because there is no hydrogen directly bonded to the carbon, and so cross polarization cannot be efficient in this case. The maximum in C=O cannot be detected in NMR.

C ls and 0 ls spectra of the samples were recorded; the C ls peak of the sample was near that of the high vacuum pump oil and was resolved by deconvolution. The 0 ls peak is due to oxygen chemically bonded to the sample (C=O and C-OH). Table II shows the values of binding energy and peak width of C ls, and the values of factor g obtained by EPR. It can be seen that the peak width, its position, and factor g have a maximum value. Ac cording to [START_REF] Wagner | Handbook of X-Ray Photoelectron Spectroscopy[END_REF], the shift toward higher values of binding energy of C ls is probably due to the formation of carbon-oxygen species in which the C ls peak has higher values of binding energy. The peak broadening can be due to the higher heterogeneity of the samples due to the simultaneous presence of oxidized and nonoxidized coke. The sample burnt at 673 K was the only one in which oxygen of oxidized coke was found. The calculated O/C atomic ratio was 0.02. In the other samples, the oxygen surface concentration was below the detection limit of the equipment.

Changes in binding energy and peak width of C ls (and also the presence of 0 ls in the sample regenerated at 673 K) with the burning temperature agree with IR results

(presence of C=O and C-OH passing through a maxi mum as a function of burning temperature) and the other techniques (TEM, XRD). Therefore, at a low oxidation temperature, an oxidized coke is produced; at a high tem peratue, the period of life of C=O species is reduced, and the burning rapidly produces C02 and H20.

Taking into account the peak appearing around 1715 cm-1 in the IR spectrum (C=O groups) and the aromatic ratio O/C of the sample burnt at 673 K in the ESCA spectrum, it is possible to calculate oxygen surface con centration for each sample. The values in weight percent are 1.6, 2.6, 1.3, 0.5, and 0.2 for samples partially burnt at 623, 673, 723, 773, and 823 K, respectively. From ESR data, the sample burnt at 673 K has the highest value of g, probably due to the presence of oxy genated groups in the proximity of aromatic groups (Marchand and Canard, 1980).

In Figure 5, the H/C atomic ratio of the fraction of coke eliminated is plotted as a function of the percentage of this coke. The combustion was performed at 0.34 MPa with a flow rate of 42 mL min-1 of oxidative mixture containing 5 vol % 02 in N2 and at two temperatures, 698 and 823 K. To calculate the H/C ratio, average concentrations of C02 and 02 were chromatographically determined in the period at which the coke was eliminated. At 698 K, the H/C value decreases when the amount of coke eliminated increases, but it rems.ins constant at 823 K. This agrees with concepts already quoted: at a low temperature, the burning is selective, the more hydrogenated coke being eliminated first, while, at a high temperature, the burning is nonselective following the shrinking core model.

Conclusions

Different techniques agree regarding the phenomenon of burning and the characteristics of coke during its burning. At low temperatures, the burning is selective, the more hydrogenated fraction with smaller aggregate sizes being first eliminated, and surface C=O and C-OH groups are formed as intermediates. At higher tempera tures, these groups rapidly complete the oxidation to C02, and the burning is nonselective: all coke components are burnt from the external layers to the inner part (shrinking core model).

  Figure 1. (A) Temperature-programmed oxidation (TPO) of the original commercially coked catalyst. (B) TPO of a sample at which the TPO was previously performed up to 698 K. (C) TPO of a sample with previous TPO up to 773 K.

Figure 3 .Figure

 3 Figure 3. Electron-diffraction pattern of coke from the sample partially burnt at 673 K.

Figure 5 .

 5 Figure 5. H/C atomic ratio of the fraction of coke elimination as a function of the percentage of this coke burnt at ( e) 698 and (•) 823 K.

  Registry No. C, 7440-44-0; Pt, 7440-06-4; Re, 7440-15-5.Literature CitedAnderson, J. R. Measurement Techniques: Surface Area, Particle Size and Pore Structure. Structure of Metallic Catalysts; Aca-

Table I .

 I Characteristics of Coke Remaining over the

	Partially Burnt Coked Catalysts
		residual	atomic	interlayer	mean
	burning temp, K sample without	carbon, 3 12.2	ratio (H/C) 0.56	spacing, thickness, A A 3.41 25.6
	burning
	623 673 723 773	10.0 8.5 7.3 6.8	0.52 0.40 0.48 0.56	3.40 3.42 3.42 3.42	28.9 29.1 30.2 28.1
	823	5.2	0.57	3.43	21.4

Figure 2. Electron-diffraction pattern of coke from the sample without burning treatment.

Table II .

 II Binding Energy of C ls, Peak Width, and g Values of Coke from the Partially Burnt Coked Catalysts

		binding	full width at
	burning temp, K	energy, eV	half-maximum, eV	g
	sample without burning	286.9	3.6	2.0025
	623	286.8	3.3	2.0027
	673 723	287.2 287.2	4.0 3.0	2.0031 2.0025
	773	286.8	2.7	2.0025
	823	286.6	2.5	2.0024


Pierre Haessig
email: pierre.haessig@centralesupelec.fr

Convex Storage Loss Modeling for Optimal Energy Management

Keywords: convex optimization, energy management, energy storage, storage losses, storage models

Models of energy storage systems used for optimal energy management commonly feature unsophisticated loss expressions (lossless or constant efficiency). Preserving the tractability (i.e. convexity) of the optimization motivates this simpleness, at the expense of physical realism. Relaxing the storage loss expression from an equality constraint to an inequality allows using more complex and physically realistic models while preserving convexity. Although this relaxation has been used quite casually for a decade in varied energy management applications (power systems dispatch, electric vehicles), only few articles studied this technique for itself. To increase the awareness of loss relaxation, along with its limitations, we depict a unified view of relaxed storage loss models. Also, we propose the "convex monomial loss model", one continuously parametrized family of nonlinear convex expressions which contains all classical loss models and unlocks new possibilities for loss dependency in storage power and energy level. This is one step to reunite piecewise linear and nonlinear models which are generally studied separately. Finally, we compare the effect of these varied loss models on the optimal charge/discharge profile of a storage in an elementary energy management application.

I. INTRODUCTION

Energy storage systems (ESSs), in particular batteries, are being used in a large and growing number of applications, down from small portable devices up to large power systems. Energy management is the specific name for the control of the power flows in a system which includes one or more energy storages, to achieve desired objectives like energy shifting, power shaving, backup... Energy management can sometimes be achieved by simple rules, but many storage applications have complex enough requirements to justify using mathematical optimization tools instead. A large subset of optimization-based control methods, like Model Predictive Control, require solving optimization problems online. This necessitates solving those optimizations quickly and reliably and this reliability depends essentially on the class of optimization problem being solved.

This need pushes control designers to formulate their problems as Convex Programs (CPs), because it is a wide class of problems for which reliable solvers exist, in particular for the common Linear and Quadratic subclasses (LPs and QPs) [START_REF] Boyd | Convex Optimization[END_REF]. True, there also exist popular and useful classes of nonconvex problems, mainly the previous ones expanded with some integrality constraints. The resulting problems (MILP, MIQP. . . ) benefit from powerful solvers, but fundamentally don't share the nice convergence properties of convex ones. This justifies the search for convex formulations of energy management problems, if physically possible.

As detailed in section II-A, for a model-based energy management problem to be convex, storage loss expression must be linear. This is very restrictive (see II-C) and means that in many applications, storage losses are simplified (e.g. [START_REF] Chen | Aggregate Power Flexibility in Unbalanced Distribution Systems[END_REF]) or fully neglected just to satisfy convexity and not for physical reasons.

To overcome this limitation, befitting constraint relaxation schemes have been used. Indeed, relaxation is a typical method to "convexify" a variety of problems, for example Optimal Power Flows (e.g. [START_REF] Gemine | Active network management for electrical distribution systems: problem formulation, benchmark, and approximate solution[END_REF]). For energy management, these methods all come to relaxing the storage loss expression from an equality constraint (losses = some expression) to an inequality (losses ≥ some expression) [START_REF] Murgovski | Component sizing of a plug-in hybrid electric powertrain via convex optimization[END_REF]. The rationale is that, at the optimum, no energy should be wasted, so the inequality should be tight. We call this the "positive energy price argument" and come back to it further down. Switching to an inequality allows using any nonlinear convex expressions for losses, while keeping the resulting energy management problem convex. Then, the choice of loss expression depends on both the application and the type of storage technology being modeled.

The most popular model (beyond the lossless case mentioned above) is, by far, the piecewise-linear-in-P expression (hereafter PWL-in-P , e.g. ∝ |P |), which is generally transcribed using a constant charging/discharging efficiency (η < 1) (see [START_REF] Chen | Aggregate Power Flexibility in Unbalanced Distribution Systems[END_REF], [START_REF] Almassalkhi | Model-Predictive Cascade Mitigation in Electric Power Systems With Storage and Renewables-Part I: Theory and Implementation[END_REF]- [START_REF] Mégel | Hybrid stochasticdeterministic multiperiod dc optimal power flow[END_REF] for a tiny subset of use cases). This is a physics-free and generic storage model which is typically formulated as a LP, while it is a non-convex MILP without the relaxation (see III-A). The recent work of Gonzalez-Castellanos et al. [START_REF] Gonzalez-Castellanos | Detailed li-ion battery characterization model for economic operation[END_REF], although being a much wider effort of convex battery modeling (not just losses), also yields a closely related model (PWL-in-P with a linear effect of the energy level).

Inspired by resistive heating, a quadratic-in-P model (i.e. ∝ P 2 ) is more physical for electrochemical storage technologies. For optimal energy management, these devices can be represented by the circuit model visible on the right-hand part of Fig. 1. Convex optimization of this model was proposed by Murgovski et al. [START_REF] Murgovski | Component sizing of a plug-in hybrid electric powertrain via convex optimization[END_REF] in a battery intensive application (Hybrid Electric Vehicle) along with a detailed, first of its kind, discussion on the "convexification" the circuit model. Also, it's maybe because quadratic terms often appears in power flow equations that we found a resembling model in a distribution grid control application [START_REF] Nick | Optimal allocation of dispersed energy storage systems in active distribution networks for energy balance and grid support[END_REF]. Section III-B is devoted to this model. Now, the most complex and genuinely nonlinear loss expression we found is in the follow-up work of Murgovski et al. [START_REF] Murgovski | Convex modeling of energy buffers in power control applications[END_REF] and subsequently used by Pinto et al. [START_REF] Pinto | Evaluation of advanced control for li-ion battery balancing systems using convex optimization[END_REF]. Their P 2 /E loss expression (with E being the storage energy level) allows modeling storage technologies for which the open circuit voltage (v 0 on Fig. 1) varies with the energy level E. We discuss and generalize this in section IV.

This literature review shows that the practice of relaxing storage loss expressions exist. Still, it is not ubiquitous (see e.g. [START_REF] Li | Microgrid sizing with combined evolutionary algorithm and milp unit commitment[END_REF] for a MILP energy management where the integrality constraints associated to the storage model may have been relaxed, although it may not be the case for other components).

Also, this relaxation is considered as natural enough that often the "positive price argument" is only lightly mentioned, if not omitted, although it cannot apply sometimes. Indeed the relaxation can fail in the sense that the relaxed loss inequality can be not tight so that more energy is wasted than prescribed by the loss expression. A notable exception to this casualness is Almassalkhi et al. [5, App. A] which not only explain the argument, but also give a worst-case bound for when the relaxation is not tight. Also, they propose a heuristic algorithm to fix this. Alternatively, Z. Li et al. [START_REF] Li | Extended sufficient conditions for exact relaxation of the complementarity constraints in storage-concerned economic dispatch[END_REF] and Duan et al. [START_REF] Duan | Improved sufficient conditions for exact convex relaxation of storage-concerned[END_REF] identified mathematically sufficient conditions for the relaxation to be exact. Still, these conditions do not always apply [START_REF] Mégel | Hybrid stochasticdeterministic multiperiod dc optimal power flow[END_REF]. Recently, Arroyo et al. [START_REF] Arroyo | On the use of a convex model for bulk storage in mip-based power system operation and planning[END_REF] refuted two Transactions articles which used too simple conditions! We thus feel the need to shed one more light on these useful storage loss relaxation techniques, while carefully referencing the limitations.

Finally, the literature shows a complete split between the popular physics-free PWL-in-P loss model and the more elaborate P 2 or P 2 /E models. However, when searching for a convex model fitted to experimental data, it may be useful to have a continuously parametrized family of loss expressions. Thus, it is our core contribution in this paper to:

1) Describe the relaxation of storage losses in a unified way which handles the various existing loss models 2) Provide a simple continuous family of nonlinear convex loss expressions which can depend on both storage power and energy level. This family contains existing models (PWL, quadratic, P 2 /E) as special cases. The unified description of loss models is covered in sections II and III, while our new "convex monomial loss model" is described in IV. Finally, in section V we create an elementary energy management application with one energy storage. We use this example to compare the effect of the varied loss models on the optimal charge/discharge profile of the storage.

II. CONVEX STORAGE MODELING WITH LOSS RELAXATION A. Recall of Convexity

We recall the conditions for an optimization problem to be convex (see [4, §2] or textbook [START_REF] Boyd | Convex Optimization[END_REF] for more background). Optimization problems, as often used in energy management, can be formalized as:

min x∈R n f (x), such that g i (x) ≤ 0, h j (x) = 0 (i, j = 1, 2 . . . )
For this problem to be convex, the cost function f must be convex and so must be the set of feasible solutions {x such that g i (x) ≤ 0 and h j (x) = 0, ∀i, j}. The latter implies that each inequality function g i must be convex and each equality function h j must be linear. This will be crucial to our discussion, because few practical storage loss expressions are linear.

B. Generic Storage Model

We start with the following generic energy storage model:

E b (k + 1) = E b (k) + (P b (k) -P loss )∆ t (1)
It is expressed in discrete time (instant k, timestep ∆ t ) and describes how the stored energy E b evolves with the storage power P b (in receptor convention) including losses P loss . Because all variables are discrete signals, we drop below the time index k for relations applying for all k. This is the most common type of model for energy management optimization, while electrochemical models are more often expressed with charge, current and voltages (see Fig. 1). The conversion between the two is discussed further (section III-B).

The stored energy is bounded, which is typically modeled with:

0 ≤ E b ≤ E rated , (2)
where E rated is the maximum usable energy of the system. From it, one can define the state of energy

SoE = E b /E rated ∈ [0, 1].
Storage power is also bounded by

-P - rated ≤ P b ≤ P + rated . (3)
Constraint (3) allows asymmetrical charge and discharge rate limits. A rich extension is to consider a more general convex set for the pair of variables (P b , E b ), instead of the rectangular box defined by ( 2)-(3), see for example [START_REF] Gonzalez-Castellanos | Detailed li-ion battery characterization model for economic operation[END_REF], [START_REF] Murgovski | Convex modeling of energy buffers in power control applications[END_REF].

Equations ( 1)-(3) describe a convex set of constraints to be embedded in a larger optimization problem involving E b (k), P b (k), P loss (k) (k = 1...K) as decision variables.

However, storage dynamics (1) is generic because the storage losses P loss have yet to be specified. So the convexity of the complete problem depends on the loss expression.

C. The Limited Usability of Linear Loss Expressions

The simplest convex model is the lossless storage

P loss = 0 (4)
which, of course, in practice is used to simplify [START_REF] Boyd | Convex Optimization[END_REF].

The only convex generalization of this, assuming losses depend on power and state of energy, is the linear expression

P loss = p 0 + c P P b + c E E b . (5)
Parameter p 0 models a constant self-discharge, while coefficient c E modulates the discharge with the state of energy. Both are physically meaningful and have some popularity (e.g. [START_REF] Chen | Aggregate Power Flexibility in Unbalanced Distribution Systems[END_REF]). In contrast, we have never seen the use of the c P coefficient. Indeed, since power P b can take both signs, there would be some physical contradiction in assuming that losses may be reduced for large charge (c P < 0) or discharge (c P > 0) rates.

D. Generic Relaxation of Storage Losses

The limited physical applicability of the linear loss model yielded befitting relaxation techniques (e.g. [START_REF] Murgovski | Component sizing of a plug-in hybrid electric powertrain via convex optimization[END_REF], [START_REF] Almassalkhi | Model-Predictive Cascade Mitigation in Electric Power Systems With Storage and Renewables-Part I: Theory and Implementation[END_REF], [START_REF] Li | Extended sufficient conditions for exact relaxation of the complementarity constraints in storage-concerned economic dispatch[END_REF]). The generic formulation, expressed again here with a dependency on storage power and energy 1 , is:

P loss ≥ g(P b , E b ) (6)
where g means here any convex function. This formulation rests on three principles, sorted by decreasing order of givenness:

1) Mathematical fact: with g convex, (6) indeed defines a convex inequality constraints as in II-A 2) Physical observation: the fact that storage losses generally increase for large positive and negative values of power P b and are minimal for P b = 0 motivates using a function g convex in P b (e.g. quadratic, absolute value. . . ) 3) Application-dependent hypothesis: most optimal energy management problems embed the property that the optimum is reached by minimizing wasted energy, so that inequality [START_REF] Gill | Dynamic Optimal Power Flow for Active Distribution Networks[END_REF] should be tight at the solution. We call this the "positive price argument", which is often mentioned in various ways (e.g. [START_REF] Gill | Dynamic Optimal Power Flow for Active Distribution Networks[END_REF], [START_REF] Nick | Optimal allocation of dispersed energy storage systems in active distribution networks for energy balance and grid support[END_REF], [START_REF] Pinto | Evaluation of advanced control for li-ion battery balancing systems using convex optimization[END_REF]). Remarks: 1) We know no principle similar to 2) to justify the convexity of g in E b nor its joint convexity. Therefore, these hypotheses appear only as useful mathematical tricks which applicability needs to be proved. This applicability question would be the same if losses should depend on more variables. 2) Under principle 3) that losses should be minimized, the loss function g does not need to be an explicit expression. Instead it can be the solution of a convex

1 Two extra meaningful dependencies for losses are temperature and aging.

minimization subproblem which can be embedded in the global problem. See examples below. 3) Some precise sufficient conditions have been identified

to replace principle 3), but they are not necessary, cannot always be checked ex-ante and are application-specific (see [START_REF] Mégel | Hybrid stochasticdeterministic multiperiod dc optimal power flow[END_REF], [START_REF] Li | Extended sufficient conditions for exact relaxation of the complementarity constraints in storage-concerned economic dispatch[END_REF], [START_REF] Duan | Improved sufficient conditions for exact convex relaxation of storage-concerned[END_REF]). 4) In the literature, the relaxation we introduce in ( 6) is more commonly written by inserting the loss expression directly in the storage dynamics [START_REF] Boyd | Convex Optimization[END_REF] and relaxing the equality as E(k + 1) ≤ E(k) + (. . . ). It is equivalent and more compact, but we carve out the losses here to ease the discussion on their expression.

III. EXISTING CONVEX LOSS MODELS

There are two commonly used convex expressions that can play the role of g in [START_REF] Gill | Dynamic Optimal Power Flow for Active Distribution Networks[END_REF].

A. Piecewise-Linear-in-P Loss Model

The most common loss model is the PWL-in-P expression (e.g. used in [START_REF] Chen | Aggregate Power Flexibility in Unbalanced Distribution Systems[END_REF], [START_REF] Almassalkhi | Model-Predictive Cascade Mitigation in Electric Power Systems With Storage and Renewables-Part I: Theory and Implementation[END_REF]- [START_REF] Mégel | Hybrid stochasticdeterministic multiperiod dc optimal power flow[END_REF], but also in many other studies):

g(P b , E b ) = c + P + b + c -P - b (7)
where c + and c -are the charge and discharge loss coefficients respectively, while P + b and P - b are the positive and negative parts of P b (i.e. positive variables with P b = P + b -P - b ). This is very popular when the energy management problem is expressed with a linear program (LP). Indeed, as said in Remark 2 above, rather than using the explicit expression [START_REF] Wu | Optimal integration of a hybrid solar-battery power source into smart home nanogrid with plug-in electric vehicle[END_REF] for g, it is better formulated as a LP:

g = min P + b ,P - b c + P + b + c -P - b (8)
with constraints on the variables so that, at the optimum, those are indeed the positive and negative parts of P b :

P + b ≥ P b (9a) P - b ≥ -P b (9b) P + b , P - b ≥ 0 (9c) P + b , P - b ≤ P + rated , P - rated (9d)
Using Remark 4, this loss formulation is more commonly injected directly in the storage dynamics (1) which becomes:

E b (k + 1) ≤ E b (k) + (η + P + b + P - b /η -)∆ t (10)
Storage model ( 10) is obtained by replacing P b with P + b -P - b and introducing the charging and discharging efficiency coefficients η + = 1 -c + and η -= 1/(1 + c -). Expression [START_REF] Nick | Optimal allocation of dispersed energy storage systems in active distribution networks for energy balance and grid support[END_REF] is more compact than (1)+( 7), but the transformations performed here show how it belongs to our generic framework of convex losses functions. Coefficients η + , η -can be chosen equal to the square root of the round-trip efficiency of a given storage technology. However, this does not make (7) a physical representation of losses, because a round-trip efficiency is valid for a given cycle only, since instantaneous losses are rarely piecewise linear.

In the literature, expression ( 10) is more often written with equality rather than ≤. An equality constraint is allowed here because the expression is linear in (P + b , P - b ), but it is still a case of relaxed losses like [START_REF] Gill | Dynamic Optimal Power Flow for Active Distribution Networks[END_REF], because the variables P + b and P - b have no complementarity constraint. This means that they can be both positive, with no effect on the net flow P b , but which creates an artificial excess of losses.

Complementarity can be enforced by either a nonlinear complementarity constraint P + b .P - b = 0 [START_REF] Almassalkhi | Model-Predictive Cascade Mitigation in Electric Power Systems With Storage and Renewables-Part I: Theory and Implementation[END_REF] or by introducing a binary variable which creates a Mixed Integer LP (MILP) [START_REF] Arroyo | On the use of a convex model for bulk storage in mip-based power system operation and planning[END_REF]. We do not detail these possibilities since they are both not convex. Only, we highlight that it is the exactness of this LP relaxation that was most thoroughly studied (no comparable results for nonlinear loss models). First, sufficient mathematical conditions have been identified for the relaxation to be tight [START_REF] Li | Extended sufficient conditions for exact relaxation of the complementarity constraints in storage-concerned economic dispatch[END_REF], [START_REF] Duan | Improved sufficient conditions for exact convex relaxation of storage-concerned[END_REF] (though we mentioned in §II-D that their applicability is limited). Second, the usage of an equality sign in [START_REF] Nick | Optimal allocation of dispersed energy storage systems in active distribution networks for energy balance and grid support[END_REF] creates an upper bound for the worst-case amount of artificially wasted energy, when the relaxation fails. This worst-case bound (given in [5, App. A]) is better (lower) when the pair of constraints (9d) is replaced by the tighter limit

P + b /P + rated + P - b /P - rated ≤ 1 (11)
which can be either viewed as the tightest convex approximation of the nonlinear complementarity constraint or derived by relaxing the integrality constraint in the MILP formulation. See Fig. 2 (draft).

B. Quadratic-in-P Loss Model

The second quite popular convex loss model is quadratic in power (detailed in [START_REF] Murgovski | Component sizing of a plug-in hybrid electric powertrain via convex optimization[END_REF], with a close variant in [START_REF] Nick | Optimal allocation of dispersed energy storage systems in active distribution networks for energy balance and grid support[END_REF]):

g(P b , E b ) = ρ.P 2 b ( 12
)
with positive parameter ρ. Compared to the PWL-in-P model from section III-A, quadratic losses penalize more strongly large charge/discharge rates, while small rates are asymptotically lossless.

For optimization purpose, expression ( 12) is convex but nonlinear, so that a QP or a general Non-Linear Programming (NLP) solver must be used (or the quadratic function must be approximated by a set of a few tangents).

Physically, the loss model ( 12) mimics, inside an energybased model, the Joule losses on the series resistance of a circuit model (see Fig. 1). To understand this, we need to state the equivalence between the two. Conversion from circuit variables i, v to power is

P b = v.i = v 0 .i + r.i 2 , ( 13
)
where the last term is equal to the Joule losses. This expression can be solved for the current:

i = P b v 0 × 2 1 + 1 + 4rP b v 2 0 ( 14
)
which is defined for P b ≥ -

v 2 0
4r (i.e. the circuit model embeds a physical discharge power limit). Expression ( 14) is the product of a term linear in P b and a nonlinear correction factor which tends to one for small power (|P b | v 2 0 4r , which is the same as assuming a small voltage drop: |v -v 0 | v 0 ). As a consequence, when substituting this expression of i in the Joule losses, one gets an approximately quadratic expression

P loss ≈ r v 2 0 P 2 b ( 15
)
from which one can identify the physical definition of the quadratic loss coefficient ρ:

ρ = r/v 2 0 . ( 16
)
This analysis shows three conditions for the quadratic-in-P losses model to indeed match the Joule losses of a circuit model:

1) Charge/discharge rate must be small 2) Series resistance r must be constant 3) Open circuit voltage (OCV) v 0 must be constant The first condition can be lifted, if needed, because the loss expression without approximation is not quadratic but still convex (not shown here).

Lifting the "constant resistance & OCV" conditions is the motivation for the model introduced in the next section.

IV. THE "CONVEX MONOMIAL LOSS MODEL"

Based on the physical limitations of the quadratic loss model we search for a loss expression which can depend on the State of Energy and be separable like [START_REF] Duan | Improved sufficient conditions for exact convex relaxation of storage-concerned[END_REF]. Also, the dependency in P is assumed a power law P a b , so that our model embeds both PWL and quadratic-in-P loss models with a = 1 and a = 2 respectively. The model would be ρ(E b ) × P a b , but to handle both signs of P b when a = 2, we split the expression for the positive and negative parts of the power, like for the PWL-in-P model in III-A. Finally, to enable an easier convexity analysis, we impose that the dependency in the SoE should also be a power law, albeit with a negative exponent (for convexity). Here is our proposed "convex monomial loss model":

g(P b , E b ) = c + (P + b ) a+ |E b -e + | b+ + c - (P - b ) a- |E b -e -| b- (17)
with power exponents a + , a -≥ 1 (likely equal), energy exponents b + , b -≥ 0, scaling parameters c + , c -≥ 0 (like in [START_REF] Wu | Optimal integration of a hybrid solar-battery power source into smart home nanogrid with plug-in electric vehicle[END_REF], but not with the same physical dimension) and parameters e + , e -which are the energy value for which losses tend to 3. Production shifting optimization for different storage loss models. Scenario is 1 kW production while low electricity price during 1 st hour and no production while high electricity price during 2 nd hour. All lossy models are approximately calibrated to show a 80% round trip storage efficiency, except the last one which is a variant at 75%. With such efficiency, all models leverage the price difference by shifting the production to the 2 nd hour. However, the discharge profile varies with the loss model. Source: [START_REF] Haessig | pierre-haessig/convex-storage-loss v1.1[END_REF].

+∞, so we must have e + , e -< 0 or > E rated . In case the behavior in charge and discharge is symmetrical, expression [START_REF] Haessig | pierre-haessig/convex-storage-loss v1.1[END_REF] simplifies to:

g(P b , E b ) = c |E b -e| b ρ(E b ) ×|P b | a (18)
Example: For a capacitor C, we have a quadratic relation between energy and open circuit voltage:

E b = C(v 2 0 -v 2 ε )/2
, where v ε is the lowest useful OCV value for which the capacitor is defined as empty (useable energy E b = 0). With it, the Joule losses expression [START_REF] Duan | Improved sufficient conditions for exact convex relaxation of storage-concerned[END_REF] becomes

P loss ≈ rC/2 E b +Cv 2 ε /2 P 2
b . Thus we identify the following coefficients in model ( 18): a = 2, b = 1, c = rC/2, e = -Cv 2 ε /2. Notice that a direct derivation of this P 2 /E model was first proposed in [START_REF] Murgovski | Convex modeling of energy buffers in power control applications[END_REF].

A. Convexity Conditions for the Monomial Loss Model

We analyze the convexity of the first term of expression [START_REF] Haessig | pierre-haessig/convex-storage-loss v1.1[END_REF] since it is the sum of two similar expressions and since addition preserves convexity. We recognize that this term is the positive scaling of the 2D monomial f a,b : x, y → x a y b , with arguments being P + b and an affine expression of E b . Indeed the first term is c + .f a+,-b+ (P + b , ±(E b -e + )). The ± is "+" if e + < 0 and "-" if e + > E rated . This transformation preserves the convexity of f (see [1, §3.2.2]). We refer to the Appendix and Fig. 4 for the convexity analysis of f a,b . (noticing that f uses a flipped sign convention for exponent b to make it symmetric between a and b). The main result is that the loss model ( 17) is convex if:

b + ≤ a + -1, b -≤ a --1 (19)
As an illustration, for the quadratic-in-P case (a ± = 2), we must have b ± ∈ [0, 1]. We can thus observe that the P 2 /E loss model of a capacitor (b ± = 1) is on the "edge of convexity". Also, the monomial model shows that making the loss coefficients c ± in the PWL-in-P model [START_REF] Wu | Optimal integration of a hybrid solar-battery power source into smart home nanogrid with plug-in electric vehicle[END_REF] depend on SoE (i.e. introduce an SoE dependent efficiency) would be not convex (a = 1 implies b = 0).

V. ILLUSTRATION OF LOSS MODELS

We now introduce an elementary energy storage application to show the effect of varied loss models on the optimal charge/discharge profile. One storage unit is used as a buffer between a power production P prod and a grid injection P grid . Conservation of power is: P grid = P prod -P b . The optimization objective is to maximize the profit from selling the production to the grid at a time varying price:

max C grid = K k=1 c grid (k).P grid (k)∆ t (20)
The price will be positive so that the loss relaxation always works (see "positive price argument" in §II-D). The scenario is a production shifting on two hours (K = 20, ∆ t = 0.1 h), with a E rated = 1 kWh storage. For the first hour, P prod = 1 kW and price is low (c grid = 0.1 C/kWh). Then for the second hour, zero production and price is high (c grid = 0.2 C/kWh). The storage, described by (1) and varying loss models, can shift the production to the high price hour, if losses are small enough.

Optimal trajectories are shown on Fig. 3 for the lossless storage and four loss models: PWL-in-P b (7) (c ± = 0.111), linear-in-E b (5) (c E = 0.29 W/Wh), quadratic-in-P b (12) (ρ = 0.122 W/W²) and the symmetric monomial loss model [START_REF] Dunning | JuMP: A Modeling Language for Mathematical Optimization[END_REF], tuned for a capacitor (a = 2, b = 1, c = 0.0685, e = -0.25 kWh), that is losses are higher at low SoE for a given power. All models are tuned to yield an 80% round-trip efficiency on this scenario to make them comparable (except the last plot, see below). We checked that for all models the relaxed loss inequality ( 6) is tight, that is the relaxation works.

It can be observed that the charging phase is the same for all, with the storage absorbing all the production (P b = +1 kW). Then, despite the same round-trip efficiency, there is a great variety in discharge profiles, depending on the loss model: ∝ P2 b losses smooth out the discharge on the entire second hour to reduce power peaks, while the lossless and PWL models allow faster discharge 2 since efficiency is constant. The ∝ E b losses push for the earliest discharge. Finally, the ρ(E b )P 2 b model shows a decaying discharge rate, to prevent loss increase at low SoE. That is P loss stays quite constant despite the rising loss coefficient ρ(E b ).

The last plot of Fig. 3 is a variant of this model, with slightly higher losses (c = 0.094, round-trip efficiency is 75%). This shows that the early charging is done at a reduced rate, with some of the production sold directly at a low price, because of the high initial value of ρ. Such time-varying partial arbitrage during a constant price period can only be obtained with an SoE-dependent loss model.

Optimization results were obtained by describing optimization problems in Julia using the JuMP package [START_REF] Dunning | JuMP: A Modeling Language for Mathematical Optimization[END_REF]. Solvers were Ipopt [START_REF] Wächter | On the implementation of an interiorpoint filter line-search algorithm for large-scale nonlinear programming[END_REF] and ECOS [START_REF] Domahidi | ECOS: An SOCP solver for embedded systems[END_REF]. The complete code is available in a Jupyter notebook [START_REF] Haessig | pierre-haessig/convex-storage-loss v1.1[END_REF].

VI. CONCLUSION

The relaxation of storage losses enable more realistic loss models while preserving convexity. We described this relaxation in a generic way, unifying the existing works which are often split by model types (e.g. piecewise linear vs. others), and we proposed one flexible model which contains all existing loss expressions.

Future work includes studying the case when the relaxation fails, e.g. generalizing the worst case bound which exists for the PWL-in-P model. Also, including other variables influencing losses (temperature, aging), would be useful. Finally, fitting our model to actual battery loss data would confirm its usefulness. + , depending on exponents a and b. Red area 0 < a < 1 is excluded because f is concave in the x-direction (same for b). Blue regions are excluded because the determinant of the Hessian is negative. The three convex areas are in green. Source: [START_REF] Haessig | pierre-haessig/convex-storage-loss v1.1[END_REF]. defined for positive arguments x, y and real exponents a, b.

The convexity of f is equivalent to the positive definiteness of its Hessian:

H(x, y) = a(a -1)x y b abx a-1 y b-1 abx a-1 b-1 b(b -1)x a y b-2 (22)
H positive definite if and only if:

• each diagonal term is positive (i.e. f is convex the x and y directions) • the determinant is positive The areas these two conditions are not fulfilled are highlighted in red and blue respectively on the (a, b) plane in Fig. 4. This leaves three convex regions in green. In this analysis, we assume that variable x will be mapped to the storage power P b (or -P b ). Since losses are generally increasing with the charging/discharging power, only the region with positive exponent a makes physical sense. Thus the single region with interesting exponent values, highlighted on Fig. 4, is:

{(a, b) such that a + b ≥ 1, b ≤ 0} (23)

Fig. 1 .

 1 Fig.1. Energy storage models. On the left: energy based model which is ubiquitous in energy management optimization. On the right: circuit based model which is more physical for electrochemical devices (v 0 and r are often state dependent). Conversion between the two is given Section III-B.

Fig. 2 .

 2 Fig. 2. Complementarity constraint (bold red lines) + convex relaxation (dark blue triangle), vs box constraint (light blue square)

5 ( 5 (

 55 Fig.3. Production shifting optimization for different storage loss models. Scenario is 1 kW production while low electricity price during 1 st hour and no production while high electricity price during 2 nd hour. All lossy models are approximately calibrated to show a 80% round trip storage efficiency, except the last one which is a variant at 75%. With such efficiency, all models leverage the price difference by shifting the production to the 2 nd hour. However, the discharge profile varies with the loss model. Source:[START_REF] Haessig | pierre-haessig/convex-storage-loss v1.1[END_REF].

Fig. 4 .

 4 Fig. 4. Convexity of the losses helper function f : x, y → x a y b on R 2+ , depending on exponents a and b. Red area 0 < a < 1 is excluded because f is concave in the x-direction (same for b). Blue regions are excluded because the determinant of the Hessian is negative. The three convex areas are in green. Source:[START_REF] Haessig | pierre-haessig/convex-storage-loss v1.1[END_REF].

There is no unique solution to the optimization problem[START_REF] Domahidi | ECOS: An SOCP solver for embedded systems[END_REF] with zero or PWL-in-P losses, so we add a tiny increasing penalty coefficient to P grid to select the "earliest injection" solution among all optimal trajectories. This avoids getting a random one. A later discharge would be equally optimal.

ACKNOWLEDGMENT

We thank Olivier Mégel for a discussion back at PowerTech 2015 which planted the seed of our reflection on the convexity of losses. We thank Bertrand Cornélusse for providing us with useful references on relaxation in optimal power flow.

APPENDIX CONVEXITY OF THE LOSS HELPER FUNCTION

We define the loss helper function as the 2D monomial