The purpose of this paper is to present a Monte Carlo solution for a random objective function coefficient linear programming problem that can be executed in Excel. A solution was given in Ridley and Khan  for random constraint limits. Both present an exceptionally easy way for students to learn Monte Carlo simulation. One of the emerging topics in business management is risk analysis. Risk analysis and business analytics are emerging topics in education and preparation for the knowledge economy. The textbook “Business Analytics (Evans  )” is designed to be Excel based. It covers useful Excel functions and illustrates how they are used to execute Monte Carlo simulation for the purpose of risk analysis. However, the chapter on linear optimization collapses to traditional deterministic linear programming where the solution tool is the add-in solver. Risk analysis in linear programming as described there relies on sensitivity analysis. In sensitivity analysis one objective function coefficient value at a time is varied while all others are held constant. In the real world, the coefficients can vary randomly and simultaneously. This paper describes a randomized linear program (LP) where all objective function coefficients are varied simultaneously. The random values are selected in accordance with any specified probability distribution. The distribution can be aligned with what occurs in reality. In this paper we present a novel approach for designing a copula of random objective function coefficients according to a specified rank correlation.
Amongst the Monte Carlo simulations are some that are feasible and some that are not. The feasible solutions are kept, and the infeasible solutions are discarded and disregarded. The set of objective function values of the solutions that are kept represent the distribution of possible objective functions values. In the typical business problem, they are profits derived from random prices. The analysis of these profits represents business risk analysis. It tells management the long term expected profit to plan for. It also tells the management the probability of low to no profits for which cash flow may be disrupted. If the expected profit is positive and acceptable, arrangement can be made to tide the business over during lean times.
The remainder of the paper is organized as follows. Section 2 is a review of related literature. Section 3 introduces the example problem used for illustration. Section 4 gives the traditional graphical solution. Section 5 gives the traditional algebraic solution and illustrates the new randomized LP solution. Section 6 summarizes some conclusions and suggestions for future research.
2. Literature Review
Linear programming can be traced back to the 1940s. Dantzig  was the creator of the simplex method for solving linear programming problems. Since then, linear programming has been applied in many fields, in business, transportation, military, etc. The Simplex method is deterministic, and only can be used to solve problems subject to no uncertainty in variables or coefficients in the objective function and constraint equations. Alternative methods have been developed to solve large problems, like the interior point algorithm (Roos et al.  ).
Given that linear programming is deterministic, one way developed to deal with uncertainty is post optimality analysis, also known as sensitivity analysis. Sensitivity analysis is a “what if” analysis where the effect on the optimal solution is measured when changes in parameters are applied. A full chapter on sensitivity analysis may be found in Eiselt et al. . A more up to date version can be found in Panik .
Some degree of uncertainty is permitted in the Fuzzy Programming method where the desired objective value is assumed to be stated in an ambiguous way (Bellman et al. ; Inuiguchi et al. ; Sakawa  ). Another approach for dealing with uncertainty is Stochastic Programming (SP), defined by Dantzig and Thapa  as the “Art and Science of deciding on the best plan of action (in some expected-value sense) while hedging against the myriad of possible ways the best laid plans can go awry”. These stochastic methods consider the expected values of the objective function coefficients when they are unknown at the time of decision making. See also Boyd, Boyd and Vandenberghe .
A randomized linear program (RLP) for computing network bit prices, where itinerary demand realization sequences are simulated and solved using deterministic linear programming (DLP) was given by Talluri et al. . They demonstrate that RLP is only a little more complicated to implement than DLP. Risk analysis using DLP and Monte Carlo simulation was applied to manage municipal solid waste (Wajs et al.  ). Uncertainty was introduced into the decision variables. Adler et al.  proposed an extension of Clarkson’s  randomized algorithm for linear programming to a general scheme for solving convex optimization problems. The aim is to speed-up the simplex method (or any vertex enumeration method) for linear programming when the number of constraints is much larger than the number of variables. Mohammed and Kassem  consider a product mix problem in which several scenarios are presented for the same problem considering two important aspects, namely, resources’ utilization and productivity. But these stochastic methods apply linear programming to the expected values of the objective function coefficients. That is, these stochastic methods do not consider solutions for all values of the random coefficient. This paper does.
3. The Linear Program
The general LP may be written as
where x = (x1, x2, x3, ∙∙∙, xn) is a vector of decision variables, p = (p1, p2, p3, ∙∙∙, pn) is a vector of independent profit contributions, b = (b1, b2, b3, ∙∙∙, bm) is a vector of constraint limits, and A is an m x n matrix of constants.
Consider an example taken from Evans [2, p. 474]. The name of the business is Slenka Ski. Slenka Ski produces two type of Skis, Jordanelle and Deercrest. The per unit profit contributions are $50 and $65 respectively. There is a fabrication constraint of 84 hours and a finishing constraint of 21 hours. Jordanelle uses 3.5 hours of fabrication time and 1 hour of finishing time. Deercrest uses 4.0 hours of fabrication time and 1.5 hours of finishing time. There is also a market mix requirement that Deercrest production must be no less than twice that of Jordanelle. As always, there are the nonnegativity constraints.
Let Jordanelle be represented by x1 and Deercrest be represented by x2.
4. Graphical Solution
The deterministic graphical solution is given in Figure 1, where x1 = 5.25, x2 = 10.5, and the Objective function value (profit) = $945. The fractional amounts of x1 and x2 are interpreted as work in progress.
Figure 1. Graphical solution of deterministic linear programming program.
5. Algebraic Solutions
The first algebraic solution will be created in an Excel spreadsheet and executed by the add-in SOLVER. A spreadsheet for the production of Slenka Skis is shown in Figure 2.
5.1. Deterministic Solution
Rewriting the LP in terms of the Excel fields, the problem formulation is
Figure 2. Monte Carlo computer simulations of linear programming profit for random objective function coefficients.
Initially, all parameters of the LP, including the objective function coefficients, are assumed to be known at the time that a decision is to be made regarding the optimal values of the decision variables. The original objective coefficients (not shown to save space) are 50 for Jordanelle and 65 for Deercrest. This LP was solved by the Excel add-in SOLVER. The solution is x1 = (B14 = 5.25), x2 = (C14 = 10.5). The profit for this configuration is $945 (not shown).
5.2. Monte Carlo Simulation
The Monte Carlo computer simulation is based on 1000 sets of two objective function coefficients. The profit contribution from Jordanelle is assumed to be lognormally distributed with mean $4.5 and standard deviation $0.5. The profit contribution from Deercrest is assumed to be normally distributed with a mean of $65 and a standard deviation of $1. The random coefficient numbers are generated from the lognorm.inv(rand(),4.5,0.5) and norm.inv(rand(),65,1) Excel functions respectively. All negative prices are discarded and ignored. So long as the slope of the objective function lies on or between the slopes of the binding constraints, the decision variables will remain constant, but the profit will vary. Higher prices will yield higher profits and lower prices will yield lower profits. For each of the 1000 simulations, each time the slope objective function (H9 = −(B9 = 176.67/C9 = 65.52) = −2.697) lies on or between the slopes of the binding constraints (H7 = −(B7 = 1.0/C7 = 1.5) = −0.66667) and (H8 = −(B8 = −2/C8 = 1) = 2, the profit is calculated and saved in D22. Otherwise, a profit of null (“”) is entered in D22. The profit in D22 is then copied into the Data Table(columns M&N). A new column of profits is then copied into column Q where the nulls are filtered out with the Excel function: (=FILTER(N4:N1003,N4:N1003<>""). Hence forth the nulls are discarded and disregarded. To save space only the first 46 profits are shown. The Excel logic and calculations are as follows:
The row (7) containing the first binding constraint for which the slack = 0 is obtained from H27 = MATCH(0,$F$1:$F$8,0).
The row (8) containing the second binding constraint for which the slack =0 is obtained from H28 = VLOOKUP(0,$F$7:$H$9,3).
The profit is obtained from =IF(OR(AND($H$9<0,$H$9<$H$7),AND($H$9>0,$H$9>$H$8)),SUMPRODUCT($B$9:$C$9,$B$14:$C$14),"").
The objective function coefficients are restricted to positive values only so the slope is always negative and therefore less than 2.
The first entry in the Data Tableis N2 = D22.
5.3. Profit Distribution
As the number of simulations increases, the distributions of profit functions are expected to converge and become similar. In the interim, two sample simulations are shown here in Figure 2. If the distribution converges to the one on the left, the expected profit will be lower and there will be times when the profit is low enough to possibly require reserve rainy day funds. If the low expected profit is acceptable then all that is required is the capacity to borrow money during the lean times. That is, to hedge against this contingency. If the low expected profit is unacceptable, then this business venture should not be undertaken. If the distribution converges to the one on the right, the expected profit will be higher. In that case, it is less likely that rainy day funds will be required.
5.4. Correlated Profit Contributions
For real-world problems, the assumption of independence among inputs may not be appropriate. Many studies have been dedicated to properly incorporating dependence among input variables (see, e.g. Iman and Conover ; Chang, Tung and Yang ; Heal and Kunreuther  ). Studies in retirement planning (Benninga  ), value-at-risk (Bohdalová  ), default loss modeling (Duffie, Eckner, Horel and Saita  ), and bankruptcy prediction Altman  stress the importance of appropriately modeling correlation among model parameters, control variables, design variables, etc. During the 2008 financial crisis, for example, failure to properly model correlation among housing defaults (Herbert  ) resulted in widespread credit ratings errors on tranched mortgage-backed securities and CDOs (Duffie, Eckner, Horel and Saita ; Securities and Exchange Commission  ).
So far, this paper modifies the Evans  Selena Ski problem to assume that per unit profit across goods is random and independent. However, profit margins may be correlated with the business cycle and exhibit seasonality (Chevalier, Kashyap and Rossi ; Hoskin, Matsa and Reiffen  ), resulting in positive correlation between profit margins of similar items. Also, profit margins may be correlated with advertising expenditures (Farris and Reibstein  ), possibly resulting in negative correlation between profit margins of Jordanelle brand and Deercrest brand if advertising strategies vary inversely.
The procedure we use to model dependence between profit variables is derived from theory related to bivariate Gaussian copulas. A copula “couples” or joins marginal distributions of single random variables into a joint distribution with a specified dependence structure. In general, a copula is a d-dimensional distribution function on [0, 1]d with uniform marginal distributions.
Our goal is to generate instances of (P1, P2) in Excel of a bivariate distribution with specified marginal distributions, and rank correlation between variables equal to a desired value τ.
The procedure is outlined below:
Step 1: Generate standard normal random variables Z1, Z2 with rank correlation τ.
Step 2: Apply the standard normal CDF to Z1, Z2, which produces uniform random variables U1,U2.
Step 3: Apply inverse CDF’s , of the desired marginal distributions to U1,U2, producing P1, P2. (see Figure 3).
In Section 5.2, we assumed per unit Jordanelle profit and per unit Deercrest profit . In this section, we will generate 10,000 instances of a bivariate distribution such that the marginal distributions are lnN(3.78, 0.5)1 (Jordanelle) and N(65, 5)2 (Deercrest). Importantly, we target a dependence structure between Jordanelle and Deercrest profits. By way of example, we target a rank correlation of τ = −0.7. A convenience of using Gaussian copulas is that there is a one-to-one relationship between linear correlation measured by Pearson’s rho (ρ) and rank correlation measured by Kendall’s tau (τ)3 for the bivariate normal:
The inverse transformations (lognormal inverse) and (normal inverse) are monotonic, so rank correlation is preserved. We can therefore target rank correlation ofP1and P2 equal to −0.7 by generating
Figure 4 demonstrates the procedure in Excel.
Step 1: Values of the standard normal random variableZ1 in column B are generated by norm.s.inv(rand()). Each instance ofZ2 in column C is produced by
A simple mathematical proof can showZ2 is distributed standard normal, with linear correlation to Z1 equal to ρ'. The targeted rank correlation τ = −0.7 is hardcoded into cell (C13). The necessary linear correlation between Z1 and Z2 (ρ') to produce rank correlation τ is shown in cell (C14 = (sin(C13*pi()/2) = −0.8910).
Figure 3. Flow chart for steps 1, 2 and 3.
Figure 4. The bivariate gaussian copula procedure to generate jointly distributed random variables.
Step 2: Values of U1 and U2 in columns E and F are produced by applying the cumulative standard normal distribution function to Z1 in column B and to Z2 in column C respectively. For instance, cell (E18 = (NORM.S.DIST(B18, true)) = 0.9942).
Step 3: Values P1 and P2 that represent instances of Jordanelle and Deercrest correlated profits are shown in columns J and K. Each instance of P1 is produced by applying the lognormal distribution inverse to U1. For example, cell (J18 = (LOGNORM.INV(E18, $L$10, $M$10)) = 154.7867). Each instance of P2 is produced by applying the normal distribution inverse to U2. For example, cell (K18 = (NORM.INV(F18, $L$11, $M$11)) = 52.8526). Due to monotonicity, the rank correlation is preserved, and is hardcoded into cell M13. The linear correlation of P1 and P2 is computed using the Excel function CORREL() (M14 = (CORREL(J18:J5017, K18:K5017) = -0.8390. Observe that this value differs from the linear correlation of Z1 and Z2 (C14 = (SIN(tau*PI()/2)) = −0.8910). (see flow chart in Figure 3).
Figure 5 describes the distribution of total profit if per unit profit margin for Jordanelle and Deercrest skis are random, distributed lognormally and normally respectfully, and per unit profit margins are correlated. As in the previous section, we include only instances for which the optimal solution is 5.25 units of Jordanelle skis and 10.50 units of Deercrest skis by bounding the slope of the objective function.
The charts within Figure 5 demonstrate the importance of accounting for possible correlation between variables when designing a model. As the rank correlation between Jordanelle and Deercrest per unit profit margins increases from −0.7 to 0.7, the shape of the total profit distribution has changed. The distribution in the τ = 0.7 case is wider with a longer right tail. The range of possible
Figure 5. Impact of parameter correlation on profit distribution.
total profit outcomes increases to $1752.74 from $1315.53 (33%). The standard deviation of the distribution increases to $158.4 from $105.1 (51%). Without accounting for correlation among the random per-item profits in the model, the increase in standard deviation (a widely accepted measure of risk in financial management) of the total profit distribution would have been overlooked. The illustrations in this paper are pedagogical in that they build on the example problem in Evans . In a real-life problem, all numbers in the problem including the rank correlation would be replaced with values estimated from actual historical data. The solution is scalable to more than two variables.
Businesses that are dealing with optimal production decisions of complementary cohorts of goods need to be particularly cautious about what Rust et al.  refer to as “profitable product death spiral” that may ensue when businesses rush to discontinue production of a product with low profit margin without fully assessing potential spiraling effect of such an action on the sales of the other complementary products of the same cohort. Others (e.g., Cannon, et al.  ) echo this view and predict suboptimal long-term profit for the incumbent company in a similar situation if consumers who buy one product from the company expect the company to supply other complementary products as well. The approach to randomized objective function linear programming presented here is well suited to address the issue of uncertainty in risk management. Instead of focusing on short term profits associated with one product at a time, the long run expected profit can be determined from the distribution of profits from multiple products with multiple cross correlations. The variance in profit can also be calculated. That way, businesses can determine the recommended amount of overdraft funds to cover periods of low profit.
The profit distribution from the model in this research was generated from a single combination of fixed product quantities obtained from a single execution of the deterministic LP solver. Thereafter only the objective function coefficient values changed randomly. A suggestion for further research is to extend randomization to include both objective function and constraints. In that case, the LP solver will need to be executed for each new combination of objective function and constraints to update the production quantities as well as the profits.
1We change μ = 3.78 from μ = 4.5, so that .
2We increase the standard deviation of Deercrest profit to $5 in this section so the impact of correlation on total profit distribution is more apparent.
3Such a one-to-one relationship also exists between linear correlation Pearson’s rho (ρ) and rank correlation Spearman’s rho (ρs).
 Ridley, A.D. and Khan, A. (2002) Randomized Constraint Limit Linear Programming. Journal of Applied Mathematics and Physics, 8, 2691-2702.
 Dantzig, G.B. (1951) Application of the Simplex Method to a Transportation Problem. In: Koopmans, T.C., Ed., Activity Analysis of Production and Allocation, John Wiley and Sons, New York, 359-373.
 Inuiguchi, M. and Ramík, J. (2000) Possibilistic Linear Programming: A Brief Review of Fuzzy Mathematical Programming and a Comparison with Stochastic Programming in Portfolio Selection Problem. Fuzzy Sets and Systems, 111, 3-28.
 Sakawa, M. (2002) Fuzzy Multiobjective and Multilevel Optimization. In: Ehrgott, M. and Gandibleux, X., Eds., Multiple Criteria Optimization: State of the Art Annotated Bibliographic Survey, Kluwer, Netherlands, 171-226.
 Wajs, W., Bieda, B. and Tadeusiewicz, R. (2001) Linear Programming and Risk Analysis Methods for Municipal Solid Waste Decision Support System. IFAC Proceedings Volumes, 34, 187-192.
 Adler, I. and Shamir, R. (1993) A Randomized Scheme for Speeding up Algorithms for Linear and Convex Programming Problems with High Constraints-to-Variables Ratio. Mathematical Programming, 61, 39-52.
 Clarkson, K.L. (1988) A Las Vegas Algorithm for Linear Programming When the Dimension Is Small. Proceedings of the 29th IEEE Symposium on the Foundations of Computer Science, White Plains, 24-26 October 1988, 452-456.
 Mohammed, A.R. and Kassem, S.S. (2020) Product Mix Optimization Scenarios: A Case Study for Decision Support Using Linear Programming Approach. 2020 International Conference on Innovative Trends in Communication and Computer Engineering, Aswan, 8-9 February 2020, 50-55.
 Iman, R.L. and Conover, W. (1982) A Distribution-Free Approach to Inducing Rank Correlation among Input Variables. Communication in Statistics-Simulation and Computation, 11, 311-334.
 Chang, C.H., Tung, Y.K. and Yang, J.C. (1994) Monte Carlo Simulation for Correlated Variables with Marginal Distributions. Journal of Hydraulic Engineering, 120, 313-331.
 United States Securities and Exchange Commission (2008) Summary Report of Issues Identified in the Commission Staff’s Examinations of Select Credit Rating Agencies. United States Securities and Exchange Commission, Washington DC.
 Chevalier, J.A., Kashyap, A.K. and Rossi, P.E. (2003) Why Don’t Prices Rise during Periods of Peak Demand? Evidence from Scanner Data. American Economic Review, 93, 15-37.
 Cannon, J.N., Cannon, H.M. and Low, J.T. (2013) Modeling Tactical Product-Mix Decisions: A Theory-of-Constraints Approach. Simulation & Gaming, 44, 624-644.