{
"cells": [
{
"cell_type": "markdown",
"metadata": {
"id": "BgpJa9pthV3-"
},
"source": [
"```{index} single: solver; cbc\n",
"```\n",
"```{index} single: solver; cvxpy\n",
"```\n",
"```{index} single: solver; highs\n",
"```\n",
"\n",
"# Extra material: Refinery production and shadow pricing with CVXPY\n",
"\n",
"This is a simple linear optimization problem in six variables, but with four equality constraints it allows for a graphical explanation of some unusually large shadow prices for manufacturing capacity. The notebook presents also contrasts AMPL with CVXPY modeling."
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"id": "immZxWu_qaff",
"outputId": "4b7675bc-207d-42c7-b8e3-d3cba434910b"
},
"outputs": [
{
"output_type": "stream",
"name": "stdout",
"text": [
"\u001b[2K \u001b[90m━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━\u001b[0m \u001b[32m5.6/5.6 MB\u001b[0m \u001b[31m13.4 MB/s\u001b[0m eta \u001b[36m0:00:00\u001b[0m\n",
"\u001b[?25hUsing default Community Edition License for Colab. Get yours at: https://ampl.com/ce\n",
"Licensed to AMPL Community Edition License for the AMPL Model Colaboratory (https://colab.ampl.com).\n"
]
}
],
"source": [
"# install AMPL and solvers\n",
"%pip install -q amplpy\n",
"\n",
"SOLVER = \"cbc\"\n",
"\n",
"from amplpy import AMPL, ampl_notebook\n",
"\n",
"ampl = ampl_notebook(\n",
" modules=[\"coin\"], # modules to install\n",
" license_uuid=\"default\", # license to use\n",
") # instantiate AMPL object and register magics"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {
"id": "_BOJTa5gqafi"
},
"outputs": [],
"source": [
"%pip install -q cvxpy"
]
},
{
"cell_type": "markdown",
"metadata": {
"id": "s4U0ctEDqafj"
},
"source": [
"This example derived from Example 19.3 from Seborg, Edgar, Mellichamp, and Doyle.\n",
"\n",
"> Seborg, Dale E., Thomas F. Edgar, Duncan A. Mellichamp, and Francis J. Doyle III. Process dynamics and control. John Wiley & Sons, 2016.\n",
"\n",
"The changes include updating prices, new solutions using optimization modeling languages, adding constraints, and adjusting parameter values to demonstrate the significance of duals and their interpretation as shadow prices."
]
},
{
"cell_type": "markdown",
"metadata": {
"id": "Slbpo7Joqafk"
},
"source": [
"## Problem data"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/",
"height": 365
},
"id": "pogEPVlQgbuT",
"outputId": "2ef99c23-ef7e-4d5d-e72e-12f580d5da8b"
},
"outputs": [
{
"output_type": "display_data",
"data": {
"text/plain": [
" capacity price\n",
"gasoline 24000 108\n",
"kerosine 2000 72\n",
"fuel oil 6000 63\n",
"residual 2500 30"
],
"text/html": [
"\n",
"
\n",
"
\n",
"\n",
"
\n",
" \n",
" \n",
" | \n",
" capacity | \n",
" price | \n",
"
\n",
" \n",
" \n",
" \n",
" gasoline | \n",
" 24000 | \n",
" 108 | \n",
"
\n",
" \n",
" kerosine | \n",
" 2000 | \n",
" 72 | \n",
"
\n",
" \n",
" fuel oil | \n",
" 6000 | \n",
" 63 | \n",
"
\n",
" \n",
" residual | \n",
" 2500 | \n",
" 30 | \n",
"
\n",
" \n",
"
\n",
"
\n",
"
\n",
"
\n"
]
},
"metadata": {}
},
{
"output_type": "display_data",
"data": {
"text/plain": [
" available price process_cost\n",
"crude 1 28000.0 72.0 1.5\n",
"crude 2 15000.0 45.0 3.0"
],
"text/html": [
"\n",
" \n",
"
\n",
"\n",
"
\n",
" \n",
" \n",
" | \n",
" available | \n",
" price | \n",
" process_cost | \n",
"
\n",
" \n",
" \n",
" \n",
" crude 1 | \n",
" 28000.0 | \n",
" 72.0 | \n",
" 1.5 | \n",
"
\n",
" \n",
" crude 2 | \n",
" 15000.0 | \n",
" 45.0 | \n",
" 3.0 | \n",
"
\n",
" \n",
"
\n",
"
\n",
"
\n",
"
\n"
]
},
"metadata": {}
},
{
"output_type": "display_data",
"data": {
"text/plain": [
" gasoline kerosine fuel oil residual\n",
"crude 1 80 5 10 5\n",
"crude 2 44 10 36 10"
],
"text/html": [
"\n",
" \n",
"
\n",
"\n",
"
\n",
" \n",
" \n",
" | \n",
" gasoline | \n",
" kerosine | \n",
" fuel oil | \n",
" residual | \n",
"
\n",
" \n",
" \n",
" \n",
" crude 1 | \n",
" 80 | \n",
" 5 | \n",
" 10 | \n",
" 5 | \n",
"
\n",
" \n",
" crude 2 | \n",
" 44 | \n",
" 10 | \n",
" 36 | \n",
" 10 | \n",
"
\n",
" \n",
"
\n",
"
\n",
"
\n",
"
\n"
]
},
"metadata": {}
}
],
"source": [
"import pandas as pd\n",
"\n",
"products = pd.DataFrame(\n",
" {\n",
" \"gasoline\": {\"capacity\": 24000, \"price\": 108},\n",
" \"kerosine\": {\"capacity\": 2000, \"price\": 72},\n",
" \"fuel oil\": {\"capacity\": 6000, \"price\": 63},\n",
" \"residual\": {\"capacity\": 2500, \"price\": 30},\n",
" }\n",
").T\n",
"\n",
"crudes = pd.DataFrame(\n",
" {\n",
" \"crude 1\": {\"available\": 28000, \"price\": 72, \"process_cost\": 1.5},\n",
" \"crude 2\": {\"available\": 15000, \"price\": 45, \"process_cost\": 3},\n",
" }\n",
").T\n",
"\n",
"# note: volumetric yields may not add to 100%\n",
"yields = pd.DataFrame(\n",
" {\n",
" \"crude 1\": {\"gasoline\": 80, \"kerosine\": 5, \"fuel oil\": 10, \"residual\": 5},\n",
" \"crude 2\": {\"gasoline\": 44, \"kerosine\": 10, \"fuel oil\": 36, \"residual\": 10},\n",
" }\n",
").T\n",
"\n",
"display(products)\n",
"display(crudes)\n",
"display(yields)"
]
},
{
"cell_type": "markdown",
"metadata": {
"id": "Mx0f8UVO6wu3"
},
"source": [
"## AMPL Model"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"id": "uhFeemwyqafn",
"outputId": "31890b00-eeee-47e9-b462-4579bcdfd5ff"
},
"outputs": [
{
"output_type": "stream",
"name": "stdout",
"text": [
"Profit: 860275.86\n",
"\n"
]
}
],
"source": [
"m = AMPL()\n",
"\n",
"m.eval(\n",
" \"\"\"\n",
" set CRUDES;\n",
" set PRODUCTS;\n",
"\n",
" param price{PRODUCTS union CRUDES};\n",
" param process_cost{CRUDES};\n",
" param yields{PRODUCTS, CRUDES};\n",
" param available{CRUDES};\n",
" param capacity{PRODUCTS};\n",
"\n",
" # decision variables\n",
" var x{c in CRUDES} >= 0 <= available[c];\n",
" var y{p in PRODUCTS} >= 0 <= capacity[p];\n",
"\n",
" # objective\n",
" var revenue = sum{p in PRODUCTS} price[p]*y[p];\n",
" var feed_cost = sum{c in CRUDES} price[c]*x[c];\n",
" var total_process_cost = sum{c in CRUDES} process_cost[c]*x[c];\n",
"\n",
" maximize profit: revenue - feed_cost - total_process_cost;\n",
"\n",
" # constraints\n",
" subject to balances{p in PRODUCTS}:\n",
" y[p] = sum{c in CRUDES} yields[p,c]*x[c] / 100;\n",
"\"\"\"\n",
")\n",
"\n",
"m.set_data(crudes, \"CRUDES\")\n",
"m.set_data(products, \"PRODUCTS\")\n",
"m.param[\"yields\"] = yields.T\n",
"\n",
"# solution\n",
"m.option[\"solver\"] = SOLVER\n",
"m.solve(verbose=False)\n",
"print(f\"Profit: {m.obj['profit'].value():0.2f}\\n\")"
]
},
{
"cell_type": "markdown",
"metadata": {
"id": "0-l7-qpxZs7Y"
},
"source": [
"## CVXPY Model\n",
"\n",
"The `CVXPY` library for disciplined convex optimization is tightly integrated with `numpy`, the standard Python library for the numerical linear algebra. For example, where `AMPL` uses explicit indexing in constraints, summations, and other objects, `CVXPY` uses the implicit indexing implied when doing matrix and vector operations.\n",
"\n",
"Another sharp contrast with `AMPL` is that `CXVPY` has no specific object to describe a set,or to define a objects variables or other modeling objects over arbitrary sets. `CVXPY` instead uses the zero-based indexing familiar to Python users.\n",
"\n",
"The following cell demonstrates these differences by presenting a `CVXPY` model for the small refinery example."
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"id": "ou4iqNTogpvd",
"outputId": "babec380-5652-4ba9-d047-c2fc9438b179"
},
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": [
"860275.8615603454"
]
},
"metadata": {},
"execution_count": 6
}
],
"source": [
"import numpy as np\n",
"import cvxpy as cp\n",
"\n",
"# decision variables\n",
"x = cp.Variable(len(crudes.index), pos=True, name=\"crudes\")\n",
"y = cp.Variable(len(products.index), pos=True, name=\"products\")\n",
"\n",
"# objective\n",
"revenue = products[\"price\"].to_numpy().T @ y\n",
"feed_cost = crudes[\"price\"].to_numpy().T @ x\n",
"process_cost = crudes[\"process_cost\"].to_numpy().T @ x\n",
"profit = revenue - feed_cost - process_cost\n",
"objective = cp.Maximize(profit)\n",
"\n",
"# constraints\n",
"balances = y == yields.to_numpy().T @ x / 100\n",
"feeds = x <= crudes[\"available\"].to_numpy()\n",
"capacity = y <= products[\"capacity\"].to_numpy()\n",
"constraints = [balances, feeds, capacity]\n",
"\n",
"# solution\n",
"problem = cp.Problem(objective, constraints)\n",
"problem.solve()"
]
},
{
"cell_type": "markdown",
"metadata": {
"id": "rc7Bvxy8VCX7"
},
"source": [
"## Crude oil feed results"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/",
"height": 112
},
"id": "kgFALj_rK1b8",
"outputId": "4579b6c4-28d6-4398-cb88-7b2a5e7cc6e3"
},
"outputs": [
{
"output_type": "display_data",
"data": {
"text/plain": [
" available price process_cost consumption shadow price\n",
"crude 1 28000.0 72.0 1.5 26206.9 0.0\n",
"crude 2 15000.0 45.0 3.0 6896.6 0.0"
],
"text/html": [
"\n",
" \n",
"
\n",
"\n",
"
\n",
" \n",
" \n",
" | \n",
" available | \n",
" price | \n",
" process_cost | \n",
" consumption | \n",
" shadow price | \n",
"
\n",
" \n",
" \n",
" \n",
" crude 1 | \n",
" 28000.0 | \n",
" 72.0 | \n",
" 1.5 | \n",
" 26206.9 | \n",
" 0.0 | \n",
"
\n",
" \n",
" crude 2 | \n",
" 15000.0 | \n",
" 45.0 | \n",
" 3.0 | \n",
" 6896.6 | \n",
" 0.0 | \n",
"
\n",
" \n",
"
\n",
"
\n",
"
\n",
"
\n"
]
},
"metadata": {}
}
],
"source": [
"results_crudes = crudes\n",
"results_crudes[\"consumption\"] = x.value\n",
"results_crudes[\"shadow price\"] = feeds.dual_value\n",
"\n",
"display(results_crudes.round(1))"
]
},
{
"cell_type": "markdown",
"metadata": {
"id": "KbaNCZtRVQB5"
},
"source": [
"## Refinery production results"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/",
"height": 174
},
"id": "lem9nBSThVoj",
"outputId": "9f8aa2ef-bac8-4a2b-86ef-1ae891604292"
},
"outputs": [
{
"output_type": "display_data",
"data": {
"text/plain": [
" capacity price production unused capacity shadow price\n",
"gasoline 24000 108 24000.0 0.0 14.0\n",
"kerosine 2000 72 2000.0 0.0 262.6\n",
"fuel oil 6000 63 5103.4 896.6 0.0\n",
"residual 2500 30 2000.0 500.0 0.0"
],
"text/html": [
"\n",
" \n",
"
\n",
"\n",
"
\n",
" \n",
" \n",
" | \n",
" capacity | \n",
" price | \n",
" production | \n",
" unused capacity | \n",
" shadow price | \n",
"
\n",
" \n",
" \n",
" \n",
" gasoline | \n",
" 24000 | \n",
" 108 | \n",
" 24000.0 | \n",
" 0.0 | \n",
" 14.0 | \n",
"
\n",
" \n",
" kerosine | \n",
" 2000 | \n",
" 72 | \n",
" 2000.0 | \n",
" 0.0 | \n",
" 262.6 | \n",
"
\n",
" \n",
" fuel oil | \n",
" 6000 | \n",
" 63 | \n",
" 5103.4 | \n",
" 896.6 | \n",
" 0.0 | \n",
"
\n",
" \n",
" residual | \n",
" 2500 | \n",
" 30 | \n",
" 2000.0 | \n",
" 500.0 | \n",
" 0.0 | \n",
"
\n",
" \n",
"
\n",
"
\n",
"
\n",
"
\n"
]
},
"metadata": {}
}
],
"source": [
"results_products = products\n",
"results_products[\"production\"] = y.value\n",
"results_products[\"unused capacity\"] = products[\"capacity\"] - y.value\n",
"results_products[\"shadow price\"] = capacity.dual_value\n",
"\n",
"display(results_products.round(1))"
]
},
{
"cell_type": "markdown",
"metadata": {
"id": "qvaoi0nbYM6P"
},
"source": [
"## Why is the shadow price of kerosine so high?"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/",
"height": 582
},
"id": "ebwOwbkbYSRt",
"outputId": "e79b8ccc-7bbb-4dc0-d54a-44239cb0f820"
},
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": [
"(0.0, 24000.0)"
]
},
"metadata": {},
"execution_count": 9
},
{
"output_type": "display_data",
"data": {
"text/plain": [
"