Homework 3 – Solution
Homework 3 – Solution¶
OPIM 5641: Business Decision Modeling – University of Connecticut
Copyright By PowCoder代写 加微信 powcoder
Due Monday, February 28, 11:59pm
Please add detailed comments to all code so that I can follow your solution. If I cannot understand your code then you will not receive full credit.
Your Name Here
Your NetID Here
%matplotlib inline
from pylab import * # for graphing
import shutil
import sys
import os.path
if not shutil.which(“pyomo”):
!pip install -q pyomo
assert(shutil.which(“pyomo”))
if not (shutil.which(“cbc”) or os.path.isfile(“cbc”)):
if “google.colab” in sys.modules:
!apt-get install -y -qq coinor-cbc
!conda install -c conda-forge coincbc
assert(shutil.which(“cbc”) or os.path.isfile(“cbc”))
from pyomo.environ import *
import numpy as np
import pandas as pd
|████████████████████████████████| 9.6 MB 4.9 MB/s
|████████████████████████████████| 49 kB 2.8 MB/s
Selecting previously unselected package coinor-libcoinutils3v5.
(Reading database … 155320 files and directories currently installed.)
Preparing to unpack …/0-coinor-libcoinutils3v5_2.10.14+repack1-1_amd64.deb …
Unpacking coinor-libcoinutils3v5 (2.10.14+repack1-1) …
Selecting previously unselected package coinor-libosi1v5.
Preparing to unpack …/1-coinor-libosi1v5_0.107.9+repack1-1_amd64.deb …
Unpacking coinor-libosi1v5 (0.107.9+repack1-1) …
Selecting previously unselected package coinor-libclp1.
Preparing to unpack …/2-coinor-libclp1_1.16.11+repack1-1_amd64.deb …
Unpacking coinor-libclp1 (1.16.11+repack1-1) …
Selecting previously unselected package coinor-libcgl1.
Preparing to unpack …/3-coinor-libcgl1_0.59.10+repack1-1_amd64.deb …
Unpacking coinor-libcgl1 (0.59.10+repack1-1) …
Selecting previously unselected package coinor-libcbc3.
Preparing to unpack …/4-coinor-libcbc3_2.9.9+repack1-1_amd64.deb …
Unpacking coinor-libcbc3 (2.9.9+repack1-1) …
Selecting previously unselected package coinor-cbc.
Preparing to unpack …/5-coinor-cbc_2.9.9+repack1-1_amd64.deb …
Unpacking coinor-cbc (2.9.9+repack1-1) …
Setting up coinor-libcoinutils3v5 (2.10.14+repack1-1) …
Setting up coinor-libosi1v5 (0.107.9+repack1-1) …
Setting up coinor-libclp1 (1.16.11+repack1-1) …
Setting up coinor-libcgl1 (0.59.10+repack1-1) …
Setting up coinor-libcbc3 (2.9.9+repack1-1) …
Setting up coinor-cbc (2.9.9+repack1-1) …
Processing triggers for man-db (2.8.3-2ubuntu0.1) …
Processing triggers for libc-bin (2.27-3ubuntu1.3) …
/sbin/ldconfig.real: /usr/local/lib/python3.7/dist-packages/ideep4py/lib/libmkldnn.so.0 is not a symbolic link
Problem 1 (20 points)¶
Solve the following LP model using the graphical method.
$Max(Z) = 30X + 32Y$
subject to:
$9X + 9Y \leq 810$
$X + 2Y \geq 50$
$X \geq 10$
$Y \geq 40$
$X,Y \geq 0$
The optimal solution is $X = 10$, $Y = 80$ and $Z = 2860$.
Answer to Problem 1¶
We create the feasible region.
# this is the example from Pyomo cookbook
# pylab makes it easy to make plots
figure(figsize=(6,6))
#subplot(111, aspect=’equal’)
axis([0,100,0,100])
xlabel(‘X’)
ylabel(‘Y’)
# Constraint 1
y = array([90, 0])
x = array([0,90])
plot(x,y,’r’,lw=2)
fill_between([0,90,100], # x area
[90,0,0], # y area
[100,100,100], # upper area
color=’red’, # color
alpha=0.15) # transparency
# Constraint 2
y = array([25, 0])
x = array([0,50])
plot(x,y,’orange’,lw=2)
fill_between([0,50,100], # x area
[0,0,0], # y area
[25,0,0], # upper area
color=’orange’, # color
alpha=0.15) # transparency
# Constraint 3
y = array([0, 100])
x = array([10,10])
plot(x,y,’green’,lw=2)
fill_between([0,10], # x area
[0,0], # y area
[100,100], # upper area
color=’green’, # color
alpha=0.15) # transparency
# Constraint 4
y = array([40, 40])
x = array([0,100])
plot(x,y,’black’,lw=2)
fill_between([0,100], # x area
[0,0], # y area
[40,40], # upper area
color=’black’, # color
alpha=0.15) # transparency
# legend([‘Carpentry Constraint’,’Painting Constraint’,’ Constraint’, ‘ Constraint’])
legend([‘Constraint 1′,’Constraint 2′,’Constraint 3′,’Constraint 4′])
#,’Painting Constraint’,’ Constraint’, ‘ Constraint’])
# the area in white is the feasible region!
text(11,43,’FEASIBLE REGION’)
Text(11, 43, ‘FEASIBLE REGION’)
Based on the feasible region, we have to test 3 points, defined at:
The intersection of the black line and the green line
The intersection of the red line and the green line
The intersection of the black line and the red line
We will find those points, determine their solution quality, and pick the best.
A = np.array([[1, 0], [0, 1]])
B = np.array([10, 40])
point_1 = np.linalg.solve(A,B)
print(“Point 1: “, point_1)
A = np.array([[1, 0], [9, 9]])
B = np.array([10, 810])
point_2 = np.linalg.solve(A,B)
print(“Point 2: “, point_2)
A = np.array([[0, 1], [9, 9]])
B = np.array([40, 810])
point_3 = np.linalg.solve(A,B)
print(“Point 3: “, point_3)
# Now we evaluate:
objective_1 = 30*point_1[0] + 32*point_1[1]
print(objective_1)
objective_2 = 30*point_2[0] + 32*point_2[1]
print(objective_2)
objective_3 = 30*point_3[0] + 32*point_3[1]
print(objective_3)
Point 1: [10. 40.]
Point 2: [10. 80.]
Point 3: [50. 40.]
(10,80) is the optimal solution with optimal value 2860.
Problem 2 (20 points)¶
Solve the following LP model using the graphical method.
$Max(Z) = 15X + 10Y$
subject to:
$3X + 5Y \leq 45$
$4X + Y \leq 20$
$X,Y \geq 0$ (nonnegativity)
The answer is $X=3.23529411764706$, $Y=7.05882352941176$ and $Z = 119.11764705882351$.
Answer to Problem 2¶
We create the feasible region.
figure(figsize=(6,6))
#subplot(111, aspect=’equal’)
axis([0,30,0,30])
xlabel(‘X’)
ylabel(‘Y’)
# Constraint 1
y = array([9, 0])
x = array([0,15])
plot(x,y,’r’,lw=2)
fill_between([0,15,30], # x area
[9,0,0], # y area
[30,30,30], # upper area
color=’red’, # color
alpha=0.15) # transparency
# Constraint 2
y = array([20, 0])
x = array([0,5])
plot(x,y,’orange’,lw=2)
fill_between([0,5,30], # x area
[20,0,0], # y area
[30,30,30], # upper area
color=’orange’, # color
alpha=0.15) # transparency
# legend([‘Carpentry Constraint’,’Painting Constraint’,’ Constraint’, ‘ Constraint’])
legend([‘Constraint 1′,’Constraint 2′])
#,’Painting Constraint’,’ Constraint’, ‘ Constraint’])
# the area in white is the feasible region!
text(11,43,’FEASIBLE REGION’)
Text(11, 43, ‘FEASIBLE REGION’)
Based on the feasible region, we have to test 4 points, defined at:
The intersection of both axes
The intersection of the y axis and the red line
The intersection of the orange line and the red line
The intersection of the x axis and the orange line
We will find those points, determine their solution quality, and pick the best.
A = np.array([[1, 0], [0, 1]])
B = np.array([0, 0])
point_1 = np.linalg.solve(A,B)
print(“Point 1: “, point_1)
A = np.array([[1, 0], [3, 5]])
B = np.array([0, 45])
point_2 = np.linalg.solve(A,B)
print(“Point 2: “, point_2)
A = np.array([[4, 1], [3, 5]])
B = np.array([20, 45])
point_3 = np.linalg.solve(A,B)
print(“Point 3: “, point_3)
A = np.array([[4, 1], [0, 1]])
B = np.array([20, 0])
point_4 = np.linalg.solve(A,B)
print(“Point 4: “, point_4)
# Now we evaluate:
objective_1 = 15*point_1[0] + 10*point_1[1]
print(objective_1)
objective_2 = 15*point_2[0] + 10*point_2[1]
print(objective_2)
objective_3 = 15*point_3[0] + 10*point_3[1]
print(objective_3)
objective_4 = 15*point_4[0] + 10*point_4[1]
print(objective_4)
Point 1: [0. 0.]
Point 2: [0. 9.]
Point 3: [3.23529412 7.05882353]
Point 4: [5. 0.]
119.11764705882354
(3.23529412,7.05882353) is the optimal solution with optimal value 119.11764705882354.
Problem 3 (20 points)¶
Formulated the dual problem for Problem 1 and explain each step.
The primal problem is
$Max(Z) = 30X + 32Y$
subject to:
$9X + 9Y \leq 810$
$X + 2Y \geq 50$
$X \geq 10$
$Y \geq 40$
$X,Y \geq 0$
We are looking to prove the best bound we can by taking linear combinations of the constraints. To do this, we will assign a new variable for each constraint (there are 4), ensure that the combination will dominate our objective function, and minimize the bound
Before we do that, we have to convert the linear programming model so that it has inequalities all facing in the same direction, so that when we combine with positive weights, we end up with the inequality in the right direction!
$Max(Z) = 30X + 32Y$
subject to:
$9X + 9Y \leq 810$
$-X – 2Y \leq -50$
$-X \leq -10$
$-Y \leq -40$
$X,Y \geq 0$
Now we can write the linear program!
$MIN(Z) = 810 z_1 – 50 z_2 – 10 z_3 – 40 z_4$
subject to:
$9 z_1 – z_2 – z_3 \geq 30$
$9 z_1 – 2z_2 – z_4 \geq 32$
Note: Although you didn’t have to solve it as a part of your homework, if you do solve this linear program, you end up with an objective value of 2860 at this solution:
(3.555556,0,2,0)
This thereby PROVES that the solution found via the graphical method is the best possible solution.
Problem 4 (40 points)¶
A farmer is considering how much of each of 3 varieties of corn to produce in his plot of land which is 20 acres in size. Every acre of land dedicted to a variety produces on average:
Variety 1: 170 bushels
Variety 2: 180 bushels
Variety 3: 185 bushels
The yield is however uncertain. The uncertainty is as follows:
Variety 1: plus/minus 15\% of dedicted land
Variety 2: plus/minus 17\% of dedicted land
Variety 3: plus/minus 20\% of dedicted land
The profit per bushel is:
Variety 1: \$4.40
Variety 2: \$4.20
Variety 3: \$4.12
Due to operational considerations, the farmer would like to plan to produce exactly 2 varieties and have at least 4 acres dedicted to each of the two varieties chosen. Additionally, the farmer only wants to make the decision up to 0.1 in acres, i.e., the allocation decision should be made up to the tenth of an acre.
Problem 4a¶
Find the land usage that maximized expected profit.
Problem 4b¶
Suppose that the farmer is concerned about risk, and wants to ensure that the expected 10th quantile of profit is as high as possible. What is the optimal allocation?
We start by defining our variables, constraints, and the objective.
Decision variables
We want to decide how much of each variety to allocate.
Constraints
Must only use at most 2 varieties, and there are at least 4 acres per variety chosen.
Maximize profit
We will also have to incorporate uncertainty, but let’s try just the above first.
best_objective_so_far = 0.0
best_variety_1 = 0.0
best_variety_2 = 0.0
best_variety_3 = 0.0
# Note that range only works for intergers, so I will loop to 200, and consider the values
# as 0.1* the values
# Let’s start with step size 10 just to make sure everything is working
for i in range(0,201,10):
for j in range(0,201,10):
for k in range(0,201,10):
proposed_solution = np.array([i,j,k])/10
## this will determine how many of the varieties are not 0
## If there aren’t exactly 2, we move forward
values_notequal_to_zero = sum(proposed_solution != 0)
if values_notequal_to_zero != 2:
## check that we don’t exceed 20 acres of allocation
if proposed_solution[0] + proposed_solution[1] + proposed_solution[2] > 20:
## check that each is at least 4, if it is above 0!
if proposed_solution[0] > 0:
if proposed_solution[0] < 4:
if proposed_solution[1] > 0:
if proposed_solution[1] < 4:
if proposed_solution[2] > 0:
if proposed_solution[2] < 4:
## Now calculate the number of bushels
bushels_variety_1 = 170*proposed_solution[0]
bushels_variety_2 = 180*proposed_solution[1]
bushels_variety_3 = 185*proposed_solution[2]
## Now calculate the profits
variety_1_profit = 4.40 * bushels_variety_1
variety_2_profit = 4.20 * bushels_variety_2
variety_3_profit = 4.12 * bushels_variety_3
total_profit = variety_1_profit + variety_2_profit + variety_3_profit
## update best if better found
if total_profit > best_objective_so_far:
best_variety_1 = proposed_solution[0]
best_variety_2 = proposed_solution[1]
best_variety_3 = proposed_solution[2]
best_objective_so_far = total_profit
#print(i,j,k)
Let’s output the result and check that it make sense
print(best_variety_1)
print(best_variety_2)
print(best_variety_3)
print(best_objective_so_far)
All looks good! This would provide us with the optimal expected profit. However, when we consider the yield uncertainty, we might be a different solution. To add that consideration we need to use simulation to simulate what the yeild will be per allocated acre.
For the solution we found, this would be code for incorporating the uncertainty:
profit_simulations = []
n_simulations = 100
for simulation_index in range(n_simulations):
bushels_variety_1 = np.random.uniform(1-0.15,1+0.15)*170*best_variety_1
bushels_variety_2 = np.random.uniform(1-0.17,1+0.17)*180*best_variety_2
bushels_variety_3 = np.random.uniform(1-0.2,1+0.2)*185*best_variety_3
## Now calculate the profits
variety_1_profit = 4.40 * bushels_variety_1
variety_2_profit = 4.20 * bushels_variety_2
variety_3_profit = 4.12 * bushels_variety_3
simulated_profit = variety_1_profit + variety_2_profit + variety_3_profit
profit_simulations.append(simulated_profit)
## Here is a histogram
plt.hist(profit_simulations)
## And the average profit
np.mean(profit_simulations)
15121.497171072582
And so, if we want to incorporate the uncertainty, we can embed this simulation within our optimization code.
best_objective_so_far = 0.0
best_variety_1 = 0.0
best_variety_2 = 0.0
best_variety_3 = 0.0
n_simulations = 1000
# Note that range only works for intergers, so I will loop to 200, and consider the values
# as 0.1* the values
# Let’s start with step size 10 just to make sure everything is working
for i in range(0,201,10):
for j in range(0,201,10):
for k in range(0,201,10):
proposed_solution = np.array([i,j,k])/10
## this will determine how many of the varieties are not 0
## If there aren’t exactly 2, we move forward
values_notequal_to_zero = sum(proposed_solution != 0)
if values_notequal_to_zero != 2:
## check that we don’t exceed 20 acres of allocation
if proposed_solution[0] + proposed_solution[1] + proposed_solution[2] > 20:
## check that each is at least 4, if it is above 0!
if proposed_solution[0] > 0:
if proposed_solution[0] < 4:
if proposed_solution[1] > 0:
if proposed_solution[1] < 4:
if proposed_solution[2] > 0:
if proposed_solution[2] < 4:
## Now calculate the number of bushels
profit_simulations = []
for simulation_index in range(n_simulations):
bushels_variety_1 = np.random.uniform(1-0.15,1+0.15)*170*proposed_solution[0]
bushels_variety_2 = np.random.uniform(1-0.17,1+0.17)*180*proposed_solution[1]
bushels_variety_3 = np.random.uniform(1-0.20,1+0.20)*185*proposed_solution[2]
## Now calculate the profits
variety_1_profit = 4.40 * bushels_variety_1
variety_2_profit = 4.20 * bushels_variety_2
variety_3_profit = 4.12 * bushels_variety_3
simulated_profit = variety_1_profit + variety_2_profit + variety_3_profit
profit_simulations.append(simulated_profit)
## update best if better found
if np.mean(profit_simulations) > best_objective_so_far:
best_variety_1 = proposed_solution[0]
best_variety_2 = proposed_solution[1]
best_variety_3 = proposed_solution[2]
best_objective_so_far = np.mean(profit_simulations)
#print(i,j,k)
print(best_variety_1)
print(best_variety_2)
print(best_variety_3)
print(best_objective_so_far)
15316.221589912198
Indeed we don’t want to change! But, now let’s look at doing the quantile.
best_objective_so_far = 0.0
best_variety_1 = 0.0
best_variety_2 = 0.0
best_variety_3 = 0.0
n_simulations = 1000
# Note that range only works for intergers, so I will loop to 200, and consider the values
# as 0.1* the values
# Let’s start with step size 10 just to make sure everything is working
for i in range(0,201,10):
for j in range(0,201,10):
for k in range(0,201,10):
proposed_solution = np.array([i,j,k])/10
## this will determine how many of the varieties are not 0
## If there aren’t exactly 2, we move forward
values_notequal_to_zero = sum(proposed_solution != 0)
if values_notequal_to_zero != 2:
## check that we don’t exceed 20 acres of allocation
if proposed_solution[0] + proposed_solution[1] + proposed_solution[2] > 20:
## check that each is at least 4, if it is above 0!
if proposed_solution[0] > 0:
if proposed_solution[0] < 4:
if proposed_solution[1] > 0:
if proposed_solution[1] < 4:
if proposed_solution[2] > 0:
if proposed_solution[2] < 4:
## Now calculate the number of bushels
profit_simulations = []
for simulation_index in range(n_simulations):
bushels_variety_1 = np.random.uniform(1-0.15,1+0.15)*170*proposed_solution[0]
bushels_variety_2 = np.random.uniform(1-0.17,1+0.17)*180*proposed_solution[1]
bushels_variety_3 = np.random.uniform(1-0.20,1+0.20)*185*proposed_solution[2]
## Now calculate the profits
variety_1_profit = 4.40 * bushels_variety_1
variety_2_profit = 4.20 * bushels_variety_2
variety_3_profit = 4.12 * bushels_variety_3
simulated_profit = variety_1_profit + variety_2_profit + variety_3_profit
profit_simulations.append(simulated_profit)
## update best if better found
if np.quantile(profit_simulations,0.1) > best_objective_so_far:
best_variety_1 = proposed_solution[0]
best_variety_2 = proposed_solution[1]
best_variety_3 = proposed_solution[2]
best_objective_so_far = np.quantile(profit_simulations,0.1)
#print(i,j,k)
print(best_variety_1)
print(best_variety_2)
print(best_variety_3)
print(best_objective_so_far)
13778.338671940213
Because of the uncertainty, we want to allocate more toward different varieties in order to ensure we have high-quality solutions more frequently!
程序代写 CS代考 加微信: powcoder QQ: 1823890830 Email: powcoder@163.com