![\"Open](\"https://colab.research.google.com/assets/colab-badge.svg\" "\\"Open"){kind=icon}
![\"Download\"](\"https://img.shields.io/badge/Github-Download-blue.svg\" "\\"Download"){kind=icon}
"
]
},
{
"cell_type": "markdown",
"metadata": {
"colab_type": "text",
"id": "BsHcQn4HT0sV",
"nbpages": {
"level": 1,
"link": "[6.3 Sensitivity Analysis with Adjoint Operators](https://ndcbe.github.io/cbe67701-uncertainty-quantification/06.03-Sensitivity-Analysis-with-Adjoint-Operators.html#6.3-Sensitivity-Analysis-with-Adjoint-Operators)",
"section": "6.3 Sensitivity Analysis with Adjoint Operators"
}
},
"source": [
"# 6.3 Sensitivity Analysis with Adjoint Operators\n",
"\n",
"Following [McClarren (2018)](https://link-springer-com.proxy.library.nd.edu/book/10.1007%2F978-3-319-99525-0), let $Q(x, \\theta)$ be a quantity of interest where $\\theta$ is a parameter. We assume there are additional constraints $f(x, \\theta) = 0$ where the Jacobian $f_x(x, \\theta)$ is everywhere non-singular. We wish to compute the **parametric sensitivity** $Q_\\theta$ where $Q_\\theta$ is the total derivative of $Q$ with respect to $\\theta$. Parameteric sensitivity appears in a wide range of applications including optimization, error estimation, uncertainty quantification, and inverse problems in imaging and siesmology. \n",
"\n",
"The subject of the notebook is the particular case where $Q(x, \\theta)$ is an inner product of $x$ with a weighting function. In the case, the computation of $Q_\\theta$ can be organized into a two step procedure. The first step is to compute an adjoint $y$ which is then followed by computing an inner product of the adjoint with functions derived from the problem data.\n",
"\n",
"This notebook reviews properties of inner products, then demonstrates use of inner products to define and solve for the adjoints of several representative applications of parametric sensitivity.\n",
"\n",
"References\n",
"\n",
"* McClarren, Ryan G., McClarren, and Penrose. Uncertainty Quantification and Predictive Computational Science. Springer International Publishing, 2018. [https://link-springer-com.proxy.library.nd.edu/book/10.1007%2F978-3-319-99525-0](https://link-springer-com.proxy.library.nd.edu/book/10.1007%2F978-3-319-99525-0)\n",
"\n",
"* Bradley, Andrew. PDE-constrained optimization and the adjoint method. [https://cs.stanford.edu/~ambrad/adjoint_tutorial.pdf](https://cs.stanford.edu/~ambrad/adjoint_tutorial.pdf)\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {
"nbpages": {
"level": 1,
"link": "[6.3 Sensitivity Analysis with Adjoint Operators](https://ndcbe.github.io/cbe67701-uncertainty-quantification/06.03-Sensitivity-Analysis-with-Adjoint-Operators.html#6.3-Sensitivity-Analysis-with-Adjoint-Operators)",
"section": "6.3 Sensitivity Analysis with Adjoint Operators"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Requirement already satisfied: autograd in /anaconda3/lib/python3.7/site-packages (1.3)\n",
"Requirement already satisfied: future>=0.15.2 in /anaconda3/lib/python3.7/site-packages (from autograd) (0.18.2)\n",
"Requirement already satisfied: numpy>=1.12 in /anaconda3/lib/python3.7/site-packages (from autograd) (1.16.2)\n"
]
}
],
"source": [
"# install modules\n",
"!pip install autograd"
]
},
{
"cell_type": "markdown",
"metadata": {
"colab_type": "text",
"id": "D5NJgzdET-P-",
"nbpages": {
"level": 2,
"link": "[6.3.1 Review of inner products](https://ndcbe.github.io/cbe67701-uncertainty-quantification/06.03-Sensitivity-Analysis-with-Adjoint-Operators.html#6.3.1-Review-of-inner-products)",
"section": "6.3.1 Review of inner products"
}
},
"source": [
"## 6.3.1 Review of inner products\n",
"\n",
"An **inner product** is the 'multiplication' of two vectors, $u$ and $v$, to produce a scalar result $\\langle u, v\\rangle$ . \n",
"\n",
"For real vectors $u$, $v$, and $w$, and a real scalar $\\alpha$, an inner product must satisfy four basic properties:\n",
"\n",
"1. Additivity: $\\langle u + v, w \\rangle = \\langle u, w\\rangle + \\langle v, w\\rangle $\n",
"2. Homogeneity: $\\langle \\alpha u, w\\rangle = \\alpha \\langle v, w\\rangle$\n",
"3. Symmetry: $\\langle v, w\\rangle = \\langle w, v\\rangle$\n",
"4. Positivity: $\\langle v, v\\rangle \\geq 0$ and equal to zero if and only if $v = 0$.\n",
"\n",
"An **inner product space** is an a vector space along with an inner product defined on that space. Common examples of inner product spaces include the following cases of general interest in computational science and engineering.\n",
"\n",
"1. Product of real numbers $u$ and $v$ $$\\langle u, v\\rangle = u v$$\n",
"1. Dot product of finite dimensional real vectors. For $u, v \\in\\text{R}^N$. $$\\langle u, v\\rangle = u^Tv =\\sum_{n=1}^{N}u_nv_n$$\n",
"1. Real vector-valued functions on a finite interval $t \\in [t_a, t_b]$. Given functions $f$ and $g$ $$\\langle f, g\\rangle = \\int_{t_a}^{t_b}\\langle f(t),g(t)\\rangle\\ dt$$\n",
"1. Expected value of the product of random variables. Given random variables $X$ and $Y$ $$\\langle X, Y\\rangle = \\text{E}(XY)$$\n",
"1. Real matrices. Given $M \\times N$ real matrices $A$ and $B$ in $R^{M \\times N}$ $$\\langle A, B\\rangle = \\text{tr}(A^TB) = \\sum_{m=1}^M\\sum_{n=1}^Na_{mn}b_{mn}$$\n",
"\n",
"A **linear operator**is an operator $L$ that maps a vector into another vector and satisfies properties of additivity and homogeneity.\n",
"\n",
"1. Additivity: $L(u + v) = Lu + Lv$\n",
"2. Homogeneity: $L(\\alpha u) = \\alpha L u$\n",
"\n",
"A basic result for linear operators is that for every linear operator $L$ defined on an inner product space there exists a unique linear **adjoint operator** $L^*$ such that \n",
"\n",
"$$\\langle L^*u, v \\rangle \\equiv \\langle u, Lv \\rangle $$\n",
"\n",
"A **self-adjoint** operator is a combination of linear operator and inner product such that $L = L^*$. Finding the adjoint operator, or finding an inner product (if one exists) for which a linear operator is self-adjoint can be challenging analytical problems.\n"
]
},
{
"cell_type": "markdown",
"metadata": {
"colab_type": "text",
"id": "JLSTKAfVaFyA",
"nbpages": {
"level": 2,
"link": "[6.3.2 Sensitivity Analysis of Linear equations](https://ndcbe.github.io/cbe67701-uncertainty-quantification/06.03-Sensitivity-Analysis-with-Adjoint-Operators.html#6.3.2-Sensitivity-Analysis-of-Linear-equations)",
"section": "6.3.2 Sensitivity Analysis of Linear equations"
}
},
"source": [
"## 6.3.2 Sensitivity Analysis of Linear equations\n",
"\n",
"Suppose there is a quantity of interest (QoI) that can be computed as an inner product\n",
"\n",
"$$Q = \\langle c, x\\rangle$$\n",
"\n",
"where $x$ is the solution to linear equality constraints\n",
"\n",
"\\begin{align*}\n",
"A x & = b\n",
"\\end{align*}\n",
"\n",
"We would like to know how $Q$ changes in response to perturbations in the problem data $A$ and $b$. That is, given perturbations $\\delta A$ and $\\delta b$ in the problem data, we would like to compute $\\delta Q$ where \n",
"\n",
"$$Q + \\delta Q = \\langle c + \\delta c, x + \\delta x\\rangle$$\n",
"\n",
"where $x + \\delta x$ is the solution to \n",
"\n",
"$$(A + \\delta A)(x + \\delta x) = b + \\delta b$$\n",
"\n",
"One way to do this, of course, is simply insert specific perturbations into the problem data and solve the equations again. The issue with this approach is that there are are many potential perturbations one would like to investigate. Solving the original equations for every possible case would get expensive and require significant data handling. Is there a better way to get the desired information?\n",
"\n",
"\n"
]
},
{
"cell_type": "markdown",
"metadata": {
"colab_type": "text",
"id": "2IWO4KPlmDOR",
"nbpages": {
"level": 2,
"link": "[6.3.3 Case 1. Perturbations in the RHS](https://ndcbe.github.io/cbe67701-uncertainty-quantification/06.03-Sensitivity-Analysis-with-Adjoint-Operators.html#6.3.3-Case-1.-Perturbations-in-the-RHS)",
"section": "6.3.3 Case 1. Perturbations in the RHS"
}
},
"source": [
"## 6.3.3 Case 1. Perturbations in the RHS\n",
"\n",
"We begin with the case where perturbations to the problem are restricted to the right hand side of the linear equatiopns. For this case, using the properties of the inner product\n",
"\n",
"\\begin{align*}\n",
"Q + \\delta Q & = \\langle c, x + \\delta x\\rangle \\\\\n",
"& = \\langle c, x \\rangle + \\langle c, \\delta x \\rangle \n",
"\\end{align*}\n",
"\n",
"subject to the linear equality constraint\n",
"\n",
"$$A(x + \\delta x) = b + \\delta b$$\n",
"\n",
"Substracting the nominal value $Q = \\langle c, x \\rangle$ and the nominal linear solutions $Ax = b$ leaves \n",
"\n",
"$$\\delta Q = \\langle c, \\delta x \\rangle$$\n",
"\n",
"subject to\n",
"\n",
"$$A\\delta x = \\delta b$$\n",
"\n",
"One way we could compute $\\delta Q$ is solve the last equation for $\\delta x = A^{-1}\\delta b$ and substitute into $\\delta Q$ to give\n",
"\n",
"$$\\delta Q = \\langle c, A^{-1}\\delta b\\rangle$$\n",
"\n",
"If $A$ isn't too large, this may be a good way to proceed. But for large dimensional problems we may seek a more direct calculation. Using the properties of adjoint\n",
"\n",
"$$\\delta Q= \\langle (A^{-1})^*c, \\delta b\\rangle$$\n",
"\n",
"There is a simple proof that $(A^{-1})^* = (A^*)^{-1}$ which results in\n",
"\n",
"$$\\delta Q = \\langle (A^*)^{-1}c, \\delta b \\rangle$$\n",
"\n",
"Taking this one step further, suppose we had a solution to a second system of equations\n",
"\n",
"$$A^*y = c$$\n",
"\n",
"where $A^*$ is the adjoint matrix, then we could write\n",
"\n",
"$$\\delta Q = \\langle y, \\delta b \\rangle$$\n",
"\n",
"The advantage of this formulation is that we no longer have to invert a matrix. All we have to do is solve a second system of linear equations using any method available to us, including efficient iterative methods. \n",
"\n",
"Further, because the coefficients of $y$ describe the dependence of the quantity of interest on the uncertain parameters, $y$ can be interpreted as sensitivity coefficients. For this reason $y$ is sometimes referred to as the **adjoint sensitivity**. \n",
"\n",
"Often we wish to study the sensitivity of $Q$ to a perturbation in a specific parameter $\\theta$ on which $b(\\theta)$ depends. In this case we can model\n",
"\n",
"$$\\delta b = b_\\theta \\delta \\theta$$\n",
"\n",
"It's a simple matter to show that $\\delta Q = Q_\\theta \\delta \\theta$ where\n",
"\n",
"$$Q_\\theta = \\langle y, b_\\theta \\rangle$$\n",
"\n",
"which shows the elements of $y$ can be interpreted as the sensitivity of $Q$ to perburbations in the corresponding elements of $b$. \n",
"\n",
"To summarize, by using the properties of inner products, we can arrive at this same result without having to introduce a matrix inverse.\n",
"\n",
"\\begin{align*}\n",
"\\delta Q & = \\langle c, \\delta x \\rangle \\\\\n",
"& = \\langle A^*y, \\delta x \\rangle \\\\\n",
"& = \\langle y, A \\delta x \\rangle \\\\\n",
"& = \\langle y, \\delta b \\rangle \\\\\n",
"\\end{align*}\n",
"\n",
"For specific parametric sensitivies,\n",
"\n",
"$$\\boxed{Q_\\theta = \\langle y, b_\\theta \\rangle \\text{ where } A^*y = c}$$\n",
"\n",
"This remarkable expression tells us that we don't need to invert a matrix. To get sensititivity coefficients for a quantity of interest we only need to solve one additional system of linear equations. Let's see how this works out for a simple example.\n"
]
},
{
"cell_type": "markdown",
"metadata": {
"colab_type": "text",
"id": "Hoy6jm8o4H22",
"nbpages": {
"level": 3,
"link": "[6.3.3.1 Example: Multi-product plant](https://ndcbe.github.io/cbe67701-uncertainty-quantification/06.03-Sensitivity-Analysis-with-Adjoint-Operators.html#6.3.3.1-Example:-Multi-product-plant)",
"section": "6.3.3.1 Example: Multi-product plant"
}
},
"source": [
"### 6.3.3.1 Example: Multi-product plant\n",
"\n",
"A plant produces two products, X and Y, the produce per-unit revenue of \\$270 and \\$210, respectively. The fixed costs of production are \\$100 and \\$90, respectively. In addition, each unit of X requires 1 of labor A at \\$50 per hour, and 2 hours of labor B at \\$40 per hour. Each unit of Y requires 1 hour of each type of labor. TA total of 80 hours of labor A and 100 hours of labor B are available for production.\n",
"\n",
"Compute the sensitivity of the total plant profit to changes in the amount of available labor."
]
},
{
"cell_type": "markdown",
"metadata": {
"colab_type": "text",
"id": "3jxtOtKTNI4r",
"nbpages": {
"level": 3,
"link": "[6.3.3.2 Solution](https://ndcbe.github.io/cbe67701-uncertainty-quantification/06.03-Sensitivity-Analysis-with-Adjoint-Operators.html#6.3.3.2-Solution)",
"section": "6.3.3.2 Solution"
}
},
"source": [
"### 6.3.3.2 Solution\n",
"\n",
"Profit for this process is given by\n",
"\n",
"$$\\text{profit} = \\text{revenue} - \\text{expense}$$\n",
"\n",
"Identifying variables as\n",
"\n",
"\\begin{align*}\n",
"x_1 & & \\text{Units of X produced} \\\\\n",
"x_2 & & \\text{Units of Y produced} \\\\\n",
"x_3 & & \\text{Hours of Labor A Used} \\\\\n",
"x_4 & & \\text{Hours of Labor B Used} \\\\\n",
"\\end{align*}\n",
"\n",
"with these definitions\n",
"\n",
"\\begin{align*}\n",
"\\text{revenue} & = 270 x_1 + 210 x_2 \\\\\n",
"\\text{expense} & = 100 x_1 + 90 x_2 + 50 x_3 + 40 x_4 \\\\\n",
"\\end{align*}\n",
"\n",
"subject to equality constraints\n",
"\n",
"\\begin{align*}\n",
"x_3 & = x_1 + x_2 \\\\\n",
"x_4 & = 2 x_1 + x_2 \\\\\n",
"x_3 & = 80 \\\\\n",
"x_4 & = 100 \\\\\n",
"\\end{align*}\n",
"\n",
"Given these modeling equations, the sensitivity formulation becomes\n",
"\n",
"$$Q = \\text{profit} = \\langle c, x \\rangle$$\n",
"\n",
"where\n",
"\n",
"$$c = \\begin{bmatrix} 170 \\\\ 120 \\\\ -50 \\\\ -40 \\end{bmatrix}$$\n",
"\n",
"with equality constraints\n",
"\n",
"$$\\underbrace{\\begin{bmatrix} 1 & 1 & -1 & 0 \\\\ 2 & 1 & 0 & -1 \\\\ 0 & 0 & 1 & 0 \\\\ 0 & 0 & 0 & 1\\end{bmatrix}}_A \\underbrace{\\begin{bmatrix} x_1 \\\\ x_2 \\\\ x_3 \\\\ x_4 \\end{bmatrix}}_x = \\underbrace{\\begin{bmatrix} 0 \\\\ 0 \\\\ 80 \\\\ 100 \\end{bmatrix}}_b$$\n",
"\n",
"We will use $\\theta_A$ and $\\theta_B$ to denote additional labor of type A and B, respectively. The perturbed system is \n",
"\n",
"$$\\underbrace{\\begin{bmatrix} 1 & 1 & -1 & 0 \\\\ 2 & 1 & 0 & -1 \\\\ 0 & 0 & 1 & 0 \\\\ 0 & 0 & 0 & 1\\end{bmatrix}}_A \\underbrace{\\begin{bmatrix} x_1 \\\\ x_2 \\\\ x_3 \\\\ x_4 \\end{bmatrix}}_x = \\underbrace{\\begin{bmatrix} 0 \\\\ 0 \\\\ 80 + \\theta_A \\\\ 100 + \\theta_B \\end{bmatrix}}_b$$\n",
"\n",
"for which \n",
"\n",
"$$\\frac{\\partial b}{\\partial \\theta_A} = \\begin{bmatrix} 0 \\\\ 0 \\\\ 1 \\\\ 0 \\end{bmatrix} \\quad \\frac{\\partial b}{\\partial \\theta_B} = \\begin{bmatrix} 0 \\\\ 0 \\\\ 0 \\\\ 1 \\end{bmatrix}$$\n",
"\n",
"The desired sensitivity coefficients are\n",
"\n",
"\\begin{align*}\n",
"\\frac{\\partial Q}{\\partial \\theta_A} & = \\langle y, \\frac{\\partial b}{\\partial \\theta_A} \\rangle = \\langle \\begin{bmatrix} y_1 \\\\ y_2 \\\\ y_3 \\\\ y_4 \\end{bmatrix}, \\begin{bmatrix} 0 \\\\ 0 \\\\ 1 \\\\ 0 \\end{bmatrix}\\rangle = y_3 \\\\\n",
"\\frac{\\partial Q}{\\partial \\theta_B} & = \\langle y, \\frac{\\partial b}{\\partial \\theta_B} \\rangle = \\langle \\begin{bmatrix} y_1 \\\\ y_2 \\\\ y_3 \\\\ y_4 \\end{bmatrix}, \\begin{bmatrix} 0 \\\\ 0 \\\\ 0 \\\\ 1 \\end{bmatrix}\\rangle = y_4\n",
"\\end{align*}\n",
"\n",
"The following cell computes the adjoint sensitivity coefficients and compares the results to a direct calculation.\n"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/",
"height": 187
},
"colab_type": "code",
"executionInfo": {
"elapsed": 575,
"status": "ok",
"timestamp": 1593721244677,
"user": {
"displayName": "Jeffrey Kantor",
"photoUrl": "https://lh3.googleusercontent.com/a-/AOh14Gg_n8V7bVINy02QRuRgOoMo11Ri7NKU3OUKdC1bkQ=s64",
"userId": "09038942003589296665"
},
"user_tz": 300
},
"id": "6HpWS8Ib3REi",
"nbpages": {
"level": 3,
"link": "[6.3.3.2 Solution](https://ndcbe.github.io/cbe67701-uncertainty-quantification/06.03-Sensitivity-Analysis-with-Adjoint-Operators.html#6.3.3.2-Solution)",
"section": "6.3.3.2 Solution"
},
"outputId": "72b61647-20b5-4508-e130-865e874ea80b"
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"X produced: 20.0\n",
"Y produced: 60.0\n",
"Labor A used: 80.0\n",
"Labor B used: 100.0\n",
"\n",
"profit = 2600.0\n",
"sensitivities [70. 50. 20. 10.]\n",
"\n",
"computed change in profit with one additional hour of A = 20.0\n",
"computed change in profit with one additional hour of B = 10.0\n"
]
}
],
"source": [
"import numpy as np\n",
"\n",
"variables = (\"X produced\", \"Y produced\", \"Labor A used\", \"Labor B used\")\n",
"equations = (\"\")\n",
"\n",
"c = np.array([170, 120, -50, -40])\n",
"A = np.array([[1, 1, -1, 0], [2, 1, 0, -1], [0, 0, 1, 0 ], [0, 0, 0, 1]])\n",
"b = np.array([0, 0, 80, 100])\n",
"\n",
"x = np.linalg.solve(A, b)\n",
"for var, val in zip(variables, x):\n",
" print(f\"{var}: {val}\")\n",
"\n",
"profit = np.dot(c, x)\n",
"print('\\nprofit = ', profit)\n",
"\n",
"y = np.linalg.solve(A.T, c.T)\n",
"print('sensitivities', y)\n",
"\n",
"# Verify\n",
"dbA = np.array([0, 0, 1, 0 ])\n",
"print('\\ncomputed change in profit with one additional hour of A = ', np.dot(c, np.linalg.solve(A, b + dbA)) - profit)\n",
"\n",
"dbB = np.array([0, 0, 0, 1 ])\n",
"print('computed change in profit with one additional hour of B = ', np.dot(c, np.linalg.solve(A, b + dbB)) - profit)"
]
},
{
"cell_type": "markdown",
"metadata": {
"colab_type": "text",
"id": "_prJFftml744",
"nbpages": {
"level": 2,
"link": "[6.3.4 Case 2. Sensitivity to changes in equation coefficients](https://ndcbe.github.io/cbe67701-uncertainty-quantification/06.03-Sensitivity-Analysis-with-Adjoint-Operators.html#6.3.4-Case-2.-Sensitivity-to-changes-in-equation-coefficients)",
"section": "6.3.4 Case 2. Sensitivity to changes in equation coefficients"
}
},
"source": [
"## 6.3.4 Case 2. Sensitivity to changes in equation coefficients\n",
"\n",
"Neglecting second-order terms\n",
"\n",
"$$A\\delta x + \\delta A x = \\delta b$$\n",
"\n",
"Beginning with $\\delta Q$ and a solution to $A^*y = C$\n",
"\n",
"\\begin{align*}\n",
"\\delta Q & = \\langle c, \\delta x \\rangle \\\\\n",
"& = \\langle A^*y, \\delta x \\rangle \\\\\n",
"& = \\langle y, A \\delta x \\rangle \\\\\n",
"& = \\langle y, (\\delta b - \\delta A x) \\rangle \\\\\n",
"& = \\langle y, \\delta b \\rangle - \\langle y, \\delta A x \\rangle \\\\\n",
"\\end{align*}\n",
"\n",
"for perturbations that depend on specific parameters,\n",
"\n",
"\\begin{align*}\n",
"\\implies Q_\\theta & = \\langle y, b_\\theta \\rangle - \\langle y, A_\\theta x \\rangle \n",
"\\end{align*}\n"
]
},
{
"cell_type": "markdown",
"metadata": {
"colab_type": "text",
"id": "qk75XQ8oF6T6",
"nbpages": {
"level": 3,
"link": "[6.3.4.1 Example](https://ndcbe.github.io/cbe67701-uncertainty-quantification/06.03-Sensitivity-Analysis-with-Adjoint-Operators.html#6.3.4.1-Example)",
"section": "6.3.4.1 Example"
}
},
"source": [
"### 6.3.4.1 Example\n",
"\n",
"As a continuation of the previous example, consider modificatins to the manufacturing process modifications that would reduce the labor required to produce each unit of product. \n",
"\n",
"1. Using a sensitivity analysis, estimate the increased profit that would result from reducing the amount of labor B required to produce each unit of product X. The solve for the actual profit as a function of the labor reduction and compare to the sensitivity result.\n",
"\n",
"2. Extend the analysis by computing the sensitivity of profit to changes in the labor required for each type of labor in each product.\n"
]
},
{
"cell_type": "markdown",
"metadata": {
"colab_type": "text",
"id": "MhAuQ3AFtdWf",
"nbpages": {
"level": 3,
"link": "[6.3.4.2 Solution](https://ndcbe.github.io/cbe67701-uncertainty-quantification/06.03-Sensitivity-Analysis-with-Adjoint-Operators.html#6.3.4.2-Solution)",
"section": "6.3.4.2 Solution"
}
},
"source": [
"### 6.3.4.2 Solution\n",
"\n",
"Refering the model equations developed above, the equality constraints read\n",
"\n",
"$$\\underbrace{\\begin{bmatrix} 1 & 1 & -1 & 0 \\\\ \\boxed{2 - \\theta} & 1 & 0 & -1 \\\\ 0 & 0 & 1 & 0 \\\\ 0 & 0 & 0 & 1\\end{bmatrix}}_A \\begin{bmatrix} x_1 \\\\ x_2 \\\\ x_3 \\\\ x_4 \\end{bmatrix} = \\underbrace{\\begin{bmatrix} 0 \\\\ 0 \\\\ 80 \\\\ 100 \\end{bmatrix}}_b$$\n",
"\n",
"where the required labor of type B required to produce a unit of product X has been highlighted. The parameter $\\theta$ has been added to this entry to indicate the reduction in the amount of labor required to produce the first product. Taking the derivative,\n",
"\n",
"$$A_\\theta = \\begin{bmatrix} 0 & 0 & 0 & 0 \\\\ -1 & 0 & 0 & 0 \\\\ 0 & 0 & 0 & 0 \\\\ 0 & 0 & 0 & 0\\end{bmatrix}$$\n",
"\n",
"The desired sensitivity coefficient is computed as\n",
"\n",
"$$ Q_\\theta = - \\langle y, A_\\theta x \\rangle $$\n",
"\n",
"where\n",
"\n",
"\\begin{align*}\n",
"A x & = b \\\\\n",
"A^T y & = c\n",
"\\end{align*}\n",
"\n",
"In this case, $Q_\\theta$ is the sensitivity of profit to reduction in the labor of type B required to produce the first product.\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/",
"height": 533
},
"colab_type": "code",
"executionInfo": {
"elapsed": 534,
"status": "ok",
"timestamp": 1593721473396,
"user": {
"displayName": "Jeffrey Kantor",
"photoUrl": "https://lh3.googleusercontent.com/a-/AOh14Gg_n8V7bVINy02QRuRgOoMo11Ri7NKU3OUKdC1bkQ=s64",
"userId": "09038942003589296665"
},
"user_tz": 300
},
"id": "ChFgIKfaF5Le",
"nbpages": {
"level": 3,
"link": "[6.3.4.2 Solution](https://ndcbe.github.io/cbe67701-uncertainty-quantification/06.03-Sensitivity-Analysis-with-Adjoint-Operators.html#6.3.4.2-Solution)",
"section": "6.3.4.2 Solution"
},
"outputId": "c42fea8e-debd-4e7b-83e4-173d0c230b0e"
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"forward solution = [ 20. 60. 80. 100.]\n",
"profit = 2600.0\n",
"sensitivities [70. 50. 20. 10.]\n",
"[[ 0. 0. 0. 0.]\n",
" [-1. 0. 0. 0.]\n",
" [ 0. 0. 0. 0.]\n",
" [ 0. 0. 0. 0.]]\n",
"sensitivity = 1000.0\n",
"\n",
"matrix of sensitivity coefficients\n"
]
},
{
"data": {
"text/plain": [
"array([[1400., 4200., 5600., 7000.],\n",
" [1000., 3000., 4000., 5000.],\n",
" [ 400., 1200., 1600., 2000.],\n",
" [ 200., 600., 800., 1000.]])"
]
},
"execution_count": 3,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"import numpy as np\n",
"import matplotlib.pyplot as plt\n",
"\n",
"c = np.array([170, 120, -50, -40])\n",
"A = np.array([[1, 1, -1, 0], [2, 1, 0, -1], [0, 0, 1, 0 ], [0, 0, 0, 1]])\n",
"b = np.array([0, 0, 80, 100])\n",
"\n",
"x = np.linalg.solve(A, b)\n",
"print('forward solution = ', x)\n",
"\n",
"profit = np.dot(c, x)\n",
"print('profit = ', profit)\n",
"\n",
"y = np.linalg.solve(A.T, c.T)\n",
"print('sensitivities', y)\n",
"\n",
"dA = np.zeros([4, 4])\n",
"dA[1, 0] = -1\n",
"print(dA)\n",
"dP = -np.dot(y, np.dot(dA, x))\n",
"print('sensitivity =', dP)\n",
"\n",
"# compare sensitivity projection with computational solutions\n",
"def profit(c, A, b):\n",
" return np.dot(c, np.linalg.solve(A, b))\n",
"\n",
"theta = np.linspace(0, 0.5, 100)\n",
"plt.plot(theta, [profit(c, A + t*dA, b) for t in theta], label=\"Sensitivity estimate\")\n",
"plt.plot(theta, [profit(c, A, b) + t*dP for t in theta], label=\"Computed solution\")\n",
"plt.xlabel('Decrease in labor B hours required per unit of A')\n",
"plt.title('Profit')\n",
"plt.legend()\n",
"\n",
"print(\"\\nmatrix of sensitivity coefficients\")\n",
"np.array([[y[i]*x[j] for j in range(len(x))] for i in range(len(y))])"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/",
"height": 312
},
"colab_type": "code",
"executionInfo": {
"elapsed": 570,
"status": "ok",
"timestamp": 1593721569250,
"user": {
"displayName": "Jeffrey Kantor",
"photoUrl": "https://lh3.googleusercontent.com/a-/AOh14Gg_n8V7bVINy02QRuRgOoMo11Ri7NKU3OUKdC1bkQ=s64",
"userId": "09038942003589296665"
},
"user_tz": 300
},
"id": "HwP_87eFxjqz",
"nbpages": {
"level": 3,
"link": "[6.3.4.2 Solution](https://ndcbe.github.io/cbe67701-uncertainty-quantification/06.03-Sensitivity-Analysis-with-Adjoint-Operators.html#6.3.4.2-Solution)",
"section": "6.3.4.2 Solution"
},
"outputId": "7f14b670-707d-4ae2-ee99-64b1a32c2f83"
},
"outputs": [
{
"data": {
"text/plain": [
""
]
}
],
"metadata": {
"colab": {
"authorship_tag": "ABX9TyNdioMTgA7TMKPIBKhcWZ74",
"collapsed_sections": [],
"name": "06.03-Sensitivity-Analysis-with-Adjoint-Operators.ipynb",
"provenance": []
},
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.7.3"
}
},
"nbformat": 4,
"nbformat_minor": 1
}