Benjamin Heymann
email: benjamin.heymann@polytechnique.edu

Pierre Martinon
email: pierre.martinon@inria.fr

Cmap

Optimal Battery Aging: an Adaptive Weights Dynamic Programming Algorithm

Keywords: Control, Aging, Dynamic Programming, Energy Management System Mathematics Subject Classification (2000) 93A13, 93C15, 90C39, 49L20, 49M27, 49M29

We present an algorithm to handle the optimization over a long horizon of an electric microgrid including a battery energy storage system. While the battery is an important and costly component of the microgrid, its aging process is often not taken into account by the Energy Management System, mostly because of modeling and computing challenges. We address the computing aspect by a new approach combining dynamic programming, decomposition and relaxation techniques. We illustrate this 'adaptive weight' method with numerical simulations for a toy microgrid model. Compared to a straightforward resolution by dynamic programming, our algorithm decreases the computing time by more than one order of magnitude, can be parallelized, and allows for online implementations. We believe that this approach can be used for other applications presenting fast and slow variables.

Introduction

A microgrid is an electric system that includes some electricity generation units and a battery to store the energy for later use. The energy management system typically enforces the optimum balance between the different energy sources, in order to minimize some global operating cost. This leads to solving an optimization problem with a short time step (say 15 minutes) over a very long time-frame (such as 10 years). The long time-frame means that the aging of the battery becomes non-negligible, both due to the gradual loss of efficiency and the high cost of the battery. However, battery aging is hard to take into account in the optimization, leading to sub-optimal policies. Indeed, there is some trade-off between the short and long term costs: a policy making heavy use of the battery may be better in terms of daily operating costs, but causes a faster aging of the battery, thus leading to a higher global cost. The two main difficulties are the profusion and complexity of aging models, as well as the numerical cost of solving a long horizon optimization problem.

The modeling and optimization of the battery aging is an active interdisciplinary field of research. In [START_REF] Haessig | Aging-aware NaS battery model in a stochastic wind-storage simulation framework[END_REF], Haessig et al. propose a simulation with an aging model to perform a cost analysis, but without optimization. In [START_REF] Riffonneau | Optimal power flow management for grid connected PV systems with batteries[END_REF], Riffonneau et al. use a discrete time dynamic programming approach to solve an optimal power flow problem. The battery aging is proportional to the discharge of the battery and linearly decreases the capacity of the battery. In [START_REF] Palma-Behnke | A microgrid energy management system based on the rolling horizon strategy. Smart Grid[END_REF], Palma et al.

integrate the battery aging in a rolling horizon strategy model. The aging is taken into account using a working zones approach as proposed in [START_REF] Guasch | Dynamic battery model for photovoltaic applications[END_REF], with the battery cost as the penalty parameter.

To our knowledge there are basically three approaches in the literature to handle battery aging in an optimization model:

-Some constraints on the control and state variables (for instance, to avoid extreme State of Charge values): The main difficulty consists in choosing the constraints, which can be too conservative or not enough. Also, the aging process is not directly taken into account in the optimization program.

-A penalization of the aging: This requires choosing an aging model and a penalty parameter.

-An aging constraint: This requires choosing an aging model and the aging level.

In the present work, we introduce a variation on the penalty approach. We focus on an efficient way to solve the optimization problem over a long horizon with a short time step, as we would like the optimization to be fast enough for an actual online implementation. We propose a new 'adaptive weight' method based on dynamic programming and a coupling of decomposition and relaxation techniques. The approach consists in two parts, the first being an offline optimization that first solves a family of short term problems associated to a discretized set of battery ages and aging penalty weights, and then store the resulting costs and age increments. The second part is an online optimization that uses this data to compute the daily optimal policy. This decomposition makes both the online and offline problems one-dimensional in state, allowing to solve them efficiently with dynamic programming.

We illustrate the adaptive weight algorithm (AWA) on a simple microgrid model. Section 2 recalls the optimal control problem for the energy management of the microgrid. Section 3 introduces the adaptive weight algorithm to solve the offline and online optimizations. In Section 4, we compare the numerical results of AWA with a brute force solving of the original problem by dynamic programming.

Motivating Application

Microgrid

We consider a simpler version of the microgrid problem studied in [START_REF] Palma-Behnke | A microgrid energy management system based on the rolling horizon strategy. Smart Grid[END_REF][START_REF] Heymann | Continuous Optimal Control Approaches to Microgrid Energy Management[END_REF][START_REF] Heymann | A Stochastic Continuous Time Model for Microgrid Energy Management[END_REF]. The electricity is produced by some non dispatchable units (solar panels) and a dispatchable unit (diesel generator).

At anytime, there is an instantaneous (uncontrollable) demand for electricity. A battery storage system can either store excess energy or provide part (or all) of the demand. The battery is not a perfect storage and it has an efficiency strictly lower than 1. The aim is to satisfy the demand at all times while minimizing the overall cost associated with the dispatchable unit (diesel).

We assume the power load (demand) P L (t) and the solar production P S (t) to be given and T -periodic, where T is the duration of a day. The diesel generator can produce a power u ∈ [u min , u max ]. We assume the resulting cost to follow a quadratic relation:

(u) = βu 2 , β > 0.

(

At all times, the power equilibrium is satisfied and allows one to express the battery power P b (with the convention P b > 0 for discharge)

P b (u, t) = -u -P s (t) + P l (t). (2)
The battery properties can be summarized by three parameters: total capacity C, charge efficiency ρ i and discharge efficiency ρ o . The state of charge (SOC) of the battery c ∈ [0, 1] is defined as the quantity of energy stored in the battery divided by the capacity, and follows the dynamics

ċ(t) = F c (u, t) = 1 C (ρ i (P b (u, t)) --(P b (u, t)) + /ρ o ). (3)
In the rest of the following, we impose c to have the same value at the beginning and at the end of the day. This is often done in practice to avoid the 'end of the world effect': with no final constraint, the optimization would typically empty the battery at the end of the day. We point out that the methodology presented thereafter also applies without this periodicity condition.

Battery Aging

We define the aging of the battery over time with a scalar a ∈ [0, 1], with the convention a = 0 for a brand new battery, and a = 1 when the battery is dead. We choose a severity factor model, adapted from [START_REF] Svoboda | Operating conditions of batteries in off-grid renewable energy systems[END_REF], in which the aging evolves depending on the state of charge and charging power:

ȧ(t) = F a (c, u, t) = (-4c 2 + 5) 5
(P b (u, t)) - K , K > 0. ( 4
)
The main idea behind this model is that charging a nearly full battery will damage it more than at a low state of charge. We model the aging impact on the performance by a simple linear decrease of the charge efficiency ratio

ρ i (a) = (1 -a)ρ 0 i
where ρ 0 i is the value for a brand new battery. The state of charge dynamics (3) thus becomes:

F c (a, u, t) = 1 C (ρ i (a)(P b (u, t)) --(P b (u, t)) + /ρ o ). (5)

The Long Term Problem

We want to minimize the operating cost (u(t))dt over a long horizon N T , and we add a penalty function ψ to account for the value of the battery at the end of the time horizon. The optimization problem thus writes

inf u N T t0 (u(t))dt + ψ(a N T ), (6)
subject to ( 5) and [START_REF] Guasch | Dynamic battery model for photovoltaic applications[END_REF]. In what follows, we denote by V (a 0 , t 0 ) the value of this problem for an initial time t 0 and an initial condition a(t 0 ) = a 0 . Note that (6) has a 2-dimensional state: state of charge and age. Combined with the long time horizon and short time step, this makes the use of global optimization such as dynamic programming rather costly.

3 The Adaptive Weight Algorithm (AWA)

Key Ideas

A first observation is that we can use a bi-level formulation to decompose the problem, using an aging constraint at the short term scale. We consider an age increment δ, and introduce the family of constraints

a(T ) ≤ a i + δ, (7)
and the corresponding 'micro' problems (P µ ) over one period T :

(P µ ) inf u T 0 (u(t))dt, (8)
subject to (4), ( 5) and [START_REF] Svoboda | Operating conditions of batteries in off-grid renewable energy systems[END_REF]. We denote by V µ (a i , δ) the associated short term value. The long term value V can be computed via dynamic programming according to

V (a i , kT ) = inf δ V µ (a i , δ) + V (a i + δ, (k + 1)T ). ( 9
)
Another key point is that aging is a slow process, therefore we can neglect the age variation over each micro problem. This reduces the number of state variables to 1 and allows for fast solving. However we still need to keep track of the age increment δ over the long time-frame.

We tackle this issue with a relaxation trick: given some weight α, we introduce the relaxed micro problems (P µ r )

(P µ r ) inf u T 0 [ (u(t)) + αF a (c(t), u(t), t)]dt, (10)
and denote V µ r (a i , α) the associated relaxed short term value. Now we present the link between the original micro problem and the relaxed micro problem, introducing some notations. Let us fix an initial age a i . For a control u over a period T , we denote by ∆(a i , u) the associated age increment and L(a i , u) the operational cost. Then for a

given weight α, we note ∆ r (a i , α) the age increment cost and L r (a i , α k ) the generator-associated in the relaxed micro problem [START_REF] Bonnans | Bocop -a collection of examples[END_REF]. Let us consider an optimal control ū for the micro problem [START_REF] Heymann | Mathematical contributions for the optimization and regulation of electricity production[END_REF], with the associated age increment ∆(a i , ū). Let us assume there exists a weight ᾱ leading to the same age increment in the relaxed micro problem, i.e. ∆ r (a i , ᾱ) = ∆(a i , ū). Then the generator-associated cost in the relaxed micro problem is actually the same as the cost of the original micro problem:

V µ (a i , ∆(a i , u)) = V µ r (a i , α) -α∆(a i , u) = L r (a i , α). (11)
This can be proven by noting that an optimal control for one problem is admissible for the other. We refer the reader to [START_REF] Heymann | Mathematical contributions for the optimization and regulation of electricity production[END_REF] for more technical details. Following this idea, our proposal consists in reducing the long term problem [START_REF] Heymann | A Stochastic Continuous Time Model for Microgrid Energy Management[END_REF] to

inf α k L r (a k , α k ) + ψ(a N T ), (12)
with (a k ) the sequence of ages defined by a k+1 = a k + ∆ r (a k , α k ). We can see (3.2) as a change of variable in (9): instead of optimizing over the one-period age increments, we optimize over the corresponding weights. Such a change of variable can be done safely in the absence of duality jumps. To our knowledge, there is no general condition for the absence of duality jumps, and we refer readers interested in a technical discussion on this decomposition (and the generality of the approach) to [START_REF] Heymann | Mathematical contributions for the optimization and regulation of electricity production[END_REF].

Unlike [START_REF] Bonnans | BocopHJB 1.0.1 -User Guide[END_REF], the new formulation (3.2) does not need to keep track of the age increment δ within the day. Therefore, the relaxed micro problems have only the state of charge as state variable, and can be solved efficiently by dynamic programming. The daily age variation can be computed simply by integrating its dynamics with the optimal control of the relaxed micro problems.

Algorithm

The adaptive weight algorithm thus solves an approximation of (6) by combining (a) slow/fast state variable separation, (b) constraint relaxation and (c) dynamic programming. The full resolution scheme is composed of an offline part (Algorithm 1) and an online part (Algorithm 2).

The offline algorithm computes and saves the values of the generator-associated cost L r and age increment ∆ r for a family of relaxed micro problems, over a discretized set of the initial age a i and weight α. For each relaxed micro problem, the dynamics use a constant age

a(t) = a i , ∀t ∈ [0, T ].
Note that all problems are independent, so this offline step could be massively parallelized.

The online algorithm uses these stored values to solve the long term problem by dynamic programming. Here the age a is the state and the weight α acts as the control. This part uses an interpolation of V , noted V .

Result: ∆ r , L r for (a i , α) in the discretization set do Compute an optimal control u for the relaxed micro problem with parameters (a i , α) ; ∆ r a i ,α ← ∆(u) and L r a i ,α ← L(u) ; end Algorithm 1: AWA algorithm: offline component Data: initial age a i , day index k 0 , ∆ r and L r (output of offline part) Result: optimal control for the next day u * V a,N ← ψ(a) for all a ; (we rescale from [0, 1] to [0, 1000]), u ∈ [-5, 15] and α ∈ [0, 250] {10 9 }. We test three levels of discretization: coarse, medium and fine, detailed in Table 1. The range for α was set by simple trial and error, since simulations only take a few seconds for the coarse discretization. We also set the parameters β = 0.5, ρ i = 0.95, ρ o = 0.95, K = 12.5, and the functions for the solar power P S (t) = max(0, 13 -0.3(4t -48) 2 ) and power load P L (t) = 3 + 3e -0.1(4t-32) 2 + 12e -0.03(4t-74) 2 .

for k ← N -1 to k 0 + 1 do for a in the discretization set do V a,k ← min α {L r α,a + V (a + ∆ r α,a , k + 1)} ; end end α * ← argmin α {L r α,a -α∆ r α,a i + V (a i + ∆ r α,a i , k 0 +
We also set ψ(a) = +∞ for a > 0.5.

Both the brute-force approach and the micro problems of the AWA method are solved via dynamic programming with the optimal control toolbox BocopHJB 1 ([9], [START_REF] Bonnans | Bocop -a collection of examples[END_REF]). The AWA algorithm itself is implemented in the R scripting language, keeping in mind that most of the

Results

We perform two batches of simulations, for a 0 = 0 and a 0 = 150. The objective values, recomputed along the optimal trajectories, as well as the CPU times (in seconds) are summarized in Tables 2 and3. The CPU time ratio is indicated, as well as the objective 'gap', i.e the relative error in objective compared to the brute force solution. As expected, the online part of the AWA method is extremely fast. It should be noted that the offline part is highly parallelizable, since each micro problem is independent from the others.

The increase in CPU time is consistent with the increased number in discretization steps.

Fig. 1 shows the age evolution for both methods with the fine discretization, with the limit a ≤ 500 touched for a 0 = 150. The values of the weights α are displayed on Fig. 2, and we see that the optimal weight is indeed not constant over time. Observe that the weights are higher when we start with an older battery, and the weight decreases for the a 0 = 0 trajectory, as the age limit becomes less stringent.

Finally, Fig. 3 illustrates the change of variable of the AWA: instead of optimizing over the one-period age increments ∆a, we optimize over the corresponding weights α. The graphs indicate that ∆a as a function of α is non-increasing, as expected. We observe that the age increments are lower for an old battery, indicating that the older and less efficient the battery becomes, the less it is used. On a side note, these graphs could be used to estimate some upper bound for the error in the method in case of duality jumps (see [START_REF] Heymann | Mathematical contributions for the optimization and regulation of electricity production[END_REF]). We take a larger initial age in order to illustrate the effect of the state constraint a ≤ 500. An infeasible state is in practice penalized by a very large value function. Due to the interpolations performed at each time step, some numerical diffusion of these large values occurs. This phenomenon tends to deform optimal trajectories that would come close to the boundary, as seen on Fig. 1 for the brute force (BF) method. The online part of AWA method uses a much larger time step (1 day versus 15 minute), which strongly reduces the diffusion and allows to compute a trajectory close to the boundary. Concerning Fig. 2, at the five last days the age reaches the limit a ≤ 500, and the corresponding weight α takes the fixed value 10 9 , which is not drawn on the graph. -Replacement of the Battery: In many applications, it is possible to buy spare parts to replace worn components. Hence we may want to introduce the possibility to buy a replacement battery in the optimization problem, probably with an impulse control.

-Periodicity: We can include seasonality by having different kinds of periods. For instance, we could model winter and summer days. In this case, one needs to perform an offline preprocessing for each kind of day. Likewise we could have week days and weekend, or even sunny and cloudy days (using a Markov Chains).

-Short Term Randomness: With Markovian dynamics and final constraints on the age increment we may want to see whether the same arguments apply on the daily average age increment.

-Infinite Horizon: If we add a discount rate, the arguments should apply for infinite horizon. Then one needs to replace the macro dynamic programing algorithm by either a policy iteration algorithm or a value iteration algorithm.

We also hope to see this approach tried on other application domains.

Conclusions

We present the adaptive weight dynamic programming algorithm (AWA), a decomposition technique for problems with periodic data, a fast state and a slow state. We illustrate the method on a toy micro grid problem, and observe that the trajectories and value functions from AWA are close to those from a brute-force approach (straightforward dynamic programming). The computing times are significantly smaller, the online part taking less than 1 second. Moreover, the offline part of the algorithm can be parallelized. We think this would make AWA suitable for a real-time optimization of the microgrid, as well as other dynamical systems with a slow state.

  1)}; Compute an optimal control u * for the relaxed micro problem with parameters (α * , a i ) ; Algorithm 2: AWA Algorithm: online component 4 Numerical Simulations 4.1 Parameters and Implementation We solve the problem with a time horizon of N = 600 days. The dynamic programming in dimension 1 and 2 both use a time step of 15 minutes, while the online part of AWA uses a time step of 1 day. For the state, control and weight we take the ranges SOC ∈ [0.1, 1], a ∈ [0, 500]

Fig. 1 Fig. 2

 12 Fig.1Age evolution for both methods -initial age: 0 and 150 -fine discretization.

4. 3

 3 Initial age a 0 = 0 For the fine discretization, the relative error in objective for AWA is only 0.5%, a factor 15 decrease of CPU time. We observe on Fig.1that the age profiles are almost identical, which confirms that AWA computes an accurate approximation of the solution.

Fig. 3

 3 Fig. 3 Daily age increments as function of weight α.

Table 1

 1 Discretizations for the brute force and AWA method.

			Brute force		AWA (offline) AWA (online)
	Grid	SOC Age Control SOC Control Age Weight
	coarse	25	25	20	25	20	25	10
	medium	50	50	40	50	40	50	25
	fine	100 100	100	100	100 100	50
	computational cost lies in solving the micro problems. All simulations are run on a standard
	laptop.

Table 2

 2 Objective values and CPU times (for AWA: offline+online) -initial age a 0 = 0

			Objective			CPU time (s)
	Grid	BF	AWA	gap	BF AWA (offline+online) BF/AWA
	coarse	121408 127271 4.8%	343	7.8+0.11	43.4
	medium 111326 107394 3.5%	2362	94+0.33	25.0
	fine	106179 106670 0.5% 22063	1432+0.87	15.4
	4.4 Initial age a 0 = 150

Table 3

 3 Objective values and CPU times (for AWA: offline+online) -initial age a 0 = 150.

			Objective			CPU time (s)
	Grid	BF	AWA	gap	BF	AWA BF/AWA
	coarse	156474 165655 5.9%	330	8.0+0.12	40.6
	medium 136942 125234 8.5%	2371	90+0.29	26.3
	fine	126980 122520 3.5% 22289 1477+0.78	15.1
	5 Perspectives and Open Problems
	We have identified several follow-up for this work. Increasing the degree of sophistication of the
	model raises modeling and mathematical questions. Possible extensions include:

http://www.bocop.org

Acknowledgements The first author BH acknowledges support from the Siebel Scholar Program. The authors wish to thank Frédéric Bonnans for his advice and comments on the paper. Thanks also to Hanyun Huang for her English suggestions.


Yoann Pascal
email: yoann.pascal@satie.ens-cachan.fr

Denis Labrousse

Mickaël Petit

Stéphane Lefebvre

François Costa

PCB-Embedding of Power Dies Using Pressed Metal Foam

This paper presents a new technique to manufacture the top-side contact of a PCBembedded power die, using a pressed piece of metal foam. The manufacturing process is detailed, prototypes are manufactured and electrically characterized, and the technique is used to embed a switching cell with a 1.7 nH switching loop. This process is easy to make and highly cost-effective; the electric resistance of the assembly is close to that of a state-of-the-art industrial package using bond wires, moreover the stray inductance is minimized.

Introduction

Aluminum bond wires are used in most power discrete components and modules to connect the top-side of the chips with the leads. This process is cost-effective and mature but could still be improved. As a matter of fact, the stray inductance of a wire bonding is rather high (about 10 nH in a TO247 package) and although paralleling bond wires can decrease the overall interconnect inductance, it can hardly be decreased to less than a few nanohenries due to magnetic couplings between close bond wires. Furthermore, the limitations in term of reliability of wire bonds often make them the weak links of packages and modules [START_REF] Ménager | Industrial and lab-scale power module technologies: A review[END_REF]. Additionally, when wire bonds are used, the heat generated within the die must be extracted nearly exclusively by the bottom-side of the die, resulting in poor thermal performances. Top-side die interconnection is therefore a highly topical field of research, investigated by many industrial and academic researchers [START_REF]Current Developments in 3D Power Packaging with Focus on Embedded Substrate Technologies[END_REF][START_REF] Wang | Review of State-of-the-Art Integration Technologies in Power Electronic Systems[END_REF] willing to develop processes with better performances than wire bonds. The resistance of aluminum wire bonds can be decreased using copper wire [START_REF] Siepe | The Future of Wire Bonding is? Wire Bonding[END_REF], though the process is more complex and still suffers from reliability issues. Replacing wire bonds by copper ribbons or claps can decrease the resistance and the inductance of the connection whilst improving the thermal performances [START_REF] Rittner | Robust top side contact technology on power semiconductors --results from the public funded project 'propower[END_REF]. In [START_REF] Stockmeier | SKiN: double side sintering technology for new packages[END_REF], a flexible Printed Circuit Board (PCB) is sintered on top of the dies, thereby improving the electrical performances and allowing the control circuit to be located at a few millimeters of the power die. Press-pack-like packages can be manufactured using either spring-like conductive wool [START_REF] Zhu | A Solution to Press-Pack Packaging of SiC MOSFETS[END_REF] or plains blocks of copper or molybdenum [START_REF] Vagnon | Study and Realization of a Low Force 3D press-pack power module[END_REF][START_REF] Vagnon | A bus-bar-like power module based on three-dimensional power-chip-onchip hybrid integration[END_REF]. Embedding the dies within a PCB is a solution expected to solve most of bond wires' drawbacks [START_REF] Wang | Review of State-of-the-Art Integration Technologies in Power Electronic Systems[END_REF]. Stray inductances should be greatly reduced thanks to smaller current loops and wide conductive tracks while the resistance should be decreased due to the possibility to use thick copper layers on both sides of the die. Doubleside cooling, enabled by this technique, should greatly improve thermal performances. In the most common embedding technique [START_REF] Gottwald | P2 Packthe paradigm shift in interconnect technology[END_REF][START_REF] Kearney | PCB-embedded semiconductors for low-voltage power electronic applications[END_REF][START_REF] Regnat | Optimized power modules for silicon carbide MOSFET[END_REF][START_REF] Manessis | Development of advanced power modules for electric vehicle applications[END_REF][START_REF] Wyss | analysis of pcb embedded power semiconductors for a 30 kw boost pfc converter[END_REF][START_REF] Sharma | PCB embedded power package with reinforced top-side chip contacts[END_REF], a layer of copper is electrodeposited on the top-side of an embedded power-die. However, dies top-side metallizations are usually made of aluminum and therefore not directly compatible with copper electrodeposition; an interface layer must be deposited beforehand in a sub-process that is not available in standard PCBmanufacturing facilities. The technique is thus rather complex and costly. It is, nonetheless, used in mass production [START_REF]Current Developments in 3D Power Packaging with Focus on Embedded Substrate Technologies[END_REF]. Some investigations on the reliability of embedded power dies have been published [START_REF] Kearney | PCB-embedded semiconductors for low-voltage power electronic applications[END_REF][START_REF] Manessis | Technical Understanding of Resin-Coated-Copper (RCC) Lamination Processes for Realization of Reliable Chip Embedding Technologies[END_REF]. In this paper, an original process for die embedding using a pressed contact made of metal-foam is under study. A simplified version of the embedding process, suited to diodes, was detailed in [START_REF] Pascal | Using Laminated Metal Foam as the Top-Side Contact of a PCB-Embedded Power Die[END_REF] and a few samples were electrically characterized. This paper considers an improved version of the process, fit for transistor embedding; it is detailed in section 2. Samples of embedded diodes are electrically characterized in section 3. A switching cell is embedded and tested in section 4. Section 5 concludes the paper.

Embedding Technique

The manufacturing process, quite straightforward, consists of only three main steps. First, the bottom side of the die is soldered on a "bottom PCB" (Fig. 1.a) using a standard soldering paste. A stack is then prepared; it is composed of (Fig. 1.b): -the bottom PCB (and the die attached to it); -layers of prepreg; -a kapton mask, with a hole above each pad of the die top-side pads; -a piece of metal foam for each top-side pad; -the top PCB. This stack is then laminated (Fig. 1.c) using the temperature and pressure profiles required by the prepreg. The kapton mask must be carefully positioned above the die, i.e. the holes in the mask must be precisely aligned with the pads. Then, pieces of foam far larger than the pads can be used without risking short-circuiting adjacent pads or guard rings since the mask is interposed between the foam and those pads (cf. Fig. 1.c). This simplifies the process insofar as, until lamination, the pieces of metal foam are not attached to the rest of the circuit -although they are roughly located above the kapton mask by the top prepreg layer (see Fig. 1.b). Although using this mask allows connecting rather small (< 400 µm) and close pads, such as the gates of transistors, the process is a priori to be used with large power dies. The copper planes located on each side of the resulting assembly can then be used as baseplate for a heat-sink or as a regular PCB -used for instance to build the control circuit of the embedded power dies.

The solder melting temperature should be high enough so that it does not liquefy during lamination. This is yet to be tested but the die could likely be sintered on to the bottom PCB.

Electrical Characterization of the Contact

Introduction

As a first step, vertical power diodes with dies dimensions of 4.2 mm × 6.5 mm × 254 µm were used. These diodes were available both as bare dies and in a state-of-the art PQFN (Power Quad Flat No-lead package) package using bond wires. The dies were embedded using our technique -a metallographic cut of one of our assemblies is given in Fig. 2. Prototypes were electrically characterized and compared with those obtained using the PQFN package, used as reference.

Nickel foam Ni5763, from Recemat BV was used [START_REF] Recemat | Data sheet Nickel Foam[END_REF]. The mean pore diameter is 350 µm. The thickness and porosity, respectively equal to 1.4 mm and 96 % before pressing, fall to a few hundreds of micro-meters and a few tens of percent after lamination -depending upon the applied pressure. Panasonics' prepreg R-1650V/7629 [START_REF] Panasonic | R1755-V & R1655-V -Datasheet[END_REF] was used; its low CTE (11 to 15 ppm/K in the x-and yaxes) is similar to that of aluminum (24 ppm/K), copper (17 ppm/K), nickel (13 ppm/K) and not very different from that of silicon (2.5 ppm/K). This should maximize the reliability of the assemblies. The top and bottom PCBs were manufactured using a standard 400 µm-tick, FR4, epoxy-based, substrate.

Forward & Reverse Characteristics

The static characteristics of our prototypes and the wire-bonded reference were measured using a Kelvin-connected Keysight B1505a power device analyzer. The diode characteristics were modeled by a constant voltage source and a series resistance (Fig. 3) accounting for the resistance of the package and the die itself. The distribution of the resistances of 59 prototypes and that of the reference are shown in Fig. 4. The PCB-embedded dies are from three different wafers from the same batch; the wirebonded die is from another batch.

The mean resistance of the prototypes of the main mode of the histogram (i.e. with resistances between 2 mΩ and 6 mΩ) was 4.1 mΩ, which is very close to the resistance of the reference (3.3 mΩ at 𝐼 𝑓 = 10 A). This mode incorporated 92 % of the samples and its standard deviation was 470 µΩ. These results show that, despite a few failed devices, the process is highly repeatable and has a high production yield.

Most samples had a breakdown voltage (measured at 𝐼 𝑓 = -1 mA) greater than 75 V (the datasheet specifies 65 V) and leakage currents similar to those measured on the reference. This proves that the blocking performances of the diodes were not noticeably affected by the embedding process.

AC-Dynamic Impedance

Nickel foam was used in our prototypes. Because of the magnetic properties of this metal, increased eddy-currents, and therefore higher losses, might be expected at high frequencies. The AC-dynamic impedance 𝑍 𝐷 (impedance vs. frequency) of one prototype and the reference were measured and compared, using the set-up illustrated in Fig. 5; the impedance is then given by (1):

𝑍 𝐷 = 𝑍 0 2 ⋅ Λ 1 -Λ , Λ = 𝑉 𝑇 𝑉 𝑅 , 𝑍 0 = 50 Ω (1)
The results (Fig. 6) show that the curves are flat up to at least 1 MHz. Then, despite the magnetic properties of nickel, no high frequency loss is to be expected in the foam.

On the other hand, the current distribution should mimic the complex and random structure of the foam which may, as a result, act as Litz wire and thus be less lossy than a plain block of metal [START_REF] Zhu | A Solution to Press-Pack Packaging of SiC MOSFETS[END_REF]. Using the technique detailed above, a switching cell composed of a PIN diode (5SLY 12M1700, 300 A, 1.7 kV, from ABB) and an IGBT (5SMY 12M1280, 150 A, 1.2 kV, from ABB) was built. Although the selected semiconductors were widely overrated vs. the intended use, they had been selected for their layout which made them easy to embed: large enough (13 mm ×13 mm) to be easily handled with bare hands; the emitter connection is not divided into several pad; the gate pad is large (1.4 mm×1.4 mm) and located in a corner of the die. A NCP81074A driver and a 4.7 Ω gate resistor were located about 15 mm away from the IGBT -though they were not PCB-embedded. The source was kelvin-connected, so as to minimize the coupling between the gate drive and the power circuits.

A 68 nF-X7R SMD capacitor (𝐶 1 ) located 5 mm away from the IGBT and the diode and not PCBembedded was used to bypass the power bus. It was supported by a 250 µF MKP bulk capacitor (𝐶 2 ) located 7 cm away from the switching cell. Fig. 7 shows the schematic of the circuit and its layout. A cross-section of the switching cell is given in Fig. 8. A photograph of the prototype is given in Fig. 9. Using pressed foam to connect the dies top-sides easily accommodates with dies having different thicknesses; the foam can be more or less pressed, and therefore exhibit different densities, depending upon its thickness. In this study, for instance, the diode was 390 µmthick whereas the IGBT was 140 µm-thick. The foam thickness was therefore 120 µm above the diode and 370 µm above the IGBT. The embedded part of the circuit (in green in Fig. 7) was 1.5 mm-thick. The overall thickness of the circuit (heatsink and distant bulk capacitor excluded) was 2.8 mm.

Switching Waveforms

The thermal properties of foam-based pressed contacts are yet to be fully characterized. The embedded switching cell was therefore submitted to a double-pulse test so as to test it at full load while avoiding over-heating the dies.

The prototype was submitted to a double-pulse test with a 10 mH-load under a supply voltage of 600 V and a load current of 40 A. In these conditions, the rise time was 170 ns (i.e. 3.5 kV/µs) and the fall-time was 580 ns (i.e. 1 kV/µs). Given the ratings of the dies, these transitions could be considered as fast.

Nonetheless, there was little ringing (only 8 % overshoot).

The voltage across the IGBT and the current in the inductive load, measured using a HDO8104 oscilloscope, a PPE2kV high-voltage passive probe and a CP031 current probe are given in Fig. 10 & Fig. 11.

Switching Loops Inductances Identification of the Loops

The electrical schematic of the cell, given in Fig. 7, shows three loops:

loop n°1, including the dies and the bypass capacitor 𝐶 1 ; -loop n°2, including both bypass capacitors; -loop n°3, including the bypass capacitor 𝐶 2 and the power source. When excited by the switching cell, each of these loops might resonate, leading to oscillations and overvoltages across capacitors and on 𝑉 𝑠𝑤 .

Usually, 𝐶 1 is in much closer proximity to the dies than 𝐶 2 whereas the power source is far away from the converter, leading to:

𝐿 1 ≪ 𝐿 2 ≪ 𝐿 𝑃𝑆𝑈 (2)
where 𝐿 1,2 represent the inductance of the layout. Furthermore, the Equivalent Series Inductance (ESL) of 𝐶 2 (represented by 𝐿 𝐶2 in Fig. 7) is generally much greater than that of 𝐶 1 (𝐿 𝐶1 ).

On the other hand, the parasitic capacitance 𝐶 𝑠𝑤 of the switching node is mainly due of the output capacitance of the IGBT and:

𝐶 𝑠𝑤 ≪ 𝐶 1 ≪ 𝐶2 (3)
As a result of this, the parallel resonant pulsation of each loop can be estimated as follow:

{ 𝜔 1 ≈ (𝐶 𝑠𝑤 (𝐿 1 + 𝐿 𝐶1 )) -0.5 𝑓𝑜𝑟 𝑙𝑜𝑜𝑝 𝑛°1 𝜔 2 ≈ (𝐶 1 (𝐿 2 + 𝐿 𝐶1 + 𝐿 𝐶2 )) -0.5 ≈ (𝐶 1 (𝐿 2 + 𝐿 𝐶2 )) -0.5 𝑓𝑜𝑟 𝑙𝑜𝑜𝑝 𝑛°2 𝜔 3 ≈ (𝐶 2 (𝐿 𝑃𝑆𝑈 + 𝐿 𝐶2 )) -0.5 𝑓𝑜𝑟 𝑙𝑜𝑜𝑝 𝑛°3 (4)
Estimation of the of the Loop inductances using Q3D Simulation

The primary loop inductance was simulated using the Q3D-Ansys software, based on the PEEC method. A 3D-view of the model is given in Fig. 12 (the top-PCB is removed for clarity). The pieces of metal foam were modeled by homogenous blocks of material having the mean conductivity and permeability of the pressed foam. The bypass capacitor was modeled by a plain block of metal (𝜎 = 26 kS/m, so that the resistance had the same value as that of the real device). The loop was excited by a current source located in the symmetry plane of the capacitor (cf. Fig. 12). The mesh was automatically refined by the software.

Simulation results are given in Fig. 13; the lowand high-frequency inductances are, respectively, 23 nH and 1.7 nH.

Measurement of the Loop Inductances

The inductance of loop n°1 was measured using the set-up schematized in Fig. 14. The IGBT was turned on and both dies were forward biased using 3 A current sources, each built using a voltage source series-connected with a 10 Ωresistor, a 60 µH-inductance (with resonant frequency 6 MHz) and a common-mode choke. The bypass capacitor 𝐶 1 was removed and a VNA E5061b was connected to its footprints; the scattering parameter 𝑆 11 was then used to obtain the impedance of the loop. It is shown in Fig. 15, together with a model fitted on the measurements. Then:

𝐿 1 + 𝐿 𝐶1 = 1.7 nH (5)
This measured value is very close to the simulated one: 1.8 nH at 100 MHz (Fig. 13). It was assumed that a power source far away from the converter was used; then: 𝐿 𝑃𝑆𝑈 ~∞.

The model (𝐶 2 , 𝐿 𝐶2 ) of the bulk capacitor was measured using an E4990a impedance analyzer:

𝐿 𝐶2 = 19 nH, 𝐶 2 = 271 µF (6)
Similarly, we measured 𝐶 1 = 34 nF (at 600 VDC biais). Last, we measured 𝐶 𝑠𝑤 = 200 pF (at 𝑉 𝑔𝑒 = 0, 𝑉 𝑐𝑒 = 600 V).

The overall secondary loop inductance 𝐿 2 + 𝐿 𝐶2 + 𝐿 𝐶1 ≈ 𝐿 2 + 𝐿 𝐶2 was measured using an E4990a:

the IGBT was turned-off and the impedance analyser was connected on the footprints of 𝐶 1 (which had been removed). From this measurement and (6) yields: 𝐿 2 = 16 nH.

The results of the measurements are gathered in Table 1.

Device

Ringing Frequencies

Using (4) and the results from the previous section, we estimate the resonant frequencies of the loops:

𝑓 1 = 𝜔 1 2𝜋 ≈ 273 MHz, 𝑓 2 = 𝜔 2 2𝜋 = 4.6 MHz (7)
The rise time during the double pulse test (cf. Fig. 11) is 170 ns -corresponding to a 𝑓 𝑒𝑞,𝑡𝑟 = 1.9 MHz-bandwidth. Since 𝑓 2 ≈ 𝑓 𝑒𝑞,𝑡𝑟 ≪ 𝑓 1 , the step is too slow to excite loop n°1 but excites loop n°2. This is confirmed by the measured ringing frequency (4.5 MHz, Fig. 11) being close to the self-resonant frequency of loop n°2 (4.6 MHz, eq. ( 7)) [START_REF] Meisser | Parasitics in Power Electronic Modules: How parasitic inductance influences switching and how it can be minimized[END_REF].

Conclusion on the Switching Loops

In the converter presented in this paper, switching loop n°1 is not excited by the steps in 𝑉 𝑠𝑤 but loop n°2 is. This shows that 𝐶 1 has limited impact in absorbing the 𝑉 𝑠𝑤 -overvoltage. This capacitor would, however, be most useful in converters using faster switches.

Discussion

The small distance between the driver and the IGBT, along with the planar layout of the gate circuit accounted for the very small gate loop inductance, resulting in very small ringing on the gate voltage 𝑉 𝑔𝑒 -even for total gate resistances as low as 1.4 Ω.

The proper separation of the power and gate drive emitter circuits, operated on the layout, accounted for the low coupling between the two circuits. This is of outmost importance to avoid unwanted 𝑉 𝑔𝑒 overvoltages and IGBT turn-on due to fast transients.

The switching loop overall inductance was measured to be about 1.7 nH. This very low value could be attained mostly thanks to the very small height of the pressed contact, resulting in a quasiplanar circuit. The absence of bond wires made possible placing the components very close to one another, resulting in a low loop inductance and a high power density.

Conclusion & Perspectives

This paper presents an original packaging method for dies in the field of power electronics using pressed metal foam as the top-side contact.

Measurements showed that the electrical performances of the realized prototypes were very similar to those of a conventional state-ofthe-art package with wire-bonded chip.

As a proof of concept, a switching cell was also embedded using our process, resulting in a 1.7 nH switching-loop inductance. It was successfully tested in a 24 kVA-double-pulse test.

Ongoing work is investigating the reliability of the process and its stability with respect to temperature, humidity and time. The thermal properties of the pressed contact are also under study.

Fig. 1 .

 1 Fig. 1. Embedding of a MOSFET using pressed metal foam. (a): the die is soldered on the bottom PCB; (b): prepreg layers, a kapton mask and two pieces of metal foam (one for each top-side pad) are positioned; (c): the stack is laminated, the epoxy resin polymerizes and floods the foam.

Fig. 3 .

 3 Fig. 3. Diode current versus voltage drop: measured values of the reference and of some embedded prototypes and steady-state models (calculated at 𝐼 𝑓 = 20 A).

Fig. 4 .

 4 Fig. 4. Distribution of the series resistances of our 59 prototypes compared with that of the reference (calculated at 𝐼 𝑓 = 10 A).

Fig. 2 .

 2 Fig. 2. Metalurgical cut of an embedded diode. The die and the foam are, respectively, 254 and 220 µmthick.

Fig. 5 .Fig. 6 .

 56 Fig. 5. Setup used to measure the dynamic impedance of a diode

Fig. 8 .Fig. 9 .

 89 Fig. 8. Schematic of a cross-section of the switching cell. Not to scale.

Fig. 7 .

 7 Fig. 7. Schematic (left) and layout (right) of the converter showing some parasitic inductances and which components are PCB-embedded. On the layout, blue zones indicates copper above the dies, red zones are for copper located bellow them.

Fig. 10 .

 10 Fig. 10. Waveforms measured in the converter during a double-pulse test: Load current 𝐼 𝑙𝑑 , voltage at the switching node 𝑉 𝑠𝑤 and gate to emetter voltage 𝑉 𝑔𝑒 .

Fig. 11 .

 11 Fig. 11. Waveforms measured in the converter: Load current 𝐼 𝑙𝑑 , voltage att the switching node 𝑉 𝑠𝑤 and gate to emetter voltage 𝑉 𝑔𝑒 . Focus on the commutations.

Fig. 12 .

 12 Fig. 12. 3D-model of loop n°1 (top-PCB removed).

Fig. 13 .

 13 Fig. 13. Simulation results: inductance of loop n°1.

Fig. 14 .Fig. 15 .

 1415 Fig. 14. Set-up used to measure the inductance of loop n°1.

Table 1 :

 1 Measured values of the elements of the model of Fig.7.

		Value	Device	Value
	𝐿 1 + 𝐿 𝐶1	1.7 nH	𝐶 𝑠𝑤	200 pF
	𝐿 2	16 nH	𝐶 1	34 nF
	𝐿 𝐶2	19 nH	𝐶 2	271 µF
	𝐿 𝑃𝑆𝑈	∞

Acknowledgements

This work was supported by a grant overseen by the French National Research Agency (ANR-15-CE05-0010).


Vasco M Manquinho

João Marques-Silva

Effective Lower Bounding Techniques for Pseudo-Boolean Optimization

Linear Pseudo-Boolean Optimization (PBO) is a widely used modeling framework in Electronic Design Automation (EDA). Due to significant advances in Boolean Satisfiability (SAT), new algorithms for PBO have emerged, which are effective on highly constrained instances. However, these algorithms fail to handle effectively the information provided by the cost function of PBO. This paper addresses the integration of lower bound estimation methods with SAT-related techniques in PBO solvers. Moreover, the paper shows that the utilization of lower bound estimates can dramatically improve the overall performance of PBO solvers for most existing benchmarks from EDA.

Introduction

Recent advances in Boolean Satisfiability (SAT) have resulted in new practical algorithms for solving the Linear Pseudo-Boolean Optimization problem [START_REF] Aloul | Generic ILP versus specialized 0-1 ILP: An update[END_REF][START_REF] Chai | A Fast Pseudo-Boolean Constraint Solver[END_REF][START_REF] Dixon | Inference Methods for a Pseudo-Boolean Satisfiability Solver[END_REF]. These algorithms perform a linear search on the possible values of the cost function, starting from the highest, at each step requiring the next computed solution to have a cost lower than the previous one. If the resulting instance is not satisfiable, then the optimal value is given by the last computed solution. By incorporating important features from SAT solvers, like non-chronological backtracking in the search tree, conflict-based learning mechanisms and lazy data structures, these solvers have been able to solve with success several classes of highly constrained pseudo-boolean instances. However, these solvers also fail to handle effectively the information provided by the cost function, which can be crucial as the experimental results in this paper clearly demonstrate. In order to prune the search due to the value of the cost function we propose the use of methods to estimate a lower bound on the value of the cost function. Whenever the lower bound estimation is higher or equal to the best solution found so far, we are able to prune the search tree. Moreover, we also establish conditions for backtracking non-chronologically in the search tree when the search backtracks due to the lower bound estimate. This paper extends [START_REF] Manquinho | Search pruning techniques in sat-based branch-and-bound algorithms for the binate covering problem[END_REF] by considering alternative lower bounding procedures.

The paper is organized as follows. The first part concentrates on describing how different lower bound estimation methods can be applied in pseudo-boolean optimization problem, focusing on linear-programming relaxation and Lagrangian relaxation. Section 4 proposes techniques for obtaining explanations that allow backtracking non-chronologically when the search is bound due to the lower bound estimate. In addition, the paper presents a number of additional techniques that have proved useful in practical PBO solving. Preliminary experimental results are analyzed in section 6 for most existing PBO benchmarks for EDA problems. Finally, the paper concludes in section 7.

Preliminaries

In a propositional formula, a literal l j denotes either a variable x j or its complement xj . If a literal l j = x j and x j is assigned value 1 or l j = xj and x j is assigned value 0, then the literal is said to be true. Otherwise, the literal is said to be false.

An instance P of a Linear Pseudo-Boolean Optimization problem can be defined as follows,

minimize j∈N c j • x j subject to j∈N a ij l j ≥ b i , x j ∈ {0, 1}, a ij , b i ∈ N + 0 , i ∈ M, N = {1, . . . , n}, M = {1, . . . , m} (1)
where c j is a non-negative integer cost associated with variable x j , j ∈ N and a ij denote the coefficients of the literals l j in the set of m linear constraints. Every pseudo-boolean formulation can be rewritten such that all coefficients a ij and right-hand side b i be non-negative.

In a given constraint, if all a ij coefficients have the same value k, then it is called a cardinality constraint, since it only requires that b i /k literals be true. A pseudo-boolean constraint where any literal set to true is enough to satisfy the constraint, can be interpreted as a propositional clause. This occurs when the value of all a ij coefficients are greater than or equal to b i .

If every constraint can be interpreted as a propositional clause then P is an instance of the binate covering problem (BCP). Covering formulations have been the subject of thorough research work that can be found in [START_REF] Coudert | On Solving Covering Problems[END_REF][START_REF] Manquinho | Search pruning techniques in sat-based branch-and-bound algorithms for the binate covering problem[END_REF][START_REF] Villa | Explicit and Implicit Algorithms for Binate Covering Problems[END_REF].

Notice that a linear pseudo-boolean optimization problem can also be viewed as a special case of linear integer programming problem. The linear integer programming formulation for the constraints can be obtained if we replace literals xj by 1 -x j . In section 3 we will use this latter formulation.

Pseudo-Boolean Optimization Algorithms

In [START_REF] Barth | A Davis-Putnam Enumeration Algorithm for Linear Pseudo-Boolean Optimization[END_REF], P. Barth first proposed an approach based on Boolean Satisfiability (SAT) techniques for solving Pseudo-Boolean Optimization (PBO). This approach consists of performing a linear search on the possible values of the cost function, starting from the highest, at each step requiring the next computed solution to have a cost lower than the previous one. If the resulting instance is not satisfiable, then the solution is given by the last recorded solution. The generalization of recent advances in SAT resulted in new successful algorithms [START_REF] Aloul | Generic ILP versus specialized 0-1 ILP: An update[END_REF][START_REF] Chai | A Fast Pseudo-Boolean Constraint Solver[END_REF][START_REF] Dixon | Inference Methods for a Pseudo-Boolean Satisfiability Solver[END_REF] for several sets of PBO instances, namely the incorporation of non-chronological backtracking in the search tree, conflict-based learning mechanisms and lazy data structures have been applied with success. The SAT-based approach focuses primarily on finding solutions for the problem constraints. Therefore, for highly constrained problems these techniques are very effective. However, these algorithms find it difficult to deal with the information from the cost function.

Unlike the SAT-based approach, branch-and-bound algorithms [START_REF] Coudert | On Solving Covering Problems[END_REF][START_REF] Liao | Solving Covering Problems Using LPR-Based Lower Bounds[END_REF] have proved to be very effective when the instances to be solved are not highly constrained since they are able to prune the search tree earlier due to estimate of the value of the cost function. In branch-and-bound algorithms upper bounds on the value of the cost function are identified for each solution to the constraints, and lower bounds on the value of the cost function are estimated considering the current set of variable assignments. For a given instance P of a pseudo-boolean optimization problem, let P.upper denote the upper bound on the value of the cost function. The search is pruned whenever the lower bound estimation is higher than or equal to P.upper. In this case it is guaranteed that a better solution cannot be found with the current variable assignments and therefore the search can be pruned. The algorithms described in [START_REF] Coudert | On Solving Covering Problems[END_REF][START_REF] Liao | Solving Covering Problems Using LPR-Based Lower Bounds[END_REF][START_REF] Manquinho | Search pruning techniques in sat-based branch-and-bound algorithms for the binate covering problem[END_REF][START_REF] Villa | Explicit and Implicit Algorithms for Binate Covering Problems[END_REF] for the binate covering problem follow this approach as well as several general integer programming solvers.

For several instances, specially for low constrained instances, the tightness of the lower bounding procedure is crucial for the algorithm's efficiency, because with higher estimates of the lower bound, the search can be pruned earlier. Several procedures can be used for lower bound estimation, namely the approximation of a maximum independent set of constraints (MIS) [START_REF] Coudert | On Solving Covering Problems[END_REF][START_REF] Manquinho | Search pruning techniques in sat-based branch-and-bound algorithms for the binate covering problem[END_REF], linear-programming relaxations [START_REF] Liao | Solving Covering Problems Using LPR-Based Lower Bounds[END_REF] or Lagrangian relaxations [START_REF] Ahuja | Network Flows: Theory, Algorithms, and Applications[END_REF].

Linear Programming Relaxations

Although the approximation of a maximum independent set of constraints (MIS) is the most widely used lower bound procedure for the binate covering problem (a particular case of PBO) [START_REF] Coudert | On Solving Covering Problems[END_REF][START_REF] Villa | Explicit and Implicit Algorithms for Binate Covering Problems[END_REF], linear programming relaxation (LPR) has also been used with success [START_REF] Liao | Solving Covering Problems Using LPR-Based Lower Bounds[END_REF] . It is also often the case that the linear programming relaxation bound is higher than the one obtained with the MIS approach. Moreover, linear programming relaxations have long been used as a lower bound estimation procedure in branch-and-bound algorithms for solving integer programming problems [START_REF] Nemhauser | Integer and Combinatorial Optimization[END_REF].

The general formulation of the LPR for a pseudoboolean problem is obtained from (1) as follows:

minimize z lpr = cx subject to Ax ≥ b 0 ≤ x ≤ 1 (2)
where vector c defines the non-negative integer cost associated with every decision variable in vector x. Entries of matrix A defines the constraint coefficients and vector b the right-hand side of every constraint. The solution of ( 1) is referred to as z * cp , whereas the solution of ( 2) is referred to as z * lpr .

It is well-known that the solution z * lpr of ( 2) is a lower bound on the solution z * cp of (1) [START_REF] Nemhauser | Integer and Combinatorial Optimization[END_REF]. Basically, any solution of ( 1) is also a feasible solution of ( 2), but the converse is not true. Moreover, for a given solution of (2) where x ∈ {0, 1} n , we necessarily have z * cp = z * lpr . Hence, the result follows. Furthermore, different linear programming algorithms can be used for solving [START_REF] Aloul | Generic ILP versus specialized 0-1 ILP: An update[END_REF], some with guaranteed worst-case polynomial run time [START_REF] Nemhauser | Integer and Combinatorial Optimization[END_REF].

Lagrangian Relaxations

Lagrangian relaxation (LGR) is a widely used method for computing bounds on the optimal value of the cost function from network optimization to nonlinear programming [START_REF] Ahuja | Network Flows: Theory, Algorithms, and Applications[END_REF][START_REF] Nash | Linear and Nonlinear Programming[END_REF]. It also known that for some instances, the bound provided by the Lagrangian relaxation method is tighter than the one obtained by the linear programming relaxation [START_REF] Ahuja | Network Flows: Theory, Algorithms, and Applications[END_REF]. Therefore, Lagrangian relaxation can be used to provide a quick and tight lower bound on the value of the cost function for pseudo-boolean optimization problems.

While in linear programming relaxations we are able to find a lower bound estimate by solving the problem constraints and relaxing the possible variable values, in Lagrangian relaxations we relax the problem constraints and incorporate them in the objective function with associated Lagrangian multipliers.

Given a generic linear optimization problem formulated as:

minimize

z * = cx subject to Ax = b x ∈ X (3)
we can define the Lagrangian function L(µ) as:

L(µ) = min{cx + µ(Ax -b) : x ∈ X} (4)
where vector µ defines the Lagrangian multiplier associated with each constraint. The Lagrangian Bounding Principle [START_REF] Ahuja | Network Flows: Theory, Algorithms, and Applications[END_REF] states that for any vector µ of the Lagrangian multipliers, the value of L(µ) is a lower bound on the optimal solution of the original optimization problem. In (3) all constraints are formulated as equalities, while in the pseudo-boolean optimization problem [START_REF]CPLEX MILP solver[END_REF] we have inequality constraints. Therefore, in that case the Lagrangian relaxation problem is formulated as:

L * = max{L(µ) : µ ≥ 0} (5)
where L * is the optimum value of the Lagrangian relaxation. The most tight lower bound estimate we can obtain using this method is given by L * . Before trying to solve the Lagrangian relaxation problem in order to obtain L * , we must determine the value of L(µ) for a given value of µ. Notice that by expanding (4) we have:

L(µ) = min{ j∈N cjxj + i∈M µi(( j∈N aijxj) -bi)} L(µ) = min{ j∈N αjxj - i∈M µibi} where αj = (cj + i∈M µiaij) (6)
In order to obtain the value of L(µ), for a given µ, we must determine the value of the decision variables x j to be able to minimize the expression. Therefore, we must have x j = 0 whenever α j ≥ 0 and x j = 1 when α j < 0.

Literature from nonlinear programming [START_REF] Nash | Linear and Nonlinear Programming[END_REF] and network optimization [START_REF] Ahuja | Network Flows: Theory, Algorithms, and Applications[END_REF] provide methods to solve [START_REF] Coudert | On Solving Covering Problems[END_REF], namely gradient methods to approximate the value of L * . For the results presented in this paper the approach outline in [START_REF] Ahuja | Network Flows: Theory, Algorithms, and Applications[END_REF] is used.

Bound-based Conflicts

Conditions for backtracking non-chronologically due to the lower bound estimate on the value of the cost function first proposed in [START_REF] Manquinho | Search pruning techniques in sat-based branch-and-bound algorithms for the binate covering problem[END_REF] assume that the lower bounding procedure is the MIS approximation. Next, we review the main ideas about pruning the search tree based on the estimated value of the cost function and describe the conditions to apply when using linear-programming relaxation or Lagrangian relaxation as a lower bound estimation procedure.

Backtracking on Bound-based Conflicts

A bound conflict in an instance of the pseudo-boolean optimization problem (PBO) P arises when the lower bound is equal to or higher than the upper bound. This condition can be written as P.path + P.lower ≥ P.upper [START_REF] Hooker | Logic-Based Methods for Optimization[END_REF] where P.path is the cost of the assignments already made, P.lower is a lower bound estimate on the cost of satisfying the constraints not yet satisfied (as given for example using Lagrangian relaxation), and P.upper is the best solution found so far.

In this situation, our approach is to identify a set of assignments responsible for the bound conflict and build a new propositional clause ω bc such that it prevents those assignments of being repeated during the search process. When ω bc is added, a conflict analysis procedure must be carried out to determine to which level of the search tree to backtrack to.

A straightforward approach to build ω bc would be to consider the decision variable assignments from all levels of the search tree, but in that case the resulting backtrack would necessarily be chronological. In [START_REF] Manquinho | Search pruning techniques in sat-based branch-and-bound algorithms for the binate covering problem[END_REF] it was already shown that the assignments responsible for the bound conflict might not be associated with all levels of the search tree.

From [START_REF] Hooker | Logic-Based Methods for Optimization[END_REF], we can readily conclude that P.path and P.lower are the unique components involved in each bound conflict. Therefore, we will analyze both the P.path and P.lower components in order to establish the assignments responsible for a given bound conflict. Our goal is to define two sets of literals ω pp and ω pl containing the explanation for P.path and P.lower, respectively. Our bound conflict clause ω bc is defined by the set union of the literals in ω pp and ω pl .

We start by studying P.path. Clearly, the variable assignments that cause the value of P.path to grow are solely those assignments with a value of 1. Hence, we can define ω pp such that each variable in ω pp has positive cost and is assigned value 1:

ω pp = {l = xj : Cost(x j ) > 0 ∧ x j = 1} (8)
which basically states that in order to decrease the value of the cost function (i.e.P.path) at least one variable that is assigned value 1 has instead to be assigned value 0. We now consider P.lower. Since different lower bound estimation procedures can be used, we will describe in the remainder of this section how to identify an explanation for the bound conflict when using either linear-programming relaxation or Lagrangian relaxation.

Lower Bound Conflicts from Linear-Programming Relaxation

When using linear-programming relaxations as a lower bound estimation procedure, the value of P.lower is obtained according to the formulation described in section 3.1. In order to determine the set of assignments we can deem responsible for P.lower, we must define S as the set of constraints with slack 1 variables assigned value 0 in the linear program solution. These are the constraints which actually limit the value of P.lower. If the literals that assume value 0 in these constraints were to have a different value, some constraints might be satisfied and the value of P.lower would be lower. Therefore, we can consider the assignments to those literals as the responsible for P.lower and define ω pl as:

ω pl = {l : l = 0 ∧ l ∈ ω i ∧ ω i ∈ S} (9)
Clearly, ω pl does not necessarily depend on all decision levels in the search tree. Hence, non-chronological backtracking can result from the conflict analysis procedure.

Lower Bound Conflicts from Lagrangian Relaxation

In order to determine ω pl using the Lagrangian relaxation lower bound estimation procedure as described in section 3.2, we can follow a similar approach to the one described for linear-programming relaxation. Let S be the set of constraints used in obtaining the value of P.lower whose Lagrangian multiplier is different from 0. We can clearly notice from (6) that the constraints with Lagrangian multiplier equal to 0 are irrelevant for computing P.lower. In this case, ω pl can be determined as formulated in [START_REF] Manquinho | Search pruning techniques in sat-based branch-and-bound algorithms for the binate covering problem[END_REF].

Another approach to determine ω pl is to consider the value of α j for each assigned variable from S. If a given variable x j is assigned value 0 and α j > 0, then by changing its value to 1 we would increase the value of P.lower. Or if variable x j is assigned value 1 and α j < 0, if we were to change the value of x j , P.lower would raise. Hence, these assignments cannot be deemed responsible for the value of P.lower and should not be considered in ω pl .

Additional Techniques

In this section we describe additional techniques that have proved useful in implementing a branch-and-bound based PBO tool.

When using LPR for computing lower bounds, one can use the information provided by the LP solution for branching purposes. Essentially, branching is restricted to vari-1 See [START_REF] Nash | Linear and Nonlinear Programming[END_REF] for a definition of slack and artificial variables.

ables for which the LP solution is not integer. Of these variables, the one closest to 0.5 is selected. In the case more than one variable has been assigned value 0.5, then the VSIDS heuristic of Chaff is applied [START_REF] Moskewicz | Chaff: Engineering an Efficient SAT Solver[END_REF].

Another additional technique is the generation of a new constraint when a new better solution is found. When this occurs, P.upper is updated, and one can generate a new constraint:

j∈N c j x j ≤ P.upper -1 (10)
which requires that a new solution must have cost lower that the lowest cost already identified for a solution. This type of constraint is referred to as a knapsack constraint [START_REF] Nemhauser | Integer and Combinatorial Optimization[END_REF]. Suppose also the existence of a cardinality constraint of the form:

j∈K x j ≥ U, K ⊆ N (11)
In this case we can infer a new constraint, that can be added to the set of constraints. Let V be the sum of the U smallest coefficients c j of variables denoted by set K. Thus we can conclude:

j∈K c j x j ≥ V ( 12
)
As a result, the following constraint can be inferred:

j∈N -K c j x j ≤ P.upper -1 -V (13)

Experimental Results

In this section we present empirical results for the techniques described in the paper using our pseudo-boolean optimizer (bsolo) which incorporates the branch-andbound techniques described in the paper and SAT-based techniques, namely boolean constraint propagation, non-chronological backtracking in the search tree and conflict-based learning mechanisms.

The CPU times presented are for a AMD Athlon processor at 1.9GHz with 1GB of physical memory. The time limit for each instance was set to one hour. If the time limit was reached, we provide an indication of which was the best upper bound (ub) value found when the search was stopped.

In order to empirically evaluate the lower bound procedures described in the paper, we ran our solver on most PBO instances available, most of which from EDA applications [START_REF] Aloul | Generic ILP versus specialized 0-1 ILP: An update[END_REF][START_REF] Zhu | Synthesis for mixed ptl/cmos circuit[END_REF][START_REF] Yang | Logic Synthesis and Optimization Benchmarks User Guide[END_REF][START_REF] Walser | 0-1 integer programming benchmarks[END_REF] 2 . Besides bsolo, we also ran PBS [START_REF] Aloul | Generic ILP versus specialized 0-1 ILP: An update[END_REF], Galena [START_REF] Chai | A Fast Pseudo-Boolean Constraint Solver[END_REF] and the commercial MILP solver CPLEX (version 7.5) [START_REF]CPLEX MILP solver[END_REF]. For bsolo, we present results when using different lower bound methods (MIS, LGR and LPR) and also when no lower bound estimation procedure is used (plain).

The bsolo solver was configured to use the constraint strengthening technique described in [START_REF] Dixon | Inference Methods for a Pseudo-Boolean Satisfiability Solver[END_REF] and widely used in mixed integer programming [START_REF] Savelsbergh | Preprocessing and probing for mixed integer programming problems[END_REF]. The probing used in the constraint strengthening is also used to detect necessary assignments during preprocessing. We also used simplification techniques described in [START_REF] Hooker | Logic-Based Methods for Optimization[END_REF][START_REF] Villa | Explicit and Implicit Algorithms for Binate Covering Problems[END_REF] in the synthesis benchmark set.

The results are presented in Table 1. As can be readily concluded, different solvers perform better for different sets of instances. In any case, there are a few clear trends. bsolo (with LPR) is by far the most regular solver, and overall the best performing solver, being able to solve 35 instances out of a total of 40 instances. The other publicly available PBO solvers, PBS and Galena, are solely able to solve 17 and 26 instances, respectively. It should be noted that the general purpose MILP solver CPLEX is extremely competitive, being clearly better that both PBS and Galena over all instances. It is also clear that CPLEX is better suited for optimizing a cost function than for solving a set of pseudoboolean constraints, as the results for the last benchmarks illustrate.

Despite being a more regular and overall better performing algorithm, bsolo performs worse than both PBS and Galena for the last class of instances [START_REF] Walser | 0-1 integer programming benchmarks[END_REF]. This is essentially due to a significantly larger time per decision required by bsolo, motivated by a less optimized implementation of pseudo-boolean constraint propagation. Finally, bsolo with LPR is significantly more efficient than bsolo with LGR. This is motivated by the slow convergence observed for the Lagrangian relaxation on most instances.

Conclusions

The paper describes the integration of lower bound estimation procedures with SAT-based techniques in linear pseudo-boolean optimization, emphasizing two lower bounding methods: linear-programming relaxation and Lagrangian relaxation. The paper also outlines a procedure for enabling non-chronological backtracking in the search tree when bound conflicts occur. In addition, the paper proposes a number of additional techniques, including an informed branching heuristic, and a technique for inferring constraints from knapsack and cardinality constraints.

Preliminary experimental results, obtained on most of existing PBO benchmarks for EDA problems, are clear: bsolo is the most effective PBO solver, far more effective than either PBS [START_REF] Aloul | Generic ILP versus specialized 0-1 ILP: An update[END_REF] or Galena [START_REF] Chai | A Fast Pseudo-Boolean Constraint Solver[END_REF], and is the only dedicated PBO solver that is competitive with (and overall significantly more efficient than) the commercial MILP solver CPLEX [START_REF]CPLEX MILP solver[END_REF].

The very promising experimental results clearly demonstrate that the integration of lower bounding techniques is a crucial aspect in the development of PBO solvers, allowing bsolo to be dramatically more efficient than the best publicly available PBO solvers, PBS and Galena, which do not utilize lower bounding techniques. As the experimental results also suggest, the time per decision in bsolo is larger than in either PBS and Galena. This motivates fine-tuning the implementation and data structures of the bsolo prototype. a For the instances from [START_REF] Walser | 0-1 integer programming benchmarks[END_REF] the CPU times for the different versions of bsolo are the same. No lower bound estimation is used, since there is no cost function.

bsolo

Table 1. Experimental results

Proceedings of the Design, Automation and Test in Europe Conference and Exhibition (DATE'05) 1530-1591/05 $ 20.00 IEEE

Proceedings of the Design, Automation and Test in Europe Conference and Exhibition (DATE'05)

Observe that the ones in[START_REF] Zhu | Synthesis for mixed ptl/cmos circuit[END_REF] represent synthesis problems for mixed PTL/CMOS circuits, and the ones in[START_REF] Walser | 0-1 integer programming benchmarks[END_REF] represent PB satisfaction problems.Proceedings of the Design, Automation and Test in Europe Conference and Exhibition (DATE'05)