## 3 Optimization Method  
The mathematical model presented above is clearly based on simulation. Machines states are simulated, random failures are also simulated according to an exponential distribution law, and time flow is also simulated. From all of this, values of stock levels, number of new and reconditioned products, lost demands, return products quantities, production rates of both machines, limit of carbon emission constraint, production plan, and total profit function are calculated, according to Equations (5) to (20). It is obviously not an integer linear mathematical model, but we have developed an optimization method to cope with the non-linearity and random variables challenges, and managed to obtain a total profit function maximization. The total profit function was programmed in the language C++ with the free software Dev-C++. The optimization method was also developed in this language. In order to facilitate the method handling, we have put the decision variables in a vector. Each value of the vector can be between a minimum and a maximum value. The optimization process is illustrated in Figure 3.  
The model is stochastic by construction. Therefore, we calculate the average of three simulations to obtain a more precise value of the objective function, in order to be able to have the best surface possible to perform optimization. The simulation time horizon is fixed to \(10^{7}\) time units. The optimization method is inspired from evolutionary algorithms and returns the optimal values of the vector which maximizes the total profit function. The quality of the results obtained by this method depends on the number of tests. The larger the number of tests, the better the results. However, the simulation time can also increase too much. Thus, we developed an improvement method based on local search that can give better solutions by testing the function the least amount possible. The steps are the following:  
Initialization: random individuals are created, with each vector of values respecting the boundaries of each of its parameters. Evaluation: In addition to the pre-determined constraints on variables, the following constraint is written: \(S_{A}>S_{Ar}\). This equation imposes that the storage of new products should be higher than that of reconditioned ones in order to ensure the return of worn products. Consequently, unadapted individuals are rejected. Selection: We retain the best adapted individuals, those who give satisfying total profit function results. Mutation: Among the best adapted individuals, we try some neighborhood tests, thanks to a deviation of some values of the vector describing an individual, to see if its characteristics are better adapted. Evaluation: If a new generated individual meets the constraints, it is evaluated with its total profit value. Replacement: If the new total profit function is better than the previous, it is retained. In the end, the final solution is represented by the following vector (\(S_{A}\), \(S_{Ar}\), \(p\)). The second part of the proposed optimization method is composed of local improvement procedures, which are very useful to perform a better optimization [32]. We developed three main procedures. The first one changes the value of a parameter, chosen at random, in one sense or in the other, here again chosen at random (positive or negative). It can go up until the appropriate bound, but the variation is also at random (it may not \"hit\" the bound). The set of parameters is then inside each variable bounds, by construction, but the constraints may be violated anyway. Then, after a constraint checking, its objective function is evaluated. If it improves the current solution, it is kept, and if not, the counter of tested solutions is incremented. We allow to test up to 100 valid neighborhoods. The second neighborhood iteratively repeats the procedure described below, until there is no improvement between two successive iterations.