Ioannis Dimitriou

Sara Alouf
email: sara.alouf@inria.fr

Alain Jean-Marie
email: alain.jean-marie@inria.fr

A Markovian queueing system for modeling a smart green base station

Keywords: Coupled queues, QBD processes, Green base station

We investigate a model to assess the performance of a base station (BS) fully powered by renewable energy sources. The BS is modeled as a three-queue system where two of them are coupled. One represents accumulated energy, the second is the data queue and the third one serves as a reserve energy queue. This smart BS is able to dynamically adjust its coverage area (thereby controlling the traffic intensity) and to generate signals to the reserve energy queue that trigger the movement of energy units to the main energy buffer. Given the randomness of renewable energy supply and the internal traffic intensity control, our queueing model is operated in a finite state random environment. Using the matrix analytic formalism we construct a five-dimensional Markovian model to study the performance of the BS. The stationary distribution of the system state is obtained and key performance metrics are calculated. A small numerical example illustrates the model and a simplified product-form approximation is proposed.

Introduction

The architectural design of cellular networks has evolved in recent years to better satisfy users needs. The traffic load generated by users exhibits a night-day pattern, having a peak of traffic during the day and almost no traffic during the night. A geographical pattern is also observed as offices areas witness a peak traffic during the day while residential areas witness a smaller peak late in the evening. Heterogeneous cellular networks are an attractive deployment solution: large powerful base stations (BSs) are used to ensure coverage and connectivity whereas smaller coverage-limited BSs are used to accommodate the peak load where needed. Those BSs can be analyzed in isolation.

In this paper, we consider a single small BS and study the question of powering it using solely renewable energy. A key factor is that the solar radiation exhibits a night-day pattern that makes solar panels fit to power a small BS. We assume that the BS is "smart" in the sense that it is able to dynamically adjust its coverage area, controlling thereby the number of mobiles with which it communicates, and consequently its offered traffic rate and its energy consumption. The harvested energy is stored in batteries that are used to power the BS.

We propose to model this BS as a queueing system that operates in a finitestate Markovian random environment (RE). The behavior of such a system is described by a five-dimensional Markov process, which is a homogeneous finite Quasi Birth-Death (QBD) process.

The literature on QBD processes is abundant and these processes have been used to study many different applications. De Cuypere et al. [START_REF] De Cuypere | Performance analysis of a kitting process as a paired queue[END_REF], have studied sensor nodes [START_REF] De Cuypere | Stochastic modelling of energy harvesting for low power sensor nodes[END_REF] and kitting processes through a paired queueing system [START_REF] Latouche | Queues with paired customers[END_REF]. In [START_REF] De Cuypere | Stochastic modelling of energy harvesting for low power sensor nodes[END_REF], a finite energy queue is paired with an infinite data queue, where customers arrive at both queues according to Markovian Arrival Processes. In [START_REF] De Cuypere | Performance analysis of a kitting process as a paired queue[END_REF], both queues are finite, but due to the sparsity of the generator matrix of the underlying Markov process, the size of the state space does not cause serious issues. Closely to our work, Takahashi et al. [START_REF] Takahashi | On a synchronization queue with two finite buffers[END_REF] consider a synchronization queue consisting of two buffers with finite capacities, the arrival processes at both buffers are Poisson for one and phase type renewal for the other. In our work we assume that energy is discretized (as in e.g. [START_REF] De Cuypere | Stochastic modelling of energy harvesting for low power sensor nodes[END_REF][START_REF] Gelenbe | A sensor node with energy harvesting[END_REF][START_REF] Gelenbe | Interconnected wireless sensors with energy harvesting[END_REF]) but it is also possible to model batteries as fluid queues; see for instance [START_REF] Jones | Fluid queue models of renewable energy storage[END_REF][START_REF] Jones | Fluid queue models of battery life[END_REF].

In the following, we present our model in Section 2 and detail the infinitesimal generator of the five-dimensional continuous-time Markov chain representing the state of the system in Section 3. Section 4 discusses several algorithms that could be used to compute the stationary distribution of the system state and relevant performance metrics, exploiting the QBD structure of the generator. The feasibility of the developed algorithm is demonstrated through a numerical example. First passage times to lower and higher levels in the QBD are discussed in Section 5. In Section 6, we propose an approximate model with product-form solution. We briefly conclude in Section 7.

The model

We are interested in a wireless communications base station (BS), isolated from the electric grid and operating thanks to renewable energy sources. The model supposes a continuous functioning of the station: it does not handles issues like startup, shutdown, malfunction and other transitory phenomena.

The model we develop has a 5-dimensional state space, representing the state of three queues, a service process and an environment process. We describe now these elements and the way they interact. Figure 1 summarizes the model under study. Although we describe here a quite specific situation, we point out that all these elements and their composition can be generalized so as to build models of other device configurations.

Queues. We consider two energy queues (EQ) of finite capacity E j , j = 1, 2 that store energy extracted from the environment, and a data queue (DQ) that keeps track of packets not yet transmitted, having also a buffer of finite capacity N . Energy is assumed to be discretized, as for instance in [START_REF] Gelenbe | A sensor node with energy harvesting[END_REF][START_REF] Gelenbe | Interconnected wireless sensors with energy harvesting[END_REF]. Assume that EQ 1, is coupled with the DQ, whereas EQ 2 is deployed close to the BS and serves as a reserve. Each EQ, is fed by different renewable energy sources, and the choice of the one that will be coupled with the DQ depends on the area where the BS

1 E N 2 E Energy packets 3 h 2 h 1 h ) ( 1 i e  ) ( , i s p  ) ( 2 i e  i PH Battery leakage 1 i u 2 i u Battery leakage ) (i s   Arrivals-Departures:
Leakage: Signals:

Fig. 1. The model when the RE is in state i.

is deployed. We will refer to this EQ, say EQ 1, as the dominant energy queue (DEQ) and the EQ 2 as the secondary energy queue (SEQ).

Environment. The system behavior depends on the state of the RE. This RE is defined by means of an irreducible regular continuous-time Markov chain Y (t) with state space {1, ..., M }, and infinitesimal generator Q Y . The state of this environment may represent a variety of factors that make energy arrival or packet arrival processes, energy consumption and service times non time-stationary. For instance, the intermittent nature of wind, the variability of sunlight, can be taken into account with an environment variable. Hourly variations of data traffic can also be modeled that way. Other features can be the variation of transmission power due to global adjustments at the network level. As usual, independent environment features can be combined in a complex environment process with multiple "phases", at the cost of large values for M . Clearly, one can skip the RE and reduce the dimension of the QBD, by assuming that the model migrates quickly from the stationary situation of one environment setting to another.

Packet service duration. The service time distribution is of phase type [START_REF] Neuts | Matrix-geometric solutions in stochastic models: An Algorithmic Approach[END_REF][START_REF] Latouche | Introduction to Matrix Analytic Methods in Stochastic Modeling[END_REF] of order ν and depends on the state of the RE. For a fixed RE state i, this distribution can be interpreted as a time until some underlying Markov process η t with finite state space {1, ..., ν} reaches the absorbing state 0 with initial probability vector (τ

(i) 0 , τ (i) )
. Transition rates of the process η t within the set {1, ..., ν} are defined by the sub-generator T (i) and transition rates into the absorbing state (which lead to a service completion) are given by the entities of the column vector t (i) = -T (i) 1. Observe that the value ν does not depend on i. A typical situation allowed by this model is to have a packet size with some given distribution, and a service speed given by the environment: T (i) kl = T kl × v i , v i begin the "velocity" typical of state i and T kl some environment-independent transition rate. If the packet length distribution itself can also depend on the environment, the model may not be accurate for the packet in service during the change of environment state.

By using a generalized PH distribution for the service process, we can incorporate several realistic concepts of the operation of a BS such as: different levels of noise in the transmission channel, hardware degradation and recovering, variations in the distance of a mobile user to the base station, etc. Packet arrival. The rate of packet arrival to the BS depends on its coverage area. We adopt a multi-threshold scheme in order for the BS to dynamically adjust the coverage area according to the available energy units in the DEQ: in this sense, we model a smart BS. More precisely, we introduce thresholds, say

h 0 = 0 < h 1 < h 2 < ... < h K < E 1 = h K+1 .
Given the state i of RE, and if h s < m 1 ≤ h s+1 , s = 0, ..., K, the users' arrival rate equals λ

(i) p,s . In practice, it is expected that λ (i) p,0 ≤ λ (i) p,1 ≤ . . . ≤ λ (i)
p,K but this is not needed for the definition of the model. Note that the packet traffic potentially includes system control packets, signaling (e.g. for handovers) etc. Superposition of different traffic sources can classically be taken care of. Energy arrival and depletion. For a fixed state i, i = 1, ..., M of the RE, energy units are stored in EQ j at rate λ (i) ej . Transmission opportunities occurs only when both the DEQ and the DQ queue are non empty.

The modeling of energy consumption by a packet transmission requires some trick. Indeed, it is unlikely that a phase of the service of some packet will consume exactly an integer amount of energy quanta. Randomness in consumption and/or non-integral energy values can be taken into account using random variables.

The number of energy units required for the transmission of a single packet depends both on the state of the RE and on the phase of the service process. Thus, the transmission of a packet may be interrupted, if the completion of a phase of service consumes the available energy. In this paper, we assume that the transmission of this packet has to restart from scratch (but with a stochastically independent value), whenever there will be available energy. If some mechanism, e.g. error-correcting codes, allows for it, transmission could resume from the same phase. Alternately, the transmission of the packet could be simply canceled in case of energy shortage. Moreover, we can also allow this cancellation to occur only after some timeout, by introducing a timer for each "impatient" packet when energy is 0. None of these alternatives would make the transition matrix much more complex. With the "restart" point of view, given the state i of RE and if there are m 1 available energy units in the DEQ, the completion of phase x, x = 1, 2, ..., ν requires k energy units with probability p

(i,m1) kx , m1 k=0 p (i,m1) kx = 1.
We allow a priori k = 0, i.e. that the completion of a phase of service may not consume even a single energy unit. However, in order to initiate a transmission we need at least one energy unit.

Energy leakage from each EQ is unavoidable. More precisely, an energy unit will be lost from EQ j, j = 1, 2 at exponential rate u ij , given that the RE is in state i. Note that if the DEQ drains during the transmission of a packet, the packet has to be retransmitted, whenever there will be available energy units.

Energy transfers. The SEQ serves as backup storage for the DEQ. The control of our smart BS generates signals to the SEQ at a rate Λ

-(i) s
, that trigger the movement of energy units to the DEQ, according to the multi-threshold rule introduced above. Given the state i of the RE, and h s < m 1 ≤ h s+1 , s = 0, ..., K, signals are generated at a rate Λ -(i) s , and trigger [START_REF] Gelenbe | G-networks with triggered customer movement[END_REF] k energy units from the SEQ towards the DEQ with probability q (m1,m2,i) ks ( m2 k=1 q (m1,m2,i) ks δ m1+k≤E1 = 1). Define also q (m1,m2) ks = diag(q (m1,m2,1) ks , . . . , q (m1,m2,M ) ks

). In case the SEQ is empty, the signal will have no effect. As above, it is expected in practice that Λ

-(i) 0 ≥ Λ -(i) 1 ≥ . . . ≥ Λ -(i)
K : the demand rate for replenishment of DEQ becomes larger as the energy in the DEQ depletes. It is also possible to reflect the urgency of replenishment by modifying the distribution q ks . Moreover, when the DEQ is full, i.e. m 1 = E 1 , there will be no signal towards the SEQ (i.e., Λ -(i) K = 0). We also adopt an overflow operation in the sense that when an EQ buffer is full, its energy traffic is rerouted to the other EQ. Clearly, the signal generation possibly consumes an amount of energy, which is assumed to be negligible in our model.

Thanks to the generality of the definition of the probabilities p

(i,m1) kx , q (m1,m2,i) ks
(which are function of the state of the RE, the phase of service process and the available energy in DEQ, SEQ), the modeler can incorporate additional realistic features regarding the energy consumption, related to the control channel data transmission, or processing and forwarding of packets. More complex packet processing architectures may however require more than one packet queue. Clearly, recent technological developments are towards smart autonomic wireless networks. The concept of signaling towards the SEQ, that triggers the instantaneous transition of energy units to the DEQ, has become an intelligent modeling tool for communication systems. Queues with signals [START_REF] Gelenbe | G-networks with triggered customer movement[END_REF] were introduced to model the behavior of control actions such as the displacement of units from one queue to another using "triggers" resulting in load balancing.

Process of the system state

The behavior of the system under consideration can be described in terms of the CTMC

X t = {Q p (t), J(t), Q e1 (t), Q e2 (t), Y (t)}, t ≥ 0, where Q p (t), J(t), Q ej (t)
and Y (t) are, respectively, the number of data packets, the phase of the service process (present only when Q p > 0 and Q e1 > 0), the number of energy units in EQ j and the state of RE at time t. The state space of X t is H = ∪ N n=0 l(n) and the "levels" l(n) are defined as:

l(0) = {(0, m 1 , m 2 , i); m j = 0, 1, ..., E j , j = 1, 2, i = 1, ..., M }, (1)
and for 1 ≤ n ≤ N :

l(n) = l(n, 0) ∪ l(n), l(n, 0) = {(n, 0, m 2 , i); m 2 = 0, 1, ..., E 2 , i = 1, ..., M }, l(n) = {(n, x, m 1 , m 2 , i); x = 1, . . . , ν; m 1 = 1, . . . , E 1 ; m 2 = 0, 1, ..., E 2 ; i = 1, ..., M }. (2)
These levels have cardinals |l(0

)| = L 0 := M (E 1 + 1)(E 2 + 1) and for n ≥ 1, |l(n)| = L := M (E 2 + 1)(νE 1 + 1)
. For convenience, define also S := M (E 2 + 1).

The state space has then cardinal

| H| = M (E 2 + 1)[N (νE 1 + 1) + (E 1 + 1)].
Using this decomposition in levels, the infinitesimal generator of X t has a quasi-birth-death structure with block matrix representation as:

Q =            B 0 C 0 0... 0 0 A 10 A 1 C 0... 0 0 0 A 21 A 1 C... 0 0 . . . . . . . . . 0 0 0 ...A 21 A 1 C 0 0 0 ...0 A 21 A 2            .
We proceed with a detailed description of the blocks. We shall need the definition of the following matrices Λ ej = diag(λ

(i) ej ), j = 1, 2, Λ p,s = diag(λ (i) p,s ), Λ - s = diag(Λ -(i) s
), s = 0, ..., K, U j = diag(u ij ), where i = 1, ..., M , and

t x = diag(t (1)
x , ..., t

(M ) x ), p (m1) lx = diag(p (1,m1) lx , ..., p (M,m1) lx
). The edge block A 10 is an L × L 0 matrix that corresponds to packet service and energy consumption events that lead to an empty data queue. It has a blocktriangular structure with A (0,m 1 ) 10 = 0 S×S , for all m 1 (this first row corresponds to the case where there is a packet waiting for service, but the transmission cannot be initiated since there is no available energy in the DEQ), A (m1,m 1 ) 10

= 0 νS×S for 1 ≤ m 1 < m 1 ≤ E 1 and for m 1 = 1, ..., E 1 , m 1 = 0, ..., m 1 : A (m1,m 1 ) 10 = A (m1-m 1 ) m1m 1 = (I (E2+1)×(E2+1) ⊗ [t (m1) 1 p (m1) m1-m 1 ,1 ], ..., I (E2+1)×(E2+1) ⊗ [t (m1) ν p (m1) m1-m 1 ,ν ]) . (3)
The sub-diagonal block A 21 is an L × L matrix that corresponds to packet service and energy consumption:

A 21 =              0 S×S 0 S×S 0 S×S ... 0 S×S 0 S×S A (1) 10 F 11 0 νS×S ... 0 νS×S 0 νS×S A (2) 20 F 21 F 22 ... 0 νS×S 0 νS×S ... ... ... . . . ... ... ... . . . A (E1) E10 F E11 F E12 ... F E1E1-1 F E1E1              , where F m1m 1 is a νS × νS matrix defined from τ x = diag(τ (1)
x , ..., τ

(M ) x
) as:

F m1m 1 = I (E2+1)×(E2+1) ⊗ (t x τ y p (m1) m1-m 1 x ), x, y = 1, ..., ν. (4)
Packet arrivals are represented in block C = diag(C 0 , C 1 , ..., C K ), which is an L 0 × L matrix where C s is of order n(s)S × n(s)Sν for s = 1, 2, ..., K, and of order n(0)S × [(n(0) -1)ν + 1]S for s = 0, with n(s) = h s+1 -h s + 1δ s=0 and

C 0 = diag(C 00 , G 0 , ..., G 0 h1 ), C s = I n(s)×n(s) ⊗ G s , s = 1, ..., K, G s = (G (1) s , ..., G (ν) s ), G (x) s = I (E2+1)×(E2+1) ⊗ (Λ p,s τ x ). (5) Matrix C = diag(C 0 , C 1 , ..., C K ) is of order L × L, where C 0 = diag(C 00 , G 0 , ..., G 0 h1 ), C s = I n(s)×n(s) ⊗ G s , s = 1, ..., K, G s = I ν×ν ⊗ C 0s , C 0s = I (E2+1)×(E2+1) ⊗ Λ p,s , s = 0, ..., K. (6) Next, B 0 = (B (0) m1,m 1 ), 0 ≤ m 1 , m 1 ≤ E 1 , is an L 0 × L 0 matrix formed of square sub-blocks B (0) m1,m 1
of order S, typical of energy movements in the DEQ that are independent of packet transmissions, such as energy leakage, energy generation and signal-triggered energy unit transition from SEQ to DEQ. Its structure is:

B 0 =          B (0) 00 B (0) 01 B (0) 02 ... B (0) 0E1 U 1 B (0) 11 B (0) 12 ... B (0) 1E1 ... ... ... ... 0 L0×L0 0 L0×L0 0 L0×L0 ... U 1 B (0) E1-1,E1-1 B (0) E1-1,E1 0 L0×L0 0 L0×L0 0 L0×L0 ...0 L0×L0 U 1 B (0) E1E1         
, where U 1 = I (E2+1)×(E2+1) ⊗ U 1 represents energy leakage from DEQ, and for 0

≤ m 1 < E 1 , h s < m 1 ≤ h s+1 , B (0)
m1,m1+1 =         Λ e1 ... Λ - s q (m1,1) 1s Λ e1 ... ... ... ... ... ...Λ - s q (m1,E2-1) 1s Λ e1 ... Λ - s q (m1,E2) 1s Λ e1 + Λ e2        
.

The remaining sub-blocks concern the energy replenishment due to signals.

Assuming E 1 ≤ E 2 , for k = m 1 + 2, ..., E 1 , h s < m 1 ≤ h s+1 , B (0) m1k =             . . . . . . . . . . . . . . . Λ - s q (m1,k-m1) k-m1s . . . Λ - s q (m1,k-m1+1) k-m1s . . . . . . . . . . . . . . . . . . Λ - s q (m1,E2) k-m1s . . .            
.

Remark: In this work we assume without loss of generality that

E 1 ≤ E 2 . For E 1 > E 2 , k = m 1 +2, ..., m 1 +E 2 , B (0)
m1k has also the same form but B (0)

m1k = 0 S×S , for k = m 1 + E 2 + 1, ..., E 1 . Furthermore, for 0 ≤ m 1 ≤ E 1 and h s < m 1 ≤ h s+1 , set F s = Q Y - 2 j=1 (Λ ej δ mj <Ej + U j δ mj >0 ) -Λ p,s -Λ - s δ m1<E1 . Then for m 2 , m 2 = 0, ..., E 2 , B (0) m1m1 = (B (0;m2,m 2 ) m1m1
) where, B (0;m2,m2-1) m1m1 = U 2 and,

B (0;m2,m2+1) m1m1 = Λ e2 + Λ e1 δ m1 =E1 , B (0;m2,m2) m1m1 = F s , h s < m 1 ≤ h s+1 . (7)
Finally, the main diagonal blocks A 1 and A 2 are square matrices of order L, also composed of sub-blocks

A l = (A (l) m1m 1 ), m 1 , m 1 = 0, 1, ..., E 1 , for l = 1, 2.
Here, A (l) 00 is of order S × S, A

(l) 0m 1 , m 1 = 1, ..., E 1 are of order S × νS, A (l) m10 , m 1 = 1, ..., E 1 , of order νS × S, and
A (l) m1m 1 , m 1 , m 1 = 1, ..., E 1 of order νS × νS.
They are given by, A ) where

Z (0,m1) x = (Z (0,m1) x;m2,m 2
), x = 1, ..., ν, are S × S matrices where, Z (0,1)

x;m2,m2 = (Λ e1 + Λ e2 δ m2=E2 )τ x , m 2 = 0, 1, ..., E 2 , Z (0,1) x;m2,m2-1 = Λ - 1 q (0,m2) 11 τ x , Z (0,1) x;m2,m 2 = 0 M ×M , elsewhere. For m 1 = 2, ..., E 1 , h s < k ≤ h s+1 , Z (0,m1) x;m2,m 2 = Λ - s q (0,m2) m1s τ x , m 2 = m 1 , ..., E 2 , m 2 = 0, 1, ..., m 2 -m 1 , (8) and Z
(0,m1) x;m2,m 2 = 0 M ×M elsewhere. Furthermore, for l = 1, 2, m 1 = 1, ..., E 1 , A (l)
Φ (l,m1,m1) xy;m2,m 2 = I (E2+1)×(E2+1) ⊗ [T xy p (m1) 0x ], m 1 = 0, 1, ..., E 1 , x = y, Φ (l,m1,m1) xx;m2,m 2 = I (E2+1)×(E2+1) ⊗ (T xx + Λ p,s δ m1<E1 ) + B (0) m1m1 + W s δ m1<E1 , W s = diag(0 M ×M , I (M E2×M E2) ⊗ Λ - s ). ( 10
)
Sparsity of the generator. The number of non-zero transition rates in the generator Q can be roughly evaluated as follows. Neglecting the "border" cases, transitions from a typical state (n, x, m 1 , m 2 , i) occur as: a) packet arrival, energy arrival, energy leakage: one transition each; b) service phase change: ν -1 possible transitions with service continuation, plus service completion, which triggers energy consumption with up to m 1 possible transitions; c) environment phase change: M -1 possible transitions (this could be much less if the phase process is itself a superposition, see Section 2); d) energy transfers: up to m 2 transitions. Summing over m 1 and m 2 , we conclude that storing the elementary blocks of

Q requires O (ν + M )M E 2 1 E 2 2 space. Storing all Q in sparse form would require O N (ν + M )M E 2 1 E 2 2 space.
The ratio of this quantity with the total number of entries in

Q: O (N M νE 1 E 2 ) 2 , is of order α = O((ν + M )/N νM ). It may even be smaller if the RE transition matrix Q Y is itself sparse, or if the distri- butions p (i,m1) kx
or q (m1,m2,i) ks have restricted supports. Note that some algorithms work with storing only the different blocks.

Stationary solution

Algorithms for stationary probabilities

Under reasonable assumptions on packet arrival rates and the representation of service distributions, the chain X t is irreducible. Since it is finite, it then admits a unique stationary distribution given by (referring to the notation in ( 1) and ( 2)):

p(n, x, m 1 , m 2 , i) = lim t→∞ P (Q p (t) = n, J(t) = x, Q e1 (t) = m 1 , Q e2 (t) = m 2 , Y (t) = i) p(n, m 1 , m 2 , i) = lim t→∞ P (Q p (t) = n, Q e1 (t) = m 1 , Q e2 (t) = m 2 , Y (t) = i).
In order to compute this stationary probability vector, several methods are available. On the one hand, the QBD structure of Q offers the possibility to use direct numerical methods, such as the Matrix Analytic ones. Those are "exact" (up to round off errors) and in general they are efficient especially for systems where their generator matrix has a special form, which is the case here. On the other hand, the sparsity of Q points at the use of iterative and other advanced techniques [START_REF] Wimmer | Symbolic partition refinement with automatic balancing of time and space[END_REF].

One way or the other, a closer investigation of the block structure of Q is fruitful. The stationary distribution vector p is represented using the level structure as p = (p 0 , p 1 , ..., p N ). Global balance equations pQ = 0 can then be written as

p 0 B 0 + p 1 A 10 = 0, ( 11
)
p 0 C + p 1 A 1 + p 2 A 21 = 0, ( 12
)
p i-1 C + p i A 1 + p i+1 A 21 = 0, i = 2, ..., N -1, ( 13
)
p N -1 C + p N A 2 = 0. ( 14
)
The matrix analytic formalism gives the opportunity to obtain the stationary probability vector p using several algorithms (see [START_REF] Latouche | Introduction to Matrix Analytic Methods in Stochastic Modeling[END_REF] and references within, and later contributions among which [START_REF] Akar | A novel computational method for solving finite QBD processes[END_REF][START_REF] Elhafsi | On the solution to QBD processes with finite state space[END_REF][START_REF] Ye | Folding algorithm: A computational method for finite QBD processes with level-dependent transitions[END_REF]). The methods of [START_REF] Gaver | Finite birth-and-death models in randomly changing environments[END_REF][START_REF] Ye | Folding algorithm: A computational method for finite QBD processes with level-dependent transitions[END_REF] are highlighted in [START_REF] Latouche | Introduction to Matrix Analytic Methods in Stochastic Modeling[END_REF]Chap. 10]. The key step in [START_REF] Gaver | Finite birth-and-death models in randomly changing environments[END_REF] resides especially in the inversion of a matrix of order L 0 and N matrices of order L. Given that L 0 ≤ L, the overall complexity is O(N L 3 ). The "folding method" in [START_REF] Ye | Folding algorithm: A computational method for finite QBD processes with level-dependent transitions[END_REF] is less expensive, reducing the "N " to log 2 N , but at a higher implementation cost. Moreover, [START_REF] Gaver | Finite birth-and-death models in randomly changing environments[END_REF] also provides an algorithm for the computation of the expected first passage times between the neighbour levels. This result is interesting for further capacity planning investigation, see Section 5 and our conclusion. In our situation, the limiting factor is clearly the value of

L ≈ M νE 1 E 2 .
The form of the block generator matrix will help us to improve the efficiency in obtaining the stationary probabilities, since the matrix C is a diagonal nonsingular matrix. Solving recursively system ( 11)-( 14) starting from [START_REF] Latouche | Introduction to Matrix Analytic Methods in Stochastic Modeling[END_REF] we get that, p n = p N R n for n = 0, ..., N , with:

R N = I, R N -1 = -A 2 C -1 , R 0 = (R 2 A 1 + R 3 A 21 )A 10 B -1
0 , and

R n = -(R n+1 A 1 + R n+2 A 21 )C -1 , n = 1, ..., N -2. ( 15
)
To obtain p N , we solve equation ( 12) along with the normalizing condition:

p N (R 0 C + R 1 A 1 + R 2 A 21 ) = 0, p N (R 0 1 0 + N n=1 R n 1) = 1, (16)
where 1, 1 0 are column vectors of order L and L 0 respectively. The computation of the R n 's involves only matrix additions and multiplications. Since C is diagonal and non singular, multiplying by its inverse is not costly. The computation of the p n 's [START_REF] Elhafsi | On the solution to QBD processes with finite state space[END_REF] requires two matrix inversions (one to compute B -1 0 , one to solve the system ( 16)), 4N + 5 matrix multiplications and 2N + 1 matrix additions. The resulting complexity is O(L 3 ). Using an iterative method to solve [START_REF] Neuts | Matrix-geometric solutions in stochastic models: An Algorithmic Approach[END_REF] and exploiting the inner structure of blocks (see Section 2) may help overcome difficulties about the computational complexity.

Performance metrics

By calculating the stationary distribution vector of the underlying Markov process we can obtain some important performance metrics, such as the depletion probability (DP ) (i.e., the probability of an empty DEQ), and the expected number of data packets and energy packets in each EQ:

DP = N n=0 E2 m2=0 M i=1 p(n, 0, m 2 , i), E(Q p ) = N n=1 E2 m2=0 M i=1 n [p(n, 0, m 2 , i) + E1 m1=1 ν x=1 p(n, x, m 1 , m 2 , i) , E(Q e1 ) = E1 m1=1 E2 m2=0 M i=1 m 1 [p(0, m 1 , m 2 , i) , + N n=1 ν x=1 p(n, x, m 1 , m 2 , i) , (17)
E(Q e2 ) = E2 m2=1 M i=1 m 2 E1 m1=0 p(0, m 1 , m 2 , i) + N n=1 p(n, 0, m 2 , i) + E1 m1=1 ν x=1 p(n, x, m 1 , m 2 , i) . ( 18
)
Using Little's law we can obtain the mean waiting time for a data packet to be served. Moreover, various optimization problems can be formulated, such as finding optimal values of λ (i) p,s that minimizes E(Q p ), asking DP ≤ a, or finding an optimal value for E 1 which minimizes DP = DP (E 1 ).

Numerical example

To demonstrate the feasibility of the developed algorithm and numerically show some features of the model, we present the results of a preliminary numerical experiment, inspired from values taken from [START_REF] Kim | Analysis of an M M AP/P H1, P H2/N/∞ queueing system operating in a random environment[END_REF]. Assume that the RE is defined by the 2-state infinitesimal generator with rates (Q Y ) 12 = 0.01 and (Q Y ) 21 = 0.1. We may interpret that the first state of the RE corresponds to a mode when the system is overloaded by the data units (peak time) and the second state to a normal mode of the system. Under the first and second state of the RE the service time distribution is characterized by the vectors τ (1) = (0.2, 0.8), τ (2) = (0.7, 0.3) and the matrices

T (1) = -0.4368 0.40768 0.42608 -1.71809 , T (2)
= -0.6402 0.56812 0.22308 -1.34325 , respectively. We assumed for m 1 + m 2 ≤ E 1 , q (m1,m2,i) ks

= 2 ki / m2 j=1 2 ji , k = 1, ..., m 2 and if m 1 + m 2 > E 1 , q (m1,m2,i) ks = 2 ki / E1-m1 j=1 2 ji , k = 1, ..., E 1 -m 1 .
Furthermore, p (i,m1) kx = ( m1 j=1 (ix) k-j ) -1 , k = 1, ..., m 1 . Table 1 contains the values of the parameters of the system.

E1 = 10, E2 = 12 M = 2, ν = 2 N = 15, h1 = 8, h2 = 12 Λp,1 = diag(1.5, 1) Λp,2 = diag(2, 1.5) U1 = diag(0.1, 0.3) U2 = diag(0.3, 0.4) Λ - 1 = diag(0.3, 0.2) Λ - 2 = diag(0.2, 0.1) Le 1 = diag(1.5, 1) Le 2 = diag(1.8, 1.2)
We focus on the depletion probability. Figure 2 (left) describes the way the depletion probability is affected by the RE for increasing values of the data arrival rate when m 1 ∈ [0, 8] (s = 0). Clearly, the depletion probability will increase, but this increase will become more apparent during the overloaded period, as expected. In Figure 2 (right) we can observe the benefits of the smart operation of the base station due to the presence of signals. By increasing the signal generation rate, the depletion probability is strongly reduced. As a result, the quality of service is thoroughly increased. Gaver et al. [START_REF] Gaver | Finite birth-and-death models in randomly changing environments[END_REF] analyzed the first passage times for general finite QBD processes. They obtained two systems of recurrence equations for the Laplace-Stieltjes transforms of passage times to higher and lower levels.

In the following we turn our attention to the computation of the average first passage times. Let the expected values of the first passage time from the states of l(n -1) to the states of l(n) be c (n) and that from the states of l(n + 1) to l(n) be v (n) . Define the following matrices:

H 0 = B 0 H 1 = A 1 + A 10 (-H -1 0 ) C, H j = A 1 + A 21 (-H -1 j-1 )C, 2 ≤ j ≤ N -1, H N = A 2 + A 21 (-H -1 N -1 )C. ( 19
)
Then we have the following recurrence relations for the expected values of first passage times to higher levels:

c (n) =      -H -1 0 1, n = 1, -H -1 1 (1 + A 10 c (1) ), n = 2, -H -1 n-1 (1 + A 21 c (m1-1) ), 3 ≤ n ≤ N.
For the determination of the expected values of the first passage time to lower levels we have first to denote by H n , n = 0, 1, ..., N the infinitesimal generator of the restriction of the process X t observed during those intervals of time spent at l(n) before the original process moves to l(n -1) for the first time. Then, all H n -processes, 1 ≤ n ≤ N are transient while H 0 is positive recurrent. The matrices H n are recursively computed as follows

H n =      B 0 + C(-H -1 1 )A 10 , n = 0, A 1 + C(-H -1 n+1 )A 21 , 1 ≤ n ≤ N -1, H N , n = N.
Then,

v (n) =      -H -1 1 (1 + Cv (1) ), n = 0, -H -1 n+1 (1 + Cv (n+1) ), 1 ≤ n ≤ N -2, -H -1 N 1, n = N -1.
These computations can be exploited, for instance, for controlling the (expected) time at which the packet queue overflows for the first time. Based on them, values for the control thresholds h s and packet arrival rates λ (i) p,s , s = 0..K can be determined numerically offline. If the computation is fast enough, online adaptation becomes possible. Another use of first passage times is to compute the average length of busy periods. From this metric, analyses about the opportunity to put the BS in sleep mode during idle periods can be envisioned.

Computing first passage times for energy levels is also of great interest, in order to determine and optimize the autonomy of the BS. However, the model we have presented is not directly a QBD for energy levels, and a more elaborate analysis is the topic of a work in progress.

An approximated product-form model

In this section, we describe a simplification of the model proposed in Section 2, strongly motivated by recent works in [START_REF] Gelenbe | A sensor node with energy harvesting[END_REF][START_REF] Gelenbe | Interconnected wireless sensors with energy harvesting[END_REF]. This new model is a G-network and therefore has product-form solution. The simplification applies to instances of our model where: a) there is no threshold control for the DEQ, therefore just environment-dependent packet as well as signal arrival rates λ (i) p , Λ -(i) , i = 1..M ; b) exactly one energy unit is required to transmit any data packet.

Consider the following simplifications: 1) infinite capacity buffers for the DEQ, SEQ, and the data queue; 2) as soon as the BS has both a data packet and enough energy to transmit that packet, the transmission is instantaneous; 3) environment-dependent parameters are averaged out using the stationary distribution q Y of the generator Q Y . Thus, defining the vector λ p = (λ ) and Λ -= (Λ - (1) , ..., Λ -(M ) ), then the arrival rates for packets, DEQ and SEQ energy units, and signals in the simplified model are as:

λ p = λ p q Y , λ -= Λ -q Y , λ ej = 1Λ ej q Y , j = 1, 2. (20)
Similarly, we will take into account the average leakage rate for each energy queue, i.e., if u j = 1U j q Y , j = 1, 2, where 1 is an 1 × M vector of 1.

The state of the approximation model can be represented by the pair (n, m), where n = 0, means that the BS has neither energy units in the DEQ, nor data packets to transmit, while n > 0, means that it currently stores n data packets but no energy units in DEQ, while n < 0, means that it stores -n energy units in DEQ but no data packets. The element m ≥ 0, counts the number of stored energy units in SEQ. Note that our system is a special type of G-networks. The underlying Markov process X = {(n, m); n ∈ Z, m ∈ Z + ∪ {0}} is ergodic provided that: λ p < λ e1 + λ -q 2 , λ e2 < u 2 + λ -, λ e1 + λ -q 2 < λ p + u 1 , (S) where q 2 = λ e2 /( u 2 + λ -), is the probability that the SEQ is not empty. These conditions come from the fact that the system must be stable with respect to the data packets (n > 0, left hand condition), with respect to the energy units in SEQ, (the condition at the middle), and with respect to the energy units in DEQ (n < 0, right hand condition). To conclude we state the following proposition.

Proposition 1. Under conditions (S) and with f (m) = q m 2 , m ≥ 0, the joint stationary distribution of X has the product form:

p(n, m) = Cg(n)f (m), (21)
g(0) = 1, g(n) = q n 1 , n > 0, g(n) = ( q1 ) -n , n < 0, ( 22)

q 1 = λ p λ e1 + λ -q 2 , q1 = λ e1 + λ -q 2 λ p + u 1 , C = (1 -q 1 )(1 -q1 )(1 -q 2 ) 1 -q 1 q1 . ( 23
)
The proof proceeds with the substitution of (21) into the global balance equations, and is omitted due to page constraint. For alternative approaches see [START_REF] Chao | Queueing Networks: Customers, Signals and Product Form Solutions[END_REF][START_REF] Harrison | Product-forms in multi-way synchronizations[END_REF]. C is obtained using the normalizing condition +∞ n=-∞ +∞ m=0 p(n, m) = 1. Extension to a network of SEQs. The approximation model is flexible enough to incorporate an arbitrary number of SEQs that are deployed close to the target BS while preserving the product-form solution. Indeed, assume that there are K reserve energy queues deployed close to our BS. In such a case our smart BS will generate signals at SEQ k, at a Poisson rate λ - k , while energy units arrive at SEQ k at a Poisson rate λ r e k , and energy leakage at a Poisson rate u 2k k = 1, ..., K. The system's state is now (n, m) = (n, m 1 , ..., m K ), n ∈ Z, m k = 0, 1, ..., where now m k , counts the number of stored energy units in SEQ k.

Then, provided that q1 = λp λe 1 + K k=1 λ - k q 2k < 1, q1 = λe 1 + K k=1 λ - k q 2k λp+ u1 < 1,

q 2k = λ r e k λ - k + u 2k
< 1, k = 1, ..., K, a modification of Proposition 1 can be proved, and the join stationary distribution of (n, m) is for f k (m k ) = q m k 2k , k = 1, ..., K, m k ≥ 0,

p(n, m) = C v(n) K k=1 f k (m k ), (24)
v(0) = 1, v(n) = qn 1 , n > 0, v(n) = ( q1 ) -n , n < 0, (25)
and the normalization constant is C = (1 -q1 )(1 -q1 ) K k=1 (1 -q 2k )/(1 -q 1 q1 ).

Conclusion

We have given the detailed specification of a versatile model of energy supply for base stations or similar devices and presented preliminary numerical results.

Our future investigations will concentrate on the resolution of models with a challenging size, the comparison with the product-form approximation, and also on ways to compute efficiently transient measures such as the mean time to depletion. With that objective, we will also investigate whether an approximate analysis through time decomposition is accurate. In a future work, we intent to specify the component interactions of our model by means of high level formalism such as a Markovian process algebra.

  m10 = ( S m11 , ..., S m1ν ) , whereS m1x = I (E2+1)×(E2+1) ⊗ [U 1 δ m1=1 +diag( m1k , k = m 1 + 1, ..., E 1 , m 1 ≥ 1.

2 )

 2 -1 = (A (l;x,y) m1,m1-1 ), k = 1, 2, ..., ν, and A(l;x,x) m1,m1-1 = U 1 , x = 1, 2, ..., ν, A (l;x,y) m1,m1-1 = I (E2+1)×(E2+1) ⊗ [T xy p (m1) 1x ], x, y = 1, 2, ..., ν, x = y. For k = m 1 -2, ..., 1, A(l) m1k = (A (l;x,y) m1k ), x, y = 1, ..., ν, where, A (l;x,x) m1k = 0 S×S , x = 1, ..., ν, and A (l;x,y) m1k = I (E2+1)×(E2+1) ⊗ T xy p (m1) m1-kx , x, y = 1, ..., ν, x = y and A m1k , k = m 1 + 1, ..., E 1 . , l = 1, 2, x, y = 1, ..., ν, m l = 1, ..., E l , where for T xy = diag(T

2 Fig. 2 .

 22 Fig. 2. Sensibility of DP w.r.t. data and RE for Λ - 0 = diag(0.35, 0.25) (left) and w.r.t. signal rate when the RE is in overload 5 First passage times

Table 1 .

 1 Overview of system's parameters


On Controlling Unknown Scalar Systems with Low Gain Feedback

Mazen Alamir, CNRS, University of Grenoble. (mazen.alamir@grenoble-inp.fr, http://www.mazenalamir.fr) I n this paper, a new feedback law if proposed to force a highly uncertain system to track a desired set-point using saturated control. The proposed feedback is not based on high gain concept, rather, a saturated integrator is used. Stability analysis is performed under realistic feasibility assumption.

Problem Statement

Let us consider the following scalar system:

ẋ = α [u -h] (α, h(•)) unknown (1)
where α ∈ [α min , α max ] is an unknown parameter with known bounds α min , α max , u is the control input, h is an unknown variable while x is the output of interest one would like to make it track some desired value x d .

Regarding the unknown term h, the following assumption is assumed to hold:

Assumption 1 (Bounded unknown h). There is a known pair of lower and upper bounds h min and h max such that for all t the following inequalties hold:

h(t) ∈ [h min , h max ] (2)
In order for the regulation problem to be feasible, it is assumed that the input control u can always dominate the term h had the latter be known. This is expressed through the following assumption which expresses also that u is constrained by lower and upper bounds:

Assumption 2 (u is authoritative enough). There exists lower and upper bound u min and u max such that for all t:

u(t) ∈ [u min , u max ] (3)
Moreover, these bounds relate to the bound of h through the following inequalities:

u max -h max ≥ + > 0 (4) h min -u min ≥ -> 0 (5)
Note that if this assumption is not satisfied, the sign of x cannot be guaranteed to be imposed by the control.

The problem we are interested in here is the following: Define a dynamic state feedback law of the form:

ż = f (z, x, x d ) (6) u = K(z, x, x d ) (7)
such that the tracking error e = x -x d is stabilized to a neighborhood of 0 that is as small as required.

Problem Statement

Discussion Regarding Alternative Solutions

There are obviously two intuitive approaches to solve the tracking part of the problem. These are, high gain solution and observer-based solutions. These are sketched in the following sections and then compared to the solution proposed in this paper.

Elementary High-Gain Solution

The idea here is to kill the unknown term h by a high gain control. This is done by using very high value of λ in the following feedback expression:

u = λ(x d -x) (8)
By doing so, the equation of the error e = x -x d becomes (assuming constant x d ):

ė = -αλe -αh (9)
which steers the errors towards the set:

E := e s.t. e ≤ max{|h min |, |h max |} λ ( 10
)
which can be made as small as required by taking high values of λ. The problem with this simple solution is that the transient of the input variable u cannot be guaranteed to meet the saturation constraints (3).

Observer Based Control

In this options, the idea is not to blindly kill the unknown term but to estimate it through a dynamic observer using the measurement of both u and x. More precisely, the observer is given by the following dynamic system. Unfortunately, the fact that α is supposed to be unknown makes the quantity h non observable. Indeed, α and h have to be both reconstructed and a simple examination of equation ( 1) leads to the fact that each pair of possibilities:

(α 1 , h 1 ) and (α 2 , h 2 ) (11)
are non distinguishable as soon as they are linked through the following relationships:

α 1 = βα 2 and h 2 = (1 -β)u + βh 1 (12)
which means that observability is rigorously lost as soon as u becomes almost constant. This is never a good news.

In what follows, a solution that avoid the explicit use of both explicit observer and high gain principle is proposed while α is still completely unknown.

The Proposed Solution

Define the saturation map according to:

S(v) :=    u max if v ≥ u max u min if v ≤ u min v otherwise (13)
The proposed feedback is defined by:

ż = λ f [S(λ(x d -x) + z) -z] (14)
u = S(λ(x d -x) + z) (15)
where λ > 0 and λ f > 0 are two design parameters to be appropriately chosen.

This control law if obviously of the form ( 6)-( 7). In the following section, the properties of the resulting closed-loop are analyzed.

Analysis of The Closed-Loop Properties

The closed-loop system equations can be written as follows:

ẋ = α [S(λ(x d -x) + z) -h] (16) ż = λ f [S(λ(x d -x) + z) -z] (17)
Let us write these equations in the three different configuration given by the definition (13), namely:

Case (u = u max ) ẋ = α [u max -h] (18) ż = λ f [u max -z] (19) Case (u = u min ) ẋ = α [u min -h] (20) ż = λ f [u min -z] (21)
Case (no saturation)

ẋ ż = -αλ α -λ f λ 0 x z + αλ -α λλ f 0 x d h (22)
or equivalently:

ẋ ż = -αλ α -λ f λ 0 x -x d z -h (23)
Note that the above three cases defines three regions in the 2D plan defined by the coordinates (x -x d , z) as it is shown in Figure 1. More precisely, these three regions are defined by:

A + := {(ξ, z) | z -λξ ≥ u max } (24) A -:= {(ξ, z) | z -λξ ≤ u min } (25) A 0 := {(ξ, z) | z -λξ ∈ (u min , u max )} (26)
Figure 1 shows also some information regarding the inclinations of the vector fields in the different regions that enable the following results to be proved.

Lemma 1. For any λ > 0 used in ( 14)-( 15), if λ f satisfies the following inequality:

λ f < + u max -u min λ (27)
then the set A 0 is attractive for all initial state such that

(x -x d , z) ∈ A + . ∆
Proof. This can be proved if one can prove that the angle θ depicted in Figure 1 is lower than arctan(λ). But Assumptions (4)-( 5) together with the fact that z necessarily belongs to [u min , u max ] enable to write (see Figure 1):

tan θ ≤ λ f (u max -u min ) + ( 28
)
which obviously gives the result.

Using the same arguments, the following Lemma follows:

Lemma 2. For any λ > 0 used in ( 14)-( 15), if λ f satisfies the following inequality:

λ f < - u max -u min λ ( 29
)
then the set A 0 is attractive for all initial state such that

(x -x d , z) ∈ A -. ∆
The previous two Lemmas obviously lead to the following corollary:

Corollary 1. For any λ > 0 used in ( 14)-( 15), if λ f satisfies the following inequality:

λ f < min{ + , -} u max -u min λ ( 30
)
then A 0 is globally attractive and invariant. ∆

x -x d z z -λ(x -x d ) = umax z -λ(x -x d ) = u min A 0 A + A - ≥ + ≤ λ f (umax -u min ) θ ≥ - ≤ λ f (umax -u min )
Figure 1: Definition of the three regions A0, A+ and A-used in the analysis of the stability of the closed-loop system with the proposed feedback law.

Proof. The attractivity is a direct consequence of Lemmas 1 and 2. The invariance results from the simple facts that when the state approaches the boundaries of A 0 with any of A + or A -, it is repulsed back before reaching the boundary whenever the requirements of these Lemmas hold. This is precisely implied by the condition (30).

The direct consequence of Corollary 1 is that the only thing that remains to be analyzed regarding the stability of the closed-loop behavior is related to the behavior of the closed-loop system inside the region A 0 . But the behavior inside A 0 is determined by the dynamic equation ( 23). This enables the following result to be established: Lemma 3. For any λ, if λ f satisfies:

λ f = r × α min λ 4 ; r ∈ (0, 1) (31)
then the matrix:

A 0 := -αλ α -λ f λ 0 (32)
possesses two real and strictly negative eigenvalues p 1,2 . More precisely:

p 1,2 = -αλ 1 ± √ 1 -r (33)
Proof. Straightforward since the discriminant of the characteristic equation is given by:

∆ = αλ [αλ -4λ f ] ( 34
)
the remaining facts directly follows. ∆

The previous results lead to the following Proposition:

Proposition 1. Take some λ > 0. If the following conditions holds:

1. λ f is such that

λ f < min min{ + , -} u max -u min , α min 4 × λ (35)
2. x d is constant and 3. the dynamic of the unknown term h satisfies:

dh dt ≤ δ (36)
then the tracking error e x = x -x d and the estimation error e h = z -h asymptotically satisfy the following inequalities:

lim t→∞ |x(t) -x d | ≤ δ λλ f ; lim t→∞ |z(t) -h(t)| ≤ δ λ f (37)
Proof. This is because condition (35) implies that the set A 0 is globally attractive. Therefore, the dynamic defined by (23) prevails after a finite time t 0 . Therefore, one has for all t ≥ t 0 :

e x (t) e z (t) = e A0(t-t0) e x (t 0 ) e z (t 0 ) - t t0 e A0(t-τ ) 0 ḣ(τ ) dτ (38)
and since (35) makes A 0 hurwitz invertible, the last expression asymptotically behaves like:

A -1 0 0 ḣ(t) = 1 λ f 0 -1 λ λ f α -1 0 ḣ(t) = -1 λλ f -1 λ f ḣ(t) (39)
This together with (36) obviously gives the result. ∆

Illustrative Simulations

Consider the unknown signal h given by:

h(t) = 5 -1.3 cos(1.2t) -1.6 sin(2t + π/6); (40)
The time evolution of h on the time interval [0, 10] is shown in Figure 2. Straightforward computation shows that δ = 4.76 satisfies the inequality (36). Moreover, h min and h max are given by:

h min = 2 ; h max = 8 (41)
Se take the unknown α = 1 which is not used in the feedback design although it changes the dynamics of the closed-loop system. Let assume however that one knows that α min = 0.2 and that the control bounds are given by:

u min = 1 ; u max = 10 (42)
so that -= 1 and + = 2 can be guaranteed in ( 4)-( 5).

It results that the condition (35) on λ and λ f becomes:

λ f < min min{1, 2} 10 -1 , 0.2 4 × λ < 0.05 × λ (43)
In the simulation, λ f = {0.04, 0.001} × λ are usedfor several values of λ are used to illustrate the properties of the closed-loop systems.

The results are shown on Figure 3, 4 and 5. More precisely: The message from these simulations are the following:

1. The theoretical predicted bounds given by (37) of Proposition 1 seems not only correct but rather tight (no over conservatism) as the effectively encountered levels get very close to the bound regularly.

2. Very small error can be achieved without uncontrollable transient on the control as the latter respects the saturation constraints anyway.

3. Low time derivative of u can be obtained low values of λ f as it can be observed by comparing Figures 4 and Figure 5 4. The better estimation of h is obtained through u and not necessarily z. This is especially true for low values of λ f . This can be confirmed by observing Figures 6, 7 and 8 that show the evolution of u and h as well as their difference |u -h| for the same scenarios.

Figure 3

 3 Figure 3 uses λ = 400 and λ f = 0.04λ = 16.

Figure 4

 4 Figure 4 uses λ = 4000 and λ f = 0.04λ = 160.

Figure 5

 5 Figure 5 uses λ = 4000 and λ f = 0.001λ = 4.

Figure 2 :Figure 3 :Figure 4 :Figure 5 :Figure 6 :

 23456 Figure 2: Illustrative simulation: Evolution of the unknown signal h(•) used in the simulation of ẋ = α(u -h) under the proposed feedback which ignores α and h.

Figure 7 :Figure 8 :

 78 Figure 7: Illustrative simulation: Comparison between u and h when λ = 4000 and λ f = 0.04λ are used.


Keywords: EMC, High Frequency, Physical Modeling of Inductor, SPICE

Multi-physics modeling in power electronics contains EMC studies, thus conception requires modeling of passive elements. These models need to be compatible with the simulation software usually used in the industry. An analysis for modeling an inductor at high-frequency is presented in this study. The method takes into account the proper inductances of the turns, the mutual inductances between them, the parasitic capacitances, and the resistance of the windings. The originality of the model is that it is based on physical geometries and the material characteristics of the inductor and covers a large frequency bandwidth. These elements are used to create an analytical model based on electrical equations in order to provide the frequency dependent impedance of the inductor. The model is validated using experimental measurements.

I. INTRODUCTION

Multi-physics modeling in power electronics contains EMC studies, thus conception requires modeling of passive elements. These models need to be physical in order to be more effective, compatible with the simulation software usually used in the industry, and predictive (before the prototyping phase).

At high frequencies, these models will have a different behavior than at low frequencies. Many factors can affect this behavior and force the components of any model to be frequency-dependent.

Previous studies have used a number of different methods to model high frequency inductors. Some are based on the finite element method (FEM) [START_REF] Zhou | Modèles de composants passifs et couplage électromagnétique pour filtres HF de puissance -Optimisation du placement[END_REF], which consumes a lot of hardware and time, and can only provide a global model that does not take into consideration the parasitic elements of the inductor.

Other models are physical analytical models based on lumped parameters that give an equivalent circuit model with global frequency-dependent L, C and R [START_REF] Kazimierczuk | High-Frequency Small-Signal Model of Ferrite Core Inductors[END_REF][2] [START_REF] Liu | High frequency model of common mode inductor for EM1 analysis based on measurements[END_REF][8] [START_REF] Ragusa | High Frequency Modelling method of EMI filters for Hybrid Si -SiC Matrix Converters in Aerospace Applications[END_REF]. These models are not optimized for use in an electric circuit simulator such as SPICE. [12] A number of methods represent the non-linear core model as a series of LR branches in order to take into consideration the hysteresis cycle [START_REF] Heinemann | Modeling of High Frequency Inductors[END_REF]. These models do not account for the physical parasitic seen at high frequencies and also replace the capacitances between coils by one equivalent capacitor. Such modifications can influence the time domain simulations as the model does not completely reproduce reality.

The purpose of this study is to give a physical model of an inductor based on its geometry and the characteristic of its core and windings while taking into account the inductances, the capacitances, and the core losses which exist. These parameters are defined using formulas from the literature [START_REF] Massarini | Lumped Parameter Models for Single and Multiple-Layer Inductors[END_REF][4] [START_REF] Grandi | Stray Capacitances of Single-Layer Solenoid Air-Core Inductors[END_REF][6] [START_REF] Hole | Stray capacitance of a twolayer air-cored inductor[END_REF].The model will represent the real aspect and behavior of the inductor at both low and high frequencies, being a physical model rather than a lumped parameter equivalent model.

The study models each turn by the electrical equation between 2 nodes. Following this, all capacitances between nodes are included in order to define a resonant model of the inductor. The frequency response of this model is compared to a real measurement in order to validate it. Finally, a time domain simulation of the model is presented.

II. MODEL DEFINITION

A. Structure

The inductor consists of n turns on the same layer as shown in fig. 1. The approach involves modeling each turn based on its physical and electromagnetic characteristics as well as modeling the coupling between the turns. Fig. 1: Inductor of n turns on the same layers Fig. 2 shows an electromagnetic representation of a winding where R i represents the resistance of the turn i (including the core losses), L i is its proper inductance , and M ij is the mutual inductance between turns i and j.

i V i
to be integrated in the system. The next subsections will present the method of calculation used for determining the parameters mentioned above.

B. Proper and Mutual Inductances

Grover [START_REF] Grover | Inductance Calculations, Chaps. 2 and 13[END_REF] presented formulas and methods for calculating the proper and mutual inductances for different coil forms.

This work focusses on coaxial circular turns. The proper inductance for turn i can be calculated using (1):

8 7 4 ( 1
)
where µ is the permeability of the core, r the radius of the conductor, and R the radius of the turn.

For the mutual inductance calculation, we consider having two turns i and j as shown in fig. 3. Based on Neumann's formula, the mutual inductance between the two turns in question is determined as follows:

4 . ( 2
)
where ds i and ds j are incremental sections of the conductors and d is the distance between them. The integral in (2) can be represented in terms of the complete elliptic functions K and E in (3) and ( 4):

1 sin (3) 1 sin (4)
Using these elliptic functions, the mutual inductance between turn i and turn j is determined as in (5):

2 √ 2 (5)
where:

2 4 2 2
Considering a solenoid consisting of n turns, for each, it is necessary to determine the proper inductance and mutual inductance between the turn in question and all other turns. Knowing all inductances in the structure, we the total inductance of the solenoid is obtained:

2 (6)
This factor gives a good image of the inductor in the first frequency range (before the main resonance frequency) where the parasitic parameters of the inductor will be neglected beside the inductive dominance.

(1) and ( 2) show that the proper and mutual inductances depend on the permeability of the core. In the case of a solenoid (rod), the electromagnetic field runs outside the core and closes its circuit through the surrounding air. This will change the effective permeability of the core. Instead of having: [START_REF] Liu | High frequency model of common mode inductor for EM1 analysis based on measurements[END_REF] where µ r is the relative permeability of the core and µ 0 the permeability of air, the relative permeability will be replaced by an effective permeability. This last depend on the core geometry as well as its relative permeability and can be obtained using the plot in fig. 4 [13]

C. Capacitances

In this section, the calculation of the capacitances between turns or between a turn and the core will be presented. The impedance of the capacitance is inversely proportional to the frequency. At low frequencies, this impedance will be neglected compared to the inductance's which is directly proportional to the frequency. In order to model an inductor at high frequency, it is important to determine the value of these capacitances and to integrate them in the model.

In [START_REF] Massarini | Self-Capacitances of Inductors[END_REF], A. Massarini proposed a method to determine these capacitances by studying the lines of the electric field. Considering 2 conductors as shown in fig. 5, the electric field between the 2 conductors is present in the insulating coating as well as in the air between the 2 conductors. The elementary capacitance of the insulating coatings can be determined as mentioned in (8):

0 2 0 (8)
where l t is the turn length, r c the inner radius and r 0 the outer radius.

As for the capacitance of the air gap, it is necessary to determine the path of the electric field:

2 1 (9)
The elementary surface of the wire is:

2 (10)
Using ( 9) and ( 10), we can determine the elementary capacitance per unit angle of the air gap: [START_REF] Ragusa | High Frequency Modelling method of EMI filters for Hybrid Si -SiC Matrix Converters in Aerospace Applications[END_REF] The series combination of the elementary capacitances mentioned above is given by:

2 (12)
And finally, to obtain the capacitance between the turns, (12) is integrated:

2 / (13)
The equation above can be used to calculate the turnto-core capacitance. Since the core is at a constantpotential and located on the vertical plane of the conductors as shown in Fig. 2, the path lengths of the electric field studied above in the air gap between the turn and the core are one half of the path lengths between two conductors on the same layer. As a result, the capacitance between a turn and the core can be determined:

2 (14)
The formulas presented in this subsection permit the determination of every parasitic capacitance in our system whether it was an adjacent turn-to-turn capacitance, non-adjacent turn-to-turn capacitance, or turn-to-core capacitance.

D. Resistances

The system resistance is due to both the winding resistances R w and the core series resistance R c . Both of these resistances are frequency-dependent. For the winding resistance, skin-effect and proximity-effect should be taken into consideration. An analytical representation of these effects is shown in (15) [START_REF] Kazimierczuk | High-Frequency Small-Signal Model of Ferrite Core Inductors[END_REF]:

2 2 2 2 (15)
where R wdc is the DC resistance and A for a round wire is:

4 . . . (16)
In ( 16), d represents the conductor diameter and t the distance between the centers of two adjacent conductors. The skin depth of the conductor δ is given by: (17) where ρ w is the wire resistivity.

On the other hand, the frequency-dependent core resistance needs the knowledge of the frequencydependent permittivity and permeability of the core.

The frequency-dependent permeability represents the core magnetic hysteresis loss:

̂ ̂ (18) r 0 r c x(θ) θ D
To include the core dielectric loss, the electric permittivity of the core is presented by a complex quantity that is also frequency-dependent:

̂ ̂ (19)
Using these definitions, the core series resistance can be calculated as shown in (20): . √2

. √2

(20)

The resistance in (20) is represented using the Bessel function of first kind,

1 2 (21)
and δ:

1 ̂ 2 ̂ ( 22
)
where δ c is the electric conductivity of the core. In this section, the determination of the frequencydependent resistance was presented, which is the series association of the conductor resistance and the core equivalent resistance.

E. SPICE Model

The purpose of this work is to develop a Spice compatible physical model. In order to do that, every section studied must be represented using Spice models. To model a solenoid with n turns in Spice, the system is considered to be composed at first of n+1 major nodes in series. Each turn will be represented between two nodes.

The electromagnetic system will be modeled by arbitrary behavioral voltage sources between two major nodes. The proper and mutual inductances in the system can be calculated as mentioned in section II-B. Using (23) for each turn i models the entire electromagnetic system of the solenoid:

(23)
where I i is the current in the turn i.

In addition to the arbitrary behavioral voltage, the wire and core resistance calculated in section II-D were added in series. After modeling all turns in such a way, it is necessary to add the electrostatic system. In section II-C, it was shown how to determine the capacitance between two turns on the same layer or between a turn and the core. The capacitance between the turn i and j will be placed between the node i and the node j. On the other hand, the capacitance between the turn i and the core will be placed between the node i and a common node for all the turns N core . Fig. 7 shows a portion of the model with the capacitances added.

III. EXPERIMENTAL VALIDATION

In order to validate our model, experimental results are necessary. Fig. 8 shows a solenoid prepared in our laboratory in order to be used in this study: The device will measure the S ij parameters as a function of frequency. These parameters will be used in order to calculate the impedance of the solenoid and thus its real part, imaginary part, module and phase.

Fig. 10 and Fig. 11 shows a comparison between the real part of the real solenoid and the physical model proposed in this article. It can be seen that the real part of the impedance of the model is close to real measurements. For low frequencies, the real part does not extend few Ohms and as soon as we approach the first resonance, it increases exponentially and hits its maximum value as the mentioned frequency. From there, increasing the frequency results in decreasing the impedance real part exponentially. It can be seen that the comparison between the proposed model and the measurements shows a good agreement until a certain frequency. Below the first resonant frequency, there is little difference between the model and the measurements. When the frequency reaches 5 MHz, slight differences can be seen but the model continue is still able to locate a second and third resonance at 12 MHz and 20 MHz. Above these frequencies, the limitations of the model prevent comparison.

IV. TIME DOMAIN SIMULATION

As previously mentioned, the model proposed in this study can used in both time and frequency domains. In the previous section, the frequency-domain validation has been studied. The efficiency of the model in the time domain must also be verified. In order to do this, the proposed model is excited by a current source of 1A and the voltage at its terminals is studied. Fig. 14 presents the result of the time simulation. In this figure, two regions can be found in the simulation. The transient state shows the effect of the parasitic elements in the solenoid and is presented by a harmonic oscillator. The inductances and capacitances in the system have an effect on the frequencies of the oscillations. On the other hand, the resistances present in the model will change the oscillations amplitude. As for the permanent state, it's a constant voltage relative to the resistive part of the inductor.

The last simulation shows the response of the inductor model to a constant current. In power electronics, currents switch quickly in the presence of active components. Fig. 15 shows the response of our model to a pulse signal. An approach providing a physical model of an inductor based on its geometry and magnetics and electrostatic characterization has been introduced here. This model gives a better understanding of the behavior at high frequencies and parasites and is SPICE compatible. Previous models which can be used in SPICE did not give this level of precision at high frequencies. This work provides an efficient model both in the frequency (low and high frequencies) and the time domains.

Fig. 2 : 1 1

 21 Fig. 2: Electromagnetic Model of a windingAs for the electrostatic system, the parasitic capacitances between the turns (consecutive or not) need

I

Fig. 3 :

 3 Fig. 3: Cross-section of 2 turns

Fig. 4 :

 4 Fig. 4: Effective permeability of rods

Fig. 5 :

 5 Fig. 5: Cross Section of 2 conductors with the

  Fig 6 represents the model for turn i which includes each of these elements:

Fig. 6 :

 6 Fig. 6: Modeling a turn of the solenoid in Spice

SHOFig. 7 :

 7 Fig. 7: A part of the model showing the 3 parts studied

Fig. 8 :

 8 Fig. 8: Solenoid used to validate the model The solenoid in question is composed of a single layer over a ferrite bobbin core. The geometric characteristics are listed as follows: • R : Radius of the solenoid: 15 mm • r : Radius of the conductor: 0.25 mm • e : Distance between 2 adjacent turns : 0.54 mm • n : Number of turns : 31 • u reff : Effective relative permeability of the core: 4.3 In order to measure the impedance of the solenoid, an Agilent E5071C Network Analyzer has been used (fig. 9).

Fig. 9 :

 9 Fig. 9: Agilent E5071C Network Analyzer

Fig. 10 :

 10 Fig. 10: Real part of the solenoid

Fig. 12 :

 12 Fig. 12: Imaginary part of the solenoid Fig. 13 compares the module of the impedance of the model presented in this work and the measurements.

Fig. 13 :

 13 Fig. 13: Impedance Module of the solenoid

Figure 14 :

 14 Figure 14: Time simulation of the model

Figure 15 :

 15 Figure 15: Voltage response of the model to a pulse current


J Taki
email: jad.taki@centralesupelec.frmohamed.bensetti@centralesupelec.frdaniel.sadarnac@centralesupelec.fr

M Bensetti

D Sadarnac

REDUCTION OF ELETROMAGNETIC PERTURBATIONS BY OPTIMIZING THE PRINTED CIRCUIT BOARD

EMC(Electromagnetic Compatibility) studies require the modeling of the printed circuit board (PCB). Such models need to be compatible with the simulation software usually used in the industry. An analysis for modeling and optimizing the PCB is presented in this study. The method of moments is used to determine a transmission line electrical model (RLGC) for traces. This model is based on the physical geometry of traces and the material characteristics of the copper and the dielectric. Given the PCB model, an optimization can be performed to minimize the conduction emissions of the board.

Introduction

Technological advances in power electronics tend to raise the switching frequency. As a result, the conducted emissions we measure at the battery may also be raised. These perturbations may interact with all components in a power electronics board. The source of these emissions includes active components, but PCB (Power Circuit Board) traces conduct these perturbations and tend to raise them at high frequencies. EMC standards are becoming more restrictive, driving an increased interest in modeling the PCB and reducing its effect on emissions. These models should be fast to generate and compatible with electric circuit simulators suck as SPICE. Previous studies have used a number of different methods to model the PCB at high frequencies. Some are based on simplified analytical equations [START_REF] Poh | Approximate Formulas for Line Capacitance and Characteristic Impedance of Microstrip Lines[END_REF] [START_REF] Sadiku | Comparison of Approximate Formulas for the Capacitance of Microstrip Line[END_REF] which sometimes neglect the mutual inductances and capacitances between traces. Others are based on the finite element method (FEM) which consumes a lot of hardware and time, and is not ideal for an optimization process. The purpose of this study is first to optimize a PCB such that the electromagnetic perturbations are reduced. In order to do this, a PCB model is needed. An analytical method has been chosen to generate a transmission line model (RLCG) that is frequency-dependent. This paper will introduce the method used to model the PCB, and give a general idea of the electrical models that we will obtain. After that, it will discuss how these models will be integrated in the electrical circuit simulator and the different types of simulations that we can use in order to determine the electromagnetic perturbations. The section that follows introduces the optimization procedure that is used in this work in order to reduce the conducted emissions. This method will be used on a DC-DC buck converter to validate the procedure.

Modeling

As the purpose of this work is optimizing the PCB, the selected model should be fast to generate and one which can be electrically simulated in SPICE tools. Fig. 1 presents a transmission line model (RLCG) that will represent a PCB trace where dz is an elementary length.

Fig. 1: Transmission line model

R represents the resistance of the trace per unit length which is frequency dependent due to the skin effect. L is the inductance per unit length of the loop formed by the trace and the ground plane. This inductance varies with frequency and depends on the distance between the trace and the ground plane. C is the capacitance per unit length between the trace and the ground plane. This element depends on the permittivity of the dielectric which can also be frequency-dependent. G is the conductance per unit length of the trace which represents the dielectric losses. Usually this term is neglected compared to the capacitance of the trace. The line transmission model presented in fig. 1 is for a single trace on top on a ground plane. In a PCB circuit, traces will interact with each other's and will create a capacitance and a mutual inductance between them.

Resistance

The resistance per unit length of the trace has been determined using analytical equations [START_REF] Lafon | Développement de techniques et de méthodologies pour la prise en compte des contraintes CEM dans la conception d'équipements du domaine automobile. Etude de l'immunité, du composant à[END_REF]. It is divided into two sections, the resistance per unit length of the trace on low frequencies and the resistance per unit length on high frequencies.

𝑅 = 𝜌 𝑆 (1 + √ 𝑓 𝑓 𝑠𝑘𝑖𝑛 ) (1)
where ρ is the copper resistivity, S the area of a section of the trace, f is the frequency and f skin is the frequency at which the skin depth is of the same order of the trace thickness. This frequency can be determined using (2) :

𝑓 𝑠𝑘𝑖𝑛 = 4ρ 𝜋𝜇 0 𝑡 2 (2)
where µ 0 is the magnetic permeability of the air and t the thickness of the trace. Using (1) and ( 2), we can get the frequency dependence of the trace resistance. Since we mentioned that the model should be SPICE-compatible, the resistance should be presented as a function of frequency too in SPICE.

There are two ways in order to model this effect. Either the frequency is presented using the s parameter where the relation between s and the frequency is given in (3):

𝑠 = 𝑗 * 2𝜋𝑓 (3)
or the resistance will be modeled by a lumped parameter circuit that will take into consideration the frequency-dependence as shown in fig. 2.

Inductance and Capacitance

A PCB circuit creates some inductances and capacitances in the system:  Proper inductance of the trace which represents the inductance of the loop formed by the trace and the ground plane.  Proper capacitance of the trace which represents the capacitance between the trace and the ground plane.  Mutual inductances between traces on the same plane.  Capacitances between traces on the same plane.

In order to determine all these inductances and capacitances, the Method of Moments [START_REF] Clayton | Analysis of Multiconductor Transmission Lines[END_REF][5] [START_REF] Tidjani | Electromagnetic Coupling Modeling of Microstrip Lines Used in Switching Power Supplies[END_REF] has been used. It is a numerical computational method of solving linear partial differential equations which have been formulated as integral equations:

∅(𝑥, 𝑦) = 1 4𝜋𝜀 ∬ 𝜌(𝑥 ′ , 𝑦 ′ ) 𝑟 𝑑𝑠 𝑠 (4)
(4) is the integral form of Poisson's equation in the transverse x-y plane where a charge distribution in the x-y plane, ρ, is distributed over some surface s. We know how to determine the potentials at every point in the region and we wish to determine the charge distribution that follows it since the capacitance is related to the potential as well as the charge distribution as we can see in [START_REF] Harrington | Field Computation by Moment Methods[END_REF] where Q is the charge distribution and V is the differential potential.

𝑄 = 𝐶 * 𝑉 (5)
In order to determine the capacitance of the system, (4) and ( 5) are used with relative permittivity of the dielectric between the traces and the ground plane. As for the inductances, [START_REF] Tidjani | Electromagnetic Coupling Modeling of Microstrip Lines Used in Switching Power Supplies[END_REF] gives a relation between the inductance L, the permittivity of the air ε 0 , the electromagnetic permeability of the air µ 0 and the capacitance of the traces if we replace the dielectric in the system by air C 0 :

𝐿𝐶 0 = 𝜇 0 𝜀 0 (6)
If the Method of Moment is used to determine the capacitances in the air of the system, (6) will give us all the inductances needed to have a proper PCB model. An important note to mention is that the Method of Moments gives us inductances and capacitances per unit length and in order to have a model that is valid and precise on high frequencies, it is important to choose the length of each section of the trace and so then number of sections of the trace.

Optimization Procedure

In order to reduce the electromagnetic perturbations, it is important to generate a PCB that does not increase these perturbations. To apply this procedure to a power electronics circuit, it is necessary to have models of every element. Passive components are modeled by lumped circuits that are valid at high frequencies using measurement from a network analyzer. After placing the components, the routing must be developed. Initially, the routing of the PCB must follow some general rules such as maintaining a distance between the traces and the edges of the board, separating the power section from the command and avoid right angles.

In the optimization procedure, it is considered at first that any routing between two pins of a component is composed of two traces as shown in fig. 3.

Fig. 3: Routing technique

By changing the position of the intersection of the two traces, the routing changes and thus its electromagnetic effect on the conductive perturbations. This would be the main variable to optimize in the procedure. The optimization algorithm used in this work is the genetic algorithm [START_REF] Mitchell | An Introduction to Genetic Algorithms[END_REF]. This algorithm should have three parts:  Selecting the right generation of variables  Producing the PCB model for the chosen individual at the current iteration  The objective function in order to find the best individual that responds to the problem which is reducing the electromagnetics conductive emisssions.

For the selection, the algortihm should at first a fixed number of individuals forming a generation. The choice of the individuals (its values) should follow the limitations giving in the description of the algorithm. For example, a limitation can be the intervals that will limit the choice of the x and y position of the intersection point. In the first iteration, the algorithm will calculate the objective function for each individual of the generation and at the end of the iteration it will choose the best individual (the lowest objective function) with other mutated ones to participate in the next iteration. The algorithm will continue throught iterations until it attends a certain precision or obtain the same lowest objective function for several iterations consecutivally. As a result, the algorithm will give the best individual as well as the value of the objective function corresponding to this individual. For each selected individual, a PCB model should be produced. The board will be meshed and the routine will begin using the technique mentioned before with only one intersection between 2 pins. After the routing is finished, the program start searching for traces that could group to form one model using the Method of Moments. The traces that the program searches for are traces who are in the same colomn and have the same angle between each others as shown in fig. 4. We can notice that there is four different regions in this example. Each region have traces with same angle between them and consist a unique PCB model. These models will be linked to either component pins or other parts of traces.

As mentioned before, the objective function is to find the best individual that reduces the electromagnetics conductive emisssions. The function chosen to get the best result is a relative comparison between the sum of the spectrum of perturbations over frequency between the case for the actual indivual and the one for the current best indivudual. If the comparaison is negative, the current individual will be the current best individual. The algorithm will keep going all over the avaible individuals until the best one is found.

Study case: A DC/DC converter

The optimization technique is being validated in this section on a DC/DC converter as shown in fig. 5.

Fig. 5: DC/DC converter

Each component in this converter has been modeled using a Vector Network Analyser (VNA) and a fitting technique to be efficient up to 1 GHz and the switching frequency is 200 KHz. The PCB optimization procedure is used to get the most efficient routing regarding the EMC conducted emissions. These perturbations are being measured by two LISN (Line Impedance Stabilization Network) that are placed between the battery and the decoupling capacitor. Fig. 6 shows an electrical model of the LISN used in this study. Fig. 7 shows a geometry comparaison between a random PCB and an optimized one created by the algorithm described in this article. The board has been meshed into 50 units through x and y axis with 5 mm distance between 2 units. It is clear that the more the number of meshing units is higher, the more choices of routing there is and a better optimized one. Components of the converter are presented by blue lines (excluding the charge) and the traces between them are in red. As we can notice, both routing are not regular and have no repetitive rules and it is clear that all right angles have been avoided even during the optimization results. The association of all models in the converter has been simulated in LTSpice [8] throught a time domain simulation where the results have been sent to Matlab in order to generate the perturbation spectrum throught the fast fourier transformation (FFT).

Spectrum comparison

Two types of results will be shown in this section. Fig. 8 compares the conductive emissions spectrum between a converter with or without PCB and Fig. 9 compares the spectrum between the converter with the random PCB shown in Fig. 7.a and the optimized one in Fig . 7.b. It is noticeable in Fig. 8 that the PCB starts taking effect around 30 MHz and that it is important to take this fact into consideration in the simulation if the prediction of the conductive emissions above this frequency is needed. Fig. 9 and 10 shows that an optimized PCB can push some resonances to higher frequencies as well as reducing the amplitude of these resonances. A difference of more than 10 dB appears around 300 MHz and few dB between 50 MHz and 100 MHz.

Conclusion and Perspective

In this study, a PCB optimization procedure has been introduced in order to reduce the electromagnetic conductive emissions. Each link between two pins has been divided into two traces and the optimization algorithm chooses the best positions for these intersections in order to satisfy the objective function. This study had the components placement as a fix variable and the only optimization was on the traces. The results shown for a DC/DC converter study case proves that optimizing the PCB reduces the perturbations on high frequencies.

In order to improve the quality of the circuit from an EMC and PCB view, multiples ameliorations can be introduced. The number of traces per link can be increase and so the number of entries to the optimization algorithm will increase. The meshing is another factor that plays a big role in the precision and choices of traces. As mentioned, the components' place are fixed in this study so adding this factor as an entry to the algorithm can improve the EMC quality of the board but can also make the study more complex and takes more time.

Fig. 2 :

 2 Fig. 2: Lumped parameter model of a frequency-dependent resistance

Fig. 4 :

 4 Fig. 4: Traces Separation for modeling

Fig. 6 :Fig. 7 .Fig. 7 .

 677 Fig. 6: LISN electrical model

Fig. 8 :Fig. 9 :

 89 Fig. 8: Comparaison between the spectrums of the converter with or without the PCB

Figure 10 :

 10 Figure 10: Envelope of the Spectrums for the converter without and with a random and optimized PCB