{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Supplying toilet papers to Walmart stores\n", "\n", "Suppose you are the distribution manager overseeing all the Walmart gorcery stores in California. There are $m$ factories producing toilet papers in the region, which need to be supplied to $n$ Walmart stores in the state.\n", "\n", "The geographic locations (latitudes and longitudes) of the factories and stores: $x_{i}^{\\text{factory}}\\in\\mathbb{R}^{2}$ and $x_{j}^{\\text{store}}\\in\\mathbb{R}^{2}$ are known, for all $i=1,...,m$, and $j=1,...,n$. At present, Walmart uses its own trucks for this delivery which incur $C_{ij}$ dollars (for fuel cost, travel time etc.) to transport **per unit mass of toilet paper** from the $i$th factory to the $j$th store. A reasonable model is that $C_{ij}$ equals to the squared Euclidean distance:\n", "\n", "$$C_{ij} = \\|x_{i}^{\\text{factory}} - x_{j}^{\\text{store}}\\|_{2}^{2}, \\quad\\forall\\;(i,j)\\in\\{1,...,m\\}\\times\\{1,...,n\\}.$$\n", "\n", "This defines a cost matrix $C \\equiv [C_{ij}]$ of size $m\\times n$. The $i$th factory has a fixed (normalized) production capacity $\\alpha_{i}>0$. The $j$th store has a fixed (normalized) storage capacity $\\beta_{j}>0$. Here \"normalized\" means that $\\sum_{i}\\alpha_{i}=\\sum_{j}\\beta_{j}=1$. Toilet papers cannot be stored elsewhere, i.e., all the toilet papers produced must be supplied to the stores." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## (a) [12 + 2 + 2 + (1 + 3) = 20 points] Primal problem\n", "\n", "As the distribution manager, your job is to decide **what is the cheapest way to transport the toilet papers from the factories to the stores**. Suppose you decide that $P_{ij}$ amount of toilet paper mass (say in pound or lb) should be transported from the $i$th factory to the $j$th store. Obviously, $P_{ij} \\geq 0$ for all $(i,j)\\in\\{1,...,m\\}\\times\\{1,...,n\\}$. \n", "\n", "(i) **Taking the matrix $P \\equiv [P_{ij}]$ of size $m\\times n$ as the decision variable**, clearly set up the optimization problem to be solved. The input parameters for this optimization problem should be the $m\\times n$ matrix $C$, the $m\\times 1$ vector $\\alpha$, and the $n\\times 1$ vector $\\beta$.\n", "\n", "**Hint:** the Frobenius inner product $\\langle C, P\\rangle = \\text{trace}(C^{\\top}P)$. Also, you may find it helpful to use the column vector of ones $\\mathbf{1}_{m}, \\mathbf{1}_{n}$.\n", "\n", "(ii) Argue **why** the optimization problem above is convex in $P$.\n", "\n", "(iii) What is the dimensionality of this problem, that is, it is a convex optimization problem in **how many variables with how many constraints**?\n", "\n", "(iv) **What type of convex optimization problem** is it? To support your answer, **explain the geometric interpretation of the objective and the constraints**." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Solution for part (a):\n", "\n", "(i) Since the total cost is $\\langle C, P\\rangle = \\sum_{i=1}^{m}\\sum_{j=1}^{n} C_{ij} P_{ij}$ dollars, the optimization problem becomes:\n", "\n", "\\begin{align*}\n", "&\\underset{P\\in\\mathbb{R}^{m\\times n}}{\\text{minimize}}\\quad \\sum_{i=1}^{m}\\sum_{j=1}^{n} C_{ij} P_{ij}\\\\\n", "&\\text{subject to} \\quad P_{ij} \\geq 0 \\quad\\forall\\:(i,j)\\in\\{1,...,m\\}\\times\\{1,...,n\\},\\\\\n", "& \\qquad\\qquad \\sum_{j=1}^{n}P_{ij} = \\alpha_{i} \\quad\\forall\\:i\\in\\{1,...,m\\},\\\\\n", "& \\qquad\\qquad \\sum_{i=1}^{n}P_{ij} = \\beta_{j} \\quad\\forall\\:j\\in\\{1,...,n\\},\n", "\\end{align*}\n", "or in matrix form:\n", "\\begin{align*}\n", "&\\underset{P\\in\\mathbb{R}^{m\\times n}}{\\text{minimize}}\\quad \\langle C, P\\rangle\\\\\n", "&\\text{subject to} \\quad P \\geq 0 \\quad\\text{(elementwise)},\\\\\n", "& \\qquad\\qquad P\\mathbf{1}_{n} = \\alpha,\\\\\n", "& \\qquad\\qquad P^{\\top}\\mathbf{1}_{m} = \\beta.\n", "\\end{align*}\n", "\n", "(ii) The objective is linear, and hence convex in $P$. All the constraints are also linear in matrix variable $P$. The set $\\mathbb{R}^{m\\times n}$ is affine. Therefore, the optimization problem in part (i) is convex.\n", "\n", "(iii) Since $P$ has $mn$ elements, this is an optimization problem in $mn$ variables with $mn$ inequality constraints, and $m+n$ equality constraints. So all in all, an optimization problem in $mn$ variables with $m+n+mn$ constraints.\n", "\n", "(iv) The problem derived in part (i) is an **LP**. \n", "\n", "This is because the objective is linear in $P$, and the constraint set is an intersection of $mn$ halfspaces with $m+n$ hyperplanes, i.e., a polyhedron." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## (b) [20 points] Numerical solution\n", "\n", "Fix $m=20, n=167$. Write a code in MATLAB/Python/Julia to load the input data from CANVAS Files section: HW problems and solutions: alpha.txt, beta.txt, x\\_stores.txt, x\\_factories.txt, and use cvx/cvxpy/Convex.jl in the same code to solve the optimization problem in part (a). **Report the numerically computed optimal value (minimized cost)** and **submit your code**. \n", "\n", "**Remark:** It is recommended (but not required) that in your code, you also check if your optimal solution obtained from using cvx/cvxpy/Convex.jl matches with linprog in MATLAB, or with scipy.optimize.linprog in Python." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Solution for part (b):\n", "\n", "For numerical computation, we need the recast our optimization problem derived in part (a) in **standard vector LP form**. To do so, let $c:=\\operatorname{vec}(C)$, $p:=\\operatorname{vec}(P)$, and\n", "$$A := \\begin{pmatrix}\n", "\\mathbf{1}_{n}^{\\top}\\otimes I_{m}\\\\\n", "I_{n} \\otimes \\mathbf{1}_{m}^{\\top}\n", "\\end{pmatrix}\\in\\mathbb{R}^{(m+n)\\times mn}, \\qquad b:= \\begin{pmatrix}\n", "\\alpha\\\\\n", "\\beta\n", "\\end{pmatrix}\\in\\mathbb{R}^{m+n},\n", "$$\n", "where $\\operatorname{vec}(\\cdot)$ denotes the column vectorization, and $\\otimes$ denotes the Kronecker product. Then we can rewrite the LP in part (a) as\n", "\n", "\\begin{align*}\n", "&\\underset{p\\in\\mathbb{R}^{mn}}{\\text{minimize}}\\quad \\langle c, p\\rangle\\\\\n", "&\\text{subject to} \\quad p \\geq 0 \\quad\\text{(elementwise)},\\\\\n", "& \\qquad\\qquad Ap = b,\n", "\\end{align*}\n", "where $\\langle c, p\\rangle = c^{\\top}p$ (the vector inner product).\n", "\n", "The numerically computed minimum cost is $$p^{*}_{\\text{cvx}} = 0.057986948821812.$$ The same via MATLAB linprog is $$p^{*}_{\\text{MATLAB}} = 0.057986941579042.$$ \n", "\n", "These numerical values were obtained by running the code $\\texttt{AM229F20HW6.m}$ in MATLAB, posted in CANVAS File section, inside the folder \"HW problems and solutions\"." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## (c) [(10 + 10 + 10) + ( 8 + 2 ) = 40 points] Duality\n", "\n", "(i) Starting from the convex optimization problem in matrix variable $P$ deduced in part (a), derive the Lagrangian $L$, the Lagrange dual function $g$, and the dual optimization problem.\n", "\n", "**Hint:** To handle the inequality constraint $P \\geq 0$ (elementwise), introduce a matricial Lagrange multiplier $\\Lambda$ of size $m\\times n$.\n", "\n", "(ii) Numerically solve the dual problem using cvx/cvxpy/Convex.jl **in the same code in part (b)**. So your submitted code should solve both the primal and dual problems separately. **Report the numerically computed optimal value of the dual problem, and explain what you observe**." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Solution for part (c):\n", "\n", "(i) Following the hint provided, the **Lagrangian** is\n", "$$L(P,\\Lambda,\\nu_{1},\\nu_{2}) = \\langle C, P\\rangle + \\langle\\Lambda, -P\\rangle + \\langle \\nu_{1}, P\\mathbf{1}_{n}-\\alpha\\rangle + \\langle \\nu_{2},P^{\\top}\\mathbf{1}_{m}-\\beta\\rangle.$$\n", "\n", "The **Lagrange dual function**, by definition, is\n", "\n", "$$g(\\Lambda,\\nu_{1},\\nu_{2}) = \\underset{P\\in\\mathbb{R}^{m\\times n}}{\\sup}L(P,\\Lambda,\\nu_{1},\\nu_{2}).$$\n", "\n", "Using that trace is invariant under cyclic permutation, we write $\\nu_{1}^{\\top}P\\mathbf{1}_{n} = \\operatorname{trace}(\\mathbf{1}_{n}\\nu_{1}^{\\top}P)$, and that $\\nu_{2}^{\\top}P^{\\top}\\mathbf{1}_{m} = \\operatorname{trace}(\\nu_{2}^{\\top}P^{\\top}\\mathbf{1}_{m}) = \\operatorname{trace}(\\mathbf{1}_{m}^{\\top}P\\nu_{2}) = \\operatorname{trace}(\\nu_{2}\\mathbf{1}_{m}^{\\top}P)$, where the last but one equality used that trace is also invariant under matrix transpose.\n", "\n", "Therefore, we obtain\n", "\\begin{align*}\n", "g(\\Lambda,\\nu_{1},\\nu_{2}) &= \\underset{P\\in\\mathbb{R}^{m\\times n}}{\\sup} \\operatorname{trace}\\left((C^{\\top}-\\Lambda^{\\top} + \\mathbf{1}_{n}\\nu_{1}^{\\top} + \\nu_{2}\\mathbf{1}_{m}^{\\top})P\\right) - \\nu_{1}^{\\top}\\alpha - \\nu_{2}^{\\top}\\beta\\\\\n", "&= \\begin{cases}\n", "- \\nu_{1}^{\\top}\\alpha - \\nu_{2}^{\\top}\\beta \\quad \\text{if} \\quad C-\\Lambda + \\nu_{1}\\mathbf{1}_{n}^{\\top} + \\mathbf{1}_{m}\\nu_{2}^{\\top} = \\mathbf{0}_{m\\times n},\\\\\n", "+\\infty \\quad \\text{otherwise}.\\end{cases}\n", "\\end{align*}\n", "\n", "Consequently, the **dual problem** is\n", "\\begin{align*}\n", "&\\underset{\\Lambda\\in\\mathbb{R}^{m\\times n},\\nu_1\\in\\mathbb{R}^{m},\\nu_2\\in\\mathbb{R}^{n}}{\\max} -\\nu_{1}^{\\top}\\alpha - \\nu_{2}^{\\top}\\beta\\\\\n", "&\\quad\\text{such that}\\quad C-\\Lambda + \\nu_{1}\\mathbf{1}_{n}^{\\top} + \\mathbf{1}_{m}\\nu_{2}^{\\top} = \\mathbf{0}_{m\\times n},\\\\\n", "& \\qquad\\qquad\\qquad \\Lambda \\geq 0 \\quad\\text{(elementwise)}.\n", "\\end{align*}\n", "We notice that the dual problem is also **LP**, as expected.\n", "\n", "(ii) The same code posted in CANVAS (as in part(b)) also solves the dual LP via cvx. The computed solution was\n", "\n", "$$d^{*}_{\\text{cvx}} = 0.057986941767268.$$\n", "\n", "We see that $d^{*}$ matches with $p^{*}$ upto 8 significant digits. This is expected since theoretically LPs have zero duality gap, i.e., $d^{*}=p^{*}$." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## (d) [20 points] To contract or not\n", "\n", "A Silicon Valley startup **Brand New Truck Company (BNTC)** has approached you, the Walmart distribution manager in California, claiming that they can lower the cost of transport for you if you give them the contract to transport the toilet papers from the $m$ factories to the $n$ Walmart stores. **Their business model is that they actually do not charge for the transportation process, but only charge variable per unit mass loading cost $\\eta_{1}\\in\\mathbb{R}^{m}$, and variable per unit mass unloading cost $\\eta_{2}\\in\\mathbb{R}^{n}$**. They also guarantee you that \n", "\n", "$$\\left(\\eta_{1}\\right)_{i} + \\left(\\eta_{2}\\right)_{j} \\leq C_{ij} \\quad\\forall\\;(i,j)\\in\\{1,...,m\\}\\times\\{1,...,n\\},$$\n", "\n", "and hence, BNTC argues that Walmart cannot loose money by awarding them the contract for toilet paper transport. \n", "\n", "To maximize revenue, BNTC determines $(\\eta_{1},\\eta_{2})$ by solving\n", "\n", "\\begin{align*}\n", "\\underset{(\\eta_{1},\\eta_{2})}{\\text{maximize}}\\quad &\\eta_{1}^{\\top}\\alpha + \\eta_{2}^{\\top}\\beta\\\\\n", "\\qquad\\qquad\\text{s.t.}\\quad &\\eta_{1}\\mathbf{1}_{n}^{\\top} + \\mathbf{1}_{m}\\eta_{2}^{\\top} \\leq C \\quad\\text{(elementwise).}\n", "\\end{align*}\n", "\n", "Is contracting the transport to BNTC profitable for Walmart? **Use the duality results in part (c) to mathematically explain your answer**." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Solution for part (c):\n", "\n", "**No, by awarding the contract to BNTC, Walmart cannot do any better than what they are doing now.** \n", "\n", "In fact, by awarding the contract, Walmart will incur the exact same cost as they are incurring now. To see why, let us start from the dual LP in part (b), and let $\\eta_1 := -\\nu_1$, $\\eta_2 := -\\nu_2$. Then the dual objective reduces to the BNTC revenue maximization objective. We can also combine the two constraints in the dual LP as the single constraint:\n", "\n", "$$C - \\eta_{1}\\mathbf{1}_{n}^{\\top} - \\mathbf{1}_{m}\\eta_{2}^{\\top} = \\Lambda \\geq 0 \\:\\text{(elementwise)}\\:\\Rightarrow \\: C \\geq \\eta_{1}\\mathbf{1}_{n}^{\\top} + \\mathbf{1}_{m}\\eta_{2}^{\\top},$$\n", "\n", "where we used $\\eta_1 := -\\nu_1$, $\\eta_2 := -\\nu_2$. But this combined inequality is exactly the fesability constraint in the revenue maximization problem of BNTC.\n", "\n", "Therefore, the revenue maximization problem of BNTC is exactly the dual LP in part (b), and by **strong duality for LPs**, its optimal value must be equal to the primal minimum: the transportation cost Walmart is currently incurring using its own trucks." ] } ], "metadata": { "kernelspec": { "display_name": "Python 2", "language": "python", "name": "python2" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 2 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython2", "version": "2.7.16" } }, "nbformat": 4, "nbformat_minor": 2 }