diff --git a/src/refinery_exercise/refinery_optimization_exercise.ipynb b/src/refinery_exercise/refinery_optimization_exercise.ipynb new file mode 100644 index 0000000..4e813db --- /dev/null +++ b/src/refinery_exercise/refinery_optimization_exercise.ipynb @@ -0,0 +1 @@ +{"cells":[{"cell_type":"markdown","id":"05b43de4","metadata":{"id":"05b43de4"},"source":["# Refinery Optimization Exercise\n","\n","\n","## Problem description\n","![refinery_schema.png](data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAARgAAADsCAYAAACv+rk3AAAABHNCSVQICAgIfAhkiAAAIABJREFUeJzt3XVcFVkfx/EPiDTSYmIgIrbYjR3YdRF77dh1rdW1u9dYFbu7uxN1MbADMCgRFBWQbnj+QO8jggrCcLl43q/Xvl4PM+fO/Ibn+mXmzMw5KklJSUkIgiBIQFXRBQiCkHuJgBEEQTIiYARBkIwIGEEQJCMCRhAEyYiAEQRBMmqKLkAQskJcXByPHz/m9evXvPTwxNjEVNEl5RhaWlqoqUKBAgWwsrLC2Ng42/atIp6DEZTd4sWL+fffFdSoUx/L8lUICwmhmIWlosvKUcJDQ/F56cbVi2do2aI5kyZNonDhwpLvVwSMoLQiIyPp27cf5NWi97CxFDEvruiScryE+Hh2blzFif3bWLd2Lba2tpLuTwSMoLRk9t0xKlSCwaMmKroUpXPt0lkWThnFmdOnKVeunGT7EQEjKKWlS5fifPs+kxY6KroUpXVkzzYe/HeeI0cOS7YPcRdJUDrx8fGsWLGSXkPHKLoUpdbBvjehEZFcvnxZsn2IgBGUzunTp6lQpRrmJSwUXYrSa9iyIxcuXJBs+yJgBKXz+PFjSperrOgycoWKNtV5+OixZNsXASMonfsPHhAUFJSl2/zzt26sWjTzm+snjOjH4hkTMrTNt/6vGShrza3rKS9BQkM+puvzj+7eZqCsNf/M/DtD+80IbW1dvLy8JNu+CBhB6RQpUpQSpbL2OZd7t5154fbkm+sf3b2N+9OHGdpmdFQUd25c48O7APmy3h0ac/zAznR9/sDOTdy5cY0dG1by7q1/hvadXjp6emjr6EqybRABIwiSKW5hyf1XYdh1sgcgPCyUx/dcSM+N26jISM4eP8igP5PPmg7t3iJlqZIRASMIP+GF2xO8PZ4T8jGY00f2sWX1Mi6dOU5EeJi8TUx08hlM4Id3REdFct/lBgCvfby4d+u/727/7LEDxMZEY9vcjqq16nFgxyYSExMlPSYpiHeRBOEnLJs7BTU1NXx9vPB6+Uy+3Lp8ZXadugbAG7/kPpjZy9ZTpnwl/ujbBYC9W9dxePdWbr388M3tH9m7jSLmxbEuX5nODv2Y+Ed/rl06S8OmraQ9sCwmzmAE4WeoqHD14hnKVrRh7xlnrjz0offgkbg9ecDF08dSNbcobc21p34AjJk677vh4v/6FQ/v3qK9rDcATVq1Q1tHhwM7NkpzLBISASMIP0lLW5upC/6ldNkK6BsaMfCPvwB46f40U9s9uHMTAO269gBAXUOT5m0789/lc7z1881c0dlMBIyQ5YKCgvjw4dt/oXOiPGrf7y1ISkokj2qeLxdQtLgF6hqa8kW6evkAiI+Py1Qtxw/soriFJa+8PLhz4xp3blyjWAlLkpKSOLRna6a2nd1EH4wgCVNTUzp37oxMJqNr166KLueHtLV1iIqK+ub69+/eUqZ8JfnPSUlJ5M2rnqqdqqoqCQkJP12H85XzvA94w/sAGChrnWr9oV1bGDJqIqqqynFuoBxVCkrp4MGDdOvWDWNjY4YNG4aTk5OiS/qmIuYl8PZ4nuY6P18fEuLjKVjYXPI6ju7bgbqGJv+5veX+q7AU/42aNJvA9wFcOXdS8jqyiggYQXJBQUGsXr0aW1tbrK2tmT59Om5uboouKwXbFna8e+vP5tVLUyz/8O4tcyf+CUBd26aZ2oe6evIZT9I3bjeHhYZw6cxxGrdsi7aOTqr1bbv0QFVVVd5HowxEwAjZyt3dnRkzZlC2bFkaNmyIo6MjgYGBii6L7v2GUrVWPf6dN5UOtlVYMHUsQxza0ayaJc5OF2jephN1GzXP1D7UNTTR1cvHrs1rmDl+RKr1pw7vJT4+jjafHsz7mqGxCfUbt+DG1YtK09krAkZQmKtXrzJ8+HBMTEzo2rUrBw4cUGg9G/adZuqClRibmnHy0B5euD2hVv3GzF62jgWOKTtXS5ctT9mKVVJto1rtBhQtVhIATS0tqtWuj0l+M/n6sdMWULykJb7engQFvk/xWV9vT+rYNvtukNn3G0LVWvW4d9s5M4eabcSAU0KWCwoK+umBpY2NjZHJZMhkMho0aJBmm1Gjx2BsbkWbzg6ZKVMg+fWFUf0643LrhiTb/2XvIkVFRXHr1i1Fl5ErhYWF/bjRNwQGBuLo6IijoyPW1tbIZDLs7e2xsrLKwgqF7PJLBoyrqysdO3bk+fO07xoIOYObmxvTp09n+vTpNGzYEHt7e2QymaLLEjLgl+yD2bdvnwgXJePk5MTQoUMpWrQo58+fw+NZzroLJaTtlwwYQfm0adOG7du3ExwcTLNmzbGwslZ0SUI6/JKXSF+TctDjX1FYWBjt2rXL9HZq1Kgh74MpVKhQFlQmZDcRMJ9IPQHVryQzw1kWL15cfhepSpXUt4EF5SICRlA4LS0teQduixYtFF3OdyXEx+Pr40l0dBQlLEqjoamVZruY6CheuLtiZGJKwcJFUVFRSbPdW//XfHgXQIlSVujoSjd0paKIgBEUpnXr1vJg+fwYfU6VlJTEv/Onsn/7RvmodXnU1Ghm15Ep81fIH+2Pj4/jz/723Lx2iYT4eCB56Mz1e09hkr+AfHu3/7vC1NFDCHjjJ182cOR4ho2ZnI1HJT3RyStkq2rVqrF48WJ8fX05efIkvXr1yvHhArB3y1q2rF5GjwHDOXLlPudcnjNq0mzOHN3PP7P+P9vAxhWLcXG+yuR5/3L1yWuWb97P+4C3TB09RN4mNjaGiX8MIH/BQuw6dY1L973o7NCP9csXcOPqRUUcnmTEGYwgOXNzc3lnrY2NjaLL+SnXLp2lbEUbho6eJF/Wo/9wnJ0u4uJ8FYDExET2bF1L8zad6CDrBUCDJi0Z8PtfLJ83hZfPXCllVZaLp44S+D6AVdsOYVWuIgCT5i3nyvlTbHZcQu0GTbL/ACUiAkaQhKampjxUWrZsqehyMm3I6ElpvuGsoaFBbEwMAB7PXPkYFEijFm1StGnUog3L503h3q3/KGVVFpcbVylQuKg8XABUVFRo2ro9h3dvJTEhAdU8ecgNRMAIkggODkZTU/PHDZVEhSrVUi3zfOHO9UvnsOuc/PbzKy8PAAwMjVK0K1CoMJA8rgyAj+dL9A0MU23PrGBhYmNjeBfwhgKFimRp/YoiAkbIckZGRj9upOQiwsMYNcAebR0dBv+ZPPNiZGQ4APpfBYyGphYamlpEfuocjoqIQN8g9e/o87LIiHApS89WImAEIYOCAz8wtGd73r19w4Z9p+RnG5+H0Exr/qLEhAT5hGtq6nlJSlK+OY5+hriLJAgZ4OP5kh5tGvLWz5d1u09QrlJV+TpTs4IAhH4MTvGZ6KhI4uJi0dPXByC/WUFCvmoDEPIx+QFFvXz6UpWf7UTACEI6Pbhzk97tGwGw5fBFKthUT7H+83MuX09u/zE4OTiKFrOQt/u6zed2GppaGJuapVqnrETACEI6PL7nwmD7NhQtbsH2Y5cobmGZqk3R4iUpXLQYzk4XUiy/4ZT8bEuNug0BqN2wKQH+r/F84S5vk5SUhPOV81SrXV9pZgxID9EHIwjpMHfSKGJjYyhXyYYDO1IOup1HTY0Bv49DVVWVXoP+YMHUsZSyKotdJ3se3rnJP7P+pn7jFhQpVgJIfjbGvIQFowbYs8BxK2YFCuO4eBYvn7kyZuo8RRyeZETACErnw/t3xObRzrb9vQ94g/vThwDs27Y+1Xp1DU0G/D4OgK49++Pr7cn8KWOYP2UMKioqNG/TiSkLVsjbq6iosGzjXob36kj3VvUAMDAyZuaSNdSq3zgbjuj/IsLDiIn+9nxQmSUCRlA65cqV411YfLbtz9SsIPdfpW8YUNU8eRg7bT4DR47H19uT4haW8hkfv1TcojTHrz3C18eLqKgISpcpr5CH66IiIyhatKhk2xcBIyidSpUqsXJN6jOJnETfwBD9ylW/20Y1Tx6KlSyVTRWlzfXBHSpVrCDZ9nNPb5Lwy2jVqhUez9x4kclJ5gU4f+IAzZtnbr6n7xEBIyilkSP/YPvqxYouQ6nt374BYwN9SQdbEwEjKKWhQ4dibKDH0tkTFV2KUjp34hD7tqxm3ry5ku5HBIygtDZv2kSe+Egm/96H566PFV2OUoiMiGDzigXsWruUXTt3YGmZ+nmerCQ6eQWltmnjBtatW8fMMYMoUqwEpctXJjw8giLmJRRdWo4SHRWBr+cznC6eo7u9jEuXLmbLS6kiYASlN2jQIAYNGsSFCxd49OgRD954E/dRQ9Fl5Sj5dLTp2KYlyxfPJ3/+/Nm2XxEwQq7RtGlTmjZtqugyhC+IPhhBECQjAkYQBMmIgBEEQTLZ1gcTHx+f5khfipCQkJDi5/j4eGJjYxVUzf/FxMSgp6en6DIEIcuoJH0exy8LeXp6cunSJZxv3sLdzQ1/fz8iIyOp26ARN65fzerdZVh4eMoX19TVNXLE3DwxMdFoaGhiampKSQsLbCpXpl69ujRpknumsRB+LVkaMMeOHWPDxs24P3tG7QaNKVelBhalrSlQqGiunBZTCvHxcXx4F4CvtyfuTx7w0MUZH88X2NvbM3BAf8zNzRVdoiCkW5YEzO3bt5k9ew7h0bG07dYn1bwwQua8fuXNmcN72L9jI7+PGM7kyblrelEh98p0wPz7778sWbqMwaMm0qqjfVbVJaThw7u3bF+zBA/3Jyz5ZzE1atRQdEmC8F2ZCpjRo0fj9tyDkZPnU6CwdIPWCCmdPX6QxdPH4bhqFV27dlV0OYLwTT99F2nYsGG8C4li3uqdWVmPkA4t2namiHlxJo0dSmxsLD169FB0SYKQpp86g/n777955uXH1EWOUtQkpNNz18eMHdSdTRs3iEfkhRwpwwGzZcsWNm3dwbLNB3PNBN3KzPnKeVbMm8xVpyuYmpoquhxBSCFDAePt7U2t2nVYsfUQltblpaxLyICNKxcT/fEt69euUXQpgpBChl4VWLhwIX0G/S7CJYfpP2Is9+8/4NKlS4ouRRBSSHfAPHnyhMtXnJD9NlzKeoSf1LXvUFauEn1iQs6S7oDZuXMnbbv2yFXTWuYmLdp25qWHJ48ePVJ0KYIgl+60OHL0GA2at5OyFiGTbFu04fjx44ouQxDk0hUwL1++xNDYhCLmxSUuR8iM6nUbc+Gi6IcRco50BYynpyeW5apIXYuQSRWqVMPNzY3IyEhFlyIIQDoDxs/fn2IlFDvFpZA+FpZWPH/+XNFlCAKQzoAJCgrGtEBBqWsRskD+AoXw9/dXdBmCAKQzYCIio9DW1pG6FiELaOvmIyQkRNFlCAKQzoBJSkpCRUVF6lqELKCimifVkKCCoCjioRZBECQjAkYQBMmIgBEEQTIiYARBkIwIGEEQJCMCRhAEyYiAEQRBMiJgBEGQjAgYQRAk89PTlghCTpWQkCAGRvuKop7EFwEj5Aq7d+/m+PHj3L13n7dv35A3r7qiS8oxHHr15djh/ZS1Lkvz5s3o3bs3RkZG2bJvETCCUrt9+zbjJ/yNgZEpDVt1pP+42Rib5Fd0WTmObOAonj19xM0rZ6hUuQozpk/jt99+k3y/ImAEpXX+/HkcevRk7PSFtGjbWdHl5Gg6unrY1KyLTc26NGtvz79zJvL6tR9Tp06RdL/iQlVQSl5eXgwcNIhZy9aLcMkgyzLlWLBmF2cvXGLbtm2S7ksEjKCU5s6dS9c+Q6hVv7GiS1FK2jo6jPh7FjNmzpJ0iFURMILS8fT0xOnqNXr0F3N0ZYZ1+crUs23KoUOHJNuHCBhB6Tg5OVHPtpmiy8gVajRoxrXr1yXbfo4NmGdPHxEU+F7RZQg5kIeHB4WKi0Hos0KpMuV4/uKlZNvPUQHj7fGcP/p2obaVGfat6tKkSkla1izDgqljiQgPT9E28MM7zp88rKBKBUWKiIxCU0sry7Z3YMdGqpjr8eju7TTXnz6yjyrmejg7XcjQdqeMGkyrWtYplu3cuOq7n/F//Yoq5nqp/uvTsQnL500hOPBDhmr4ER1dPSpWtsnSbX4px9ymDg78QP8uLdHU0sLht6FUqFKdj0GBPLhzgz1b1uLr48XKrQfl7fu0b0z9Ji1pZtdRgVULwreVtLQiLjZG/vOGFYvYuHJxuvqOqtaqR7Xa9QEI+RjMax8vtqxexvEDu9l75j+MTc2yrE7n69eybFtfyzEBc/zAToIC33PqxlMKFjaXL+9g35v8BQuzfvkCXrg/xbJMOQAiIsJJSkpSVLmC8EP9ho1O8XNkeBiJiYnf/czn77RNjToMGTUxxToX56sMsrfj5KE99B48MmuLlUiOCZjQkI8AxMfFp1rXyb4Pqqqq8vdL1i6dR0x0FE8f3mXt0nn0HPg7Orq6ABzdtwMX56t8DA6klJU1Dr8NI3+BQvJt7dm8Bkvr8jxzfcTdm9dp2MyOdl17yD/78O5NAt74Y2hkgkXpMnTrPRAdXT355/1fv+Ls8YM8unsLA0MjOtr35V2APzHRUdh16i5v5/nCnf3bN+Dr44WxiSmNWrTFtrld1v/ihGz34d1bDu7cTLfeAzl+YCf3XW6QN686ZSvaYN93EJpa2gBcOXcSP19vevQfzqUzx3lw5yaxMdGsXTqP6nUaYFOzbob2a2GVfLkVFqo809LkmD6Yeo1bADDpzwG4OF9N0edSoHBRhoyaiEXp5F9wEv8/c/nyfw9xaMf0sUOJjorAvLgFZ48fQtayDs9dH8vb7Nq8mkUzxvPv/Om88XvNreuXAZg6eghLZv1NdFQUFpZleOX9kn/nT2Pa2KHyz/r5+mDfsg5b1ywjf4FCfAwOok/HJiybM5njB3bL2928dgkHuwZcvXAa8+IWBAV+YNQAezauXJzFvzVBEd6/e8uapXOZ/OdAdm1aTeGixYmMjGD5vCksnP6XvN3F08fYsX4lQIqz7SSSfurse8+WtQDUb9wyk0eQfXLMGUzlarVYufUgs/8eySD75L/0lmXKUa9Rc6rXbUjtBk3kbYeMmsjeresoV6mq/DRy+/oV3Lp+mW1HL1OhSjUA/pw4kwGy1sz4awQ7TzjJP//s6SMu3ffC0NgEAPenDzl+YCfzV26mRbsu8nYLpo5lz5a1REZEoK2jw+p/ZqOqqsrBiy7y912unDvJqAH2FClWUv65GX+NoKJNddbtOSlfdv7kYf4a2ptGLdpQ0rJMVv/6BAXw8XzBoUt30dZJnpRw3JBenDy0l6kLVqZq26RVO54+uMPTR/dTXfp86fNbzycP7eG+yw0AIsLC8PJ4TnRUJIvX7KBi1RoSHI00cswZDEDdRs056nSf6YtXU7dRcwLe+LF59VKG9exA/y4tiI2J/uZnr5w7Sc16jeThAqCuoUnPASNwfXQPzxfu8uU2NerIwwWgTLlKHLlyn4bNWqfYptan2SyjoiI+7eMErTrKUrxMZ9vcjmIl/3/L9MbVi7z18+W34WNSbKuZXUc0tbQ5dXhvRn4lQg4m6ztYHi4AZcpXIjYmmpjoqCzdj46eHuUrV8WsYGEWz/obbw/lmXs8x5zBfKauoUn7bj1p360nAC/cnrBx1T+cPXaA7etX0n/E2DQ/5/nCPc0+jsJFiwHg7fFCfuZgVqhIqnYGRsbs27ae+y438Hr5jLf+fhgaGSevTErifcAbIsLDKVS4aKrPFilWkvi4OAD8XnkDMLRH+zTr9Hr57DtHnzsEBQVx7NgxZDIZWll4O1kqamp5gZSX21+Kj0/+//brISAKfXEzAiBPHrVP7RPQ+MlaPl862XWyZ9jYlC8iRkVGImtZm5G/deOo04Of3EP2yjEBM8ShHQ2atsTht2Eplltal2f+ys3c/u8K92470/8bn9fR0SUiLCzV8vCwUOD/p55pTYPr//oV/bu2JD4ujoo21ek/YixNWnXg2L7tLJg2jsTERPQNk8fPCA9PvY+w0I9oaiZ37PFp27OXrcOsYOogy6dv8O1fQi7Sr18/hgwZgkwmQyaT0bp16x9/SEH08ukDEBkRnub6d2/fAKCr96mz/1MIqGTzoFZa2to0btmWrWuW88bvVYq7rTlVjrlECv0YzK6Njt+8DFJBBX0DQ/nPqiqqKYKipGUZ3J6kTvUnD+4CUKKU1Tf3vWP9St76+bL71DX+WbeLtl16oK2jw4f3AQAkJiairq5BScsyqR7GCg8LxeOZq7yWosVKAMlBVq12/RT/vXB/8kvdWo+JiWHbtm3Y2dlhbm7OuHHjuHPnjqLLSuXz3RmX/5zSXP/wzk2ATPedqaiqZnpkOX/fVwDo5VOOP1Q5JmAGj/obP18fhvZonyIovD2eM35YH4IC32PXUSZfrq6pmeIvTq9Bf/Dax4tDu7bIlwW+D2DPlrVYV6hCcQvLb+5bQzP5hDY66v/Xzu5PH7Jv23oA4j5d/vQePJKb1y6xeskcggLf8/KZK+OH9Ulxx6tmvUZYlLbGcfEcvD1eABATHcWi6eNZOO0vXrg/+Zlfj9Lz9fVl8eLFVK9enWrVqrFo0SJevXql6LIAKG5Rmgo21dm9ZS0XTx0FIDYmGrcnD1gwdSzXLp2l54ARaGhm7nJPXV3jh/0z8gD6KohcH91j+7p/OX/yMC3adUFXL1+maskuOeYSqWGz1izdsIfJfw7AoXX9FOu0tLUZO20+dRs1ly+zsi7Psf07ObZ/J3tO/0f1Og34Y8IMZk34nZOH96CqqorXy+foGxiyfNP3O1ZlfQZz+sh+ureuR5nylYiLjcX/9Ss6OfRl65rl+Hg8p4h5cdp364mfrzfrls1n3bL5AFiUtqZw0WJoamrKt7ds415+79uFjo1ssCpXkVdeL4mKjMTht2G06eyQhb815XT37l3u3r3LX3/9RatWrbC3t0cmk6Gh8bM9F5k3999NDO/VgbFDeqZa16Vnf8ZMnZfpfZQuWwGAKuZ62PcbwvgZi77Zdv3yBaxfviDV8mZ2Hfl79pJM15JdVJLScc4+Y9YcilpXy/CDQT8jNiaaZ65P8HjmSlhYCOUq2lCukk2qvx4hH4M5vHsrGpoatGzXVX5X6IX7Ux7du837gDeYFShMx+59Unzu8f076OrppXnJdPH0MTyfu6Gnb0DrjjLy6Rtw58Y1ChYxl3cWA0SEh+P53I2YmGisK1SmR5sGVLSpycwla1Ich7PTRbw9XmBsakr1Og0pkEbnclZbPGMCLW1r0bNn6n8o2SUoKAhjY+MMfUZLS0veX9Oy5fef8xg1egzG5laShPWZo/vx9nyBmlpeihYvSbVa9VI9lh8ZEYHro3uUsiqLgdH/j/Ot/2te+3jJH/H39nhOWEgIFWyqy9ucO3EIH8+X1Kxrm+bt5tiYaB7dc0mztnz6BvKQyirhYaGM6tcZl1s3snS7n+W4gMnJRg2wx9DYJMVzDm/8XtG6djnGz1yMfd/BCqwumbIGzJeKFSsmDxsbm9Qv4kkZML8aqQMmx1wiKYO6ts2YM/FPEuLjqd+kJa+8PDi4cxNFipWgbZec8WWPiY5izpw5bNiwQWE1xMenft0jI3x8fFi4cCELFy6kevXq8rApUkT6M0Aha4mAyYAuPfuTmJjIycN7OHv8EDHRUdSxbca0hStTvK+kSAkJCbi7u+Pu7v7jxkrAxcUFFxcXxo4di52dHTKZ7IcvDAo5hwiYDOrWeyDdeg9UdBm/pJMnT6Krq0tEZBSmxa1//AFB4XLMbWpB+BZbW1vWrFlDcHAwe/bsoVQpMZqdshBnMLmMto4u27dvV+pOXoCyZcvK+16srL79kKSQs4mAEXIMExMTeajUr1//xx8QcrxcFTDb168gMjz1+yR2nbtTxLy4/Of7Ljc4c3Q/gR/eUc+2OS3adUFLWzvFZ2Kioziydzs3r12iiHkJmtp1oFLVmlIfwi+pa9euyGQyOnfO2ROohXwM5ti+7dy77YxZwcLUbdSc+p/GMfrS5wGp3B7fp1ylqrRo1wXzEhap2jk7XeD8ycNERURQr3FzWrbvKn/xMrfINQETHRXJkllpj7NRqVpNecCcPXaACSP6UaOuLSVKlWblopns3LiKNbuOyR+oSkxIYFjPDrxwf0qrDt14/cqLvh2bMmzMZAaOHJ9dh5SrNWzYUH62kl0TsWfGx6BAHOwaEBoSROOW7Xkf8IYxAx1oateBuf9ulLd7H/AGWcs6aGvr0LCZHdcunWHDysUsdNyaYjiQHRtW8s/Mv2nQpCUFChdl0fTx7Nu+gdU7juaYO5JZIdcEjOvj5PeX1u89JX+S8msfgwKZNeEPmrfpxALHrQD0HfInXZrVZOWimUxbmDzi+85Njty77cyWwxfkZy1rl87D8Z/ZNGzWOsufpvxVWFtby0OlTBnlGnRr1yZH3ge84dClOxQtnjy4mPOV8wzv3YlO3fvKv3MLpo0DYPuxKxgam5CYmMjvfbswaWR/Lt7zRENTCz9fH5bPnUr3fkP5a8ZCALr2GoCsZR02rljEH3/PVMxBSiDX3EV69vQhQIoBp7524dQRIsLDGPTnBPmyAoWL0qhFW5zOn5K/6Xx491ZsatRJcUlk328Iamp5uXrxjERHkDsZGRkxbNgwnJyccHV1Zdq0aUoXLgDGpmZ06z1QHi7w/3eL3B7fB5JfIbl46igduvWSv7qiqqqKQ78hRISH8+DTW9knDuwiMTGB/iP+PyhZKauy1KrXCKcLp7PrkLJFrjmDcX/ykJKWZQgOCuSG00VioqOoVL0W1uUry9s8fXQPXb188rF9PytXyYYTB3fx1t8XA0MTvF4+o+kff6Voo29gSKGi5jx9eDdbjkfZdenSBZlMRpcuXX7cWAnI+gxKtezG1YsAWJWvBMDj+8lDeVStVS9Fu3KVqgLg9vghNes1wu3xfYoUK5HqHadylWxwdrpAbGwM6uqKe/EzK+WagHnh9pSAN37Y1S2PhWUZ/Hy9iYyIoJNDX6bMXwEkz72U1oBPBp8Gk/o8sBBAvi/GnvlM39AoRRvh2/bv36/oEiSzd+s6XB/d48ShPQwbO4UadRoCyCdFy2dq1i8IAAAUaUlEQVSQsk8pn4EhqqqqvHvrD0BQ4PsUYxt99nlQs3dv/CnyaVwhZZdrLpFiY2MwNDZh5wkn9p27yeWHPnR26MehXVs4um8HkDyui14aAfN5RLPoqEhiY2IB0NXTT9Uun75hlo+3mhspQ6dtZty4epGHd29BUhKvfTz5GBQIQGxs8ndHL1/KsVpUVVXR0tYhJib5uxMXG5fmgFGf//hF56LvWK45gzlwIeVIc+rqGvw1YyFnjh3g6oVTtO/WE0MjY3zSGDD587CaWlo68tfvI8JDU7ULC/2Y6UGHBOW3bGPy+EIez93o27EpH96/Y9W2QxiZmAIQ8dWwqomJiURGhKOllTxAuKGxSao28MX3UFsn1TpllWvOYNKirqFJcYvSvHntC4BJfjP5BG9fCvkYDICRiSn6BoaoqeVNs13ox48Yf/oSCYJFaWuat+2E85XzJCYkYGKaPNvE5+/TZ2GhISQlJcn/eBmbmhHyMSjV9j5/zsAwc09B5yS5ImD8X79i6ughOF85n2J5UlIS/q99KGSePFiUTc16hIWG4OvtmaLdo3su6Ojqya97q9aqJx/L97PQkI+88vagzKcOPeHXsnLhDNb/uzDV8qSkJNTU8pKYlIildXn08unz9KvvzqN7yWfX1p++O1Vr1uWVt2eqP2KP7t7GvISFfJbS3CBXBIy+gREnD+1m65rlKZYf3LmJ4MAPtP00MFG9Rs0xNjVj65pl8jZv/Xy5cu44zdp0lC/rIOvFzWuXePnMVb5s96bVJCYk0Lxtzn7aVJCG25OHbHb8h7f+r+XLfDxfcu74QWxq1kVNLS9qanmx62TP0X3b5dO7JiYmsmfzGgyNTahRzxaAZm06oaaWl92bVsu35fHcjVv/XaGZXUdyk1zRB6Ojq8tvw8ewYcUiRvTpTNmKVXh09za3rl+mW++B8icoVVVVGTBiLAumjcPP14cSpaw4d+IQSUlJDPzj/0/oNm7ZlpKlrOjfpQWtOnQj4I0fV86dpE1nB0pZlVXUYQoKNHzsZPp3bUW/Ts3oIOuNj9dLLp46irauLlMX/n+EQ1mfwZw6vBcHu/rYNrPj4d2bPL5/h0lzl8nnVdLR1aX3oN9Zs3Qu7k8fYlawMGeO7sfE1Ixeg/7I9mMzMTH5caOflCvOYACGj5vK9MWriYmO4qHLTUzNCjJ1wcpUAyTb9xvCkvW7MclfgJfPnsqneC1U5P9zzKhraLLl8EU69/iNt36+hIeF0mfISCbNXZrdhyWkQUM9L/FxmRs1L6PKVrRh18mrlK1YhTs3rxEfF8egPydw4vrjFOM1F7ewZOeJq9Rr3ALPF+5oamkzae4yuvRMOaPX8HFTmblkDVraOng8d6NmvUas3nk0zdvXUoqPiyP009mWFMSYvLlMThiTV2qOjo48fu7N4DFTFV2K0nt09zan921g965dkmw/15zBCL+OOnXq4OJ8VdFl5Ar3b16hVk3pRgkQASMoncqVK1OooBlnjx9UdClK7WNQIIf2bKNdu3aS7UMEjKCUxv/1F2v/mcUrLw9Fl6K0Vi2cxsABAyhRQrrXEkTACEqpbt26TPx7AlNG9k9zTnLh2xITElg4dTQ66qpMnjxZ0n2JgBGU1oABAxg/bjRjBtizzXERr195K7qkHC0qMpKzR/cxsGtzipgasGXzJsn3mSuegxF+XQ4ODjRo0IA1a9YwfrA9oEqevHnTfGv+V1Whsg0uztd4/cqH5s2bsXjhfBo1apQt+xYBIyi9IkWKMHv2bGbPno2vry/Ozs5iJoKvjBo2kFKlSqGqmr0XLSJghFylaNGiyGQyRZchfCL6YARBkIwIGEEQJCMCRhAEyYiAEQRBMiJgBEGQjAgYQRAkIwJGEATJiIARBEEyImAEQZBMugImJCQEP18fqWsRskBIcCB+fn6KLkMQgHS+KlDAzAxdMR+QUtDW0cbS0lLRZQgCkM4zGAMDfQI/vJO6FiELBH94j6mp+GMg5AzpCpiCBcxSTVYm5Ey+3p6SjlAmCBmRroApWbIkL90eS12LkElv/F6hqpo8fIEg5ATpChhra2ueuT4mJDj1fLpCznH3xjVsGzZUdBmCIJeugFFVVaVtGzuuXTgpdT1CJly7cJLWrVsrugxBkEv3czD29vacOrRbylqETHjy4C7v3vhhZ2en6FIEQS7dAWNra0uB/KacPbJPynqEn3Rw+1qGDxuq6DIEIYUMPck7btxYNq5aJPpicphzJw7xMfAdAwYMUHQpgpBChgKmdu3a9O7Zg5ULxJzAOYWvtycr509l7uxZii5FEFJRSUpKSsrohxwcHDAsVIKBI/+WoiYhnSLCw5kw1IGe3WUMHjxY0eUIQio/9bLjunXr8Hx6j82r/snqeoR0Cg8LZeqf/WnRtLEIFyHH+qmA0dXVZefOnXi73WPV/ClZXZPwAx7P3Rg70J4mDetKPvWnIGTGTw/XYGRkxJEjR9BQiePPfp1xfXQvK+sSvuHovh0M69GOwQP6MXHiREWXIwjf9VN9MF/bvn07M2fNxrZZK9rK+mFewiIrahO+4HzlPId3bURLXY2ZM6ZTuXJlRZckCD+UJQEDEBYWxqpVq9i4aRPW5SpRo0FjylepSXGL0lmx+V9OdFQk7k8f8eSOM1cvniGfni5DBg+iS5cuii5NENItywLmS8eOHePixYvcuHkTX9/XFCxUGN18+ugbGBEW+jGrd5dhH9N4jsfA0EgBlaT0/l0AhoaGBAcFEhwURKWKFalbtw4tW7akWrVqii5PEDJMkoD5UlRUFP7+/oSEhJCQkCDlrtJt/fr1rF+/Xv7z6tWrqVq1qgIrSvb06VNsbGwwNjamcOHCii5HEDItXSPaZYaWlhYWFjmrT+bkyZQvbZYpU4bq1asrqJr/ywk1CEJWEoN+C4IgGREwgiBIRgSMIAiSkbwPRhCyS0BAAI8fP+bGjRuULi0ej/iSnp4epUuXplSpUtm6XxEwgtK7evUqjo6O3Lp9m3IVqxAXn4ibl5gb6jMdHV0+vPXH44Ubhvr69O//G717986WfYuAEZTavHnz2LFrN70Gj2LCog2KLifHu3fbmb3b1nLq9BlWO67C0NBQ0v2JPhhBac2cOZMr12/iuOskLdt3VXQ5SsGmRh3mrNyKWbHSyOy7ExcXJ+n+RMAISunIkSOcOHWG6UvWo28g7V/h3Kjf8LGYW5Zj/Pjxku5HBIyglBYv/ocBf05ES1tb0aUoraFjpnDm3HlcXV0l24cIGEHpXL9+HTV1DWrUtVV0KUotj5oabTp3T/Vke1YSASMoHRcXFypVq63oMnKFStXrctvljmTbF3eRBKXj5++PUYGse7/t8T0XNDQ1KV22Qprrnz68i6qqKtYVqmRou3duXKNYyVKYmhX8qbpioqO4c+Mani+foaOjRymrslSsWuOntvUtRcxL4Ofvn6Xb/JIIGEHpJCQkkidP1p18/zWsN0WKlWD93lNprp82Zijq6hrsOnUt3dt8H/CGgbLWTJyzlK69kqeTWbN0Lvn0DXD4bdgPP79hxSI2rVpMVGQk1uUr8/bNa4IDP1CgUBFWbT9MScsy6a7le/KoqREXF58l20qLuEQSBAlo6+gxeNTflKtkA0BiYiJrl84jPDT0h5+dNeF3Vi2aSb+ho/nP7S27Tl3j0n0v5q/awsfgIEb07kxMdJTUh5AlRMAIwg+oqKiAikqa6z68e4vb4/upluvo6jJk1ETKVrTJ0L7OHN3PoV1b+GPCDAaOHI+2jo58XYu2nZk4Zwlv/F5xeM+2jB2EgohLJEH4gaSkJPhiXLZxQ3qRmJiAgZExh3ZtAcDAyJheA3/nt+FjgORLpObVSzNxzlLqN25Bq9plAVi9ZA6rl8zBxTMINbW8qfZ1dN8OALr3G5JmLW279CAxMYkmrdpl5SFKRpzBCMJPuH75PA/v3GLUpNksXruTCpWrsWLBdG5dv5yqbYHCRVm3J/lWcPtuPVm/91Sa4QLJU9KUr1wVTa1vP9/TvltPdPXyZc2BSEycwQjCT4iNieafdbsoVjL57eTqdRrQsEJRbl2/Qs16jVK1r1qrHgCFihSjWu3639zu+4A3lK2Qe2aMEGcwQpYLCko9qHpOpqKimtzP8s31KqjmySP/OYkkTM0KysMFIJ++ATq6esRER2eqFl29fERnchs5iQgYQRK1a9dm2bJlvHnzRtGl/JCmlhYxMd/+Rx0cFIiGpmaKZTq6eqnaqaiokJiYuYHtCxUtxhu/V99t4/HcjciIiEztJ7uIgBEkcfPmTUaNGkWhQoVo3749u3btIjExUdFlpckkvxlBH95/c33g+wBMTM2ypZZ6jZrzyssDj+duaa6PCA/DoXV9po4elC31ZJYIGEFyx44do0ePHhgYGDBw4EAuXLig6JJSqFHXltc+Xrx8lvqlv5OHdgNQpXrmXk1QVU3+p/a9SzFIvnukraPDvMmj01w/b/JoYmNj6DFgRKbqyS6ik1fINmFhYWzYsIENGzZgYWGBTCbD3t6eChXSfkQ/u3TvN5RDu7bQt2NTmrXpSJ2GTYmLjeXxgzsc3r2VwkWL0cE+8yPA6eXT59rls5iXsKBFu7Rn6DTJX4BpixwZP6wP9q3qMuD3ceQzMMTzuTunjuzl8T0XuvTsn+nAyy4iYASF8PDwYO7cucydO5fatWvLw8bMLHsuRb6ko6vLzhNOLJk9iSN7tnHki4fY2nR2YPTkOd+9bZxe7br2ZOfGVUy4148y5Sun6CT+UvM2nYiOimLL6qWMG9JLvlxbR4ffho9h8J8TMl1LdpF8ZsecaPr06cyYMUP+8+XLl7G1tVVcQblMUFAQxsbGP/XZ9u3bY29vj0wm++blxKjRYzA2t6JNZ4fMlPlNL9yeoKKqSimrspJsPyNCgoPw9fFCVy8fxS0ss3z74WGhjOrXGZdbN7J82yD6YIQc5ujRo3Tv3h0DAwMGDx7MpUuXsr0GS+vyOSJcAPQNjShfuaok4ZIdxCUS0KhR6gejBMUKDQ1l3bp1rFu3jlKlSsnPasqXL6/o0oQMEGcwQo738uVLZs+eTYUKFahXrx4PHtwnIjxM0WUJ6fBLBkyLFi0UXYKQQfr6+gwaNIhZs2ZRuXKVNB90E3KeX/ISqXbt2pw6dYpbt24puhTJfO95C6n79aOioli4cGGWbKtDhw7yy6PPjh0/kSXbFqT3SwYMQKtWrWjVqpWiy8iVgoKCMhUwderUkd+2zp8/fxZWJmS3XzZghJylVKlS8lARHbm5hwgYQWHy5cuHTCZDJpPRpEkTRZeTIdPHDuV9wFtWbT+cat3m1UtxOncS96ePMDI2pXbDxoyfuRh1dQ15G//Xr1i1aCaP7t7i/bsASpayYujYydRvnLv6B3/JTl5BsT6//BgcHMy6deuUKlzi4mKZOX6EfOS5r508tId/502lSLGSzFuxiaZ2HTi0awsrF85M0W7ckF44O12gW59BzFq6FoCxg3vi7fFC8mPITuIMRsgWtWrVknfWFihQQNHl/JSnD+8ybcxQPJ67ffPVgQ0rFlGzXiNmL1sHQKMWbYiPi2PnhpX0GfwHxqZm3Lp+GddH91i+eT8NmrQEoE7DpjSvbsnaZfOYt2JTth2T1MQZjCCZkiVLMnHiRB4+fMiNGzcYOXKk0oYLwMzxvxMZEYHjjiNYlkn9pO9bP1+8PZ7Ttkv3FMvbdLYnMTERZ6fkt8j/u3weAyNj6jVqLm+jo6uHbfM2XL90VtqDyGbiDEaQxPnz52natKmiy8hSw8dNpVY9W9Q1NHFcPCvV+mduTwAwLVAoxfJiJZMf8/fxfAnAC/enmOYvIB/C4bPiJS05FbaXwA/vMDbJHXfPxBmMkOWMjIxyXbgANGjSEnWN5JHtkicaSPk8UXjoRyB5+Mwv6ejqoaqqSmhI8vrQkI/o5dNPtf18BoYAhH1qlxuIgBGEn5A8VVLaDzNqamqlXqalRUJC8gyKSUlJ8qD6kpZ28hxI8fHSzbSY3UTACMJPSOsMxujTZU3Ix+Cv2iYRGRGBjk7y6w0mpvnlZzNfCgkOBEA3F70GIQJGEH5CWmcw+T/1vXwdHkGByeP9mhUqDIBpgYKEhaYRMB+DyaOmhpFp7uh/AREwgpBlipe0xNDYhPu3/kux/P5tZwAqVqkBQI06DXnt48Vb/9cp2t29eZ0y5SqleCBP2YmAEYQskkdNja69BrBv+wbcnz4EkkekW/3PHEqUsqJ8lWoANGndHn1DIxZOG0fCp/6WM0f38+DOTTrIen1z+8pI3KYWlM79e/cweRcs2ZCZmeHQbyj/XT5PD7sGWFqXx9vjBQaGRqzYelB+SaWmlpfJc5cz8Y/+NLEpibGpGZ4v3GnXtQedHPpla72REeF8eP9Osu2LgBGUjq1tQ2LUFNsROnHuMtTypP7no29oxKYDZ7h26Ry+3p6UKFWaOrZNyZtXPUW7Jq3bs9eqLHdvXScqMoKa9RphWaZcdpUvl5iQQImSJSXbvggYQemUKVOGY6fPK7QG6/Lfnj9aXUOTJq3a/XAbxS0sFT7W7jPXx5QtYyXZ9kUfjKB0mjRpwuVzpwkLDVF0KUrP+eJJGjZsKNn2RcAISsfU1JQ+fXqza/1yRZei1G5dv8wL96d07dpVsn2IgBGU0pQpU7h/6zrH9u9UdClKydvjOUtmjGf2rJk/bpwJImAEpaSjo8OmjRs4snM929YsVXQ5SuX6pbOMHdSdCePH0bp1a0n39UvO7CjkHgEBAUyaNIl79x/SrE1HrCpUJTo6ikJFiim6tBwlKjKcl89cuX31PEHv3jJ9+jSaNWsm+X5FwAi5gqurK9evX8fVzR0PT08MPr2ZLEDFSlXwePmMksWLUaNGjWydJlkEjCAIkhF9MIIgSEYEjCAIkhEBIwiCZETACIIgGREwgiBIRgSMIAiSEQEjCIJkRMAIgiAZETCCIEhGBIwgCJIRASMIgmREwAiCIBkRMIIgSEYEjCAIkvkfLXm6s738hXIAAAAASUVORK5CYII=)\n","\n","\n","\n","The figure above shows a simplified schema of a refinery with three units (a, b, and c). Each unit consumes crude oil from the same storage (S).\n","The different units and storage have the following characteristics\n","\n","### Storage:\n","- Capacity (Cs): 500 barrels\n","- Initial volume (Vs): 300 barrels\n","\n","### Units:\n","- Unit A capacity (Ca): 200 barrels\n","- Unit B capacity (Cb): 300 barrels\n","- Unit C capacity (Cc): 100 barrels\n","\n","The volumes of crude oil consumed by units are:\n","- $V_a$: Volume consumed by Unit A\n","- $V_b$: Volume consumed by Unit B\n","- $V_c$: Volume consumed by Unit C\n","\n","Your goal is to model the allocation of crude oil from storage to the processing units to ensure constraints are satisfied and optimize certain objectives.\n","\n","First, let's define problem decision variables.\n","\n","---"]},{"cell_type":"code","execution_count":1,"id":"2859bcfb","metadata":{"colab":{"base_uri":"https://localhost:8080/"},"id":"2859bcfb","executionInfo":{"status":"ok","timestamp":1738056682074,"user_tz":-60,"elapsed":11186,"user":{"displayName":"Saad Balbiyad","userId":"16984689236939898284"}},"outputId":"cd0824c5-72fe-48d9-a14b-cc884a6151e8"},"outputs":[{"output_type":"stream","name":"stdout","text":["Collecting pulp\n"," Downloading PuLP-2.9.0-py3-none-any.whl.metadata (5.4 kB)\n","Downloading PuLP-2.9.0-py3-none-any.whl (17.7 MB)\n","\u001b[2K \u001b[90m━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━\u001b[0m \u001b[32m17.7/17.7 MB\u001b[0m \u001b[31m14.6 MB/s\u001b[0m eta \u001b[36m0:00:00\u001b[0m\n","\u001b[?25hInstalling collected packages: pulp\n","Successfully installed pulp-2.9.0\n"]}],"source":["# install pulp if needed\n","%pip install pulp"]},{"cell_type":"code","execution_count":25,"id":"eccc3231","metadata":{"id":"eccc3231","executionInfo":{"status":"ok","timestamp":1738057420986,"user_tz":-60,"elapsed":233,"user":{"displayName":"Saad Balbiyad","userId":"16984689236939898284"}}},"outputs":[],"source":["# Define decision variables\n","# We can use an optimization library like PuLP\n","from pulp import LpProblem, LpVariable, LpMinimize, LpStatusOptimal, LpStatus, LpConstraint, LpAffineExpression\n","\n","# Create a problem instance\n","problem = LpProblem(\"Refinery_Optimization\", LpMinimize)\n","\n","# Define decision variables\n","V_a = LpVariable(\"V_a\", 0, 200) # Volume for Unit A\n","V_b = LpVariable(\"V_b\", 0, 300) # Volume for Unit B\n","V_c = LpVariable(\"V_c\", 0, 100) # Volume for Unit C\n","V_l = LpVariable(\"V_l\", 0, 1, cat='Binary') # Volume for Unit S\n","V_S = LpVariable(\"V_S\", 0, 300) # Volume for Unit S"]},{"cell_type":"markdown","id":"864f1868","metadata":{"id":"864f1868"},"source":["## Exercise Questions\n","\n","### Step 1: Write constraints modeling the crude oil flow from Storage (S) to the different units"]},{"cell_type":"code","execution_count":26,"id":"e9604727","metadata":{"id":"e9604727","executionInfo":{"status":"ok","timestamp":1738057423767,"user_tz":-60,"elapsed":272,"user":{"displayName":"Saad Balbiyad","userId":"16984689236939898284"}}},"outputs":[],"source":["# Constraints can be added to the problem using natural syntax: problem += x + y <= 1, \"{cstr_name}\"\n","\n","problem += V_a + V_b + V_c == 300 , \"Total_Consumption_Limit\""]},{"cell_type":"markdown","id":"d44b56ef","metadata":{"id":"d44b56ef"},"source":["Refinery units (a,b and c) operate with a certain cost. For each barrel of crude oil processed, each unit incurs a cost ($c_a$, $b_c$ and $c_c$). Let's consider following operating costs:\n","- $c_a$ = 3\n","- $c_b$ = 4\n","- $c_c$=1"]},{"cell_type":"code","execution_count":19,"id":"23ee9a77","metadata":{"id":"23ee9a77","executionInfo":{"status":"ok","timestamp":1738056774230,"user_tz":-60,"elapsed":249,"user":{"displayName":"Saad Balbiyad","userId":"16984689236939898284"}}},"outputs":[],"source":["c_a, c_b, c_c = (3, 4, 1)"]},{"cell_type":"markdown","id":"b91b4e45","metadata":{"id":"b91b4e45"},"source":["### Step 2: Define an objective function minimizing total operation cost of all units?"]},{"cell_type":"code","execution_count":27,"id":"c6077884","metadata":{"scrolled":true,"id":"c6077884","executionInfo":{"status":"ok","timestamp":1738057427764,"user_tz":-60,"elapsed":5,"user":{"displayName":"Saad Balbiyad","userId":"16984689236939898284"}}},"outputs":[],"source":["# Define the objective function (fill the missing part)\n","problem += (c_a*V_a) + (c_b*V_b)+ (c_c*V_c), \"Total_Cost\""]},{"cell_type":"markdown","id":"ee65da96","metadata":{"id":"ee65da96"},"source":["Let's now solve the problem. Can you guess what would be the optimal solution?\n","\n","Use the solve function to solve the LP and confirm wether the solution is coherent."]},{"cell_type":"code","execution_count":21,"id":"c0cc1d8c","metadata":{"id":"c0cc1d8c","colab":{"base_uri":"https://localhost:8080/","height":105},"executionInfo":{"status":"ok","timestamp":1738056791049,"user_tz":-60,"elapsed":233,"user":{"displayName":"Saad Balbiyad","userId":"16984689236939898284"}},"outputId":"1edeb60a-7c41-42f6-fc34-b2ffe53bcf81"},"outputs":[{"output_type":"stream","name":"stdout","text":["Problem solver. Optimal value = 700.0\n","V_a = 200.0\n","V_b = 0.0\n","V_c = 100.0\n"]},{"output_type":"execute_result","data":{"text/plain":["'Optimal'"],"application/vnd.google.colaboratory.intrinsic+json":{"type":"string"}},"metadata":{},"execution_count":21}],"source":["def solve(problem: LpProblem) -> LpStatus:\n"," problem.solve()\n"," status = LpStatus[problem.status]\n"," if problem.status != LpStatusOptimal:\n"," print(\"Problem is infeasible\")\n"," return LpStatus[problem.status]\n"," print(f\"Problem solver. Optimal value = {problem.objective.value()}\")\n"," for var in problem.variables():\n"," print(f\"{var.name} = {var.varValue}\")\n"," return LpStatus[problem.status]\n","\n","solve(problem)"]},{"cell_type":"markdown","source":["### Step 3: Priority constraint\n","In the reality, we have some operational constraints resulting in a priority order between units c and b. Basically, we cannot use unit c if unit b is not fully used.\n","\n","- Write a constraint modelling the priotiy constraint between unit b and unit c.\n","- Solve the problem and analyse results"],"metadata":{"id":"GtLsGIfMQsWv"},"id":"GtLsGIfMQsWv"},{"cell_type":"code","source":["# Add as many constraints/variables as you need\n","problem += 300 - V_b <= V_l * 300\n","problem += V_c <= (1-V_l) * 100\n","\n","# problem += ..., \"Priority constraints\"\n","solve(problem)"],"metadata":{"colab":{"base_uri":"https://localhost:8080/","height":122},"id":"bj3X-OoHRSsf","executionInfo":{"status":"ok","timestamp":1738057432996,"user_tz":-60,"elapsed":233,"user":{"displayName":"Saad Balbiyad","userId":"16984689236939898284"}},"outputId":"d5782a51-4ad6-4bb1-af0f-1be39041ef0c"},"id":"bj3X-OoHRSsf","execution_count":28,"outputs":[{"output_type":"stream","name":"stdout","text":["Problem solver. Optimal value = 1000.0\n","V_a = 200.0\n","V_b = 100.0\n","V_c = 0.0\n","V_l = 1.0\n"]},{"output_type":"execute_result","data":{"text/plain":["'Optimal'"],"application/vnd.google.colaboratory.intrinsic+json":{"type":"string"}},"metadata":{},"execution_count":28}]},{"cell_type":"markdown","id":"71c70a76","metadata":{"id":"71c70a76"},"source":["In the real world, refieneries operate on a 1-week schedule. Each unit has a daily processing capacity:\n","- Unit A: 20 barrels/day\n","- Unit B: 30 barrels/day\n","- Unit C: 20 barrels/day\n","\n","The refinery needs to process all volume in storage unit within a week.\n","\n","Before going any futhure, let's define a function that builds and solves a model given a set of constraints and an objective function."]},{"cell_type":"code","execution_count":null,"id":"41fb0053","metadata":{"id":"41fb0053"},"outputs":[],"source":["def solve_problem(problem_name: str, constraints: list[LpConstraint], objective_function: LpAffineExpression) -> LpStatus:\n"," problem = LpProblem(problem_name, LpMinimize)\n","\n"," for cstr in constraints:\n"," problem += cstr\n"," problem += objective_function\n"," problem.solve()\n"," status = LpStatus[problem.status]\n"," if problem.status != LpStatusOptimal:\n"," print(\"Problem is infeasible\")\n"," else:\n"," print(f\"Problem solver. Optimal value = {problem.objective.value()}\")\n"," return LpStatus[problem.status]\n","\n","def display_results(V_a: list[LpVariable], V_b: list[LpVariable], V_c: list[LpVariable], V_s: list[LpVariable]):\n"," for t in range(T):\n"," print(f\"Day {t + 1}:\")\n"," print(\" V_a:\", V_a[t].varValue)\n"," print(\" V_b:\", V_b[t].varValue)\n"," print(\" V_c:\", V_c[t].varValue)\n"," print(\" V_s:\", V_s[t].varValue)"]},{"cell_type":"markdown","id":"5d4dcf54","metadata":{"id":"5d4dcf54"},"source":["\n","### Step 4:\n","instructions:\n","- extend the model to include daily time steps\n","- ensure appropriate storage volume daily-evolution\n","- ensure daily processing capacities are not exceeded\n","- ensure that volume storage $V_s$ starts at initial value and is entirely consummed at the end of time horizon T\n","- formulate objective function representing total operational cost over the tme horizon\n","\n","Costs of operating units are considered constant on the whole time horizon."]},{"cell_type":"code","execution_count":null,"id":"707df2b6","metadata":{"id":"707df2b6"},"outputs":[],"source":["# Extend model for time dimension\n","T = 7 # Days horizon\n","# Redefine variables\n","V_a_t = [... for t in range(T)]\n","V_b_t = [... for t in range(T)]\n","V_c_t = [... for t in range(T)]\n","V_s_t = [... for t in range(T)]\n","\n","# initialize constraints list\n","constraints = []\n","# add storage evolution constraints\n","for t in range(T):\n"," if t < T-1:\n"," constraints.append(V_s[t+1] == V_s[t] + ...)\n","\n","# storage volume initial value\n","constraints.append(...)\n","# storage volume end value\n","constraints.append(...)\n","\n","# define objective function\n","total_operational_cost = sum(... for t in range(T))"]},{"cell_type":"markdown","id":"6a5c4c64","metadata":{"id":"6a5c4c64"},"source":["### **Step 5: Run Optimization and Explore Solutions**\n","#### Instructions:\n","- Solve the extended problem.\n","- Display the optimal solution over the time horizon."]},{"cell_type":"code","execution_count":null,"id":"0f9f0e74","metadata":{"id":"0f9f0e74"},"outputs":[],"source":["solve_problem(\"refinery_optimization\", constraints=constraints, objective_function=total_operational_cost)\n","display_results(V_a, V_b, V_c, V_s)"]},{"cell_type":"markdown","id":"e94fd0ed","metadata":{"id":"e94fd0ed"},"source":["### **Step 5: Add priority Constraints**\n","#### Instructions:\n","- Add constraints ensuring that unit c is only used after unit a runs at full capabity.\n","- Use binary variables defined below 'MAX_CAPA_a'"]},{"cell_type":"code","execution_count":null,"id":"d5768afd","metadata":{"id":"d5768afd"},"outputs":[],"source":["MAX_CAPA_a = [LpVariable(f\"MAX_CAPA_a_{t}\", 0, 1, cat='Binary') for t in range(T)]\n","\n","# add priority evolution\n","for t in range(T):\n"," constraints.append(...)\n"," constraints.append(...)"]},{"cell_type":"markdown","id":"a0cccc84","metadata":{"id":"a0cccc84"},"source":["- Solve and display the results."]},{"cell_type":"code","execution_count":null,"id":"1fdcccf6","metadata":{"id":"1fdcccf6"},"outputs":[],"source":["solve_problem(\"refinery_optimization\", constraints=constraints, objective_function=total_operational_cost)\n","display_results(V_a, V_b, V_c, V_s)"]},{"cell_type":"markdown","id":"afa856fc","metadata":{"id":"afa856fc"},"source":["### **Step 6: Add Outage Constraints**\n","#### Instructions:\n","- Introduce necessary outages for maintenance.\n"," - Each unit must proceed in an outage at least one time per 7 days\n"," - outage lasts 1 day\n"," \n","We defined turn-on status binary variables 'Y_a, Y_b, Yc' to be used to enforce outages."]},{"cell_type":"code","execution_count":null,"id":"e2d052c3","metadata":{"id":"e2d052c3"},"outputs":[],"source":["# Define binary variables for turn-on status\n","Y_a = [LpVariable(f\"Y_a_{t}\", 0, 1, cat=\"Binary\") for t in range(T)]\n","Y_b = [LpVariable(f\"Y_b_{t}\", 0, 1, cat=\"Binary\") for t in range(T)]\n","Y_c = [LpVariable(f\"Y_c_{t}\", 0, 1, cat=\"Binary\") for t in range(T)]\n","\n","# Link turn-on status to volumes processed\n","for t in range(T):\n"," constraints.append(...)\n"," constraints.append(...)\n"," constraints.append(...)\n","\n","# ensure 1-day outage per week\n","constraints.append(...)\n","constraints.append(...)\n","constraints.append(...)"]},{"cell_type":"markdown","id":"9647966c","metadata":{"id":"9647966c"},"source":["- solve and display results"]},{"cell_type":"code","execution_count":null,"id":"8b2542a4","metadata":{"id":"8b2542a4"},"outputs":[],"source":["solve_problem(\"refinery_optimization\", constraints=constraints, objective_function=total_operational_cost)\n","display_results(V_a, V_b, V_c, V_s)"]},{"cell_type":"markdown","id":"5605d951","metadata":{"id":"5605d951"},"source":["### Step 7: furthure extend the model with refined products\n","In the reality, processing units a,b and c will produce refined products: $product_1, product_2$. Each barrel of a refined product can be sold at a correspondin price: $p_1$ and $p_2$. Furthuremore, each unit converts crude oil to refined products with a specific ratios $r_{unit, product}$. Per example, unit a, after processing $x$ barrels of crude oil will generate $r_{a,1}x$, $r_{a,2}x$ barrels of refined products. We assume that crude ol and refined products have the same density, thus $r_{a,1} + r_{a,2} = 1$.\n","\n","Convertion ratios are define for each unit:\n","- unit a: [$r_{a,1}, r_{a,2}$] = [0.7, 0.3]\n","- unit a: [$r_{b,1}, r_{b,2}$] = [0.6, 0.4]\n","- unit a: [$r_{c,1}, r_{c,2}$] = [0.8, 0.2]\n","\n","a barrel of $product_1$ generates a revenue of $p_1=4$ and $p_2=7$.\n","\n","instructions:\n","- Extend the model to include refined products\n","- Include refined products revenue in the objective function"]},{"cell_type":"code","execution_count":null,"id":"f4ebf058","metadata":{"id":"f4ebf058"},"outputs":[],"source":["# Define binary variables for starting each unit\n","r = {\"a\": [0.7, 0.3],\n"," \"b\": [0.6, 0.4],\n"," \"c\": [0.8, 0.2]}\n","p_1, p_2 = 4, 7\n","\n","# Define refined products variables\n","V_product_1 = [... for t in range(T)]\n","V_product_2 = [... for t in range(T)]\n","\n","# Add refined products production cstr (from crude oil)\n","for t in range(T):\n"," # product_1\n"," constraints.append(...)\n"," # product_2\n"," constraints.append(...)\n","\n","# Add refined products revenue to the objective function\n","total_operational_cost = sum(... for t in range(T))\n","refined_products_revenue = sum(... for t in range(T))\n","objective_function = total_operational_cost + refined_products_revenue"]},{"cell_type":"markdown","id":"1c085b99","metadata":{"id":"1c085b99"},"source":["- solve problem and display results"]},{"cell_type":"code","execution_count":null,"id":"8f7e1812","metadata":{"id":"8f7e1812"},"outputs":[],"source":["solve_problem(\"refinery_optimization\", constraints=constraints, objective_function=total_operational_cost)\n","display_results(V_a, V_b, V_c, V_s)"]},{"cell_type":"code","execution_count":null,"id":"14b5ec35","metadata":{"id":"14b5ec35"},"outputs":[],"source":[]}],"metadata":{"kernelspec":{"display_name":"Python 3 (ipykernel)","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.10.12"},"colab":{"provenance":[]}},"nbformat":4,"nbformat_minor":5} \ No newline at end of file diff --git a/src/refinery_exercise/refinery_optimization_exercise.md b/src/refinery_exercise/refinery_optimization_exercise.md new file mode 100644 index 0000000..3f86aaa --- /dev/null +++ b/src/refinery_exercise/refinery_optimization_exercise.md @@ -0,0 +1,289 @@ +# Refinery Optimization Exercise + + +## Problem description + +The figure above shows a simplified schema of a refinery with three units (a, b, and c). Each unit consumes crude oil from the same storage (S). +The different units and storage have the following characteristics + +### Storage: +- Capacity (Cs): 500 barrels +- Initial volume (Vs): 300 barrels + +### Units: +- Unit A capacity (Ca): 200 barrels +- Unit B capacity (Cb): 300 barrels +- Unit C capacity (Cc): 100 barrels + +The volumes of crude oil consumed by units are: +- $V_a$: Volume consumed by Unit A +- $V_b$: Volume consumed by Unit B +- $V_c$: Volume consumed by Unit C + +Your goal is to model the allocation of crude oil from storage to the processing units to ensure constraints are satisfied and optimize certain objectives. + +First, let's define problem decision variables. + +--- + + +```python +# install pulp if needed +# %pip install pulp +``` + + +```python +# Define decision variables +# We can use an optimization library like PuLP +from pulp import LpProblem, LpVariable, LpMinimize, LpStatusOptimal, LpStatus, LpConstraint, LpAffineExpression + +# Create a problem instance +problem = LpProblem("Refinery_Optimization", LpMinimize) + +# Define decision variables +V_a = LpVariable("V_a", 0, 200) # Volume for Unit A +V_b = LpVariable("V_b", 0, 300) # Volume for Unit B +V_c = LpVariable("V_c", 0, 100) # Volume for Unit C +``` + +## Exercise Questions + +### Step 1: Write constraints modeling the crude oil flow from Storage (S) to the different units + + +```python +# Constraints can be added to the problem using natural syntax: problem += x + y <= 1, "{cstr_name}" + +problem += ..., "Total_Consumption_Limit" +``` + +Refinery units (a,b and c) operate with a certain cost. For each barrel of crude oil processed, each unit incurs a cost ($c_a$, $b_c$ and $c_c$). Let's consider following operating costs: +- $c_a$ = 3 +- $c_b$ = 4 +- $c_c$=1 + + +```python +c_a, c_b, c_c = (3, 4, 1) +``` + +### Step 2: Define an objective function minimizing total operation cost of all units? + + +```python +# Define the objective function (fill the missing part) +problem += ..., "Total_Cost" +``` + +Let's now solve the problem. Can you guess what would be the optimal solution? + +Use the solve function to solve the LP and confirm wether the solution is coherent. + + +```python +def solve(problem: LpProblem) -> LpStatus: + problem.solve() + status = LpStatus[problem.status] + if problem.status != LpStatusOptimal: + print("Problem is infeasible") + return LpStatus[problem.status] + print(f"Problem solver. Optimal value = {problem.objective.value()}") + for var in problem.variables(): + print(f"{var.name} = {var.varValue}") + return LpStatus[problem.status] + +solve(problem) +``` + +In the real world, refieneries operate on a 1-week schedule. Each unit has a daily processing capacity: +- Unit A: 20 barrels/day +- Unit B: 30 barrels/day +- Unit C: 20 barrels/day + +The refinery needs to process all volume in storage unit within a week. + +Before going any futhure, let's define a function that builds and solves a model given a set of constraints and an objective function. + + +```python +def solve_problem(problem_name: str, constraints: list[LpConstraint], objective_function: LpAffineExpression) -> LpStatus: + problem = LpProblem(problem_name, LpMinimize) + + for cstr in constraints: + problem += cstr + problem += objective_function + problem.solve() + status = LpStatus[problem.status] + if problem.status != LpStatusOptimal: + print("Problem is infeasible") + else: + print(f"Problem solver. Optimal value = {problem.objective.value()}") + return LpStatus[problem.status] + +def display_results(V_a: list[LpVariable], V_b: list[LpVariable], V_c: list[LpVariable], V_s: list[LpVariable]): + for t in range(T): + print(f"Day {t + 1}:") + print(" V_a:", V_a[t].varValue) + print(" V_b:", V_b[t].varValue) + print(" V_c:", V_c[t].varValue) + print(" V_s:", V_s[t].varValue) +``` + + +### Step 3: +instructions: +- extend the model to include daily time steps +- ensure appropriate storage volume daily-evolution +- ensure daily processing capacities are not exceeded +- ensure that volume storage $V_s$ starts at initial value and is entirely consummed at the end of time horizon T +- formulate objective function representing total operational cost over the tme horizon + +Costs of operating units are considered constant on the whole time horizon. + + +```python +# Extend model for time dimension +T = 7 # Days horizon +# Redefine variables +V_a_t = [... for t in range(T)] +V_b_t = [... for t in range(T)] +V_c_t = [... for t in range(T)] +V_s_t = [... for t in range(T)] + +# initialize constraints list +constraints = [] +# add storage evolution constraints +for t in range(T): + if t < T-1: + constraints.append(V_s[t+1] == V_s[t] + ...) + +# storage volume initial value +constraints.append(...) +# storage volume end value +constraints.append(...) + +# define objective function +total_operational_cost = sum(... for t in range(T)) +``` + +### **Step 4: Run Optimization and Explore Solutions** +#### Instructions: +- Solve the extended problem. +- Display the optimal solution over the time horizon. + + +```python +solve_problem("refinery_optimization", constraints=constraints, objective_function=total_operational_cost) +display_results(V_a, V_b, V_c, V_s) +``` + +### **Step 5: Add priority Constraints** +#### Instructions: +- Add constraints ensuring that unit c is only used after unit a runs at full capabity. +- Use binary variables defined below 'MAX_CAPA_a' + + +```python +MAX_CAPA_a = [LpVariable(f"MAX_CAPA_a_{t}", 0, 1, cat='Binary') for t in range(T)] + +# add priority evolution +for t in range(T): + constraints.append(...) + constraints.append(...) +``` + +- Solve and display the results. + + +```python +solve_problem("refinery_optimization", constraints=constraints, objective_function=total_operational_cost) +display_results(V_a, V_b, V_c, V_s) +``` + +### **Step 6: Add Outage Constraints** +#### Instructions: +- Introduce necessary outages for maintenance. + - Each unit must proceed in an outage at least one time per 7 days + - outage lasts 1 day + +We defined turn-on status binary variables 'Y_a, Y_b, Yc' to be used to enforce outages. + + +```python +# Define binary variables for turn-on status +Y_a = [LpVariable(f"Y_a_{t}", 0, 1, cat="Binary") for t in range(T)] +Y_b = [LpVariable(f"Y_b_{t}", 0, 1, cat="Binary") for t in range(T)] +Y_c = [LpVariable(f"Y_c_{t}", 0, 1, cat="Binary") for t in range(T)] + +# Link turn-on status to volumes processed +for t in range(T): + constraints.append(...) + constraints.append(...) + constraints.append(...) + +# ensure 1-day outage per week +constraints.append(...) +constraints.append(...) +constraints.append(...) +``` + +- solve and display results + + +```python +solve_problem("refinery_optimization", constraints=constraints, objective_function=total_operational_cost) +display_results(V_a, V_b, V_c, V_s) +``` + +### Step 7: furthure extend the model with refined products +In the reality, processing units a,b and c will produce refined products: $product_1, product_2$. Each barrel of a refined product can be sold at a correspondin price: $p_1$ and $p_2$. Furthuremore, each unit converts crude oil to refined products with a specific ratios $r_{unit, product}$. Per example, unit a, after processing $x$ barrels of crude oil will generate $r_{a,1}x$, $r_{a,2}x$ barrels of refined products. We assume that crude ol and refined products have the same density, thus $r_{a,1} + r_{a,2} = 1$. + +Convertion ratios are define for each unit: +- unit a: [$r_{a,1}, r_{a,2}$] = [0.7, 0.3] +- unit a: [$r_{b,1}, r_{b,2}$] = [0.6, 0.4] +- unit a: [$r_{c,1}, r_{c,2}$] = [0.8, 0.2] + +a barrel of $product_1$ generates a revenue of $p_1=4$ and $p_2=7$. + +instructions: +- Extend the model to include refined products +- Include refined products revenue in the objective function + + +```python +# Define binary variables for starting each unit +r = {"a": [0.7, 0.3], + "b": [0.6, 0.4], + "c": [0.8, 0.2]} +p_1, p_2 = 4, 7 + +# Define refined products variables +V_product_1 = [... for t in range(T)] +V_product_2 = [... for t in range(T)] + +# Add refined products production cstr (from crude oil) +for t in range(T): + # product_1 + constraints.append(...) + # product_2 + constraints.append(...) + +# Add refined products revenue to the objective function +total_operational_cost = sum(... for t in range(T)) +refined_products_revenue = sum(... for t in range(T)) +objective_function = total_operational_cost + refined_products_revenue +``` + +- solve problem and display results + + +```python +solve_problem("refinery_optimization", constraints=constraints, objective_function=total_operational_cost) +display_results(V_a, V_b, V_c, V_s) +``` + + +```python + +```