From Quantecon, we'll start with a simple production problem, very similar to our baker problem.
We have two products we produce/sell, $x_1$ and $x_2$, which create respectively 3 and 4 of profits per unit sold. We have two inputs and two associated resource constraints.
We formulate the following problem:
$$\max_{x_1, x_2} 3x_1 + 4x_2$$$$\text{subject to}$$$$2x_1 + 5x_2 \leq 30\\ 4x_1 + 2x_2\leq 20\\ x_1 \geq0, x_2 \geq 0$$As we will see, there are many ways of formulating a linear program, but for the purposes of entering a linear program into a computer and solving numerically, it is important to have a certain standard form.
We state all objective functions as minimization problems. If we initially have a maximization problem, just take the negative.
write all resource constraints as equality constraints in the form:
Let's look at our production problem and see how we can transform it into standard form.
But, for the sake of argument, let's say that we could loan 20 of $x_1$ in our production, so that we have a constraint:
$$x_1\geq -20$$How could we transform this into standard form? We would replace $x_1$ with the difference of two new variables, which we call $x_1^-$ and $x_1^+$: $x_1^+ - x_1^- = x_1$ and a slack variable, $s_3$ so that we get a new constraint:
$$x_1^+-x_1^- + s_3 = -20$$$$s_3, x_1^+, x_1^- \geq 0$$Thus if our optimal solution would require loaning such that $x_1=-15$, then $x_1^- = 20$, $x_1^+ = 0$ and $s_3=5$
Now that we have our problem in standard form, we can do our computation
import numpy as np
from scipy.optimize import linprog
import matplotlib.pyplot as plt
# Construct parameters
c_ex1 = np.array([-3, -4])
# Inequality constraints
A_ex1 = np.array([[2, 5],
[4, 2]])
b_ex1 = np.array([30,20])
# Solve the problem
# we put a negative sign on the objective as linprog does minimization
res_ex1 = linprog(c_ex1, A_ub=A_ex1, b_ub=b_ex1, method='highs')
res_ex1
message: Optimization terminated successfully. (HiGHS Status 7: Optimal)
success: True
status: 0
fun: -27.5
x: [ 2.500e+00 5.000e+00]
nit: 2
lower: residual: [ 2.500e+00 5.000e+00]
marginals: [ 0.000e+00 0.000e+00]
upper: residual: [ inf inf]
marginals: [ 0.000e+00 0.000e+00]
eqlin: residual: []
marginals: []
ineqlin: residual: [ 0.000e+00 0.000e+00]
marginals: [-6.250e-01 -4.375e-01]
mip_node_count: 0
mip_dual_bound: 0.0
mip_gap: 0.0
We could also reformulate our problem in the dual
$$\min_{p_1, p_2} 30 p_1 + 20 p_2$$Subject to
$$2p_1 + 4p_2\geq3\\ 5p_1 + 2p_2\geq4\\ p_1, p_2\geq0$$To solve this we can just transform some of our constraints to be negative so that the inequalities are in the right direction and then reuse the arrays and matrices we defined above:
res_ex1_dual = linprog(b_ex1, A_ub=-A_ex1.T, b_ub=c_ex1, method='highs')
res_ex1_dual
message: Optimization terminated successfully. (HiGHS Status 7: Optimal)
success: True
status: 0
fun: 27.5
x: [ 6.250e-01 4.375e-01]
nit: 2
lower: residual: [ 6.250e-01 4.375e-01]
marginals: [ 0.000e+00 0.000e+00]
upper: residual: [ inf inf]
marginals: [ 0.000e+00 0.000e+00]
eqlin: residual: []
marginals: []
ineqlin: residual: [ 0.000e+00 0.000e+00]
marginals: [-2.500e+00 -5.000e+00]
mip_node_count: 0
mip_dual_bound: 0.0
mip_gap: 0.0
We'll consider the second problem from Quantecon: The investment problem. Here we'll formulate a portfolio investment problem as a linear program.
The setup is as follows:
100,000 to invest
3 year horizon
Choice of three investments:
Annuity ($x_1 \geq 0$): place the same amount every year and at the end of the third year receive payoff of 130\% the total amount put in
Bank account ($x_2, x_3, x_4 \geq -20,000$): you can place any amount at the beginning of each year (1,2,3) and receive 6\% interest at the end of the year. You can also borrow up to 20,000 each year at 6\% interest.
Corporate bond ($0 \leq x_5 \leq 50,000$): at the beginning of the second year you can place upto 50,000 in a corporate bond that pays our 130\% of your capital at the end of the third year.
Considering the timing and investment constraints we can write:
This will be our objective function
We can write out our LP as
$$\max 1.3*3x_1 + 1.06x_4 + 1.3x_5$$$$\text{Subject to}$$$$x_1 + x_2 = 100000\\ x_1 + - 1.06x_2 + x_3 + x_5 = 0\\ x_1 - 1.06x_3 + x_4 = 0\\ x_2 \geq -20,000\\ x_3 \geq -20,000\\ x_4 \geq -20,000\\ x_5 \leq 50,000\\ x_1, x_5 \geq 0$$After transforming to standard form we can calculate out our optimal solution
# Objective function parameters
c_ex2 = np.array([1.30*3, 0, 0, 1.06, 1.30])
# Inequality constraints (x1...x5)
A_ex2 = np.array([[1, 1, 0, 0, 0],
[1, -1.06, 1, 0, 1],
[1, 0, -1.06, 1, 0]])
b_ex2 = np.array([100000, 0, 0])
#With the linear program solver, we can enter our wished-for bounds on our variables:
bounds_ex2 = [( 0, None),
(-20000, None),
(-20000, None),
(-20000, None),
( 0, 50000)]
# Solve the problem
res_ex2 = linprog(-c_ex2, A_eq=A_ex2, b_eq=b_ex2,
bounds=bounds_ex2, method='highs')
res_ex2
message: Optimization terminated successfully. (HiGHS Status 7: Optimal)
success: True
status: 0
fun: -141018.24349792697
x: [ 2.493e+04 7.507e+04 4.649e+03 -2.000e+04 5.000e+04]
nit: 0
lower: residual: [ 2.493e+04 9.507e+04 2.465e+04 0.000e+00
5.000e+04]
marginals: [ 0.000e+00 0.000e+00 0.000e+00 1.650e-01
0.000e+00]
upper: residual: [ inf inf inf inf
0.000e+00]
marginals: [ 0.000e+00 0.000e+00 0.000e+00 0.000e+00
-1.470e-03]
eqlin: residual: [ 0.000e+00 0.000e+00 0.000e+00]
marginals: [-1.376e+00 -1.299e+00 -1.225e+00]
ineqlin: residual: []
marginals: []
mip_node_count: 0
mip_dual_bound: 0.0
mip_gap: 0.0
According to our results if we want to maximize our profit we should:
Place 24,927 in the annuity in each period
Place 75,072 in the bank the first period (just the rest amount after the annuity)
Place 4,648 in the bank the second period, while placing the full 50,000 possible in the corporate bond
In the third period, borrow the full 20000 that is possible in order to finance the final installment of the annuity.
# Objective function parameters
c_ex2_dual = np.array([100000, 0, 0, -20000, -20000, -20000, 50000])
# Equality constraints
A_eq_ex2_dual = np.array([[1, -1.06, 0, 1, 0, 0, 0],
[0, 1, -1.06, 0, 1, 0, 0],
[0, 0, 1, 0, 0, 1, 0]])
b_eq_ex2_dual = np.array([0, 0, 1.06])
# Inequality constraints
A_ub_ex2_dual = - np.array([[1, 1, 1, 0, 0, 0, 0],
[0, 1, 0, 0, 0, 0, 1]])
b_ub_ex2_dual = - np.array([1.30*3, 1.30])
# Bounds on decision variables
bounds_ex2_dual = [(None, None),
(None, None),
(None, None),
(None, 0),
(None, 0),
(None, 0),
( 0, None)]
# Solve the dual problem
res_ex2_dual = linprog(c_ex2_dual, A_eq=A_eq_ex2_dual, b_eq=b_eq_ex2_dual,
A_ub=A_ub_ex2_dual, b_ub=b_ub_ex2_dual, bounds=bounds_ex2_dual,
method='revised simplex')
res_ex2_dual
con: array([-2.22044605e-16, 0.00000000e+00, 0.00000000e+00])
fun: 141018.2434979269
message: 'Optimization terminated successfully.'
nit: 4
slack: array([0., 0.])
status: 0
success: True
x: array([ 1.37644176, 1.29852997, 1.22502827, 0. , 0. ,
-0.16502827, 0.00147003])
We could give the solutions a direct interpretation: something like the value (return) of loosening each constraint by a dollar. For example, $p_1 = 1.37$ says that if we had one more dollar of initial capital, we could make 1.37 in return at the end of the third year. We could interpret this as a version of return on capital.
Our last example, also from Quantecon, will be what at first appears to be a logistics problem, but actually has many applications in economics.
The full description of the problem can be found at the QuantEcon lecture but here is a short summary:
We have m factories (indexed by i) that produce goods that must be sent to n (indexed by j) locations
A planner wants to minimize the total transport costs subject to the constraints:
Thus we could write out our LP problem as:
$$\min_{x_{ij}} \sum_{i=1}^{m} \sum_{j=1}^{n} c_{ij} x_{ij}$$$$\text{subject to}$$$$\sum_{j=1}^{n} x_{ij} = p_i, \quad i = 1,2,...,m\\ \sum_{i=1}^{m} x_{ij} = q_j,\quad j = 1,2,...,n\\ x_{ij} \geq0$$The Quantecon lecture gives a numerical example with $m=3$ factories and $j=5$ locations:
m = 3
n = 5
p = np.array([50, 100, 150])
q = np.array([25, 115, 60, 30, 70])
C = np.array([[10, 15, 20, 20, 40],
[20, 40, 15, 30, 30],
[30, 35, 40, 55, 25]])
Thus the i'th row and j'th column of C represents the cost of shipping from factory i to location j. Each p represents the capacity of each of the three factories. Each of the q represents demand at each location, 5.
A somewhat technical problem is that our solver expects a vector in the objective function, but here we have a matrix of costs. It is pretty straightforward to reshape our C vector as a column vector using the reshape command
C_vec = C.reshape((m*n, 1), order='F')
C_vec
array([[10],
[20],
[30],
[15],
[40],
[35],
[20],
[15],
[40],
[20],
[30],
[55],
[40],
[30],
[25]])
Now we want to create one A matrix that includes our constraints. We split this into two parts.
First we want to create the constraints for the m=3 factories
We start by creating an identity matrix of size 3 x 3 (m=3), as well as a vector of 1s with n=5 entries
identM = np.identity(m)
print(identM)
onesN = np.ones((1, n))
print(onesN)
[[1. 0. 0.] [0. 1. 0.] [0. 0. 1.]] [[1. 1. 1. 1. 1.]]
Then we use kroniker product to create a 3x15 matrix. When we take a kroniker product of two matrices, A and B, written $A \bigotimes B$, then we take each element of A times the entire matrix of B.
So below we multiply each of the 1s in the row vector with the whole 3x3 identity matrix, combining the result into one 3x15 matrix
A1 = np.kron(onesN, identM)
A1
array([[1., 0., 0., 1., 0., 0., 1., 0., 0., 1., 0., 0., 1., 0., 0.],
[0., 1., 0., 0., 1., 0., 0., 1., 0., 0., 1., 0., 0., 1., 0.],
[0., 0., 1., 0., 0., 1., 0., 0., 1., 0., 0., 1., 0., 0., 1.]])
This can then represent the left part of the constraint $\sum_{j=1}^{n} x_{ij} = p_i$, where each row represents one of the 3 factories, where we are summing over the 5 locations, represented by a one. Basically we can imagine matrix multiplying this by a long 15 entry vector:
$[x_{11}, x_{21}, x_{31}, x_{12}, x_{22}, x_{32}, x_{13}, ...x_{35}]$
to get the left-hand side of our first constraint: $\sum_{j=1}^n x_{ij}$
Then we construct the second part of the constraints:
A2 = np.kron(np.identity(n), np.ones((1, m)))
A2
array([[1., 1., 1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
[0., 0., 0., 1., 1., 1., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
[0., 0., 0., 0., 0., 0., 1., 1., 1., 0., 0., 0., 0., 0., 0.],
[0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 1., 1., 0., 0., 0.],
[0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 1., 1.]])
Again matrix multiplying by the vector of $x_{ij}$ terms will give the left-hand side of the second constraint: $\sum_{i=1}^m x_{ij}$
We'll then stack the two A matrices on top of each other as well combining the two vectors of limits (p, q)
A = np.vstack([A1, A2])
A
array([[1., 0., 0., 1., 0., 0., 1., 0., 0., 1., 0., 0., 1., 0., 0.],
[0., 1., 0., 0., 1., 0., 0., 1., 0., 0., 1., 0., 0., 1., 0.],
[0., 0., 1., 0., 0., 1., 0., 0., 1., 0., 0., 1., 0., 0., 1.],
[1., 1., 1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
[0., 0., 0., 1., 1., 1., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
[0., 0., 0., 0., 0., 0., 1., 1., 1., 0., 0., 0., 0., 0., 0.],
[0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 1., 1., 0., 0., 0.],
[0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 1., 1.]])
b = np.hstack([p, q])
b
array([ 50, 100, 150, 25, 115, 60, 30, 70])
Just to confirm our math, we could make up a vector of x's:
x = np.random.uniform(0,10,15)
x
array([5.32819219, 4.79879356, 5.7128987 , 0.88632538, 2.93967514,
3.28784464, 0.55935787, 8.91241387, 9.37797808, 8.48694383,
7.47636113, 1.66244424, 9.30249403, 4.13165592, 2.15410871])
np.dot(A,x.T)
array([24.56331331, 28.25889963, 22.19527437, 15.83988445, 7.11384516,
18.84974982, 17.62574921, 15.58825867])
Which we can compare to the constraints, b
b
array([ 50, 100, 150, 25, 115, 60, 30, 70])
If these were inequality constraints, then we can see that this would be a feasible solution.
To find an optimal solution we can then solve using our linprog command:
res = linprog(C_vec, A_eq=A, b_eq=b, method='highs')
res
con: array([0., 0., 0., 0., 0., 0., 0., 0.])
fun: 7225.0
message: 'Optimization terminated successfully.'
nit: 12
slack: array([], dtype=float64)
status: 0
success: True
x: array([15., 10., 0., 35., 0., 80., 0., 60., 0., 0., 30., 0., 0.,
0., 70.])
We can look closer at the x variable
res.x.reshape((m,n), order='F')
array([[15., 35., 0., 0., 0.],
[10., 0., 60., 30., 0.],
[ 0., 80., 0., 0., 70.]])
We can also double-check that this satisfies our constraints
We can interpret this as saying that from the first factory we want to send 15 to location 1, 35 to location two and nothing to the 3-5th locations, and so on.
x= res.x
np.dot(A,x.T)
array([ 50., 100., 150., 25., 115., 60., 30., 70.])
b
array([ 50, 100, 150, 25, 115, 60, 30, 70])
We can also calculate out the cost:
np.dot(C_vec.T, x)
array([7225.])
Doing the above optimisation, we get a warning message saying that A "is not full row rank". Recall from our linear algebra section that this means that the rows are not all linearly independent. This means that at least one of the rows is redundant, and it also means that we don't necessarily have a unique solution.
np.linalg.matrix_rank(A)
7
The rank of the matrix is 7 (we have 7 linearly independent rows) but 8 total rows.
We could try running the optimisation, eliminating one of the constraints/rows:
res2 = linprog(C_vec, A_eq=A[:-1], b_eq=b[:-1], method='highs')
res2
con: array([0., 0., 0., 0., 0., 0., 0.])
fun: 7225.0
message: 'Optimization terminated successfully.'
nit: 13
slack: array([], dtype=float64)
status: 0
success: True
x: array([ 0., 25., 0., 35., 0., 80., 0., 60., 0., 15., 15., 0., 0.,
0., 70.])
Notice now that we get a different solution, but if we calculate out the costs:
x2= res2.x
np.dot(C_vec.T, x2)
array([7225.])
We get the same (minimized) costs.
We can also formulate a dual formulation of our transport problem:
$$\max_{u_i, v_j} \sum_{i=1}^{m}p_i u_i + \sum_{j=1}^n q_j v_j$$$$\text{subject to}$$$$u_i + v_j \leq c_{ij}$$Here we can again interpret the $u_i$ and $v_j$ values as shadow prices of the production and demand constraints in our primal problem.
We can solve the dual problem:
res_dual = linprog(-b, A_ub=A.T, b_ub=C_vec,
bounds=[(None, None)]*(m+n), method='Revised simplex')
res_dual
con: array([], dtype=float64)
fun: -7225.0
message: 'Optimization terminated successfully.'
nit: 7
slack: array([ 0., 0., 0., 0., 15., 0., 15., 0., 15., 0., 0., 15., 35.,
15., 0.])
status: 0
success: True
x: array([ 5., 15., 25., 5., 10., 0., 15., 0.])
Consider $u_1=5$: since our initial problem was one of minimizing costs, this shadow cost then represents the extra (transport) cost of shipping one more product from factory 1.
We can also consider $v_1=5$ (the 4th entry in our output array). This represents the marginal cost of shipping to the 1st location.
The constraints of the dual problem now makes sense. For example: $u_1 + v_1 \leq c_{11}$: the sum of shipping-from and shipping-to costs equals the shipping cost from a factory to a location.
There is value of using a general linear programming solver and manually configuring our optimal transport problem. But there exists a specialized Optimal Transport software for Python, which can make setting up and solving such a problem easier. You can explore this if you are interested:
As we introduced optimal transport problems with, this is a frameworks that is relevant to more than just logistics problems. A thorough (and technical) review of optimal transport proplems can be found here:
https://arxiv.org/pdf/1803.00567.pdf
Basically, whenever we have a problem that we can formulate as moving a commodity, person, or even bit of information from one place to another, then we can make use of the optimal transport frameworks.
We can, for example apply it to problems such as:
In this exercise we imagine with have a small academic department with 10 academic employees who each work 120 hours in a month doing three tasks:
The 10 employees have varying competencies which we rank between 1 and 5 which reflects the quantifiable productivity/value of work.
| Employee | Teaching | Research | Admin |
|---|---|---|---|
| 1 | 5 | 5 | 5 |
| 2 | 1 | 1 | 1 |
| 3 | 3 | 5 | 4 |
| 4 | 2 | 1 | 4 |
| 5 | 5 | 3 | 2 |
| 6 | 5 | 4 | 2 |
| 7 | 3 | 4 | 1 |
| 8 | 2 | 1 | 1 |
| 9 | 1 | 5 | 2 |
| 10 | 4 | 2 | 1 |
All employees are limited to work no more than 120 hours in a month and the total hours of each category need to be filled.
Formulate an optimal transport problem to divide the hours among the employees while maximizing the total value of work.
The Ministry of Education wants to put emphasis on teaching. The department therefor decides to increase the value on teaching by a factor of 1.5 (So a value on teaching of 1 becomes 1.5, 2 becomes 3, so on). How will this change the optimal assignment of teaching vs research and admin. Consider especially how those that are good or bad across the board.
Now consider some extra constraints. The academic staff have written into their contracts that they are guaranteed 30 hours a month of research time. How does this affect the total. What are the shadow price associated with this constrains
We could formulate our problem as follows
$$\max_{x_{ij}} \sum_{i=1}^{10} \sum_{j=1}^3 c_{ij} x_{ij}$$$$\text{subject to}\\ \sum_{j=1}^3 x_{ij} \leq 120, \quad i=1,2...10 \\ \sum_{i=1}^{10} x_{i1} = 500\\ \sum_{i=1}^{10} x_{i2} = 500\\ \sum_{i=1}^{10} x_{i3} = 200\\ x_{ij} \geq 0, \quad i=1,2....10, \quad j=1,2,3$$We can then create our matrices and vectors:
C = np.array([[5,5,5], [1,1,1], [3,5,4], [2,3,4], [5,3,2], [5,4,2], [3,4,2], [2,2,3], [3,5,3], [4,2,3]])
p= np.repeat(120, 10)
q = np.array([500, 500, 200])
C_vec = C.reshape((30, 1), order='F')
C_vec
array([[5],
[1],
[3],
[2],
[5],
[5],
[3],
[2],
[3],
[4],
[5],
[1],
[5],
[3],
[3],
[4],
[4],
[2],
[5],
[2],
[5],
[1],
[4],
[4],
[2],
[2],
[2],
[3],
[3],
[3]])
Then we create our A vector
A1 = np.kron(np.ones((1, 3)), np.identity(10))
A2 = np.kron(np.identity(3), np.ones((1, 10)))
A = np.vstack([A1, A2])
# Construct vector b
b = np.hstack([p, q])
A
array([[1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 0., 0., 0., 0., 0.,
0., 0., 0., 0., 1., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
[0., 1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 0., 0., 0., 0.,
0., 0., 0., 0., 0., 1., 0., 0., 0., 0., 0., 0., 0., 0.],
[0., 0., 1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 0., 0., 0.,
0., 0., 0., 0., 0., 0., 1., 0., 0., 0., 0., 0., 0., 0.],
[0., 0., 0., 1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 0., 0.,
0., 0., 0., 0., 0., 0., 0., 1., 0., 0., 0., 0., 0., 0.],
[0., 0., 0., 0., 1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 0.,
0., 0., 0., 0., 0., 0., 0., 0., 1., 0., 0., 0., 0., 0.],
[0., 0., 0., 0., 0., 1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1.,
0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 0., 0., 0., 0.],
[0., 0., 0., 0., 0., 0., 1., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 0., 0., 0.],
[0., 0., 0., 0., 0., 0., 0., 1., 0., 0., 0., 0., 0., 0., 0., 0.,
0., 1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 0., 0.],
[0., 0., 0., 0., 0., 0., 0., 0., 1., 0., 0., 0., 0., 0., 0., 0.,
0., 0., 1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 0.],
[0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 0., 0., 0., 0., 0., 0.,
0., 0., 0., 1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1.],
[1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 0., 0., 0., 0., 0., 0.,
0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
[0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 1., 1., 1., 1., 1.,
1., 1., 1., 1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
[0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
0., 0., 0., 0., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.]])
b
array([120, 120, 120, 120, 120, 120, 120, 120, 120, 120, 500, 500, 200])
Above, the first 10 rows have 3 "1"'s and thus represent the summation for each employee of the total hours of work divided between the three categories.
The last three rows have 10 "1"'s and represent the summation of the total "demand" for teaching, research and administration of 500, 500 and 200 hours respectively
Then we optimize our problem. Since we are maximizing, C_vec enters negatively:
resAcad = linprog(-C_vec, A_eq=A, b_eq=b, method='Revised simplex')
/var/folders/fs/kr83xdyx2jq_8gww_xfmtxfr0000gn/T/ipykernel_25525/2473598620.py:1: OptimizeWarning: A_eq does not appear to be of full row rank. To improve performance, check the problem formulation for redundant equality constraints. resAcad = linprog(-C_vec, A_eq=A, b_eq=b, method='Revised simplex')
Here again we see that we do not have full rank, thus some constraints are not necessary. We can check the rank of the A matrix:
np.linalg.matrix_rank(A)
12
We have 13 rows and rank 12, so one of the constraints is not strictly necessary
But for now let's consider the (non-unique) optimal solution found
resAcad
con: array([0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.])
fun: -4880.0
message: 'Optimization terminated successfully.'
nit: 15
slack: array([], dtype=float64)
status: 0
success: True
x: array([100., 0., 0., 0., 120., 120., 0., 40., 0., 120., 20.,
120., 120., 0., 0., 0., 120., 0., 120., 0., 0., 0.,
0., 120., 0., 0., 0., 80., 0., 0.])
AcadSchedule = resAcad.x.reshape((10,3), order='F')
AcadSchedule
array([[100., 20., 0.],
[ 0., 120., 0.],
[ 0., 120., 0.],
[ 0., 0., 120.],
[120., 0., 0.],
[120., 0., 0.],
[ 0., 120., 0.],
[ 40., 0., 80.],
[ 0., 120., 0.],
[120., 0., 0.]])
Let's compare this with our productivity score:
C
array([[5, 5, 5],
[1, 1, 1],
[3, 5, 4],
[2, 3, 4],
[5, 3, 2],
[5, 4, 2],
[3, 4, 2],
[2, 2, 3],
[3, 5, 3],
[4, 2, 3]])
A few things to notice:
- Most employees end up being specialized where they are most productive. A few (1,8) split their time between tasks.
- Employees 1 and 8 are especially interesting. Employee 1 is, in absolute terms, the best teacher but spends less time teaching than other less productive teachers. Why? Because 1 is also good at research, which is a scarcer resource. So where they end up being placed isn't always because of their absolute productivity, but rather because of their productivity relative to the scarcity of that productivity.
Now we increase the value of teaching due to a change in policy at the Ministry of Education
C = C.astype(float)
array([[5., 5., 5.],
[1., 1., 1.],
[3., 5., 4.],
[2., 3., 4.],
[5., 3., 2.],
[5., 4., 2.],
[3., 4., 2.],
[2., 2., 3.],
[3., 5., 3.],
[4., 2., 3.]])
C_new = C.copy()
C_new[:,0] = C_new[:,0]*1.5
C_new
C_new_vec = C_new.reshape((30, 1), order='F')
C_new_vec
array([[7.5],
[1.5],
[4.5],
[3. ],
[7.5],
[7.5],
[4.5],
[3. ],
[4.5],
[6. ],
[5. ],
[1. ],
[5. ],
[3. ],
[3. ],
[4. ],
[4. ],
[2. ],
[5. ],
[2. ],
[5. ],
[1. ],
[4. ],
[4. ],
[2. ],
[2. ],
[2. ],
[3. ],
[3. ],
[3. ]])
resAcad_new = linprog(-C_new_vec, A_eq=A, b_eq=b, method='Revised simplex')
/var/folders/fs/kr83xdyx2jq_8gww_xfmtxfr0000gn/T/ipykernel_25525/821962993.py:1: OptimizeWarning: A_eq does not appear to be of full row rank. To improve performance, check the problem formulation for redundant equality constraints. resAcad_new = linprog(-C_new_vec, A_eq=A, b_eq=b, method='Revised simplex')
resAcad_new
AcadSchedule_new = resAcad_new.x.reshape((10,3), order='F')
AcadSchedule_new
array([[120., 0., 0.],
[ 0., 120., 0.],
[ 0., 120., 0.],
[ 0., 20., 100.],
[120., 0., 0.],
[120., 0., 0.],
[ 0., 120., 0.],
[ 20., 0., 100.],
[ 0., 120., 0.],
[120., 0., 0.]])
We can compare this to the original:
AcadSchedule
array([[100., 20., 0.],
[ 0., 120., 0.],
[ 0., 120., 0.],
[ 0., 0., 120.],
[120., 0., 0.],
[120., 0., 0.],
[ 0., 120., 0.],
[ 40., 0., 80.],
[ 0., 120., 0.],
[120., 0., 0.]])
And the original productivity matrix
C
array([[5., 5., 5.],
[1., 1., 1.],
[3., 5., 4.],
[2., 3., 4.],
[5., 3., 2.],
[5., 4., 2.],
[3., 4., 2.],
[2., 2., 3.],
[3., 5., 3.],
[4., 2., 3.]])
We note that a few shifts:
We could reformulate our problem with new constraints
We can adjust our computation by setting bounds on our $x_{ij}$
#acadBounds = [(0, None), (30, None), (0, None)]*10
teach = [(0, None)]*10
research = [(30, None)]*10
admin = [(0, None)]*10
acadBounds = teach + research + admin
acadBounds
[(0, None), (0, None), (0, None), (0, None), (0, None), (0, None), (0, None), (0, None), (0, None), (0, None), (30, None), (30, None), (30, None), (30, None), (30, None), (30, None), (30, None), (30, None), (30, None), (30, None), (0, None), (0, None), (0, None), (0, None), (0, None), (0, None), (0, None), (0, None), (0, None), (0, None)]
resAcad_bound = linprog(-C_vec, A_eq=A, b_eq=b, bounds = acadBounds, method='Revised simplex')
/var/folders/fs/kr83xdyx2jq_8gww_xfmtxfr0000gn/T/ipykernel_25525/533781487.py:1: OptimizeWarning: A_eq does not appear to be of full row rank. To improve performance, check the problem formulation for redundant equality constraints. resAcad_bound = linprog(-C_vec, A_eq=A, b_eq=b, bounds = acadBounds, method='Revised simplex')
AcadSchedule_bound = resAcad_bound.x.reshape((10,3), order='F')
AcadSchedule_bound
array([[ 90., 30., 0.],
[ 70., 30., 20.],
[ 0., 120., 0.],
[ 0., 30., 90.],
[ 90., 30., 0.],
[ 90., 30., 0.],
[ 70., 50., 0.],
[ 0., 30., 90.],
[ 0., 120., 0.],
[ 90., 30., 0.]])