Application of Multi-Objective Evolutionary Algorithms for Planning Healthy and Balanced School Lunches

: A multi-objective formulation of the Menu Planning Problem, which is termed the Multiobjective Menu Planning Problem, is presented herein. Menu planning is of great interest in the health ﬁeld due to the importance of proper nutrition in today’s society, and particularly, in school canteens. In addition to considering the cost of the meal plan as the classic objective to be minimized, we also introduce a second objective aimed at minimizing the degree of repetition of courses and food groups that a particular meal plan consists of. The motivation behind this particular multi-objective formulation is to offer a meal plan that is not only affordable but also varied and balanced from a nutritional standpoint. The plan is designed for a given number of days and ensures that the speciﬁc nutritional requirements of school-age children are satisﬁed. The main goal of the current work is to demonstrate the multi-objective nature of the said formulation, through a comprehensive experimental assessment carried out over a set of multi-objective evolutionary algorithms applied to different instances. At the same time, we are also interested in validating the multi-objective formulation by performing quantitative and qualitative analyses of the solutions attained when solving it. Computational results show the multi-objective nature of the said formulation, as well as that it allows suitable meal plans to be obtained.


Introduction
Recent findings suggest the importance of fostering healthy nutrition habits in order to improve the quality of life of children [1][2][3]. The increased presence of fast food and snacks in our diet has led to its consumption even when there is no real need related to lack of time [4,5]. In addition, in recent years, the low frequency of physical activity have induced high levels of excess weight and obesity in the population [6]. At the same time, the economic crisis in many developed countries has increased the poverty rates of the population, with the ensuing malnutrition problems, which is especially serious when it arises in children. It is vital, therefore, to have a healthy and balanced diet [7]. Changes in nutritional habits for children require the effort and cooperation of families, educational institutions and governments [1,8,9]. We note that it is a hard task, from the regional administration standpoint, to control the nutritional habits that children have at home. Nevertheless, as children usually have lunch at schools, regional administrations could promote policies in order to improve those nutritional habits at school. In general, a balanced diet should be designed manually and reviewed by nutrition experts. However, by considering general nutritional criteria or recommendations for school-age children, as well as a predefined set of meals, it is possible to correctly formulate the problem of automatically planning school menus.
In this work, we propose a multi-objective formulation of the Menu Planning Problem (MPP) [10], which we will refer to as the Multi-objective Menu Planning Problem (MMPP). Firstly, school canteens are in use mainly at lunchtime, so the MMPP only focuses on recommending lunch menus, thus excluding the remaining daily meals. Secondly, the target audience are children ages 4 to 13. The amount of nutrients recommended for children between 4 and 13 does not differ considerably by gender, according to data from the White Book on Child Nutrition (this document is available at https://www.aeped.es/ sites/default/files/documentos/libro_blanco_de_la_nutricion_infantil.pdf) of the Spanish Association of Paediatrics and the Spanish Nutrition Foundation. In addition, no gender distinction is made in practice in school cafeterias, and the same meals are given to boys and girls. At the same time, our formulation is intended for groups of people, rather than a particular person. Consequently, gathering physical data, nutritional goals and user preferences, in an effort to design a personalized plan, is out of the scope of this work.
Some other aspects that will be considered are those related to the cost, variety, and nutritional quality of the meals. In this paper, the cost of the menu and the variety of the meals are taken as the objective functions of the problem. The goal is to minimize the cost of the menu while minimizing the degree of repetition of the courses and food groups that the menu consists of. The nutritional features of the meal plans are handled as the problem constraints, i.e., any feasible plan must satisfy a set of nutritional requirements.
Bearing the above in mind, the main contributions of this work are as follows: • A multi-objective formulation of the MPP, termed MMPP, is described. It not only considers the cost of the menu but also the degree of repetition of courses and food groups of the plan, as the two objective functions to be optimized. • A rigorous experimental evaluation is conducted in which well-known Multi-Objective Evolutionary Algorithms (MOEA) are applied to solve different instances of the MMPP. • Quantitative and qualitative analyses of the resulting meal plans are provided that show their suitability from the standpoint of their nutritional value. • The multi-objective nature of the MMPP formulation is demonstrated.
A review of the literature on the MPP is presented in Section 2. Our formulation of the MMPP is introduced in Section 3. Then, the experimental evaluation is described in Section 4, which is also devoted to presenting and exhaustively analyzing the computational results obtained. Finally, the conclusions and some lines of further research are given in Section 5.

State of the Art in Menu Planning
The MPP [10] consists of producing healthy meal plans. Those plans are usually generated by considering a day, a week or a month. Moreover, each meal involved in the menu plan has to be stated precisely. Considering only lunch, for instance, a first course, second course and a dessert could be specified. Some interesting considerations in the problem field involve the design of approaches that are not only limited to offering nutritionally adequate diet plans, but that are also able to adapt these plans to the specific needs of the user: personal preferences involving food types, cost of the menus, meal variety, aesthetic characteristics of the food, preparation methods, incompatibilities with certain foods due to allergies, intolerances, diseases or specific lifestyles, and the quality of the food based on its seasonality, among others [10]. In this sense, the problem could involve many other objectives or constraints beyond offering a healthy menu, and consequently, it is necessary to define which objectives are going to be considered or optimized, and which set of constraints is to be satisfied by the feasible solutions.
Different aspects have been considered as problem objectives, while other have been considered as constraints, something that varies depending on authors. Defining one objective function as the price or cost of the meal plan, however, is one of the most frequently approaches [11][12][13]. In other works, in opposition to considering the cost of the meal plan, objective functions are defined by considering the level of acceptance [14], food preferences [15,16], or the level of suitability [17].
In general, in the MPP, it is very important to foster a healthy and balanced meal plan from the nutritional perspective. Since recommended minimum and maximum quantities of different macronutrients and micronutrients have to be provided, the definition of constraints related to nutritional requirements is something frequently addressed [11,12,18]. Additional aspects, such as the consistency or predominant color of the meals, their variety, or cooking time, among others, are also considered as constraints [12,16,19]. Finally, depending on the level of completeness that is desired, there may be many other aspects to consider. First and foremost is the users' data or personal preferences. If a meal plan has to be designed for a specific person, additional features, such as their weight, height, physical activity level, or health goals, could be considered in order to produce a more specific plan [15,17,18]. The above information is not required when menus are planned for groups of people, which is the case in the current work. Another important aspect not considered in many of the papers reviewed in the literature is whether the user in question suffers from any type of allergy, intolerance, or if their lifestyle simply prevents the person from consuming certain foods [15,17]. Whether or not the user belongs to one of these groups should be considered as a specific constraint so as to exclude or include the corresponding food in the meal plan.
Authors in the field have proposed different single-and multi-objective formulations of the MPP and have applied various types of meta-heuristics and optimization methods to obtain a solution [10]. Following an exhaustive review of the literature, Table 1 shows a summary of the most relevant works, including information on the number of objective functions addressed and the types of optimization methods applied. As we can see, most existing works in the literature use single-objective optimization methods to address the problem. Among these methods, Evolutionary Algorithms (EA) and, more specifically, Genetic Algorithms (GA) are very popular (see Table 1). In some works, several objectives are considered, although, in most cases, multiple objectives are reduced to a single-objective function during the optimization phase. In spite of the above, we have found three references where multi-objective formulations of the MPP are addressed by multi-objective optimizers. In all of them, MOEA are applied.

Meta-Heuristics
Other Methods

Reference
In Seljak [37], the Non-dominated Sorting Genetic Algorithm II (NSGA-II) [39] was applied to solve a multi-objective formulation of the MPP. In this particular case, plans were designed for a whole week, i.e., considering seven days, as well as five meals per day. The price of the meal plan, as well as its seasonal quality, flavor, and temperature, among other factors, were considered to define the objective functions. Constraints were defined by considering personal preferences, as well as nutritional needs. Finally, meal plans were designed bearing in mind adults. Kaldrim and Köse [36] also addressed a multi-objective variant of the MPP through the application of the NSGA-II. This time, objective functions considered the cost of the meal plan, as well as food preferences. Other aspects to provide a personalized meal plan were the gender and age of the user. A more recent paper that also solves the problem through meta-heuristics and considers a multi-objective formulation is the one written by Moreira et al. [38]. Once more, the NSGA-II was used to minimize the cost of the menu while also minimizing the nutritional error in Brazilian schools, i.e., for children between four and five years old, based on the daily nutritional requirements set by the government. Anyway, none of these proposals contemplates the same set of objectives here proposed. Moreover, these previous MPP formulations have also some differences regarding the constraints and the planning features considered-for all meals a day instead of a single main meal.
Regarding previous research carried out by the authors, a single-objective memetic approach was applied to a single-objective formulation of the MPP where global and daily constraints related to the nutritional requirements of the meal plans generated were taken into account [40]. The memetic algorithm was used to minimize the cost of the meal plans. The same set of global and daily constraints were also considered in a multiobjective formulation of the MPP recently proposed [41], where the single-objective memetic algorithm introduced in Reference [40] was compared, by performing a wide experimental assessment, to a novel multi-objective memetic approach based on the well-known Multiobjective Evolutionary Algorithm by Decomposition (MOEA/D) [42]. With respect to the objective functions optimized, the cost and the variety of specific courses and food groups that the meal plans consist of were considered. Since promising results were obtained by the multi-objective memetic approach based on decomposition in comparison to those results attained by the single-objective memetic algorithm, in this work, we carry out an exhaustive experimental comparison among multi-objective optimizers applied to a multi-objective formulation of the MPP with global nutritional constraints, in which the aim is to demonstrate the multi-objective nature of the said formulation.

Multi-Objective Formulation of the Menu Planning Problem
In this paper, we propose an alternative multi-objective formulation of the MPP, which we term the MMPP [43]. The main goal is providing, in an automatic manner, healthy and balanced meal plans for large groups of children at school canteens. As a result, meal plans only consider lunch, thus excluding the remaining daily meals. With respect to the objective functions being optimized, the first function consists of the meal plan cost, which should be minimized, since in the context of school canteens, meal plans should be affordable. The second objective function consists of the degree of repetition of courses and food groups, which should be minimized, as well. Diversity in a diet is required by children, and not only from the health standpoint [7,44] but also to not abhor certain foods. Past research has demonstrated that the more diverse a diet, the healthier [45][46][47]. A diverse diet involves a larger intake of micronutrients and macronutrients and, at the same time, is much more suitable in terms of its nutritional value and quality [48]. The second objective function was designed by taking into consideration the above motivation, as well as nutrition expert advice from the Intervention Programme for the Prevention of Childhood Obesity (information on this program is available at http://www.programapipo.com/) endorsed by the Health Service of the Government of the Canary Islands. At this point, we should note that, in past research, the diversity of courses and food groups has been considered as constraints [12,31]. Considering it as an objective function is a novel approach addressed by the authors [41].
The formulation we are introducing herein considers, as constraints, the nutritional quality of the meal plan in terms of the intakes of energy, macronutrients (fats, carbohydrates, and proteins), and micronutrients, by following the suggestions given by the White Book on Child Nutrition. Other documents, recommended by the Spanish Government Ministry of Education, Culture and Sports and the Ministry of Health, Social Services and Equality, were also taken into consideration (these consensus documents are available at http://www.aecosan.msssi.gob.es/AECOSAN/docs/documentos/nutricion/ educanaos/documento_consenso.pdf and https://www.aepnaa.org/).
Finally, it is important to remark that a database of ingredients and courses was designed in order to provide proper meal plans. The database not only contains information on their nutritional value but also on their allergens. The aforementioned documents were also considered to this end.

Solution Encoding
A solution, which in the Evolutionary Computation field is also known as an individual, consists of a set of meals-represented by a starter, a main course, and a dessert-for a given number of days. From the computational perspective, a solution can be represented by means of a vector of integer values I, in which the length is given by |I| = 3 · n, being n the number of days of the meal plan. Each integer value in the vector, i.e., i q=1,...,3·n ∈ I, matches with the identification number (id) of a particular course in the database, as Figure 1 shows. Moreover, the database also stores additional data for each course: amount in grams, nutritional value, food incompatibilities, allergens, food groups, name, and price. The size of the solution space S is given by Equation (1), where l st , l mc , and l ds represents the number of starters, main courses, and desserts stored in the database, respectively. As it can be observed, the size of the solution space exponentially increases with the number of days considered to design the meal plan. For example, by taking into consideration a meal plan for n = 10 days and l st = l mc = l ds = 10, there exist |S| = 10 30 potential solutions. |S| = (l st · l mc · l ds ) n . (1)
Formally, a feasible solution I should comply with the following constraints: where the global intake of nutrient h in a meal plan I is given by in(I, h), r h denotes the daily intake at lunch time of element h, and the minimum and maximum amounts of element h that the meal plan can include considering n days is given by r min h and r max h , respectively.

Objective Functions
The objective function modeling the degree of repetition of courses and food groups is calculated as: where v MC j represents the compatibility, in terms of food groups, among courses st, mc, and ds for day j; p st , p mc , and p ds are the penalty constants, one per course type; d st j , d mc j , and d ds j are the number of days since the starter, the main course, or the dessert last appeared before day j; and v FG j is the penalty value for repeating food groups in the previous five days before day j. The food groups considered for the meals available in this work are G = {other, meat, cereal, fruit, dairy, legume, shellfish, pasta, fish, vegetable}, in keeping with the suggestions given by the Intervention Programme for the Prevention of Childhood Obesity.
Penalties are the mechanism used by the decision maker to set preferences in terms of the particular way a meal plan is automatically generated. As shown in Table 2, penalties take constant float values in the range [0, 10], and the higher the values the greater the penalties [43]. As it can be observed, penalties are determined by the repetition of food groups (p 1 -p 10 ), the repetition of the same food group from one to five days prior to the current day (p 11 -p 15 ), and the repetition of specific courses (p 16 = p st , p 17 = p mc , and p 18 = p ds ). A preliminary analysis was carried out in order to study different solutions in terms of the repetition of food groups and specific courses. The above allowed proper values to be determined for penalties [43].
Considering those penalties for repeating food groups (p 1 -p 10 ), if the penalty value of a given food group is too large in comparison to the penalty values of the remaining food groups, then a plan with a lower number of courses belonging to that food group will be provided. For example, preference has been given to those courses consisting mainly of vegetables (p 10 = 0.1) over courses consisting mainly of meat (p 2 = 3). Penalties for repeating food groups are considered during the calculation of v FG and v MC (see Equation (3)).
The higher the number of days without repeating a specific food group, the lower the value assigned to those penalties for repeating food groups in previous days (p 11 -p 15 ). Since school meal plans are usually designed for five days (a school week goes from Monday to Friday), a time window T was set to five days . Repeating specific courses is more restrictive in comparison to the repetition of food groups. A meal plan should not include the same course several times, at least in a short period of time. Nevertheless, food groups will likely be repeated because a specific food group usually is included in several courses. Penalties for repeating food groups are used in the computation of v FG (see Equation (3)). Finally, the values of those penalties for repeating a specific course (p 16 = p st , p 17 = p mc , and p 18 = p ds ) are assigned depending on the number of starters, main courses, and desserts that the database includes. In our case, the number of desserts included in the database is lower in comparison to starters and main courses and, consequently, desserts will be repeated more frequently. Bearing the above in mind, the penalty value associated to the repetition of desserts will be lower in comparison to the penalty values associated to the repetition of starters and main courses.
The value of v MC j is computed by means of Equation (4). The term x g ∈ {0, 1, 2, 3} refers to the number of repetitions of a specific food group in the three courses, i.e., starter, main course, and dessert, of a particular day j, while the penalty value for repeating the food group g is given by p g . As it can be observed, there exist |G| different food groups.
The value of v FG j is calculated through Equation (5). The reader should recall that a time window T = 5 was fixed and that there exist |G| different food groups. Considering day j, y g d = 1 means that the food group g was repeated on day j − d, while y g d = 0 means that it was not. At the same time, z d = 1 means that any food group was repeated d day(s) before day j, while z d = 0 means that no food group was repeated d day(s) before day j. Finally, penalty values are given by p g and p |G|+d .
The total cost of a plan, i.e., the other objective function considered in this formulation, is computed as follows: where c st j , c mc j , and c ds j are the costs for the starter, main course, and dessert, respectively, for day j. In solution I, st j , mc j , and ds j correspond to elements i 3·j−2 , i 3·j−1 , and i 3·j , respectively. We should note that the database stores the price per kilogram for each ingredient, as well as the amount of grams of each ingredient that a specific course requires to be prepared. Bearing the above in mind, the cost of a specific course is given by the sum of cost of the ingredients that the said course consists of.

Experiments and Evaluation of Results
To deal with the particular formulation of the MMPP proposed in this paper, we have considered the following algorithms: The adaptive version of the Indicator-based Evolutionary Algorithm (IBEA) [50].
It is true that only a few contributions have been made considering the application of multi-objective optimizers to multi-objective formulations of the MPP, as we stated in Section 2. Nevertheless, as far as we know, NSGA-II has been selected as one of the optimizers to be applied in all those previous contributions. Bearing the above in mind, it would be interesting to determine how it performs with our MMPP formulation. NSGA-II is a Pareto-based MOEA, meaning it applies Pareto optimality concepts to guide the search. The SPEA2, which is another Pareto-based MOEA, was also selected in order to compare its performance with respect to the NSGA-II. Finally, the IBEA was also selected so as to include in the comparison a MOEA that applies a quality indicator, specifically, the binary additive -indicator [51], to assign a fitness value to solutions.
In order to complete the definition of the MOEA selected for comparison, the crossover, mutation and repair operators must be described. The crossover method aims to uniformly combine a pair of solutions or solutions I and I to create a new offspring. Given a crossover rate p c , which is the probability of applying the operator, and q = {1, ..., |I|}, all the elements of both solutions will be traversed at the same time, and, for each element i q ∈ I and i q ∈ I , there will be a 50% probability that said elements will be exchanged. The mutation operator is used to randomly replace every gene of a particular solution, i.e., a daily meal consisting of a starter, a main course, and a dessert, with another daily meal with probability p m . The repair method is responsible for assessing each solution according to the nutritional value constraints and for determining if the solution is feasible, based on the definition of feasibility given in Section 3.2. In particular, the same procedure used in the mutation operator will be applied to change a daily meal of a non-feasible solution. The repair operation, however, will be applied always, i.e., no rate is taken into consideration. Each daily meal will be modified by following the aforementioned steps, until a feasible solution is provided.
Experimental method: The Meta-heuristic-based Extensible Tool for Cooperative Optimization (METCO), proposed by León et al. [52], was used to implement all the algorithms, as well as the MMPP. Further information about the source code of the MMPP, the specific instances defined for the current work, and our results can be found at the following Github repository. The experiments were executed on a Debian GNU/Linux computer with four AMD ® Opteron™ processors (model number 6348 HE) at 2.8 GHz with 64 GB of RAM. Unless otherwise specified, each experiment was run 30 times, since all the approaches considered are stochastic. With the aim of statistically comparing the results achieved, the following procedure [53] was applied. In order to check for the normality of the results, a Shapiro-Wilk test was applied. Then, if results are normal, a Levene test is applied to check for homogeneous variances. If results present homogeneous variances, then the ANOVA test is considered. In the opposite case, the Welch test is applied. In the case that results do not follow a Gaussian distribution, the Kruskal-Wallis test is considered. A significance level α = 0.05 is taken into consideration when applying all the above tests.
The hypervolume [54] was selected as the metric for comparing the different approaches. The hypervolume indicator I hyp (A) computes the volume of the region H delimited by a given set of points A and a set of reference points N, as it is shown by Equation (7). Therefore, larger values of the indicator will correspond to better solutions.
Ideally, the set of reference points should be the nadir points, i.e., the worst points in the Pareto Front or, in other words, the elements of the Pareto Front that do not dominate any other point. In our case, the worst-maximum-values attained for the two objective functions of the MMPP, considering all the runs performed in each of the experiments, were used to set the reference points. Additionally, the hypervolume was normalized in the range [0, 1]. As we stated before, the higher its value, the better the performance of the algorithm in question.
Instances: In order to define an instance of the problem, the number of days for which the meal plan will be designed has to be specified. Furthermore, a starter, a main course, and a dessert have to be selected for each day. As a result, all potential starters, main courses, and desserts that can be chosen also have to be described somehow. In the particular case of our implementation, a total number of 67 courses consisting of l st = 19 starters, l mc = 34 main courses, and l ds = 14 desserts are provided in three different files containing information on each type of course. Assuming a meal plan for 5 days, for instance, the size of the search space would be |S| = (19 × 34 × 14) 5 = 6.051 × 10 19 potential solutions (the reader is referred to Equation (1)). The format designed for those files, and for each of their lines in particular, is depicted in Figure 2, where each field must be filled in with the following information:  (1), cereal (2), fruit (3), dairy (4), legume (5), shellfish (6), pasta (7), fish (8), vegetable (9). Only those food groups that the course in question contains have to be specified through their corresponding integer numbers.

First Experiment: Parameter Setting
In this first experiment, we conduct a study involving the values that parameters belonging to the MOEA selected for comparison should take. Hence, the goal in this first experiment is not focused on comparing the performance of those MOEA but analyzing their corresponding parameters separately, in order to properly set them for subsequent experiments.
To this end, we tested different values for the parameters of each algorithm, as Table 3 shows. Specifically, a total of 36 different configurations for each MOEA were applied in order to yield a 20-day meal plan. Those configurations were obtained by combining four values for the population size (n), three values for the mutation rate (p m ), and three values for the crossover rate (p c ). In the specific case of parameter p m , test values were selected based on the operation of the mutation operator. The reader should recall that a mutation operator specifically designed for the MMPP is considered herein. Specifically, for each day, the operator modifies the three courses with probability p m . Since this experiment requires a 20-day meal plan, p m = 0.05 is the minimum probability that allows, at least, one day's courses to be changed. Bearing the above in mind, the values 0.05, 0.1, and 0.15 were selected; therefore, the courses for one, two, and three days were altered by the mutation operator, respectively. So as to avoid the overly disruptive behavior of the operator, larger values for p m were not considered. In the case of the crossover rate, high values are usually set. That is the reason why the test values for parameter p c were 0.8, 0.9 and 1. In order to complement the information given in Table 3, we should note that all configurations based on the IBEA applied a scaling factor equal to 0.002, while, for each configuration of the SPEA2, the archive size was set to the same value considered for the population size. Finally, note that a sufficiently long stopping criterion was set, consisting of 4 × 10 8 function evaluations. The idea was to study whether the MOEA selected converged prematurely or not, and thus determine a proper value for the stopping criterion to be used in subsequent experiments. Furthermore, for this first analysis, and due to the large number of configurations tested, the executions were only repeated 10 times. With respect to the running time of the different approaches, we note that the mean time, considering 10 repetitions of the executions, required by NSGA-II to perform 4 × 10 8 function evaluations, was 31,438 seconds, while, in the case of IBEA and SPEA2, the mean time was equal to 129,601 and 50,783 seconds, respectively. It can be observed how NSGA-II is the fastest optimizer, while IBEA is the slowest. The above is probably due to the fact that IBEA needs to calculate a performance metric in order to guide the search, something that is not required by NSGA-II and SPEA2. Anyway, performing a running time analysis of the approaches selected for comparison is out of the scope of the current work. Table 4 shows some statistical data related to the hypervolume values attained by some configurations of the NSGA-II at the end of the executions. The reader should recall that 36 different configurations of each MOEA were run, but in order to simplify tables, information is only shown for five configurations. Tables including information for the whole set of configurations executed are provided as Supplementary Material. The last three columns can be used to rank the different configurations tested. In particular, the number of times a given configuration was able to statistically outperform other configurations (W), and the number of times a particular configuration was statistically outperformed by other configurations (L) are shown. Configuration A statistically outperforms configuration B if the p-value obtained is lower than the significance level α , and, if at the same time, A provides a higher mean and median for the hypervolume (when the stopping criteria are reached). The ranking (R) was calculated as the difference between the number of times a configuration outperformed others, and the number of times a configuration was outperformed by others (R = W − L). Note that the configurations are shown in descending order based on their position in the ranking as the first criterion, and the mean and the median of the hypervolume as the second and third criteria. With the aim of differentiating the configurations, the nomenclature Algorithm_ps_n_pm_p m _pc_p c is used. Finally, Tables 5 and 6 show the same information but for the IBEA and SPEA2 approaches, respectively. The same nomenclature is used to differentiate the configurations, but adding some slight modifications to include information on the scaling factor (sFactor) in the case of the IBEA, and the archive size (as) in the case of the SPEA2.   With respect to the NSGA-II, we observed that, in general, larger population sizes yielded better performance in terms of hypervolume values, in comparison to smaller populations. For instance, the values for the first-ranked configuration were p m = 0.1, p c = 0.8, with n = 250 solutions, which was the largest population tested. If the same values are set for parameters p m and p c , the worst-ranked configuration contained 25 solutions, which was the smallest population considered. In order to better analyze the above fact, Figure 3 shows the mean hypervolume values attained by the different configurations of the NSGA-II (first row), the IBEA (second row), and the SPEA2 (third row) at the end of the executions with respect to the population size (first column), the mutation rate (second column), and the crossover rate (third column). In the case of the NSGA-II, note that, for almost every possible combination of values for parameters p m and p c (first row, first column), the best results were achieved by applying the largest population considered (n = 250), while the worst results were provided by using the smallest population tested (n = 25).
Regarding the crossover rate, we should note that, generally speaking, lower values yielded better hypervolume values, and therefore better performance, while higher values of the crossover rate decreased the performance of the NSGA-II significantly. As proof of the above, the last five configurations in the ranking were executed with p c = 1, while, in the first five configurations, p c = 0.8.
Moreover, as Figure 3 shows (first row, third column), in the general case, the lowest hypervolume values were achieved with p c = 1, while the highest hypervolume values were attained with p c = 0.8.  In regard to the mutation rate, however, a clear conclusion cannot be extracted. It is true that, for instance, the first three configurations of the NSGA-II in the ranking had the same population size and crossover rate, with a different value for the mutation rate. Bearing the above in mind, we may claim that the mutation rate p m does not seem to have a significant effect on the performance of the NSGA-II, at least considering the values selected for this study. However, as Figure 3 shows (first row, second column), in some cases, lower values for the mutation rate provided better results, while higher values of the mutation rate yielded better performance over other cases. As a result, more experiments should be conducted to shed more light on the behavior of the NSGA-II with respect to the mutation rate, which seems to be its most sensitive parameter when dealing with this particular problem.
Finally, it is important to remark that the first configuration of the NSGA-II in the ranking, which was run with n = 250, p m = 0.1, and p c = 0.8, was able to statistically outperform 28 out of 35 configurations (80%), while it was not statistically outperformed by any other configuration. Furthermore, that configuration achieved the highest mean and median of the hypervolume at the end of the executions, taking into account the 36 different parameterizations of the NSGA-II tested.
Similar conclusions to those given for the NSGA-II can be extracted for the IBEA (Table 5 and Figure 3). In the general case, the application of larger population sizes and lower crossover rates yielded better performance, while the mutation rate seems to be its most sensitive parameter, as in the case of the NSGA-II. In fact, the configuration of the IBEA that yielded the best performance, i.e., the first configuration in the ranking, was executed with n = 250, p m = 0.1, and p c = 0.8, which was also the best-performing parameterization of the NSGA-II. The best-performing configuration of the IBEA was able to statistically outperform 32 out of 35 configurations (91.4%), while it was not statistically outperformed by any other configuration. Moreover, we should note that the best-performing configuration of the IBEA achieved the highest median and mean of the hypervolume, considering the 36 different configurations tested, as was the case with the NSGA-II.
The behavior of the SPEA2, however, was slightly different in terms of its parameters (Table 6 and Figure 3), and particularly of the population size, with respect to the other two approaches. The best-performing configuration of the SPEA2 had n = 100 solutions, in contrast to the best-performing configurations of the NSGA-II and the IBEA, which were run with 250 solutions. The above may be explained by the fact that the SPEA2 makes use of an archive in which the size is usually set equal to the value selected for the population size, which allows the non-dominated solutions considered so far to be stored during the executions. The use of that archive may avoid the need to employ larger population sizes. With regard to the remaining parameters (p m and p c ), similar conclusions to those drawn previously for the NSGA-II and IBEA can be also extracted for the SPEA2. Finally, we note that the first configuration of the SPEA2 in the ranking statistically outperformed 34 out of 35 configurations (97.1%), while it was not outperformed by any other configuration. As was the case with the NSGA-II and the IBEA, that first-ranked configuration yielded the highest mean and also median for the hypervolume indicator, considering the 36 different parameterizations of the SPEA2.
The previous analysis was carried out by considering the results obtained at the end of the executions. However, it would be interesting to study the behavior of the different approaches during the runs. As we said at the beginning of this section, a long stopping criterion was set with the aim of analyzing if any of the MOEA selected for comparison converged prematurely. Figure 4 shows the evolution of the mean hypervolume achieved by some configurations of the NSGA-II, the IBEA and the SPEA2. Specifically, for each of the approaches tested, the first and last configurations, i.e., the best-and worst-performing parameterizations, respectively, given their corresponding rankings, were selected for comparison purposes. The configurations ranked 9, 18, and 27 were also included in the comparison. In fact, the above five configurations are those included in Tables 4-6. As we can see, for the three approaches, significant differences arose in terms of the performance yielded during the runs by the parameterizations applied, something that we had already stated previously. In addition, despite having set a long stopping criterion consisting of 4 × 10 8 function evaluations, the mean hypervolume was improved even at the end of the runs, allowing us to conclude that none of the MOEA converged prematurely. Bearing the above in mind, runs may be executed for even longer in order to study the behavior of the approaches in the long term.

Second Experiment: Designing Meal Plans for a Variable Number of Days
The main goal of the second experiment is to analyze the performance of the different MOEA when they are applied to produce meal plans for a variable number of days. Considering the conclusions extracted from the first experiment, we decided to apply the NSGA-II, IBEA, and SPEA2 with the same parameter values used by their corresponding first-ranked configurations. Hence, the size of the population was set to n = 250 solutions for NSGA-II and IBEA, and to n = 100 solutions for SPEA2. The crossover rate was set to p c = 0.8 for all of the approaches. In addition, the scaling factor of the IBEA was set to 0.002, and the archive size of the SPEA2 was set to 100 solutions. Finally, since the first experiment provided no clear conclusions on the mutation rate p m , three different values were tested in the second experiment for this parameter. Since for this experiment we decided to design meal plans for 5, 10, 20, and 40 days, p m was set such that the meals for one, two or three days are changed by the mutation operator. For instance, values of 0.2, 0.4, and 0.6 were selected when dealing with meal plans for five days, and therefore, the meals for one, two and three days were altered by the mutation operator, respectively. The reader should recall that, in this case, the chromosome would consist of five genes, i.e., the meal plan length, with each gene encoding a starter, a main course, and a dessert. When designing meal plans for 40 days, however, the values 0.025, 0.05, and 0.075 were tested.
Bearing the above in mind, three different configurations for each MOEA were applied to design meal plans for 5, 10, 20, and 40 days. The stopping criterion was also set to 4 × 10 8 function evaluations, as in the first experiment, but this time every run was repeated 30 times rather than 10 times, in an effort to provide much more statistically significant conclusions. Table 7 shows, for each meal plan length, statistical information related to the hypervolume values provided by each configuration of the different MOEA at the end of the runs. The last three columns show data on the statistical ranking achieved by each configuration, which was calculated by following the same procedure explained in the first experiment. Note how, for a small number of days, i.e., 5 and 10 days, the MOEA that yielded the best performance at the end of the executions was the NSGA-II. When designing a meal plan for 5 days, the three highest-ranking configurations were the three parameterizations of the NSGA-II. In fact, each of those three configurations statistically outperformed every IBEA and SPEA2 configuration. When designing a meal plan for 10 days, the three NSGA-II configurations also obtained the first three positions in the ranking. Although each of its configurations statistically outperformed every SPEA2 configuration, some of its configurations did not exhibit statistically significant differences in comparison to some IBEA parameterizations.
For longer meal plans, however, i.e., 20 and 40 days, the best performance at the end of the runs was attained by the SPEA2. When producing meal plans for 20 days, two configurations of the SPEA2 where the highest ranked. Those two configurations were able to statistically outperform, or did not exhibit statistically significant differences with respect to, all of the NSGA-II and IBEA configurations. Finally, when dealing with 40-day meal plans, the three SPEA2 configurations obtained the top three positions in the ranking, with each of them statistically outperforming each NSGA-II and IBEA configuration.
Regarding the mutation rate, note how small values yielded better performance in comparison to larger ones. In fact, when designing meal plans for 5 and 10 days, where the NSGA-II was the best-performing MOEA, those configurations that used the smallest mutation rate tested obtained the first position in their corresponding rankings. When dealing with longer meal plans, i.e., for 20 and 40 days, the SPEA2 behaved similarly, with lower mutation rates providing better performance in comparison to higher values of this parameter.
As was done in the first experiment, it is interesting to study the behavior of the MOEA both during and at the end of the executions. Figure 5 shows, for each meal plan length considered in this experiment, the evolution of the mean hypervolume attained by the NSGA-II. Figures 6 and 7 show the same information for the IBEA and the SPEA2, respectively. Note how the three MOEA exhibited differences during the runs depending on the specific value selected for the mutation rate. Although a few exceptions arose for some instances, particularly in the case of the IBEA, in general, every approach that used lower values for the mutation rate was able to attain a higher mean for the hypervolume, and therefore better performance, during almost every execution, in comparison to using larger mutation rates. Table 7. Performance of the different Multi-Objective Evolutionary Algorithms (MOEA) in terms of the hypervolume obtained at the end of 30 runs. For each menu plan duration, configurations are sorted in descending order by considering their ranking as the first criterion, and the mean and the median of the hypervolume as the second and third criteria, respectively.

Configuration
Min

Third Experiment: Qualitative and Quantitative Analysis of Solutions
In this third experiment, the main goal is to analyze the solutions attained by the MOEA considered in this paper from both a quantitative and qualitative point of view. Additionally, we would like to validate the multi-objective formulation of the MPP we are providing herein in terms of the features of the solutions achieved. In order to do this, we considered, for each meal plan length, the solutions obtained in the second experiment by the best-and worst-performing configurations of each MOEA, i.e., the configurations obtaining the first and last positions in their corresponding ranking. Figure 8 shows, for each meal plan length, examples of front approximations obtained after the tests by those best-and worst-performing parameterizations. First of all, we should note that the shape of the different front approximations obtained demonstrates that the two objective functions are in conflict with each other. A decrease in the cost of a meal plan entails an increase in the frequency of course repetitions, and vice-versa. Bearing the above in mind, we can state that the multi-objective formulation of the MPP we are proposing herein is valid from the standpoint of its multi-objective nature.
Note also that differences among the front approximations attained by the bestperforming and worst-performing configurations arose for each plan length. Generally speaking, for each case, the best-performing configuration was able to provide a better front approximation not only in terms of its convergence but also in terms of its spread and uniformity, features that are closely related to each other and that are usually referred to as the diversity of a front. For instance, when dealing with meal plans for 10 and 20 days, it holds that the corresponding best-performing configuration is able to provide a significant number of solutions that dominate a large number of solutions attained by the corresponding worst-performing configuration. The improvement in performance provided by the best configuration with respect to the worst one was even more noticeable when designing 40-day meal plans. The best-performing configuration was able to provide solutions located in regions of the front that are not covered by those solutions attained by the worst-performing configuration. It can be observed that in the case of smaller instances, the best-performing approach was NSGA-II, while, in the case of larger instances, the bestperforming scheme was SPEA2. It is important to recall that one of the main differences of SPEA2 in comparison to NSGA-II and IBEA is that the former makes use of an archive, which allows the non-dominated solutions achieved so far to be stored during the executions. The use of an archive could have a significant impact over the performance of an optimizer when dealing with larger instances. Nevertheless, more experiments are required to clarify the above. Finally, it would also be interesting to analyze one of those solutions in depth, not only from the point of view of the different meals that it offers but also from the point of view of its nutritional value. In order to do that, we selected one of the 10-day meal plans obtained by the corresponding best-performing configuration for which the values for the cost and the degree of repetition were 8.16 and 16.145, respectively. Table 8 enumerates the different meals in said meal plan, while Table 9 describes the total amount of nutrients it provides. According to the White Book on Child Nutrition, the main nutrients to consider in the case of lunch are carbohydrates, fats and proteins. The total number of kilocalories recommended for a child is 2000 kcal per day. Lunch accounts for approximately 35% of the total daily nutrient intake. As a result, the recommended number of kilocalories in a lunch should be approximately 700 kcal. The meal plan depicted in Table 8 contains a total of 7269.34 kcal, which provides the recommended kilocalorie intake for 10.4 days. Regarding the total recommended amount of fats, they should contribute no more than 35% of the total amount of energy, meaning that, for a lunch of 700 kcal, 245 kcal should come from fats. One gram of fat provides 9 kcal, so one lunch should contain 27.2 g of fat. The 10-day meal plan in our example provides 282.05 g of fat, which corresponds to the total amount of fats for 10.4 days. The recommended amount of carbohydrates is 50% of the total energy intake provided by one lunch; consequently, 350 kcal should come from carbohydrates.
One gram of carbohydrate provides 4 kcal, and, as a result, a lunch should contain 87.5 g of carbohydrates. Our meal plan provides 769.24 g of carbohydrates, which corresponds to the total amount for 8.8 days. Finally, the recommended amount of proteins is 15% of the total energy provided by one lunch, which means that 105 kcal should come from proteins. One gram of protein provides 4 kcal; therefore, a lunch should contain 26.25 g of protein.
Our meal plan example provides 291.08 g of proteins, which is equivalent to 11.1 days in terms of the intake recommendation. In conclusion, although the recommended amounts of nutrients may vary slightly from one recommendation to the next, and in practical terms there is a margin of excess and deficiency of nutrients above and below the recommended amounts, the meal plan selected for this study provides, in general terms, an adequate and balanced amount of nutrients. Chicken salad with mayonnaise Potato omelette Tangerine   10 Mashed vegetables Marinera potatoes Watermelon Table 9. Example of the total amount of nutrients for a 10-day menu plan. For each nutrient h, its corresponding amount r h must be between the minimum (r min h ) and maximum (r max h ) amounts allowed. All amounts are for 10 days.

Conclusions and Further Research
Most approaches found in the literature on menu planning deal with single-objective formulations of the problem. However, in this paper, we propose a multi-objective formulation, termed MMPP, which considers two conflicting goals: minimizing the menu cost and minimizing the degree of repetition of the courses consumed. Our meal plans are designed for school cafeterias; as a result, the menu prices should not exceed reasonable amounts. At the same time, the courses should be varied to keep the students interested in the food. The computational results attained by a set of well-known MOEA show that suitable meal plans-in terms of their nutritional values-are provided. For short plans-5 and 10 days-the algorithm achieving the highest hypervolume values at the end of the executions is NSGA-II. For longer plans-20 and 40 days-the best-performing algorithm is SPEA2. In addition, the experimental analysis also demonstrates the multi-objective nature of the introduced formulation. The cheapest meal plans are those with the lowest variety of courses and food groups, while the largest variety of courses and food groups is given by the most expensive plans.
This work presents an initial approach to the MMPP. During and after its development, we identified some aspects worthy of further research, and other features to consider for inclusion. A larger number of courses could be added to the database in order to provide a greater variety in the production of meal plans. This would expand the search space significantly, however. Seasonal products could also be considered as a determining factor in the price and quality of the food. Regarding the problem formulation, it would be interesting to carry out more tests by including additional tailored approaches, such as local searches and variation operators specifically designed for the said formulation. Furthermore, other MOEA could be applied, including aggregation-based multi-objective optimizers, such as Global WASF-GA. Additionally, the repair method could be improved in such a way that it can obtain a feasible solution in a much smarter manner.
Since we are designing menus for a group of people, the greatest limitation of this approach is that it precisely does not allow to independently optimize menus for certain people in the group who involve specific restrictions in relation to allergies, intolerances, or to lifestyle. In these cases, we could provide a tool that helps nutritionists to automatically evaluate pre-designed meal plans or alternative plans to those initially obtained by the optimizer. This tool could automatically calculate costs, the degree of repetition of specific courses, and food groups and nutritional values, among others. Since different algorithms yielded the best results depending on the particular instance of the problem being solved, it would be interesting to apply machine learning techniques in order to classify instances in terms of their features. Hence, when solving a novel instance, the best-performing approach for that particular instance could be recommended in light of its characteristics. Finally, we note that the proposed formulation may be adapted, with little effort, to deal with the specific nutritional requirements of other groups of people, including adults, such as patients in hospitals and convicts in prisons, among others. Other interesting analyses could be performed by considering personalized nutritional plans that model, for instance, intermittent fasting, one meal a day, and ketogenic diet, among others.
Supplementary Materials: The following are available at https://www.mdpi.com/2227-7390/9/ 1/80/s1, Table S1: Performance of the NSGA-II in terms of the hypervolume obtained at the end of the runs.
Funding: This work was partially funded by the Spanish Ministry of Economy, Industry and Competitiveness as part of the "I+D+i Orientada a los Retos de la Sociedad" program (contract number TIN2016-78410-R). It was also partially funded by the Spanish Ministry of Science, Innovation and Universities, as well as by the University of La Laguna, as part of the program "Nuevos Proyectos de Investigación: Iniciación a la Actividad Investigadora" (contract number 1203_2020).
Institutional Review Board Statement: Not applicable.

Informed Consent Statement: Not applicable.
Data Availability Statement: Publicly available datasets were analyzed in this study. This data can be found here: https://github.com/Tomas-Morph/On-the-Design-of-Healthy-Menus-through-Evolutionary-Computation.