Mehrnoosh Mohammadi
email: m.mohammadi94@ut.ac.ir

Simon Thevenin
email: simon.thevenin@imt-atlantique.fr

Alexandre Dolgui
email: alexandre.dolgui@imt-atlantique.fr

-Morocco "New advances and challenges for sustainable and smart industries" Multi-period Multi-item Lot sizing under Random Yield and Stochastic Capacity

Keywords: Lot-sizing, Uncertainties, Random yield, Random capacity, Stochastic production planning

This study addresses a multi-product multi-period lot-sizing problem with stochastic yield and capacity. Typically, random yield is the result of defective items, whereas the capacity is uncertain due to machine breakdowns. Considering both stochastic yield and capacity leads to a complex problem, since these two uncertainty types are in conflict. Indeed, the planners tend to produce large batches to overcome random yields, but the resulting plan might not be feasible due to random capacity. We provide an analytical expression for the service level function. A dynamic programming algorithm is also proposed to find the optimal released quantities. We validate the proposed methodology with a numerical example.

INTRODUCTION

To respond to today's competitive business world, companies tend to increase the size of their product assortments, launch new products frequently, and adjust their production processes constantly (acquire new equipment). In these conditions, production systems can become very difficult to manage due to the lack of regularity and consistency in the production process, making it very difficult for manufacturers to ensure product-process quality. To face these challenges, inventory control systems must be able to cope with uncertainties from both process (e.g., equipment failure) and product quality. The goal of this work is to provide a methodology to size the production lots in a multi-item multi-period production line subject to break-down and defective items. In inventory control, machine breakdowns leads to stochastic resource capacity, whereas bad product quality leads to random yield.

One of the first studies on inventory control with imperfect item quality and its effect on production run time and economic production quantity (EPQ) was [START_REF] Karlin | Studies in the mathematical theory of inventory and production[END_REF]). However, despite the large amount of works following this seminal work, there is a lake of study considering both random yield and capacity. The joint consideration of these uncertainties is critical since they are conflicting. Indeed, companies tends to produce large lot-size to hedge against yield uncertainty. This large lot sizes lead to a tight usage of the production capacity, which is not desirable when the capacity is unknown.

The only works considering jointly these two types of uncertainties are limited to single period models [START_REF] Schemeleva | Evaluation of solution approaches for a stochastic lot-sizing and sequencing problem[END_REF]). To avoid nervousness, companies manage their production plan with a frozen horizon. In this context, companies multi-periods problem, where the uncertainties in a period results in over and under productions that affects future periods.

This study, we consider a single machine production line which operates with a lot for lot policy to satisfy the demand of each item in each period. Multi-item production planning requires the known sequence of products; there are some works in literature that propose a method for schedule the production sequence (e.g., [START_REF] Dolgui | Decomposition approach for a problem of lot-sizing and sequencing under uncertainties[END_REF][START_REF] Schemeleva | Evaluation of solution approaches for a stochastic lot-sizing and sequencing problem[END_REF][START_REF] Schemeleva | Multi-product sequencing and lot-sizing under uncertainties: A memetic algorithm[END_REF], regarding to focus on the lot sizing in this paper, we assume that the sequence is predefined and known. Based on our reviews, lot sizing under yield and capacity uncertainties is consider as a single period problem; it means that at the end of the period there are maybe backlog or inventory which are ignored totally in the probable next period. In this regard, in this study we consider backlog/inventory in all periods. In other word, if there is a shortage in a period we add it to the next period demand and if there is inventory we reduce the next period demand. The goal of the study is to maximize the total service level which is probability of meet the demand of all item in all period at the end of the horizon. We also assume that in each stage of planning the released quantity is frozen and the next stage is planned based on non-changeable decided released quantities of previous stages.

The contributions of this paper include a mathemati-cal analysis of the multi-item lot-sizing under stochastic yield and capacity. This analysis leads to the definition of the service level function, that incorporates the probabilistic backlog/inventory from one period to the next. In addition, we propose a dynamic programming algorithm that find the optimal released quantities, along with some efficient speed up strategies. We evaluate the impact of the multi-period model compare to single period. Finally, we provide some managerial insights.

The rest of this study is as follows: Section2 reviews the literature, section3 introduces the considered problem and gives the proposed mathematical model, section5 presents the employed optimization approach, section6 presents the empirical validation of proposed methodologies and section7 discusses about conclusion of work.

LITERATURE REVIEW

Inventory control plays a crucial role in the management of a production system since it impacts the customer satisfaction and the supply chain costs. Indeed, an improper inventory policy may lead to shortage of products or needless stocks, and both are costly for companies. Therefore, there is a significant motivation for developing efficient and cost-effective inventory policy/methods.

There is a vast literature on lot-sizing problem under uncertainty, and various uncertain parameters have been considered such as demand, yield, lead time or capacity. Review papers on this topic include [START_REF] Aloulou | A bibliography of non-deterministic lotsizing models[END_REF][START_REF] Dolgui | Supply planning under uncertainties in MRP environments: MOSIM'20 -November 12-14, 2020 -Agadir -Morocco A state of the art[END_REF][START_REF] Yano | Lot sizing with random yields: A review[END_REF]. In the recent review of [START_REF] Aloulou | A bibliography of non-deterministic lotsizing models[END_REF]) on non-deterministic lot-sizing models, the authors classified such problems based on: the number of periods, the number of products, the number of machines, the uncertain parameters, and the modelling approaches. The present work focuses on multi-period, multi product, single machine, under uncertain yield and capacity. Therefore, the review below focuses on these two types of uncertainties.

Since the seminal paper of (E. [START_REF] Silver | Establishing the order quantity when the amount received is uncertain[END_REF] where the traditional EOQ model is extended with yield uncertainty, an extensive amount of research has been done on lot sizing problems with random yields (e.g., [START_REF] Gerchak | Multiple lot-sizing, and value of probabilistic information, in production to order of an uncertain size[END_REF][START_REF] Gerchak | Lot sizing in assembly systems with random component yields[END_REF][START_REF] Inderfurth | How to protect against demand and yield risks in MRP systems[END_REF][START_REF] Khan | Economic order quantity model for items with imperfect quality with learning in inspection[END_REF][START_REF] Yano | Lot sizing with random yields: A review[END_REF]. There exists three approaches to represent random yield in a production context. A first stream of research (e.g., [START_REF] Agnihothri | Lot sizing with random yields and tardiness costs[END_REF][START_REF] Dolgui | Supply chain engineering: useful methods and techniques[END_REF][START_REF] Haji | Lot sizing with non-zero setup times for rework[END_REF][START_REF] Li | Periodic-review inventory systems with random yield and demand: Bounds and heuristics[END_REF][START_REF] Papachristos | Economic ordering quantity models for items with imperfect quality[END_REF]Wang and Gerchak 2000) considers that the fraction of defective items is known, and the problem is deterministic. A second stream of research (e.g., [START_REF] Singh | Planning for production of a set of components when yield is random[END_REF][START_REF] Teunter | Lotsizing for a single-stage single-product production system with rework of perishable production defectives[END_REF] represents random yield as Bernoulli process, where a certain fraction of goods is defective, and the number of defective items depends on the lot size. The third research stream (e.g., [START_REF] Maddah | Economic order quantity for items with imperfect quality: revisited[END_REF][START_REF] Salameh | Economic production quantity model for items with imperfect quality[END_REF] models random yield with a geometric distribution, where the manufacturing system becomes out of control with a given probability, and all produced items are defective. However, yield uncertainty does not necessarily result from bad product quality, but it can also arise due to shortage from the suppliers (e.g., [START_REF] Güllü | Analysis of an inventory system under supply uncertainty[END_REF][START_REF] Moon | Inventory systems with variable capacity[END_REF].

Random capacity is a practical issue in different types of industries such as high-tech industries like semiconductor industry [START_REF] Geng | Capacity planning for semiconductor wafer fabrication with uncertain demand and capacity[END_REF]. (Wang and Gerchak 1996a) illustrate the possible impacts of variable capacity on the optimal solution in the continuous lot sizing problem. The uncertainty on the capacity has various causes, such as unexpected breakdowns and unplanned maintenance, material shortages; an uncertain duration of repair, even when the repair is planned and strikes [START_REF] Ciarallo | A periodic review, production planning model with uncertain capacity and uncertain demandâoptimality of extended myopic policies[END_REF]Wang and Gerchak 1996a). Among these causes, machine breakdowns attracted the most attention in the literature (e.g., [START_REF] Chakraborty | Production lot sizing with process deterioration and machine breakdown under inspection schedule[END_REF][START_REF] Chiu | Robust planning in optimization for production system subject to random machine breakdown and failure in rework[END_REF][START_REF] Giri | Optimal lot sizing for an unreliable production system based on net present value approach[END_REF][START_REF] Lin | On a productioninventory system of deteriorating items subject to random machine breakdowns with a fixed repair time[END_REF]. The mathematical modelling of machine breakdown is usually based on two factors, namely, the mean time to/between failure (MTTF/MTBF) and mean time to repair. For instance, [START_REF] Lin | On a productioninventory system of deteriorating items subject to random machine breakdowns with a fixed repair time[END_REF] considered the EPQ model with random breakdownsand constant repair times. Similarly, [START_REF] Giri | Optimal lot sizing for an unreliable production system based on net present value approach[END_REF][START_REF] Giri | Optimal lot sizing for an unreliable production system under partial backlogging and at most two failures in a production cycle[END_REF] examined the optimal lot-sizing problem with at most two failures during the planning period (both failures and repair times follow a general (arbitrary) distribution.

Very few papers [START_REF] Dolgui | Decomposition approach for a problem of lot-sizing and sequencing under uncertainties[END_REF][START_REF] Jain | The single period procurement problem where dedicated supplier capacity can be reserved[END_REF][START_REF] Schemeleva | Evaluation of solution approaches for a stochastic lot-sizing and sequencing problem[END_REF][START_REF] Schemeleva | Multi-product sequencing and lot-sizing under uncertainties: A memetic algorithm[END_REF]Wang and Gerchak 1996b) consider simultaneously random yields and capacity. (Wang and Gerchak 1996b) proposed a stochastic dynamic programming for the single item case, to minimize the expected total cost (i.e., production, holding, and shortage cost). [START_REF] Jain | The single period procurement problem where dedicated supplier capacity can be reserved[END_REF] investigated the case where suppliers allow to reserve a dedicated capacity level by paying a premium charge. These assumptions illustrated the independence of the optimal order size and the random capacity. [START_REF] Dolgui | Decomposition approach for a problem of lot-sizing and sequencing under uncertainties[END_REF]) study a sequencing and lot sizing problem for a multi-item single-period production system under uncertainties. The goal of their problem is maximizing service level, and they use dynamic programming for finding the solution of lot sizing part of the model. This work was extended in [START_REF] Schemeleva | Multi-product sequencing and lot-sizing under uncertainties: A memetic algorithm[END_REF]) and [START_REF] Schemeleva | Evaluation of solution approaches for a stochastic lot-sizing and sequencing problem[END_REF], where the author propose a memetic algorithm for lot sizing (production planning) part and a heuristic method without decomposition the two sub-problems.

To conclude this literature review, despite its relevance in practice, very few studies considered random yield and capacity simultaneously. As these two uncertainty types are in conflicts, it is crucial to manage them simultaneously. Indeed, companies tend to produce pore than required to overcome random yield issues, but a random capacity should not be used tightly. The few studies on lot-sizing with random yield and capacity are limited to one period. To the best of our knowledge, the multi-period multi-item lot sizing problem with random yield and capacity uncertainties has never been study in the literature.

PROBLEM STATEMENT

This section gives a formulation of the considered problem. Table 1 summarizes the notations used in the rest of the paper. Note that we adopt the following convention. M 1...m,1...n denotes a m × n matrix, M i,1...n (resp. M 1...m,i ) denotes the i th lines(resp. column) denotes the i th lines, and M i,j denotes an elements of the matrix. Finally, M i,1→n denotes the sum of the components in the vector M i,1...n .

We consider a multi-item multi-period lot sizing problem with stochastic yield and stochastic capacity. Given the demand D it for each item i in period t of the planning horizon, the considered problem is to determine the production quantity x it to released for each items i in each period t, in order to maximize the service level.

The production must respect the capacity of the bottleneck resource. More precisely, there is C time capacity in each period, and the production of one unit of item i consumes k i time of the capacity. However, the machine can breakdown during operation, and it becomes unusable during a random repair time. Similarly to [START_REF] Dolgui | Decomposition approach for a problem of lot-sizing and sequencing under uncertainties[END_REF], we assume that the failure and repair incidents follow independent exponential distributions, with mean time between failures 1/U and mean time to repair 1/ Ū .

In addition, the machine can produce defective items. More precisely, we assume that each item i is of good quality with probability p i . Therefore, the actual produced quantities y it is lower than or equal to the released quantity x it .

Because of these uncertainties, the produced quantity might not meet the demand, and there can be backlog and inventory. In other words, the unmet demand in period t can be met in period t + 1.

The decision framework assumes that the production quantities are decided for the entire horizon, and they are frozen. Probability to yield at least z it good quality items, given the processed quantities

p c= it (y it |x 1...it )
Probability to process exactly y it item i in period t given the released quantities.

p c> it (y it |x 1...it )
Probability to process at least y it item i in period t given the released quantities.

p cy= it (z it |x it )
Probability to yield exactly z it good quality items, given the released quantities

p cy> it (z it |x it )
Probability to yield at least z it good quality items, given the released quantities

p r (R ≤ K|K)
Probability to have a repair time R lower than K, with a required operation time K.

Table 1 -Notations

MOSIM'20 -November 12-14, 2020 -Agadir -Morocco

In the rest of this paper, we make a distinction between the quantity released, processed and produced. The released quantities (x 1...n,1...t ) is the amount of item the planner asked to produced. The quantity processed (y 1...n,1...t ) is the quantity that can be process by the machine with its actual capacity. The quantity produced (z 1...n,1...t ) is the amount of item process that are of good quality.

The objective is to maximize the service level S, which is defined as the probability to meet the demand for all items in all period. That is

S = n i=1 T t=1 p(z i1→t ≥ D i1→t |x 1...i1...t ), (1)
where

p(z i1→t ≥ D i1→t |x 1...i1...t
) is the probability that the total produce quantities z i1→t from period 1 to t is larger than the total demand of these periods.

MATHEMATICAL ANALYSIS

The service level can be decomposed based on three stochastic parameters

p(z i1→t ≥ D i1→t |x 1...i1...t ) = Di1→t Dit=Di1→t-xi1→t-1 (p b it ( D it |x 1...i1...t-1 ) ( xit yit=0 p y> it ( D it |y it ) × p c= it (y it |x 1...it ))), (2)
where p(z it ≥ D it |y it ) is the probability to produce at least D it good items of type i given that y it items of type i are processed in period t, p c= it (y it |x 1...it ) is the probability to process y it given the quantity released in period t, and p b ( D it |x 1...i1...t-1 ) is the probability to have a total demand D it (including backlog and inventory) in period t . This function tries to maximize the probability of meeting demand of item i from period one to the period t when all i -1 released quantities are planned. The total demand D it can vary from subtraction of sum of the demands of item i from first period to period t (D i1→t ) from the total released quantity of item i until end of period t-1 (x i1→t-1 ), which means all released items have good quality and there maybe inventory from the previous periods, to sum of the demands of item i from first period to period t, which indicates a large backlog. The probability of this value depends on the previous periods released quantities, then base on this value, the function calculate the probability of the meeting this demand with the preassumed released quantity of the item i in period t.

It should be noted that for the first period, there is no variation for the total demand and it is fixed as the initial demands of each item, so we have to delete the external summation and replace D it with D it .

Basically, the larger are the sizes of the lots, the larger is the value of p y it (z it ≥ D it |y it ), and the smaller is p c= it (y it |x it ) since there may be less capacity to produce the last lots.

We detail below the computation of these three probability functions.

Yield uncertainty

The probability to p y it (z it |y it ) to yield exactly z it good quality items when y it items i are processed can be computed with Bernoulliâs formula:

p y= it (z it |y it ) = C zit yit (1 -p i ) (yit-zit) p zit i , (3)
where p i is the given perfection probability.

The probability to yield more than z it item i in period t if y it items are processes is

p y> it (z it |y it ) = yit a=zit C a yit (1 -p i ) (yit-a) p a i (4)
And the probability to yield more than z it item i in period t if x 1...it items are released is

p cy> it (z it |x 1...it ) = xit yit=0 p y> it (z it |y it ) × p c= it (y it |x 1...it ) (5)
Similarly,

p cy= it (z it |x 1...it ) = xit yit=0 p y= it (z it |y it ) × p c= it (y it |x 1...it ) (6)

Capacity uncertainty

The processing of items i starts only once the first i -1 lots are completed. More precisely, y it items i are processed in period t if the repairing time R t respects the capacity constraint

j=i-1 j=1 k j x jt + y it k i + R t ≤ C (7)
where k i is the capacity consumption per unit of item i, and C is the capacity per period.

For sake of clarity, we define:

K = K(x 1...i-1t , y it )) = j∈1...i-1 k j x jt + k i y it (8)
MOSIM'20 -November 12-14, 2020 -Agadir -Morocco and

K = K(x 1...i-1t , y it ) = C -K(x 1...i-1t , y it )) (9)
The probability p c> it (y it |x 1...it ) to process more than y it items i in periods t coincides with the probability p r (R t < K|K) to have less than K units of repair time in period t if the required processing time is K. That is,

p c> it (y it |x 1...it ) = p r (R t < K)|K) (10)
If the time to failure and repair time follow an exponential distribution, [START_REF] Dolgui | Decomposition approach for a problem of lot-sizing and sequencing under uncertainties[END_REF]) showed that such a repair time probability can be computed as follows:

p r (R t ≤ K|K) = 1 -e -(KU + K Ū ) ∞ ν=1 [ (KU ) ν ν! ν-1 j=0 ( K Ū ) j j! ] ( 11
)
where 1 U is the mean time between failures, and 1 Ū is the mean time to repair. As mentioned in [START_REF] Dolgui | Performance analysis model for systems described by renewal process[END_REF], the infinite sum in ( 9) can be approximated with a large enough number, since the function converge quickly.

Finally, the probability p c= it (y it |x 1...it ) to process precisely y it items i in period t. can be inferred from p c> it (y it |x 1...it ) as follows:

p c= it (y it |x 1...it ) =p c> it (y it |x 1...it ) -p c> it (y it + 1|x 1...it ) (12)
In the case of y it = x it , if p c> it (y it |x 1...it ) would be equal to 1, we set p c= it (y it |x 1...it ) for y it = x it equal to 1 and for all y it ≤ x it the value of the probability will ignore and consider it as 0.

Backlog and inventory uncertainty

The backlog/inventory level in period t can be computed from the backlog/inventory level in period t-1. More precisely, the flow conservation constraint can be expressed as follows

D it = D it-1 -z it-1 + D it (13)
Equation ( 14) gives a recursive formula to compute the actual requirements (including backlog/inventory) in period t, where s 1 refers to the production in period t -1, and s 2 refers to the actual requirements in period t -1.

p b it ( D it |x 1...i1...t-1 ) = s1,s2|s2-s1+Dit= Dit p cy= it-1 (s 1 |x 1...it-1 )p b it-1 (s 2 |x 1...i1...t-2 ) (14)
For the first period equation 14 must be changed as equations 15 and 16.

p b i1 ( D i1 = D i1 |∅) = 1 (15) p b i1 ( D i1 = D i1 |∅) = 0 (16)
The computation can be simplified as follows.

p b it ( D it |x 1...i1...t-1 ) = s1=xit-1 s1=0 p cy= it-1 (s 1 |x 1...it-1 )p b it-1 ( D it -D it + s 1 |x 1...i1...t-2 ) (17)
The values of x it are limited to predefined intervals [x it , x + it ] which are computed based on a user defined minimum service level β i .

x + it = min{x it , x it + (C - n j=1 k j x j,t )/k i } (18)
where x it is the minimum lot size of item i such that p cy> it (z it ≥ D it |x it ) ≥ β i , and xit is the minimum lot size of item i such that p cy> it (z it ≥ D it |x it ) ≥ 1 -for a sufficiently small .

OPTIMIZATION APPROACH

We propose below a dynamic programming algorithm to compute the released quantities. In each stage, the algorithm computes the cost H i (V i 1...T ) of the optimal released quantities to complete the first i lots before the vector of times V i 1...T in each period.

H i (V i 1...T ) = max n i=1 p(z i1→T ≥ D i1→T |x 1...i1...T ) | i j=1 k j x jτ = V τ , x jt ∈ [x jt , xjt ], τ ∈ 1 . . . t (19)
with

p(z i1→T ≥ D i1→T |x 1...i1...T ) = T t=1 p(z i1→t ≥ D i1→t |x 1...i1...t ) (20)
The cost H i (V i 1...T ) of the optimal lots sizes up to item i with completion times vector V i 1...T , can be computed from the optimal cost of the previous stage MOSIM'20 -November 12-14, 2020 -Agadir -Morocco

H (i-1) (V (i-1) 1...T ) : H i (V i 1...T ) = max{H (i-1) (V (i-1) 1...T )× p(z i1→t ≥ D i1→t |x 1...i1...t ) |k j x jτ = (V jτ -V (j-1)τ ), x jτ ∈ [x jτ , xjτ ], V (i-1)τ ∈ [V (i-1)τ , V (i-1)τ ], (V jτ -V (j-1)τ ) ∈ [w jτ , wjτ ], τ ∈ 1 . . . T , j ∈ 1 . . . i} (21)
where

V it = i j=1 w jt (22) Vit = min{T t - n j=i+1 w jt , i j=1
wjt } (23)

w jt = k j x jt (24) wjt = k j x + jt ( 25
)
The initial state V 0 1...T = 0 . . . 0, and H 0 (V 0 1...T ) = 1.

Speedup strategy

This section provide the complexity of the proposed dynamic program, and it shows that the number of states grows exponentially with the number of period. To circumvent this issue, we propose to prune the states based on a lower bound and dominance properties.

A lower bound on the service level: Such a lower bound can be computed with a heuristic. To do so, we use the method proposed in ... for a the single period lot-sizing with random yield and capacity. This method is applied successively on each period, and this yield a set of released quantity x h 1...n,1...T . The service level associated with x h 1...n,1...T is computed with equation ( 6). The resulting service level is sub-optimal, but this lower bound allows to prune a large number of node.

Dominated states:

A state V i 1...T dominates the state W i1...T , if V i 1...T ≤ W i1...T and H i (V i 1...T ) ≥ H i (W i1...T ).
Algorithm 1 Dynamic programming algorithms Set V 0 1...T = 0 . . . 0, and H 0 (V 0 1...T ) = 1 for all items i ∈ {1 . . . n} do for all The optimization procedure is applied on the following small example.

V i 1...T such that V it ∈ V it , Vit do H i (V i 1...T ) = 0 for all x i1...T such that x it ∈ [x it , xit ] do for all states V (i-1) 1...T such that V (i-1) 1...T + K(x i1...T ) = V i 1...T do if H (i-1) (V (i-1) 1...T ) × p(z i1→T ≥ D i1→T |x 1...i1...T ) ≥ H i (V i 1...T ) then H i (V i 1...T ) = H (i-1) (V (i-1) 1...T ) × p(z i1→T ≥ D i1→T |x 1...i1
We assume a single-machine line that produces two types of item (n=2) during two periods (T=2).Each period duration is equal to 1.2 hours (C = 1.2).The other feature of the problem is also known as follows(the times are given in hours).

• The initial demand of the items (D it )

D 11 = 2, D 12 = 1, D 21 = 1, D 22 = 3
• The processing time for each type of items k 1 = 0.17, k 2 = 0.09

• The mean time between failures (MTBF) U = 0.6667

• The mean time to repair (MTTR) Ū = 4

• Probability of defectiveness for each type of items p 1 = 0.85, p 2 = 0.83

• Sequence of producing lots γ = {γ 1 = 1, γ 2 = 2}
It is considered that probability of minimum service for all lots in all periods is equal to 0.93 (β 1 = β 2 = β 3 = β 4 = 0.93). As it is obvious to maximize the service level we should obtain at least the produced quantities (z it ) as same as demand quantities for each item at the end of each period.

The objective is to determine the optimal lot sizes x it when we know the sequence in order to maximize the probability of equation 13.

In the optimization algorithm, the items appears as the stages of dynamic programming and number of stages varies by the number of items. The number of periods has impact on the number of states in each stage and also the complexity of computation of each loop in the algorithm.

EXPERIMENTAL RESULTS

With considering probability of minimum service level as 0.93 for items and using equation ( 18) we can obtain x it ∈ [x it , x + it ] interval as:

x 11 ∈ [3, 6], x 12 ∈ [2, 4], x 21 ∈ [2, 5], x 22 ∈ [5, 8]
To be able to apply optimization procedure the components of dynamic programming (equation 19 to 25) are calculated in the table 2.

This instance has two stages. The first stage is producing item 1 and the second is producing item 2. Below, all possible released quantities and states of each stage based on the x it and V it bounds are shown.

Stage 1: ) for all possible states of producing item 1 is in the vector below: .9143; 0.9337; 0.9381; 0.9822; 0.9868; 0.9878; 0.9966; 0.9975; 0.9977; 0.9994; 0.9996; 0.9996] If we assume that periods are independent (single period model) the released quantities for lots are:

x 11...2 ∈ {[3, 2], [3, 3], [3, 4], [4, 2], [4, 3], [4, 4], [5, 2], [5, 3], , [6, 2], [6, 3], [6, 4]} V 1 1...2 ∈ {[0.
x 21...2 ∈ {[2, 5], [2, 6], [2, 7], [2, 8], [3, 5], [3, 6], [3, 7], [3, 8], [4, 5], [4, 6], [4, 7], [4, 8], [5, 5], [5, 6], [5, 7], [5, 8]} V 2 1...2 ∈ {[0.
p(z 11→1 ≥ D 11→1 |[3]) = 1 × (0 × 0 + 0 × 0 + 0.7225 × 0 + 0.9393 × 1) = 0.9393 The next is calculation of p(z 11→2 ≥ D 11→2 |[3, 2]) value. p(z 11→2 ≥ D 11→2 |[3, 2]) = 3
p(z 11→2 ≥ D 11→2 |[3, 2]) = p(z 11→1 ≥ D 11→1 |[3]) × p(z 11→2 ≥ D 11→2 |[3, 2]) = 0.
H 1 (V 1 1..2 ) = [0
x 11 = 5, x 12 = 3, x 21 = 3, x 22 = 6
When we consider this amount as production quantities in multi-period model the service level will be: 0.9913. We set this value as a lower bound for H i (V i1..T ) and delete the states with H lower than this bound for the next item in algorithm. So the new set of states of stage 1 and their related H value are: The optimal answer is max of H 2 (V 2 1..2 ), so: customer service =H 2 ([1.12, 1.06]) = 0.9916

V 1 1...2 ∈ {[0.
And the released quantities are :

x 11 = 5, x 12 = 2, x 21 = 3, x 22 = 8

Analysis and comparison

In this study we develop a multi-item multi-period stochastic lot sizing under yield and capacity uncertainties,in literature the models proposed for lot sizing under yield and capacity uncertainties are singleperiod horizon such as [START_REF] Dolgui | Decomposition approach for a problem of lot-sizing and sequencing under uncertainties[END_REF][START_REF] Schemeleva | Evaluation of solution approaches for a stochastic lot-sizing and sequencing problem[END_REF][START_REF] Schemeleva | Multi-product sequencing and lot-sizing under uncertainties: A memetic algorithm[END_REF],regarding to the literature which indicate single-period method is popular in the problem, here,we compare our model which consider the backlog/inventory in multi-period horizon production with a single-period horizon one in which the conditions and assumptions are the same. For comparison we use the example which is given in section 6.1.

The considered example include two periods. First we solve the problem as a multi-period horizon and then for single-period method we divided the horizon to the two independent period and ignore the backlog/inventory assumption.

The table 3 illustrates the results. In the table the multi period method refers to the method in this study assuming backlog/inventory and single period method refers to model which ignore this assumption.

The table 3 indicates that the results for the first period are precisely the same in the both methods,it is found out easily, the reason is that there is no backlog/inventory in the beginning of the first period.The difference is in the second period. In multi-period horizon, the probability of maximizing service level is higher than single period. It should be noted, since in the multi-period method the periods are depended it's not possible to compute the each period service level.

If we produce as amount as optimal quantities in the single-period horizon, in the multi-period horizon the As it is mentioned before, here backlog/inventory means that the production quantity is less or more than enough, in other words if we may obtain more good items than the demand quantity we keep remaining to the next period or if we may get less than the demand quantity we add the backlogs to the demand of the same item in the next period. That is why we could get higher probability of service level in comparison with single period horizon.

Overall, this work in comparison with the lot sizing problem in literature results in a higher probability of satisfying demands, but it needs more time to solve the problem. It worth to mention the instance in this study is very small, if this method is applied on a larger problem the efficiency will be more obvious.

CONCLUSIONS

In this research, a multi-item multi-period lot sizing problem for a single-machine production line considering random yield and random capacity was studied.In which the uncertainties rise from item imperfection and machine breakdowns. A predefined item sequences and backlog/inventory allowance are assumed in the production line. The problem is to estimate the optimal lot size of each item to maximize the probability of having a sufficient number of good quality items by the end of a given planning horizon, i.e. maximize the overall customer service level.

Regarding the similarity of the problem to the wellknown knapsack problem, dynamic programming was employed to estimate the optimal lot sizes to be produced. As we went deep into the arts this problem has not been studied in the literature sufficiently specially for the multi-period horizon.

In order to validate the model and the proposed optimization approach, a numerical example was illustrated. Later, a comparison between this study, multi-period stochastic lot sizing and single-period stochastic lot sizing was drawn. The obtained results demonstrate the better efficiency of considering multi-period than single-period setting which is effect of considering backlog/inventory. This research can be extended in some direction as follow: Extending the approach to solve similar problems considering other objective functions such as minimizing production cost and also Implementing other optimization approaches such as (meta)heuristics to empower the model to solve the large problems.

  12 ( D 12 |y 12 ) × p c= 12 (y 12 |[2]) To compute the probabilities of actual requirements (p b 12 ( D 12 |[3])) of the period to equation 15, 16 and 17 are used. Sop(z 11→2 ≥ D 11→2 |[3, 2]) = 0.6141 × (1 × 0 + 1 × 0 + 1 × 1)+0.3251×(0×0+0.85×0+0.9775×1)+0.0574×(0× 0+0×0+0.7225×1)+0.0034×(1×0+1×0+0×1) = 0.9734

  Then it continues by considering first stage or item, i.e i = 1. The first selected state of first stage in the set of possible state is V 1 1...2 = [0.51, 0.34] and H 1 ([0.51, 0.34]) = 0 is set for this state.

	The dynamic programming algorithm calculation be-
	gins with setting V 0 1...2 = [0, 0], and H 0 (V 0 1...2 ) = 1;
	Procedure is started with x 11...2 = [3, 2], since this is
	the first stage of dynamic programming there is just
	one state for the V 0 1...2 which is equal to [0,0], so it
	should be checked if K([3, 2] = [0.51, 0.34] with using
	equation 8.K[3, 2] = 0.17 * [3, 2] = [0.51, 0.34] which
	is equal to V 1 1...2 = [0.51, 0.34].
	The next step is to check if H 0 ([0, 0]) × p(z 11→2 ≥
	D 11→2 |[3, 2]) ≥ H 1 ([0.51, 0.34]) is satisfied or not;
	where H 0 ([0, 0]) is equal to 1. According to equa-
	tion 20, p(z 11→2 ≥ D 11→2 |[3, 2]) has 2 component as
	follows:
	p(z 11→2 ≥ D 11→2 |[3, 2]) =		2 t=1 p(z 11→t ≥
	D 11→t |x 1...11...t )
	where for t=1, p(z 11→1 ≥ D 11 |[3]) should be calcu-
	lated.
	p(z 11→1	≥	D 11 |[3])	=	p b 11 (D 11 |∅)×
	3 y11=0 p y> 11 (D 11 |y 11 ) × p c= 11 (y 11 |[3])
	At last:
	69, 0.79], [0.69, 0.88], [0.69, 0.97], [0.69, 1.06],
	[0.69, 1.05], [0.69, 1.14], [0.69, 1.13], [0.78, 0.79], [0.78, 0.88],
	[0.78, 0.97], [0.78, 1.06], [0.78, 0.96], [0.78, 1.05], [0.78, 1.14],
	[0.78, 1.13], [0.87, 0.79], [0.87, 0.88], [0.87, 0.97], [0.87, 1.06],
	[0.87, 0.96], [0.87, 1.05], [0.87, 1.14], [0.87, 1.13], [0.96, 0.79],
	[0.96, 0.88], [0.96, 0.97], [0.96, 1.06], [0.96, 0.96], [0.96, 1.05],
	[0.96, 1.14], [0.96, 1.13], [0.86, 0.79], [0.86, 0.88], [0.86, 0.97],
	[0.86, 0.96], [0.86, 1.05], [0.86, 1.14], [0.86, 1.13], [0.95, 0.79],
	[0.95, 0.88], [0.95, 0.97], [0.95, 1.06], [0.95, 0.96], [0.95, 1.05],
	[0.95, 1.14], [0.95, 1.13], [1.04, 0.79], [1.04, 0.88], [1.04, 0.97],
	[1.04, 1.06], [1.04, 0.96], [1.04, 1.05], [1.04, 1.14], [1.04, 1.13],
	[1.13, 0.79], [[1.13, 0.88], [1.13, 0.97], [1.13, 1.06], e[1.13, 0.96],
	[1.13, 1.05], [1.13, 1.14], [1.13, 1.13], [1.03, 0.79], [1.03, 0.88],
	[1.03, 0.97], [1.03, 1.06], [1.03, 0.96], [1.03, 1.05], [1.03, 1.14],
	[1.03, 1.13], [1.12, 0.79], [1.12, 0.88], [1.12, 0.97], [1.12, 1.06],
	[1.12, 0.96], [1.12, 1.05], [1.12, 1.14], [1.12, 1.13], [1.2, 0.79],
	[1.2, 0.88], [1.2, 0.97], [1.2, 1.06], [1.2, 0.96], [1.2, 1.05],
	[1.2, 1.14], [1.2, 1.13]}

Table 3 -

 3 Results of comparison of two methods probability of satisfying the demand, objective function value would be equal to 0.9913. The result shows that the objective value (service level) wouldn't be optimal.

		multi period method	single period method
	periods t	1		2		1		2
	items i
	1	x 11	5	x 12	2	x 11	5	x 12	3
	2	x 21	3	x 22	8	x 21	3	x 22	6
	period service level	-		-		0/9926	0/9873
	total service level			0/9913				0/9800


Sakna Bazzi

Emmanuelle Schulz

Mohamed Mellah

Electrogenerated Sm(II)-Catalyzed CO 2 Activation for Carboxylation of Benzyl Halides

Sm(II)-catalyzed carboxylation of benzyl halides is reported through the electrochemical reduction of CO 2 . The transformation proceeds under mild reaction conditions to afford the corresponding phenylacetic acids in good to excellent yields. This user-friendly and operationally simple protocol represents an alternative to traditional strategies, which usually proceeds through C(sp 3 )-halide activation pathway.

Scheme 1. Carboxylation strategies for the preparation of phenylacetic acids using CO 2 as a C-1 building block

Carboxylic acids, which are high-added-value industrial compounds are in general the most targeted products for this type of transformation. More specifically, phenylacetic acids have received increasing attention as their scaffolds can be found in drugs that are used on a daily basis (e.g., ibuprofen, naproxen). [START_REF] Maag | Prodrugs of Carboxylic Acids[END_REF] In recent years, elegant metal-catalyzed reactions were developed to produce phenylacetic acids from styrene derivatives 5 or organozinc reagents. 6 In 2013, Martin and co-workers described the first catalytic carboxylation of primary benzyl chlorides, featuring a nickel-based catalytic system, to generate phenylacetic acids. [START_REF] Correa | Direct Carboxylation of Benzyl Halides with CO2[END_REF] Although remarkable, this reaction still required the use of over-stoichiometric amounts of Zn and MgCl 2 . [START_REF] Sayyed | The Crucial Roles of MgCl2 as a Non-Innocent Additive in the Ni-catalyzed Carboxylation of Benzyl Halide with CO2[END_REF] Electrochemistry represents a viable alternative by eliminating the use of toxic and hazardous reducing agents. [START_REF] Senboku | Electrochemical Fixation of Carbon Dioxide[END_REF]10 Major advances have been made with the use of sacrificial anodes (e.g. Mg, Al), which enabled the in situ formation of an organometallic species. 11 In general, the electrocarboxylation of benzyl halides involves an initial reduction, forming either the benzyl radical or the anion, which can then react with CO 2 to yield phenylacetic acid. The CO 2 activation offers probably an alternative to reach more efficient carboxylation by avoiding the side products resulting from dehalogenation or dimerization, which are common drawbacks in all these transformations. [START_REF] Correa | Direct Carboxylation of Benzyl Halides with CO2[END_REF] To the best of our knowledge, the direct carboxylation of alkyl halides starting with CO 2 reduction has not yet been described in the literature.

Recently, we reported a procedure for the carboxylation of aryl halides with the use of a "soluble" Sm anode, affording functionalized benzoic acids from aryl chlorides. We demonstrated that the electrolysis enabled the generation in situ of a powerful Sm(II) species, which directly reduced CO 2 . [START_REF] Bazzi | CO2 Activation by Electrogenerated Divalent Samarium for Aryl Halides Carboxylation[END_REF] We depict herein our findings on the catalytic carboxylation of benzyl halides via CO 2 activation promoted by electrochemically generated Sm(II) species to access highly sought-after phenylacetic acids, including Naproxen. As a first step, we explored the carboxylation of benzyl chloride (1a) as a benchmark reaction, using a Sm anode to produce Sm(II) reagent under 1 atm CO 2 (Figure 1). We performed the electrolysis in an undivided cell under standard stoichiometric conditions. [START_REF] Bazzi | CO2 Activation by Electrogenerated Divalent Samarium for Aryl Halides Carboxylation[END_REF] The reaction was implemented by introducing 1a in a solution of nBu 4 NBF 4 in DMF with a continuous bubbling of CO 2 .

Figure 1. Electrocarboxylation of 1a using Sm "soluble" anode a

a Standard stochiometric conditions: Undivided cell, fitted with two electrodes a samarium anode (diameter rod = 1 cm) and a stainless-steel cathode grid (20 cm 2 ), 1a (5 mmol), solvent (50 mL), 4 hours at 0.1 A. Isolated yields; b the temperature was monitored during the addition of dry ice, and it remained between 20-23°C, c 1a (1 mmol) after 2 h, d nBu4NBF4 (0.04 M).

The reaction occurred with low conversion and led to the formation of acid 2a in a low yield (Figure 1, entry 1), which could be improved to 29% by replacing DMF with CH 3 CN (Figure 1, entry 2). This result is likely due to the higher solubility of CO 2 in CH 3 CN (0.59 M versus 0.38 M in DMF). [START_REF] Gennaro | Solubility and Electrochemical Determination of CO2 in Some Dipolar Aprotic Solvents[END_REF] In order to increase the availability of CO 2 in solution, dry ice was employed and, to our delight, 2a was obtained in quantitative yield (Figure 1, entry 4). Then, we evaluated the scope of the carboxylation with a series of benzyl chlorides. The corresponding phenylacetic acids were obtained in yields ranging from 60 to 99% (Figure 1). It is noteworthy that this carboxylation was not limited to primary benzyl chlorides but could be also extended to secondary benzyl chlorides (see 2e in Figure 1).

We then turned our attention to the establishement of a catalytic version (Table 1). Building on our precedent reports regarding the implementation of catalytic procedures, 14 we conducted a preelectrolysis using a solution of nBu 4 NBF 4 [0.02 M] in CH 3 CN at a constant current intensity of 0.1 A. This first step was carried out to generate 20 mol % of Sm(II) with respect to the substrate (1.0 mmol). The benzyl chloride 1a and an oxophilic reagent (3.0 equiv.) were then added to the solution. From there, the polarity of the electrode was inverted to act as a cathode [START_REF]When the samarium electrode is used as the cathode, this later is not consumed allowing its use for a very long time. Of note, before each electrolysis[END_REF] throughout the electrolysis, achieved at a constant current intensity of 0.1 A. It is important to stress that the use of an oxophilic reagent is critical in catalytic processes involving SmI 2 . 16,[START_REF] Huang | Of note, Procter and his group reported very recently the first catalytic use of SmI2, avoiding the need of stoichiometric co-reductant and additives[END_REF] Indeed, this additive promotes the cleavage of the Sm III -O bonds. In this way, the released trivalent Sm(III) species can be reduced onto the cathode to regenerate the Sm(II) active species. The use of TMSOTf resulted in the formation of the targeted acid 2a in 45% yield. Interestingly, switching from TMSOTf to TMSCl increased the yield to 60% (Table 1, entry 1 and 2). Various Brønsted acids were also screened as a proton source (Table 1, entry 3-5) in the hope that the resulting Sm salts could be reduced without having to rely on a Si-based additive. However, the reactions were completely shut down, and 1a remained intact. The effect of Sm ligand was also evaluated by replacing nBu 4 NBF 4 by nBu 4 NOTf, but it did not have a significant impact (Table 1, entry 6). In contrast, the use of nBu 4 NI improved drastically the reactivity, providing 2a in 98% yield and the amount of TMSCl could be lowered to 1.5 equivalent without being detrimental to the outcome of the electrocarboxylation (Table 1,entry 7,8). Eventually, we found out that using a catalytic loading of 20 mol% was essential for a smooth functioning of the reaction. Control experiments were carried out to ascertain the role of each parameter. Thus, the use of a glassy carbon cathode instead of a samarium one furnished 2a in only 30% yield, which confirmed the need for a Sm cathode in order to efficiently regenerate the Sm(II) catalytic species. 14a The direct reduction onto the Sm cathode without catalyst left 1a intact (Table 1, entry 11). As anticipated, the catalytic turnover cannot be achieved without TMSCl as only 8% of 2a was isolated (Table 1,entry 12). With these optimized conditions in hand, we explored the scope of the reaction by using a series of commercially available benzyl chlorides and bromides (Scheme 2). We were pleased to find that the transformation tolerated a wide range of substituted benzyl halides bearing both electron-donating (1b-1d and 1f-g) and withdrawing groups (1h-1m). Besides, our electrocatalytic reaction demonstrated excellent chemoselectivity. Accordingly, additional functionalities (halides (2h, 2i, and 2k), esters (2m), and alkenes (2n)) remained intact without any dehalogenation or dicarboxylation side reactions. The reaction was not precluded by the presence of ortho substituents (2c, 2i, and 2k). Of note, benzyl bromides could be also subjected to the carboxylation conditions to access similar compounds in higher yields (see 2b, 2d, 2f, and 2n). Another interesting trademark is the excellent reactivity of secondary benzyl chlorides, which gave the relevant products with no yield drop (2e and 2o-2t). a Standard catalytic conditions: Undivided cell, fitted with two electrodes: cathode of samarium rod (diameter 1 cm) and anode of stainless-steel grid (20 cm 2 ), i = 100 mA, (1a) (1 mmol), nBu4NI (1 mmol), TMSCl (1.5 mmol), solvent (50 mL), 2 h. Dry ice was added in small pieces each 15 min. b Isolated yields obtained from benzyl bromides.

Scheme 2. Sm(II)-catalyzed electrocarboxylation of benzyl chlorides a

It is even more appealing when one considers the low propensity of secondary alkyl chlorides to undergo oxidative addition in organometallic-based strategies, making their carboxylation still challenging. [START_REF] Correa | Direct Carboxylation of Benzyl Halides with CO2[END_REF]18 Finally, as a proof of the robustness of our electrocatalytic process, we succeeded to prepare (±)Naproxen 2t in 43% yield on a 5 mmol scale. Then, we conducted complementary experiments to gain some insight into the mechanism (Figure 2 and SI for details). Therefore, the electrolysis was performed in the absence of benzyl chloride, which produced a white-grey precipitate. [START_REF] Gennaro | Solubility and Electrochemical Determination of CO2 in Some Dipolar Aprotic Solvents[END_REF] C-NMR of the crude reaction mixture revealed the formation of oxalic acid, demonstrating that the samarium can reduce efficiently CO 2 .

Furthermore, running the electrolysis without CO 2 let 1a intact. Both experiments clearly confirmed that the carboxylation was initiated by the reaction between the CO 2 radical anion and the benzyl chloride. The influence of the samarium ligand was also evaluated by carrying out the electrocatalytic carboxylation with chloride as the only possible source of ligand for samarium, using benzyl chloride, TMSCl and 20 mol% of SmCl 3 as a catalyst (see SI). The reaction delivered the phenylacetic acid in 75% yield, indicating that the ligand exchange was not the key step of the electrocarboxylation. In addition, radical trapping experiments were performed by adding TEMPO (1 equiv.) to the reaction mixture (Figure 2). Even though the TEMPO was completely consumed, the formation of product 3 was not observed, which suggests that the generation of the benzyl radical did not occur. It is also noteworthy that 2a was not obtained either. On the other hand, prolonging the electrolysis time to 4 hours afforded 2a in a 90% yield. Thus, we assumed that TEMPO had quenched the CO 2 radical anion in the first place but, after the complete consumption of the radical scavenger, the reactivity could be restored to yield 2a. Moreover, a radical clock experiment was performed with the aim of isolating a compound resulting from the insertion of CO 2 . Indeed, the reaction led to product 4 in 23% yield, providing further evidence of the CO 2 radical anion formation. Finally, cyclic voltammetry measurements were conducted to assess the redox behavior of the catalyst. In the absence of CO 2 , the quasireversible Sm(III)/Sm(II) system was observed at -1.4 V/ SCE. Once CO 2 was introduced, a weak reduction wave was noticed with a potential shift around -1.8 V/SCE and the oxidation wave disappeared. This behavior confirmed the chemical interaction between Sm(II) and CO 2 . Then, the addition of a mixture of substrate 1a and TMSCl (ratio 1:1.5) triggered a significant reduction wave, which started at -0.95 V/SCE and increased significantly upon further addition of this mixture. Overall, this electrochemical behavior of the catalyst proves the appearance of a catalytic current.

Based on these experiments, the following mechanism is proposed for this electrochemical process (Scheme 3). Initially, Sm(II) could reduce the carbon dioxide into the CO 2 radical anion A, before being engaged in a radical substitution with benzyl chloride to produce the samarium carboxylate intermediate B.

Then, TMSCl would furnish the silyl ester C, which leads to the carboxylic acid 2a upon treatment. Finally, the catalyst would be recovered by reduction onto the samarium cathode.

Scheme 3. Proposed mechanism.

In conclusion, we devised a novel catalytic carboxylation of primary and secondary benzyl halides relying on electrogenerated SmI 2 .

When compared to previous reports, this electrocarboxylation process proceeds through different pathway involving a CO 2 activation instead of C(sp 3

)-halide one. This represents an alternative to other existing 2 fixation methodologies that require the preparation of welldefined and sensitive reagents. This catalytic process provides streamline to valuable phenylacetic including Naproxen. Further developments of an enantioselective carboxylation reaction and an extension to unactivated alkyl halides are currently under investigation.

ASSOCIATED CONTENT Supporting Information

The Supporting Information is available free of charge on the ACS Publications. Experimental procedures, compounds data, and spectra.

Table 1 .a

 1 Electrogenerated Sm(II)-catalyzed carboxylation of benzyl chloride 1a. Reaction conditions: anode of samarium (rod) and cathode of stainless-steel grid (20 cm2), with a solution of nBu4NX (1 equiv) in CH3CN (50 mL). The electrolysis was performed for 386 seconds at i= 100mA, before switching the polarity of the electrodes and adding dry ice, 1a (1 mmol) and the additive (3 mmol). Dry ice was added in small pieces each 15 min; b the temperature was monitored during the addition of dry ice, and it remained between 20-23°C; c Isolated yields; d TMSCl 1.5 mmol. e Glassy carbon cathode instead of Sm cathode.

Figure 2 .

 2 Figure 2. Preliminary mechanistic experiments Cyclic voltammetry conditions: Glassy carbon working electrode (surface 20 mm 2 ), Pt wire as counter electrode and SCE reference, scanning potential between -0.5 and -2 V in CH3CN with nBu4NPF6 [0.1 M]. Scan rate: 100 mV/s. (a) : nBu4NPF6 [0.1 M] in CH3CN; (b) : after introducing 0.2 M of electrogenerated SmI2; (c) : After addition of dry ice to the solution (b); (d): Addition of 0.5 mL of CH3CN solution containing BnCl [1 M] and TMSCl [1.5 M]; (e): Addition of 0.5 mL of CH3CN solution containing BnCl [1 M] and TMSCl [1.5 M] to solution (d).

ACKNOWLEDGMENT

We thank the CNRS, the Ministère de l'Enseignement Supérieur et de la Recherche and the Charm3at LABEX(ANR-11-Labex-0039). We also thank Dr. David Leboeuf (ICMMO-UMR CNRS 8182) for fruitful discussions.

AUTHOR INFORMATION

Corresponding Author * mohamed.mellah@u-psud.fr

Author Contributions

Notes

The authors declare no competing financial interests.


W Bekir

O Messal

L Lambourg

F Dubas

S Harmand

A Kedous-Lebouc

F Gillon

on Ecological Vehicles and Renewable Energies (EVER) Loss Model Impact on Multi-Physical Coupling Designed for Sizing Optimization of Electrical Machines

Keywords: Electrical machines, electric traction vehicle, iron losses, multi-physical model

This article analyzes and highlights the impact of an iron loss model on temperature estimation within a multiphysical model of a radial-flux permanent-magnet (PM) synchronous machine (PMSM) for traction application. Three iron loss models have been compared in terms of computing time and accuracy within the machine torque-speed plane. This study underlines the iron loss calculation impact on the temperature estimation in the PMSM in relation with the required performances.

I. INTRODUCTION

The energetic efficiency improvement in the transport sector and particularly in the automotive field has become a major challenge of crucial importance. How to reduce the energy consumption? One possible way is by enhancing the performances of traction systems. In order to improve or optimize these systems, we should be able to correctly model them considering all the main phenomena involved including losses. Thus, a multi-physical model is required. This model should be sufficiently accuracy but at the same time, its computing time should be short due to the high number of model evaluations [START_REF] Caillard | Mono and bi-Level optimization architectures for powertrain design[END_REF].

One of the difficulties faced during electrical traction chain design is to be able to estimate enough accurately the iron loss [START_REF] Frias | Loss Minimization of an Electrical Vehicle Machine Considering Its Control and Iron Losses[END_REF]. In practice, the electric traction machines operate in a wide range of torque and speed (Torque/Speed plane), leading to variable working frequency. In addition, these machines have weakening mode control. For a relevant iron loss prediction, it is mandatory that the model takes into consideration all these behaviors.

During the last years, the machine iron loss evaluation has become a challenge of many researches. Even though the iron loss relation to the total loss varies according to the machine under consideration, thy can be dominant when the frequency is high [START_REF] Xue | Thermal-Loss Coupling Analysis of an Electrical Machine Using the Improved Temperature-Dependent Iron Loss Model[END_REF]. Indeed, at high-speed, the iron loss can reach 80% of the total losses [START_REF] Bui | Core losses measurement technique for high frequency and flux density of switched reluctance machines[END_REF]. This is prejudicial, because, the loss downgrades the energy efficiency of the machine on the one hand, and causes the overheating of its magnetic parts on the other hand. In this way, the machine temperature rise depends partially on the iron loss estimation. Studying this relation requires a multi-physical model [START_REF] Xue | A New Iron Loss Model for Temperature Dependencies of Hysteresis and Eddy Current Losses in Electrical Machines[END_REF]. The use of such model coupled to a thermal one is not common. The developed model will allow estimating not just the traction machine yield but also the temperature inside the machine, which means that it will be possible to analyze the impact of iron loss estimation error as a temperature variation.

Three iron loss models with different complexity levels have been investigated. The Loss Surface (LS) hysteresis model will be taken as a reference [START_REF] Bui | Core losses measurement technique for high frequency and flux density of switched reluctance machines[END_REF]. This model is considered to be the most accurate but is time-consuming. An intermediate model fast enough and with few coefficients is later on selected [START_REF] Ployard | Iron loss computation of a hybrid excitation synchronous generator[END_REF]. Finally, the classical Steinmetz model [START_REF] Steinmetz | On the law of hysteresis[END_REF] is used for its simplicity and high speed.

The paper is organized as follows: the section Ⅱ introduces the multi-physical coupling as well as the iron loss models and their calculation times. The section Ⅲ describes first the coupling developed for to the PMSM, then the iron loss results are presented and discussed from a local and global point of view is studied. The section Ⅳ shows the iron loss calculation impact on the machine performances for eight operating points, and then the loss and efficiency cartographies are presented. Finally, a comparison between iron loss models is made by considering the root-mean-square deviation (RMSD) between the loss models in terms of power (W) and temperature (°C).

II. MULTIPHYSICAL MODEL FOR SIZING OPTIMIZATION

Traction machine designing requires a multi-physical model. The magnetic and thermal model are based on lumped constant managing a compromise that includes computing time and accuracy. Fig. 1 represents the multi-physical model structure with three related layers.

The first layer integrates a cost model and a generalized nonlinear adaptive magnetic equivalent circuit (MEC). This model is mesh-based reluctance network. It takes into consideration the saturation effect and uses the "air-gap sliding-line technic" [START_REF] Benlamine | Nonlinear Adaptive Magnetic Equivalent Circuit of a Radial-Flux Interior Permanent-Magnet Machine Using Air-Gap Sliding-Line Technic[END_REF] in order to incorporate the machine motion. The method allows a different discretization between the stator and the rotor, so a possible adaptation of the mesh density. The second layer contains three iron loss models with different level of accuracy and calculation times. The third layer includes an aero-thermal model (SAME) allowing fast and accurate determination of the thermal machine behavior by taking into account the convective exchanges, thermal conduction and the power loss distribution. This model uses the nodal method, which consists in discretizing the geometry into many isothermal volumes wherein the centers are linked by thermal conductances [START_REF] Jandaud | Aero-Thermal Optimization of a Heat Sink using Variable Neighbourhood Search[END_REF]. Knowing the local iron and Joule losses within these volumes, the model provides temperature cartography of the machine. This coupling is built as a single function. The inputs are: the geometrical machine parameters [START_REF] Messal | Iron Losses in Electromagnetic Devices: Nonlinear Adaptive MEC & Dynamic Hysteresis Model[END_REF], the 3 control parameters (stator current Ismax (A), the phase angle between the stator current and voltage ѱ (°) and frequency f (Hz)), as well as the initial PMs and coil temperatures. The outputs are: the machine electromagnetic performances (torque, etc.), the Joule loss, the iron loss according to the used model (LS, M4+, St) and the temperatures obtained from the whole calculated losses. The function includes also some options that allow the model configuration (models choice, parallelization, display, etc.). The generated quantities are spatio-temporal. For example, the radial and tangential components of magnetic flux density (viz., Br and Bθ) are calculated in each cell of the MEC as a function of the rotor position (i.e., time dependence). The model structuring as well as its data are of prime importance to rigorously study the machine behaviors.

A. Zoom on the iron loss models

Many models are used in literature to estimate the iron losses in magnetic circuits. Three models have been selected herein. The first model is based on a global iron loss calculation according to the classic Steinmetz (St) equation [START_REF] Steinmetz | On the law of hysteresis[END_REF]:

max Losses se Ir k f B   (1)
Eq. ( 1) is an empirical formulation requiring three coefficients to quickly and roughly estimate the iron loss. The second model called M4+ uses loss decomposition into three terms. The first term represents hysteresis losses (major loops of ∆Bpp pic to pic flux density variation), the second one represents eddy currents and the third one is introduced to consider minor hysteresis loops with  of B variation [START_REF] Ployard | Iron loss computation of a hybrid excitation synchronous generator[END_REF]. As the previous model, the M4+ model includes three adjustable coefficients (kh1, kh2, αp) to characterize the magnetic loss under variable flux density waveform.

    2 2 1 2 0 2 1 2 1 1 T Losses h pp h pp p n h i h i i dB Ir k B k B f dt T dt k B k B f                 (2)
The third model is the so-called LS (Loss Surface) model [START_REF] Messal | Iron Losses in Electromagnetic Devices: Nonlinear Adaptive MEC & Dynamic Hysteresis Model[END_REF], it is considered herein as the reference model. It is a dynamic H(B) model relating the applied excitation field H to the average magnetic flux density B in the cross section of the material sheet. This model precisely evaluates the iron loss but requires some specific measurements (under a controlled triangular waveform B(t) and variable frequency) introduced in the form of a dynamic surface depending on B and dB/dt. In order to describe the macroscopic behavior of the material, a static hysteretic model [START_REF] Messal | Advanced approach for static part of loss-surface iron loss model[END_REF] is introduced in the main formulation [START_REF] Xue | Thermal-Loss Coupling Analysis of an Electrical Machine Using the Improved Temperature-Dependent Iron Loss Model[END_REF].

      , , , stat dyn H B dB dt H B history H B dB dt   (3)
To predict the 2D vector excitation behavior observed in some parts of the machine magnetic circuit, the LS and M4+ models use the magnetic flux density decomposition according to two axes Bθ and Br that are projected on a new orthonormal coordinate [START_REF] Messal | Iron Losses in Electromagnetic Devices: Nonlinear Adaptive MEC & Dynamic Hysteresis Model[END_REF]

B. Principle of the iron loss machine computation

The MEC discretization principle is shown in Fig. 3. The computational domain is discretized according to the tangential and radial directions [START_REF] Benlamine | Nonlinear Adaptive Magnetic Equivalent Circuit of a Radial-Flux Interior Permanent-Magnet Machine Using Air-Gap Sliding-Line Technic[END_REF][START_REF] Messal | Iron Losses in Electromagnetic Devices: Nonlinear Adaptive MEC & Dynamic Hysteresis Model[END_REF]. Each Block of the MEC contains four branches and thus gives four cells, having two flux density components. The iron loss is computed in each spatial discretization cell. The principle is well detailed in [START_REF] Messal | Iron Losses in Electromagnetic Devices: Nonlinear Adaptive MEC & Dynamic Hysteresis Model[END_REF].

C. Computing time

The main advantage of the developed multi-physical model is the reduced calculation time for both magnetic and thermal models. Table I shows time calculation for the different models.

The nonlinear adaptive MEC is relatively fast with the used spatio-temporal discretization level. The thermal model SAME without thermal transition has its own spatial discretization. The model is so fast. Calculation times of iron losses models depend on the spatial discretization of MEC. In this case, and with the used LS model version, loss calculation time is much higher that of the MEC, while it is insignificant for the Steinmetz model.

III. APPLICATION TO A TRACTION MOTOR

The study has been done on a radial-flux PMSM with 24slots/16-poles. The stator has 3 phases with double-layer concentrated winding distribution supplied by sinusoidal current waveform with a maximum amplitude Ismax. The topology of the machine is shown in Fig. 7. The Iron losses within the machine stator are studied from both local and global points of view. The rotor losses are neglected.

A. Local point of view

A study has been conducted on iron loss calculation in the mesh cells with LS, M4+ and Steinmetz models. Fig. 5 shows the iron loss calculation inside the stator yoke and for the same operating point using the three models for the 56 cells in this part of the teeth. The models M4+ and Steinmetz lightly underestimate the loss calculation in comparison to the LS model. The calculation error remains low. However, the behavior is different for other cases. For example, Fig. 6 illustrates the loss calculation using the three models in the case of the 75 Nm@6,000 rpm operating point (high speed, low torque) and in the 8 cells of tooth-tips. The result shows that the M4+ model overestimates the loss and while the Steinmetz model underestimates it in this same case. For this operating point, the error becomes significant.

It is possible to study the reasons behind these calculation differences by focusing on the flux densities time variations in the cells as illustrated in Fig. 6 (right side), and by analyzing more closely the calculation principle of each model. However, this analysis should be achieved at a local, the aim being to examine not the errors sources but their effects on the coupling and the final temperature distribution in the machine.

B. Global point of view

Table Ⅱ gives the local loss computed, at 160 Nm@1,000 rpm operating condition, in each part of the stator for one pole as well the global one reported for the whole machine (8 poles). The total loss difference obtained is 16 W and 6 W between the Steinmetz and the LS model and between LS and M4+ models respectively.

Table II. Iron losses calculation at 160 Nm@1,000 rpm. Fig. 7 shows the mass specific iron loss distribution within different parts of the stator. This result is obtained using the LS model for 160 Nm@1,000 rpm. Iron loss is focused in the low part of the tooth-tips.

Iron losses [W]

IV. IMPACT OF IRON LOSS MODELS ON THE MOTOR PERFORMENCES

We should not conclude on the accuracy of the loss models based on a single operating point. For electric/hybrid/fuel cell

Bθ [T]
vehicle application, extended the whole torque-speed variation should be investigated. It is necessary to evaluate the iron loss proportion against the Joule loss to bring to light how relevant it is to consider these iron losses for a correct evaluation of the machines thermal behaviour.

A. Results on chosen operating points

In order to study the machine within the torque-speed plane, it is necessary to calculate the machine control parameters for each operating point. Considering the case of a maximum torque per ampere (MTPA) classical control, the optimal control consists in identifying the maximum current amplitude Ismax and the control angle ѱ minimizing the current for a given speed Ωref and torque Temref. This is should be performed under a voltage constrant Veff and bounded ranges for both Ismax and ѱ. It is therefore a mono-objectif non-linear constrained optimization problem, that could be written as in [START_REF] Bui | Core losses measurement technique for high frequency and flux density of switched reluctance machines[END_REF]. A gradient method is an effective solution for this kind of problems, which are relatively smooth. In this example, the current minimization Ismax is equivalent to the Joule loss minimization.

I I I V V T T with                
Fig. 8 shows the ratio in % of the iron losses to the Joule losses for eight operating points of the studied motor numbered from 1 to 8 and distributed on the torque-speed plane. These results are obtained using MTPA control. The iron losses are calculated using the LS model. The iron losses can be neglected in the case of point no.1 i.e., at low speed and high torque. However, it becomes important for the other points, especially at low speed and low torque, and in the weakening regime at high speed. The advantage of the developed multi-physical model is its ability to evaluate the temperature distribution in the machine for the steady state. Indeed, the model inputs are the initial temperatures of the active winding part T0(1), the winding heads T0(2) and the magnets T0 [START_REF] Xue | Thermal-Loss Coupling Analysis of an Electrical Machine Using the Improved Temperature-Dependent Iron Loss Model[END_REF]. They are used in the resolution of the MEC magnetic model. So to evaluate the steady state final temperatures, an iterative procedure or an optimization process should be added. An iterative loop has been introduced in the MTPA problem, minimizing the difference between the initial (T0) and final temperatures (Tf) in the thermal model. An overview of the optimization problem MTPA_T is shown in (5): (5) explicitly shows the two nested loops. Nesting these loops increases the running time.

                  max max 0 0 0 max 0 , , 0
Fig. 9 illustrates iron and Joule losses rates, under MTPA_T control and for the same eight operating points studied before. The point no.1 becomes unattainable. Iron and Joule loss proportions have also changed. For instance, for the point no. 8 the iron loss percentage is decreased from 32% to 16%. The evolution of the temperature conditions change, thus the values and the distribution of global losses. Fig. 10 and Fig. 11 show respectively iron and Joule losses for the eigth operating points and the two control cases. Therefore, the control does not impact the iron loss calculation but rather the Joule losses, especially for the points no. 4 and 8. The point no.1 has the highest torque and current, which leads to an important increase of the temperature and thus a reduction of the PMs flux. Given the imposed current limit, the desired torque could not be reached.

B. Results in term of temperature

Starting from the MTPA control previously defined, loss models have been compared in terms of loss and temperature prediction an observed difference in temperature being directly linked to an error on the loss. An important error might lead to an important increase of the temperature. Conversely, a low error might increase the temperature beyond the tolerated limit. However, the iron loss becomes overestimated in some regions of the torque-speed plane and still remain underestimated in others. Fig. 16.a shows the impact of the relative error on the computed temperatures when both Steinmetz and LS models are used. The loss model influence is analyzed through the temperature error. The NRMSD of Fig. 16.a is 3% and the RMSD of Fig. 16.b is 15°C. However, temperature discrepancies are significant and reach -42°C at high speed while they remain acceptable in the constant torque area. These discrepancies significantly increase in the weakening mode zone where iron loss prevails. Notice that the same comparison has been achieved with the M4+ model. The results show that the M4+ is more accurate than the Steinmetz one and provides an RMSD of 117 W for loss leading to 7°C deviation in temperature.

V. CONCLUSION This article presents the impact of the iron loss model on a multi-physical coupling developed for sizing by optimization of electrical machines. Section Ⅱ describes the multi-physical approach and three iron loss models. Section Ⅲ presents a comparison from a local and global point of view in the case of a radial-flux PMSM. These results provide an idea on the iron loss distribution within the machine and show that it is not sufficient to compare models for just one operating point. Section Ⅳ develops the iron loss calculation impact on the machine performance evaluation. Firstly, a study was conducted at serval operating point about the control and the temperature impact on the computation of the iron loss as well as their relative importance in relation to Joule losses. Then a comparative study has been carried out in the whole torquespeed plane regarding the predicted iron losses and their influence on the evaluated thermal behavior of the machine.

Fig. 1 .

 1 Fig. 1. Multi-physical model for sizing by optimization.

  . The Steinmetz model, calculates the loss based only on the flux density module. The M4+ and Steinmetz model coefficients were identified referring to the LS model under sinusoidal waveform for different frequencies and flux density amplitudes. Fig. 2 shows the iron loss IrLosses(Bmax) curves for different frequencies. The relative error is mainly low for important flux density but might become more important for low flux density values.

Fig. 2 .

 2 Fig. 2. IrLosses(Bmax)curves for different frequencies under sinusoidal waveforms.

Fig. 3 .

 3 Fig. 3. Spatial discretization principle.

  Fig. 4.a and Fig.4.b show respectively the results obtained one stator teeth within the cells with Steinmetz and M4+ models for the operating point 160 Nm@1,000 rpm. The models are quite similar with few differences in amplitudes.

Fig. 4 .Fig. 4 .

 44 Fig. 4. a. Steinmetz Fig. 4. b. M4+ Fig. 4. Iron loss distribution in a stator tooth.

Fig. 5 .

 5 Fig. 5. Iron loss calculation in the 56 cells of the stator yoke part of the tooth shown in Fig.4.

Fig. 6 .

 6 Fig. 6. Loss calculation in a stator tooth-tip cells.

Fig. 7 .

 7 Fig. 7. Specific iron losses in (W/kg) distribution at 160Nm@1000rpm.

Fig. 8 .

 8 Fig. 8. Bubblepie of loss percentage for MTPA.

Fig. 9 .

 9 Fig. 9. Bubblepie of loss percentage for MTPA_T.

Fig. 10 .

 10 Fig. 10. Iron losses for MTPA and MTPA_T.

Fig. 11 .

 11 Fig. 11. Joule losses for MTPA and MTPA_T.

Fig. 12 .

 12 Fig. 12. Iron losses cartography.

Fig. 13 .

 13 Fig. 13. Joule losses cartography.

Fig. 14 .

 14 Fig. 14. Efficiency cartography.

Fig. 15

 15 Fig. 15.a shows the relative error on iron losses made by using the Steinmetz model in comparison with the LS model. Fig.15.b illustrates the absolute error, which means the difference (in Watt) between the two model losses. Errors are important for this model with the realized identification. To effectively give a ruling about the models in a global point of view, the normalized root mean square derivation (NRMSD) and the root mean square derivation (RMSD) have been calculated respectively for the relative and the absolute errors. The obtained values are respectively 58% and 248 W. In this example, Steinmetz model underestimates the iron loss. A purely mathematic global correction of the Steinmetz model has been tested by multiplying the iron loss by a correction factor of 2.02 obtained by optimization procedure. The obtained NRMSD is then reduced from 58 to 28%.

Fig. 15

 15 Fig. 15. a. NRMSD[%]: Steinmetz Vs LS. Fig. 15. b. RMSD[W] : Steinmetz Vs LS. Fig. 15. Error of the Steinmetz model in comparison with LS model for iron losses.

Fig. 16

 16 Fig. 16. a. NRMSD[%]: Steinmetz Vs LS.

TABLE I .

 I Calculation times.

			Computing time [s]
	MEC	LS	M4+	Steinmetz	SAME
	3.8	56	0.1	0.003	0.8

ACKNOWLEDGMENT

This work was supported by RENAULT-SAS, Guyancourt, France, through the program COCTEL (Conception Optimale des Chaines de Traction ELectrique) financed by the French Agency for Environment and Energy Management (ADEME).