The application of nonlinear programming on ration formulation for dairy cattle

The objective of this study was to compare the application of iterative linear programming (iteLP), sequential quadratic programming (SQP), and mixedinteger nonlinear programming-based deterministic global optimization (MINLP_DGO) on ration formulation for dairy cattle based on Nutrient Requirements of Dairy Cattle (NRC, 2001). Least-cost diets were formulated for lactating cows, dry cows, and heifers. Nutrient requirements including energy, protein, and minerals, along with other limitations on dry matter intake, neutral detergent fiber, and fat were considered as constraints. Five hundred simulations were conducted, with each simulation randomly selecting 3 roughages and 5 concentrates from the feed table in NRC (2001) as the feed resource for each of 3 animal groups. Among the 500 simulations for lactating cows, 57, 45, and 21 simulations did not yield a feasible solution when using iteLP, SQP, and MINLP_DGO, respectively. All the simulations for dry cows and heifers were feasible when using SQP and MINLP_DGO, but 49 and 11 infeasible simulations occurred when using iteLP for dry cows and heifers, respectively. The average ration costs per animal per day of the feasible solutions obtained by iteLP, SQP, and MINLP_DGO were $4.78 (±0.71), $4.45 (±0.65), and $4.44 (±0.65) for lactating cows; $2.39 (±0.52), $1.48 (±0.26), and $1.48 (±0.26) for dry cows; and $0.98 (±0.72), $0.97 (±0.15), and $0.91 (±0.14) for heifers, respectively. The average computation time of iteLP, SQP, and MINLP_DGO were 0.59 (±1.87) s, 1.15 (±0.62) s, and 58.69 (±68.45) s for lactating cows; 0.041 (±0.070) s, 0.76 (±0.37) s, and 14.84 (±39.09) s for dry cows; and 1.60 (±2.90) s, 0.51 (±0.19) s, and 16.45 (±45.56) s for heifers, respectively. In conclusion, iteLP had limited capability of formulating least-cost diets when nonlinearity existed in the constraints. Both SQP and MINLP_DGO handled the nonlinear constraints well, with SQP being faster, whereas MINLP_DGO was able to return a feasible solution under some situations where SQP could not.


INTRODUCTION
Feed costs account for around 50 to 70% of the expenses of operating a dairy farm (Bozic et al., 2012), which highlights the importance of minimizing this cost when feasible. Diet formulation relies on the nutrient requirements of the animal, nutrient composition of available feeds, and nutrient interactions, all of which require systems for estimation, such as Nutrient Requirements of Dairy Cattle (NRC, 2001) and the Cornell Net Carbohydrates and Protein System (Fox et al., 2004). One method to formulate least-cost diets that fulfill an animal's nutrient requirements is linear programming (LP), which optimizes a linear objective function subject to a set of linear constraints (Chandler and Walker, 1972;O'Connor et al., 1989). A linear objective function (or a linear constraint) is a function with a polynomial of degree equal to 1 or 0 (e.g., a + b 1 x 1 + b 2 x 2 + . . . + b n x n ) so that the value to be optimized has a linear relationship with the variables. Several studies have implemented LP for ration formulation to explore the effects of feed composition variability in the constraints (Tozer, 2000), feed price variability in the objective function (Alqaisi and Schlecht, 2021), and objective functions other than least cost (Qu, et al., 2019). One limitation of LP is that it allows only linear objective functions and constraints, while some equations in the dairy nutrition model are nonlinear when adapted to an LP structure. Nonlinear functions include a broad range of functions, in which the change of the output is not proportional to the change of the input [e.g., x 2 , exp(x), log(x), and piecewise function].
In NRC (2001), nonlinearity mainly exists in the cal-culation of energy and protein contents of feeds. For energy, the TDN values of feeds are corrected by the animal intake level, calculated as follows: where Kp i = protein passage rate of feed i, %/h; PercentConc = dietary concentrate percentage = the amount of concentrate (kg)/total amount of feed (kg) × 100%, percent of DM. These rates are then used to calculate rumen degradable and undegradable protein (RDP and RUP) contents. The intake level and dietary concentrate percentage are unknown before formulating the diet, so adapting these equations to an optimization programming creates the nonlinearity. Not many studies have investigated how to handle nonlinearity in ration optimization, which is important because increasingly more nonlinear equations may appear with advances in dairy nutrition. Moraes et al. (2012) used an iterative linear programming (iteLP) method to deal with the nonlinearity in the NRC (2001) for ration formulation. However, iteLP has certain limitations, which will be discussed in this paper. A nonlinear programming optimization strategy, sequential quadratic programming (SQP; Boggs and Tolle, 1995), was employed in the Cornell-Penn-Miner (CPM) dairy model to formulate rations based on the Cornell Net Carbohydrates and Protein System (Boston et al., 2000), but the performance of using SQP based on NRC (2001) is unknown. Mixed-integer, nonlinear programming-based, deterministic global optimization (MINLP_DGO) is another method to solve nonlinear programming problems and is capable of solving a broader range of problems involving continuous, binary, and integer variables compared with SQP.
Mathematical modeling has been an important technique to evaluate production and environmental impacts of dairy systems (France and Kebreab, 2008). Whole-farm models such as the Integrated Farm Systems Model can provide holistic estimates of production and environmental outcomes in response to changes in weather and management (Rotz et al., 2016;Veltman et al., 2018) and require a method to predict feed use and delivery over the timespan of the model simulation. The Integrated Farm Systems Model currently uses an iteLP approach for ration formulation with requirements modified from a previous version of the NRC dairy nutrient requirement model (NRC, 1989), the NRC beef nutrient requirement model (NRC, 2000), and early versions of the Cornell Net Carbohydrates and Protein System model (Rotz et al., 1999(Rotz et al., , 2016. Limitations of the ration formulation method in the Integrated Farm Systems Model are among the factors that have led our group to develop a new whole-farm model, the Ruminant Farm Systems Model, to evaluate connections between dairy system components, including animal husbandry and feeding, manure management, field and crop management, and feed storage (Kebreab et al., 2019). The nonlinear programming framework for ration formulation developed in this study will automate simulation of feed use and production within the Ruminant Farm Systems Model and represent an advancement over extant whole-farm models (Kebreab et al., 2019). The objectives of this study are (1) to introduce a methodology to formulate least-cost rations using MINLP_DGO and (2) to compare the performance of iteLP, SQP, and MINLP_DGO on ration formulation based on NRC (2001).

METHODS
Least-cost diets were formulated to meet the nutrient requirements of dairy cattle in different life stages, with nutrient compositions of feeds as described by the NRC (2001). Ration formulation constraints were determined based on NRC (2001) and included the requirements for net energy for maintenance (NE M , Mcal/d), lactation (NE L , Mcal/d), and growth (NE G , Mcal/d), and metabolizable protein (MP, g/d), calcium (g/d), and phosphorus (g/d). We introduced additional ration formulation constraints to guide formulation of diets that meet requirements for rumen function. These included a fat constraint of less than 7% of diet DM (NRC, 2001), an NDF constraint of greater than 25% and less than 40% of diet DM (NRC, 2001;Moraes et al., 2012), and a constraint of forage NDF greater than 19% of diet DM (NRC, 2001). Finally, DMI was limited to be less than or equal to the predicted DMI Li et al.: NONLINEAR PROGRAMMING FOR RATION FORMULATION in NRC (2001), so that low-quality and low-cost feeds would not be overfed. We provide detailed information of all the constraints in Supplemental File S1 (Li, 2021). Three optimization strategies, iteLP, SQP, and MINLP_DGO, use different approaches to deal with the nonlinearity existing in the energy and protein constraints as will be described.

Optimization Strategies
iteLP. Formulating least-cost diets using LP can be written as follows: Min Z c x subject to a x b for i 1, 2, . . . , n, where Z ($) is the diet cost; c j ($/kg of DM) is the feed price of feed j; x j (kg of DM) is the amount of feed j; a ij is the coefficient for x j in constraint i; and b i is the lower bound of constraint i. For example, if constraint i represents NE M requirement, then b i is the minimum NE M requirement (Mcal) and a ij is the NRC (2001) predicted concentration of NE M in feed j (Mcal/kg of DM). The iteLP method uses the following procedures to deal with the nonlinearity existing in Eq.
[1] and [2]. The feed ingredients (a ij ) in iteLP starts at certain initial values and then are updated according to the DMI at the solution. The iteration process is repeated until DMI and dietary concentrate percentage at the solution are relatively constant, to ensure that energy and protein requirements are met (Moraes et al., 2012). However, we found that iterating based on DMI and dietary concentrate percentage does not always result in a satisfactory solution, depending on the feed ingredients. During the iteration process, 2 solutions may have similar DMI and concentrate percentage but very different values of dietary TDN concentration and RDP intake, leading to a large discrepancy between nutrient requirement and nutrient supply at the final solution. Therefore, the iteration was based on intake level and microbial crude protein (MCP) production in this study. The algorithm stopped when the differences of intake level and MCP production between 2 iterations were both lower than 0.1%. To prevent infinite iterations when intake level or MCP does not converge, we set the maximum number of iterations to be 1,000. SQP. Sequential quadratic programming allows nonlinear constraints, so the nutrient requirement constraints can be built according to NRC (2001) directly. The basic structure of an SQP problem is as follows: Min Z c x subject to g 0, for i 1, 2, . . . , n, where x is a vector of feed inclusion rates (x j defined as above) and g i (x) is constraint i, which must be twice continuously differentiable with respect to all x j in x. However, several equations in NRC (2001) where TDN intake (kg) is the discounted TDN intake; RDP intake (kg) is the RDP intake. These constraints present challenges to formulating diets with SQP, which we discuss in the following section. Additionally, the SQP algorithm converges to a local optimum (the optimum within a neighboring set of solutions), which could be far away from the global optimum (the optimum among all possible solutions) solution in some complicated nonlinear problems (Boggs and Tolle, 1995). MINLP_DGO. Another optimization strategy, MINLP_DGO, allows the use of a mixture of continuous and binary decision variables. Note that mixed-integer nonlinear programming (MINLP) includes a wide range of nonlinear optimization problems containing continuous, binary, and integer variables. Solving an MINLP problem usually involves multiple algorithms and techniques (Kronqvist et al., 2019). Several solvers, including Couenne (Burer, 2009), BARON (Tawarmalani and Sahinidis, 2013), and Gurobi (Gurobi Optimization LLC, 2021), find deterministic and global solutions in MINLP problems. Formulating least-cost diets using MINLP can be written thus: where y is a binary variable vector. Using binary variables enables the conversion of Eq.
[1], the intake level calculation can be written into 3 constraints: where M is a large positive number (suppose M = 100,000) and y 1 is a binary variable. When y 1 is equal to 0, Eq.
[7] is ineffective, because it models an intake level greater than a very small number and smaller than a very large number; Eq.

Simulation
A simulation study was conducted to compare the 3 optimization strategies. We simulated 3 animal groups: lactating cows, dry cows, and heifers. The animal inputs needed to define the nutrient requirement constraints are summarized in Table 1. A total of 500 simulations were conducted by randomly selecting a set of feeds containing (arbitrarily) 5 concentrates and 3 forages from 100 feeds (feeds and prices shown in Supplemental File S2; Li, 2021) in the feed table provided by NRC (2001). Calcium phosphate (monobasic) was kept for all 500 sets of feeds, to minimize excess mineral feeding. The 2020 annual average feed prices were taken from the Pennsylvania State University feed price list (https: / / extension .psu .edu/ files/ feed -price -lists; State University Cooperative Extension, 2020). For each animal group, 3 least-cost diets were formulated based on each set of the randomly selected feeds using iteLP, SQP, and MINLP_DGO, respectively, which gave 4,500 (3 × 3 × 500) diets in total.

Evaluation
To demonstrate the validity of each optimizer, all the formulated diets were evaluated with the NRC (2001) equations to examine the feed nutrient supply. The difference between requirement and supply was calculated for all the nutrient constraints. In addition, we compared the diet cost and computation time for 3 optimizers implemented with an i7-5500U computer processor (2.40 GHz). We report the results of diet cost and computation time as mean (±SD).

Software
All simulations and computations were conducted in Python 3.7.4 (Van Rossum and Drake, 2009). We used the SciPy package (Virtanen et al., 2020) as the solvers for iteLP and SQP. For MINLP_DGO, Gurobi (Gurobi Optimization LLC, 2021) was used because it is free to academic users and can be easily implemented in Python.

An Example Simulation
A set of feeds in one simulation is shown in Table 2. The diets formulated based on the feeds by iteLP, SQP, and MINLP_DGO are shown in Table 3. For lactating cows and dry cows, the rations obtained by SQP and MINLP_DGO were the same, but those obtained by iteLP were more expensive. For heifers, the rations obtained by all 3 optimizers were very close.

Feasibility and Nutrient Balance
In formulating diets for lactating cows, 57 simulations for iteLP were infeasible, 6 of which were infeasible because the SciPy package could not find a feasible solution in the first several iterations. The other 51 simulations were infeasible because the maximum number of iterations was reached before the convergence of intake level and MCP. In these cases, the value of intake level or MCP oscillated between 2 or several values and failed to converge during the iteration process ( Figure  1), which highlights a flaw of iteLP. For MINLP_DGO, 21 simulations had no feasible solutions. To investigate the reasons for infeasibility, the 21 simulations were rerun after adding 2 dummy feed variables, including a "super protein" feed with 100% of protein and a "super energy" feed with 30 Mcal of digestible energy. Both feeds were set to be very expensive (i.e., $100/kg of DM), so that they would not be used when other feeds can fulfill the protein and energy requirements. Their  appearance in the solution indicates a lack of protein or energy from feeds. The results showed that all 21 simulations were infeasible due to protein deficiency. The same 21 simulations were also infeasible for SQP, and 24 additional simulations were infeasible due to failure to fulfill the protein requirement. These 24 simulations were feasible when using MINLP_DGO, which means SQP did not fully explore the feed potential in these cases. Optimization with SQP relies on the gradient of the objective function and constraints, but the MP constraint is not continuous due to the discrete choice between prediction of MCP based on energy or protein availability represented in Eq.
[5]. The SciPy package calculates the gradient numerically instead of analytically, which makes fitting a discrete constraint into SQP possible. In the MP constraint, the MCP production is limited by either TDN intake or RDP intake, which creates 2 discrete value domains. In most cases, the optimum exists in one of the domains and is far away from the other, so the SQP solver is able to search for the optimum within one domain without influence from discreteness. However, sometimes when the optimum is close to the boundaries of the 2 domains, SQP may not be able to find a proper solution based on the gradient, due to the discreteness on the optimum direction. However, by including binary variables, MINLP_DGO is able to find the optimum through the branch-andbound approach (Land and Doig, 1960), which is able to evaluate the suboptimal solutions in each domain Li et al.: NONLINEAR PROGRAMMING FOR RATION FORMULATION Figure 1. An example of microbial crude protein (MCP) production oscillating between several values and failing to converge when using iterative linear programming (iteLP) to design feed formulation for lactating cows (using corn grain, cotton seed meal, corn distillers, corn cob, canola meal, mixed grass-legume hay, alfalfa meal, mixed grass-legume silage, and calcium phosphate monobasic). and find the best one (Taylor, 2009). Heuristically, MINLP_DGO will consider both TDN-limited and RDP-limited cases for MCP production and select the better solution. The requirement of NE L did not result in any infeasible solutions for SQP, even though the NE L constraint contained a nondifferentiable function (Eq. [1]). Because IntakeLevel represents the incremental intake above maintenance, the IntakeLevel value of diets fulfilling the requirements of NE M , NE L , and NE G at the same time would always be greater than 1, and the nondifferentiable point (IntakeLevel = 1) would not disrupt the solution. All the infeasible simulations for SQP and MINLP_DGO were also infeasible for iteLP. We also investigate the feasibility of different lactating groups by changing the milk production level (results not shown). The comparison among the 3 optimization strategies was similar, but more infeasible simulations were observed for higher milk production groups, who demand more protein, so it is more difficult to obtain a solution with the same feeds. When formulating diets for dry cows and heifers, a feasible solution was obtained for all the simulations using SQP and MINLP_DGO, but 49 and 11 infeasible simulations, respectively, occurred using iteLP. For dry cows and heifers, the protein requirement did not cause infeasibility in SQP, probably due to much smaller nutrient requirements of dry cows and heifers. All the infeasible simulations were removed for the following analysis.
Boxplots of the difference between requirements and estimated nutrient delivery of the ration solutions for lactating cows from all feasible simulations are shown in Figure 2. Only NE L and MP balance values are shown, because the other constraints were linear and did not differ greatly between the 3 optimizers. The NE L and MP balance values were positive for the solutions obtained by iteLP, SQP, and MINLP_DGO in all feasible simulations, but the NE L balance values were higher for iteLP, which suggests a tendency of overfeeding energy for rations formulated with iteLP. Differences between the rations simulated for dry cows and heifers ( Figure  3 and 4) by iteLP, SQP, and MINLP_DGO were similar to those for lactating cows, except that iteLP also showed a tendency of overfeeding protein.

Ration Cost
Average ration costs per animal obtained by iteLP, SQP, and MINLP_DGO were $4.78 (±0.71), $4.45 for heifers, respectively (Figure 2, 3, and 4). Ration costs were similar for SQP and MINLP_DGO, which indicates that the solution obtained through SQP was very close to the global optimum in ration optimization. However, the costs obtained by iteLP were greater, especially for dry cows. The iteration process in iteLP does not push the solution along the optimum direction, but blindly replaces the initial value with the solution from last iteration; therefore, the ration cost is minimized only within each iteration instead of the whole process and is likely to be greater than the one obtained by SQP or MINLP_DGO.
The least-cost diets for a given level of milk production estimated by the algorithms we present here have practical applications for the dairy industry and can be implemented in whole-farm simulation studies. However, other methods of diet optimization, such as maximization of income over feed costs or income minus feed costs, are important to consider as well. In our study, milk production serves as an input and thus is not estimated in the optimization process. To use maximization of milk income minus feed cost as the objective function, instead of minimization of feed cost, we would need to consider the interactions among milk production, milk composition, ration formulation, and DMI, the complexity of which requires future studies to investigate.

Computation Time
The  (Figure 2, 3, and 4). The computation time of iteLP and SQP was much shorter than that of MINLP_DGO, especially when formulating diets for lactation cows. Because MINLP_DGO solves the problem involving binary variables, the resulting problem becomes computationally expensive. When the solver Gurobi solves MIMLP problems, the best bound and the suboptimal solution are updated throughout the computation. The best bound is the upper bound of the objective function value for a maximization problem  (or the lower bound of a minimization problem), and the suboptimal solution is the best feasible solution found so far. By default, the computation stops when the gap between the 2 values decreases below 0.01%, but it may take hours to obtain such a solution. The computation time was limited to 3 min in this study. In total, 98, 21, and 33 simulations for lactating cows, dry cows, and heifers, respectively, exceeded the time limit (computation time of 3 min), but the gaps for them were all less than 0.1%, meaning the solution found was close to the true optimum.

CONCLUSIONS
The linear optimization iteLP had limited capability to formulate least-cost diets when nonlinearity existed in the constraints. Both SQP and MINLP_DGO were able to handle nonlinear constraints well, with SQP being faster, but MINLP_DGO was able to return a feasible solution under some situations where SQP could not. Thus, both nonlinear programming frameworks for least-cost ration formulation represent advancement over previous linear programming techniques. In addition, the nonlinearities in the upcoming NRC revision can likely be addressed using the techniques described in this paper.