Production planning of a furniture manufacturing company with random demand and production capacity using stochastic programming

In this article two multi-stage stochastic linear programming models are developed, one applying the stochastic programming solver integrated by Lingo 17.0 optimization software that utilizes an approximation using an identical conditional sampling and Latin-hyper-square techniques to reduce the sample variance, associating the probability distributions to normal distributions with defined mean and standard deviation; and a second proposed model with a discrete distribution with 3 values and their respective probabilities of occurrence. In both cases, a scenario tree is generated. The models developed are applied to an aggregate production plan (APP) for a furniture manufacturing company located in the state of Hidalgo, Mexico, which has important clients throughout the country. Production capacity and demand are defined as random variables of the model. The main purpose of this research is to determine a feasible solution to the aggregate production plan in a reasonable computational time. The developed models were compared and analyzed. Moreover, this work was complemented with a sensitivity analysis; varying the percentage of service level, also, varying the stochastic parameters (mean and standard deviation) to test how these variations impact in the solution and decision variables.


Introduction
Carrying out an aggregate plan is important in manufacturing industries, especially those where it is planned in periods of 3 to 18 months, or medium term. The production plan seeks to determine the optimal levels of production, hiring, layoffs, inventories, subcontracting, etc. This work presents an aggregate plan that was made for a company that manufactures furniture in the State of Hidalgo. Initially, a first approach to the solution of the problem was made in [1]. In this work, only the production capacity was considered as a random variable using two models, one with a continuous probability distribution and the other with a discrete one. The contribution of this work is a real problem where uncertainty affects the production system, generally, the models used in the literature consider demand as a random variable with a discrete approximation (one model), in this work, in addition, the human factor is considered as a stochastic parameter that can be modeled and 2 models are compared. The rest of the article consists of the following sections. Section 2 presents the literature review, section 3 describes the model under study, first as deterministic and then as stochastic. The methodology is introduced in section 4. The results are given in section 5. An extensive sensitivity analysis is performed in section 6. Section 7 concludes the article.

Literature review
Considering uncertainty within optimization problems remains a trending topic to be investigated because organizations face to stochastic variables when making decisions. A search was carried out in Scopus and in Web of Science written during 2020 and 2021 that used stochastic programming to solve optimization problems. Huang et. al [5] develop a multistage stochastic optimization model for system operators to efficiently schedule power-generation assets to cooptimize power generation and regulation reserve service under uncertainty. Ghayour et al. [6] present an approach called MLPR with linear programming used as its core in order to solve the influence maximization problem in the linear threshold model, that is one of two classic stochastic propagation models that describe the spread of influence in a network. Robust Multi-product Newsvendor Model with Substitution, where the demand and the substitution rates are stochastic and are subject to cardinality-constrained uncertainty sets that is an NP hard problem is presented in [7]. Also, Basciftci et. al [8] reformulate the robust facility location problem, in which they interpret the moments of stochastic demand as functions of facility-location decisions. In Shone et. al [9], stochastic modeling applications within aviation are presented, with a particular focus on problems involving demand and capacity management and the mitigation of air traffic congestion; using operations research perspective, including analytical queueing theory, stochastic optimal control, robust optimization and stochastic integer programming. Ghasemi et. al [10] present an Evolutionary Learning Based Simulation Optimization (ELBSO) method embedded within Ordinal Optimization. In ELBSO a Machine Learning (ML) based simulation metamodel is created using Genetic Programming (GP) to replace simulation experiments aimed at reducing computation; ELBSO is evaluated on a Stochastic Job Shop Scheduling Problem (SJSSP). Zhang et. al [11], consider a stochastic vehicle routing problem with probability constraints; the probability that customers are served before their (uncertain) deadlines must be higher than a pre-specified target. Wang et. al [12] propose a model to solve a project scheduling problem where resource assignments and activity schedules need to be determined to achieve a set of due-date requirements as well as possible. Torres et. al [13] present multistage stochastic program for the design and management of flexible infrastructure networks with stochastic demands.
In the methods for solving stochastic programming, Dowson and Kapelevich [14] develop the Julia package for multistage stochastic and dual programming and Gangammanavar et. al [15] work with stochastic decomposition for two-stage stochastic linear programs with random cost coefficients.
Talking specifically of Aggregate planning, there are some works related with it, Jamalnia et al. [16] mention that the methodologies applied to deal with aggregate production plans (APP) under uncertainty can be classified into six main categories: stochastic mathematical programming, possibilistic programming, fuzzy mathematical programming, simulation modelling, metaheuristics, and evidential reasoning. Here are some important works.
Using Multi objective stochastic optimization, Nowak [17] introduced a work that combines linear mathematical programming with multiple objectives (multi-objective), simulation and an interactive approach with uncertain demand. In contrast, Jamalnia et al. [16] presented a nonlinear stochastic optimization model with multiple objectives for an aggregate production plan under uncertainty. The WWW-NIMBUS software was used, using more than 500 decision variables and 1,000 restrictions with 7 objectives. Zhao et al. [18], showed a case study for the pharmaceutical industry where they optimize the quantity of production, minimizing the duration of clinical trials and operating costs. The problem is formulated by the multistage stochastic programming using C# programming language as well as CPLEX as a solver, comparing the Progressive Hedging Algorithm (PHA) and the Sample Average Approximation (SAA) to test the optimality GAP, the time to solve a scenario and period, the memory used, among others. Rakes, Franz, and Wynne [19] and Chen and Liao [20] also submitted works where there is uncertainty in production plans with multiple objectives. Recently, Tirkolaee et al. [21] investigates a novel fuzzy multi-objective multi-period Aggregate Production Planning (APP) problem under seasonal demand. As two of the main real-world assumptions, the options of workforce overtime and outsourcing are studied in the proposed Mixed-Integer Linear Programming (MILP) model. Stochastic optimization has been applied in various production plans where the environment tends to have uncertainty as mentioned by Birge and Louveaux [22]. Wagner and Whitin [23] are considered the first to study production plans under uncertainty, solving a one-state problem through dynamic forward programming for a single product. Kazemi et al. [24] perform a multi-stage stochastic mixed integer linear programming model for sawmill production (lumber industry) in Canada; giving a solution through an approximation by decomposition of scenarios. CPLEX software was used for the solution. It was shown that problems with more than 22,000 restrictions and 44,000 decision variables can be solved within a reasonable time, this was a single objective problem. Huang [25] also presented models for production plans through multi-stage stochastic mathematical programming.
Using, nonlinear stochastic optimizatiom, Ning et al. [26], Mirzapour Al-e-hashem et al. [27] and Lieckens and Vandaele [28] developed mixed integer nonlinear mathematical programming methodologies to study the decision problem of an aggregate production plan and they consider demand and lead time as variables with uncertainty in their proposals. Nasiri et al. [29] suggested a non-linear stochastic model for a production and distribution plan in a three-level supply chain (suppliers, production centers and customers). This model was solved by using commercial packages such as Lingo in addition to a heuristic programmed in Matlab. Similarly, Nasiri et al. [30] proposed a nonlinear stochastic integer mixed mathematical programming model for a supply chain to later extend their work in a multi-stage system [31]. Ning et al. [26] presented a multi-product nonlinear application model where market demand and production cost are uncertain.
A key concept in stochastic programming is scenario tree. Kazemi et al. [24] mentioned that a scenario tree is a computationally viable way of discretizing the dynamic stochastic data underlying a problem over time, this allows to reach viable solutions in reasonable times. In a competitive environment it is required to have results in faster times. Hu and Hu [37] used a scenario tree to optimize the job shop scheduling problem (JSSP) and the economic order quantity (EOQ) under uncertain demand, Körpeoglu, Yaman and Aktürk [38] used it to solve the master production scheduling problem (MPS). In Table 1 the studies that closer to our proposal are compared in different characteristics.
Reviewing the works found that are related to stochastic programming, it was observed that all of them use mixed integer programming, Jamalnia et al. [16], Zhao et al. [18] and Tirkolaee et al. [22] with nonlinear multiobjective, while Kazemi et al. [24] uses multistage stochastic programming which is the same as that used in our proposal. Tirkolaee et al. [22] use demand and costs as stochastic variables, while Kazemi et al. [24] demand and yield; the others, only use a single variable as a stochastic. In their approaches, a single approximation is used to explain the uncertainty, which is discrete, that is, a single stochastic model, and in our case two models, that allow us to compare between the normal distribution and its discretization in order to offer to the company a good solution in a reasonable computational time. The proposals found use the scenerio tree, except Tirkolaee et al. [22] that deal with the problem with weighted goal using GAMS. The level of service is only handled in Zhao et al. [18], but the impact of the service level on the solution is not considered. The sensitivity analysis varying stochastic parameters was not used in the approaches found. After analyzing the characteristics of the studies found, the gap was in considering the level of service, which was a very important restriction for the company in the case study, considering demand and labor as stochastic variables, which were two variables that generate a lot of uncertainty in the company and also varying the stochastic parameters. Finally, generate an efficient model, that is to say, a model that obtains a good response in a reasonable computational time, for that reason 2 models were tested for comparison. In the next section the development of mathematical models will be described.

Development of mathematical models
To develop the stochastic programming models a proposed methodology was used, where before developing the stochastic programming models directly, a first deterministic model was created, in order to observe the behavior of the model and after extended it by incorporating the stochastic parameters.

Deterministic model
The deterministic base model for the aggregate production plan for optimizing the total cost of production, determines the number of workers that must be hired and fired per month, the number of parts that will be produced, the parts that will be sent to the finished product warehouse (inventory) and the backlog per month with a service level of 90% each period. The model considers the following assumptions that are considered in the literature review and some others are specific for the company: ◾ The demand is known for all periods.
◾ The production capacity per worker is the same for all months.
◾ Backlog may exist in the company, but it is penalized for missing units, this is because there is no possibility of buying parts (outsourcing) and overtime is not allowed in the company.
◾ Backlogged parts must be satisfied in the next period.
◾ The costs associated with the production of a part, inventory and backlogs are linear.
◾ There is enough raw material for production of parts (X t ).
The notation used for the model is shown in the following Table 2. Where the indices, parameters and decision variables are detailed, later the objective function and the model restrictions are explained. The deterministic base model used was based by the one in Ravindran and Warsing [39] adapting it to the needs of the company by incorporating the service level and minimum inventory restrictions. The objective function of the deterministic model in Eq (1), minimizes the total cost of the APP, considering the costs associated with the workers assigned to production, fired, hiring, as well as inventory costs, backlogs and production, as defined in the following equation: The model constraints are the following: Constraint (2) focuses on the size of the workforce in the company, and it indicates that the total number of workers in period t, it must be equal to those existing in period t−1, plus those hired in period t−1, minus those laid off in period t−1. Constraint (3) is about the allocation of the workforce; it simply defines how many workers will be assigned to production and the number of workers that will be laid off in period t. Constraint (3) complements the balance over the number of workers denoted in Eq (2) and specifies the assignment of workforce in production and the number of workers to be fired each month. Constraint (4) refers to the balance of demand and inventory in the company, where what is produced in period t plus the inventory of period t−1, it must be equal to the demand of period t, plus the backlog of period t−1 plus the inventory in period t minus the backlog of period t. Constraint (5) addresses the production capacity; it ensures that the workers assigned to production can manufacture the units required in period t. Constraint (6) defines that the inventory of the period is greater than the safety inventory defined by the company. Constraint (7) indicates that the service level is greater than or equal to 90% per month, this restriction is a company's policy. The expression (8) defines the constraint of non-negativity decision variables. The expression (9) specifies that the decision variables must be integers, which implies a mixed integer linear stochastic programming model. The deterministic model is the basis for developing Model-I and Model-II, both are of the multi-stage stochastic programming type, where production capacity and demand are random parameters. The Model-I associates these random parameters to a normal distribution, so the solution strategy is to use the stochastic programming solver integrated in Lingo 17.0 for Model-I, which internally creates a scenario tree with the same probability of occurrence for each scenario. Subsequently, a Model-II is proposed through a discretization of the probability distributions to create a scenario tree, each with different probabilities.
Model-I has the same assumptions as the base model, in addition to the following: ◾ Demand follows a normal distribution with mean 353 and standard deviation 29 and production capacity a normal distribution with mean 12 and standard deviation 2. These values were obtained through historical data of the company during 3 years prior to the start of the study using 40 data. The Arena ™ Input Analyzer was used to determ|wine these values. The statistical tests made to the data are shown complete in https://doi.org/10.6084/ m9.figshare.14444609.v1.
◾ When considering a random demand, it is not possible to know how much to produce before the realization of this random event, but it must be determined how many workers are needed in each period. Therefore, the decisions of the same period can be made in different states. In Table 3 there are the parameters associated with this model: Model-II has the same assumptions as the base model, in addition to the following: ◾ The discrete probability distribution is obtained with the maximum likelihood values calculating the probability of some defined intervals, then the Gaussian quadrature method is applied to calculate the probabilities (areas) of the intervals. ◾ Similar to the Model-I, some decision variables of the same period could be done in different state due the stochastic process.
In both cases, because the models use a scenario tree, the equivalent deterministic model for each problem consist of the construction of a large mixed integer linear programming Table 3. Parameters in Model I. Own elaboration based in the organization information. problem, where at the end, the probability of the sum of each scenario should be one as indicated Eq (10). The notation for Model-II is detailed in Table 4, Model-I is similar in structure to Model-II, the difference is that Model-I, the probability of occurrence for each scenario is the same, in addition the occurrence values are obtained by Lingo using identical conditional sampling and for this, it uses variance reduction techniques to improve the quality of the solution.

Model I-Parameters:
The objective function of the stochastic Model-II (11), as in (1) minimizes the total cost of the APP, considering the costs associated with the workers assigned to production, fired

P ij
Probability of occurrence of random events {ω,δ}, where the following equation is obtained: Which indicate that the probability sum of all scenarios across the scenario tree must be one as indicated Eq (10). workers, hired workers, in addition to inventory costs, backlogs and production.
The model constraints of the stochastic model are the following: Constraints (12) and (14) are zero state (or first state) constraints, these are considered before the occurrence of the random event and are equivalent to the constraints (3) and (5) of the deterministic model. Constraint (13) is equivalent to constraint (4) only for the first period. Constraints (15) to (20) are the same of (2) to (7) of the deterministic model, considering the random events ω, δ the necessary corrections are made in order to deal with the uncertainty, these variables are called recourse variables.
Constraints (21) and (22) indicate the nonnegativity constraint and assigning some variables as integers for the zero state, it should also be noted that not all variables are forced to be integers, a computational advantage is obtained and makes the model of the mixed integerlinear class. Constraints (23) and (24) indicate the non-negativity constraint and assigning some variables as integers for the next states. The nonanticipativity constraints (25)  where ; and are the column vectors for the decision variables in the scenarios for the production capacity production low, medium and high respectively, and it is indicated if the demand is low when they do not have a quote, medium demand if they have one quote ( 0 h ; 0 h ; 0 h ) and double quote if the demand is high (, @ h ; @ h y @ h ). In Fig 1 it can be observed how for a period nine scenarios are generated considering 3 possible scenarios for each random variable, each one with its respective probability. Then for the first node another nine scenarios are generated, that is, for two periods the scenario tree grew up to a total of 81 scenarios. The sum of all the probabilities of the nine scenarios must be equal to the probability of the predecessor node to those scenarios, and in the same way the sum of all the 81 scenarios must be 1 and the filtering process of the sigma algebras is respected, which implies respecting the nonanticipativity constraints.

Methodology
The methodology used in this article was divided into five stages. Each stage is briefly described below: 1.-Data collection and analysis: In order to prepare the aggregate production plan, historical company data was used to obtain production parameters and associated costs. An ABC classification of products was carried out, where the main product was determined, this product generates more income for the company. In this case it turned out to be a rustic chair.
2.-Deterministic mathematical modeling: a deterministic model was developed that minimized the production costs of the product. The company's policies to maintain a minimum inventory every month and have a service level of at least 90%, both were considered in the model. The solution strategy for this problem was to implement the Monte Carlo technique of identical conditional sampling. Lingo uses identical conditional sampling and calculates the expectation for each state in a similar way to SAA, associating the random variable "Production capacity" and "demand" to a normal distribution with known distribution parameters. Additionally, the variance reduction technique is carried out using Latin-hyper-cubes in order to improve the solution.
The sample size is obtained in such a way that it does not generate an excessive computational effort in time and iterations and that in the solution obtained is still of quality. However, the use of variance reduction techniques in numerous studies (Kleywegt et al. [40]; Verweij et. al [41]; Shapiro & Homem-de-Mello [42]; Ruszczynski & Shapiro [43]) has shown that the number of samples to require a solution with the same quality as using a simple sampling is 10 to 100 times less, so the use of 3 samples is justified in this case. 5.-Resolution and comparison of the models: Both stochastic models were processed to solve from two periods to the maximum number of periods that the computer was able to process. An optimization gap between the models was calculated. Also, a sensitivity analysis was performed, the impact of the percentage of service level and the parameters of probability distributions were analyzed.
Once the robustness and structure of the solution obtained has been confirmed, finally the optimal values of the related aggregate plan production will be prepared to be implemented in the manufacturing system. It is worth mentioning that this phase will be outside the domain of this study and remains as a possible future research direction. The lingo models are in https:// doi.org/10.6084/m9.figshare.14450430. The work methodology is summarized in the Fig 2, in a flow chart.

Results
For the solution of the problem, a Dell Inspiron 5570 computer with an Intel1 Core™ i5-8250U processor with 1.6GHz, with a 4Gb RAM memory with Windows 10 operating system was used. The Lingo 17.0 software was run using the Lingo stochastic solver, this solver uses an identical Monte Carlo sampling and Latin-hyper-cubes techniques to reduce the variance of the samples when a continuous distribution is used.
When the discrete distribution was implemented, only the empirical distribution and its probabilities were indicated to Lingo. It should be mentioned that in both cases the advantage that Lingo has is that by indicating the model as stochastic, it is unnecessary to indicate nonanticipativity restrictions, Lingo performs them internally.
For the solution of the equivalent determinists, five of the six available cores of the parallel processor were used (this option can be registered in the Lingo options), for the prior relaxation of the problem, it was specified in the linear solver to use 4 cores, each one with different strategies to solve the problem (two primaries using the simplex method, one barrier and the other dual).
For the mixed integer pre-solver, the maximum amount of heuristics (100) that Lingo has enabled were used to find optimal pre-solutions, the other capabilities were left as they come by default in Lingo. Finally, the Branch-and-Bound (B-and-B) algorithm was applied in the integer solver with Lingo's default criteria. In the case where a discrete distribution is used, a matrix decomposition was tested to perform a faster solution time.
Approximate results were found regarding the objective function of Model-I with respect to Model-II, which shows that the approximation of the continuous model with the discrete model demonstrated to be efficient when the problem becomes larger. The following comparative table shows how the problem grows when more periods are added in the aggregate production plan. This growth is because the model must generate the equivalent deterministic models for each branch of the scenario tree. Table 5 shows how the problem grows as the number of periods increases, making it more difficult to find a solution to the problem if it is assumed that now the random variables are The number of iterations continues to be significantly greater in Model-I with respect to Model-II, therefore the time required to solve the problem is greater in Model-I as observed in Table 6.
From Table 6, as the problem increases in periods, the time taken to find the solution of Model-I is longer. In all cases, a global optimum was found for Model-II, for Model-I only a feasible solution was found in period 4. Moreover, an advantage can also be observed from the perspective of the number of iterations and CPU Time after 3 periods. number of iterations and CPU Time after 3 periods. Fig 3 reports the Expected Value (EV) indicator, which is the expectation of all scenarios, Lingo reports this value also as the value of the objective function. The values are close (see Fig  2) and this causes a low gap between both models, so from a computational perspective the Model-II is convenient, even more if the interest is to seek quick decisions in organizations. Fig 4 reports the wait-and-see (WS) indicator, this reports the expected value of the objective function by removing the nonanticipativity constraints. In practice, this is impossible due to the randomness cannot be anticipated. The values tend to be very similar; this is because the approximation using the maximum likelihood values of the normal distribution, minimize the error of the approximation.
Generally, the use of heuristic and computational techniques to approximate the solution is a motive for research, seeking speed with the cost of a solution. From Fig 3 it can be seen that from the perspective of the WS solution the proposed model (Model-II) closely approximates to the real problem unlike other studies (see [29,45]). Fig 5 reports the indicator perfect Information of expected value (EVPI) is the absolute value of the difference between EV and WS. It is the maximum amount that would be paid to obtain information [22] and reduce randomness. It should be mentioned that the existence of a high EVPI justifies the use of a stochastic model, rather than using a deterministic model. Table 7 summarizes what was mentioned about the EV, WS and EVPI stochastic optimization indicators, adding the EV GAP. Note that the GAP obtained is low, so both models are similar in their result of the objective value as already mentioned above. These results show the efficiency of the proposal Model-II even when there are two random variables, it can be noted that when the random variables are associated with a normal distribution with known parameters and in a certain way compact (i.e. standard deviation and mean close), it allows the discretization and generation of the scenario tree to be efficient computationally.
When continuous distributions are approximated with discrete distributions, it can be observed that it affects the quality of the Model-II solution, despite this drawback, the quality of the solution is still acceptable in terms of the EV GAP, in this case about 5%.

Sensitivity analysis
In this section a sensitivity analysis is carried out to observe how robust the objective function is, and the decision variables are when costs are changing. The problem with 3 periods was taken as a reference for both models. When performing this sensitivity analysis (Tables 8 and 9), it can be observed that there are parameters that greatly affect the objective function, such as production costs C p and product production costs C X , if it varies towards positive and negative percentages, and the firing cost C F , when it varies negatively, however, there is a small difference in terms of the number of workers assigned to production in the first period P 1 (first state decision).
In general, the variation of the objective function is similar between both models except for the variation in C p which implies a greater increase in Model-I with respect to Model-II, and  in C F which results have greater reduction of Model-II with respect to Model-I. One consideration is the hiring policy, Model-I considers hiring in 11 of the 12 sensitivity results, this because it must be considered that the demand, from a statistical perspective, has infinite achievements. In order to deal with this, the resource variable evaluates how many workers to hire.
On the other hand, Model-II considers that there are only three possible realizations, therefore its evaluation from a computational perspective is less complex. As there are not infinite realizations, the corrective actions by the resource variables are minimal. Another consideration is when C F is null, it should be noted that the model prefers to hire and assign many workers P 1 and then carry out the fired workers, this suggests that an outsourcing policy would mean a viable strategy for the company in a situation where the behavior of demand is known.
For the cost C I , C R and C S there were no large significant changes in the objective function, but reducing C R , it is possible to have a smaller workforce in the first periods, to later carry out hiring without having an additional cost, thus indirectly reducing inventory costs and production.
It is clear, that the proposed Model-II turns out to be in general very similar as mentioned previously, perhaps it would be worthwhile to analyze in detail the hiring policies obtained through the variation of parameters, particularly in economic situations where companies see their capital reduced and they are forced to decrease payments or costs could increase.

Impact of service level on Model-I and Model-II
In this section, an analysis on the impact of the service level in Model-I and Model-II is developed, considering the same scenarios and number of periods previously used, starting from 86% to 98% of service level. Because the service level restriction is important in the model, the  Table 7. Comparison between Model-I and Model-II with respect to the importance indicators of stochastic optimization, " � " indicates approximation. Own elaboration. sensitivity should be analyzed in the same way as another variation parameters of the model. From Tables 10 and 11 it can be observed that, at a lower service level of Model-I and Model-II, the objective function is lower (Expected Value). The EV solution gradually increases as the service level increases. For the indicator EVPI, this increase does not occur, in the case of Model-I a lower EVPI implies that the model is not so affected by random issues, therefore, managing a service level of 90% is viable for the Model-I. However, in Model-II we find that with a service level of 92%, the model does not show positive increases, hence, managing this level of service allows reducing the EVPI, but this company policy implies a higher cost due to the increment in the EV.

EV
It can be noticed that the best service level policies are obtained with low service levels, which can improve EV and EVPI, but this would imply too many backlogs that could leave a not unfavorable impression of the company. The other option is to manage a service level policy of between 90% and 92%, being able to improve EVPI and without a drastic increase in EV, however, the cost is higher, even though the level of service is good, and that is easy to agree between company and clients. Tables 12 and 13 show that the level of service significantly affects the decision variables. It is evident that as the level of service increases, more workers are assigned to production, and with this measure the company produces more units. This also increases the inventory level; with these countermeasures the backlogs decrease.

Impact of variation in the distribution probability parameters of the random variables
Finally, an analysis on the impact of variation in the distribution probability in Model-I and Model-II is developed, considering the same scenarios and number of periods previously used in the las two sections, two different cases were analyzed, the first where the mean has different values, but with the same standard deviation and the second case, where the mean is the same, but the standard deviation has different values. Figs 6 and 7 explain what has been done. Table 14 presents the sensitivity for the normal distribution parameters (mean and standard deviation) of the two random variables (production capacity and demand). Cases 1 to 3 and 7 to 9 show how the behavior is when the mean is fixed, and the standard deviation varies for both random variables.
Cases 4 to 6 and 10 to 12 show the sensitivity achieved when the mean of both probability distributions is varied, leaving the standard deviation corresponding to each random variable fixed. It is evident that the best results regarding EV are obtained in cases 1, 6, 9 and 10, which are when the standard deviation decreases for production capacity or its mean increases, as well as when there is a high variation in demand or its average low. For the decision variables, the greatest changes occur when the parameters of the random variable demand are varied, in case 12, particularly for the variables P 1 , S 3 y X 1 . In the case of the random variable production capacity, its behavior is like that found in Model-I, that a lower variation in production capacity or that its average is higher leads to better results.
Comparing both random variables, varying the standard deviation of the production capacity has more impact on the decision variables and on the model in general, but if the mean is varied, the demand has more impact, this because the cost of the production plan will be proportional to the demand. Table 15 is the counterpart of Table 14, where the parameters of the discrete distributions (their values and probabilities of occurrence) were obtained using the Gaussian quadrature method proposed in the methodology. It can be observed in Table 15 similar behaviors in terms of the variation of the parameters of the probability distribution for production capacity, than those observed for Model-I.
It should be noted that the best results with respect to the EV indicator occur again in cases 1, 6, 10 and now in case 7. For case 10 the changes in the EV indicator They are explained because less products are produced and the inventory reduces compared to the base case, this is explained by the decrease in demand. Again, the higher the demand, it is clearly observed in the results how the decision variables are affected. These results show the quality of Model-II, using the proposed methodology because the solutions are similar.

Conclusions
This research considers an aggregate production plan applied to a company. Aggregate production plans play an important role in SMEs because they allow them to manage their resources and operations efficiently. This plays an important role in planning as it reduces the risk of SMEs disappearing, particularly in Mexico. Due to various policies such as hiring and firing, the production capacity in the company was not constant, which is why it was initially considered as a random variable, later demand is incorporated as a second random variable. Aggregate production plans under uncertainty have been studied in the literature for different types of problem structures such as linear optimization models, mixed linear-integers, non-linear-mixed integers, among others, in addition to whether they are single-objective or multipleobjective. Although some studies report solutions under the here and now or wait and see solution, both solutions have never been considered in the same study, in addition we incorporated the result of the expected value of perfect information in this work. Model-I could not solve the problem for more than three periods, for which a second model was proposed (Model-II), following a proposed methodology, achieving good results when comparing both models that is, there is no significant difference between the objective  function, with model 2 being the most effective at the computational level. The advantages offered by the methodology proposed is its flexibility, being able to use it in other problems where there is uncertainty. Some of its limitations, are the need for goodness of fit tests that ensure that the data have a certain probability distribution; if they are not done, the results will be far from optimal. Also, the sample size or number of possible realizations of a random variable when it is discretized can make the problem difficult to solve, because the equivalent  deterministic grow exponentially as the number of states increases, for this study; solving a fifth state would imply solving a deterministic equivalent of more than two million decision variables and four million constraints because more than 40,000 scenarios are required, which Approximation to the variation of the mean μ of the capacity production is the maximum number of scenarios that the Lingo software can process without having a computational memory deficit.

Approximation to the variation of the standard deviation σ of the capacity production
The study was complemented with sensitivity analysis, in the literature few studies report these analyzes, generally only perform it by varying parameters associated with costs, but in this study is carried out to see the impact of varying the percentage of the policy of level of service, also, a sensitivity analysis is also carried out, varying the parameters of the probability distributions or stochastic parameters (mean and standard deviation) and evaluate the impacts in the solutions and decision variables, so that the company has sufficient information for correct planning in case these parameters could change in the future, or, as an area of opportunity to improve productivity, for example, it could be observed that reducing the variability of the random variable production capacity, that is, reducing its standard deviation, reduces the total cost of the APP.
Since many of the APPs occupy neither linear functions, a direction for future research is to make a model considering some non-linear functions, such as the inventory cost, also reformulate the problem, removing some restrictions that will allow to have a lower inventory level, allowing to solve the problem in a more efficient way, it is also interested in the use of another algorithm to solve the equivalent determinists, in this work the algorithm B-and-B was used, being able also to use algorithms of cut of plan, consider using multiple kernels in parallel, using multiple heuristics to pre-solve the problem, and using robust algorithms for relaxation of the problem.