{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "197a2e0f",
   "metadata": {},
   "outputs": [],
   "source": [
    "# a function of two variables, x and z, that we are minimizing in this example\n",
    "#\n",
    "def objective(x, z):\n",
    "    return x**4 + z**4 + 2*x*x*z*z - 10*x**3 - 3*z**3 - 5*x*x*z - 17*x*z*z \\\n",
    "                + 32*x*x + 36*z*z + 34*x*z - 39*x - 63*z + 68\n",
    "\n",
    "# wrapper around the function above where the arguments are passed in the form of a list;\n",
    "# this format is required by the minimizers from the scipy library\n",
    "#\n",
    "def objlistf(x):\n",
    "    return objective(x[0], x[1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c03511ba",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "step = 1\n",
    "xlow, xhigh = -1, 6\n",
    "zlow, zhigh = 0, 4\n",
    "\n",
    "# create a two-dimensional array using numpy\n",
    "#\n",
    "x_param_values = np.arange(xlow, xhigh + step/2, step)\n",
    "z_param_values = np.arange(zhigh, zlow - step/2, -step)\n",
    "gridx, gridz = np.meshgrid(x_param_values, z_param_values)\n",
    "\n",
    "# evaluate the objective function for all x, z combinations\n",
    "#\n",
    "objective_values = objective(gridx, gridz)\n",
    "\n",
    "# the output already looks like a map of the parameter space\n",
    "#\n",
    "print(\"parameter values for x:\", gridx, sep=\"\\n\")\n",
    "print(\"parameter values for z:\", gridz, sep=\"\\n\")\n",
    "\n",
    "# objective function values over this map, indicating that there are two local minima\n",
    "#\n",
    "print(\"objective function values y = f(x, z):\", objective_values, sep=\"\\n\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e51aac14",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "# same as above, but we are refining the step size in order to draw diagrams\n",
    "#\n",
    "step = 0.05\n",
    "xlow, xhigh = -0.8, 7.2\n",
    "zlow, zhigh = -2.4, 4.2\n",
    "\n",
    "x_param_values = np.arange(xlow, xhigh + step/2, step)\n",
    "z_param_values = np.arange(zhigh, zlow - step/2, -step)\n",
    "gridx, gridz = np.meshgrid(x_param_values, z_param_values)\n",
    "\n",
    "objective_values = objective(gridx, gridz)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bb9c11cb",
   "metadata": {},
   "outputs": [],
   "source": [
    "# first variant of visualizing a function over two variables as a colour map\n",
    "#\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "fig, ax = plt.subplots()\n",
    "fig.set_size_inches(15, 9)\n",
    "plt.xticks(fontsize=18, color=\"#322300\")\n",
    "plt.yticks(fontsize=18, color=\"#322300\")\n",
    "ax.set_xlabel(\"parameter x\", fontsize=24, color=\"#322300\")\n",
    "ax.set_ylabel(\"parameter z\", fontsize=24, color=\"#322300\")\n",
    "\n",
    "plt.imshow(objective_values, cmap = plt.cm.twilight_shifted, interpolation='bilinear', \\\n",
    "           extent=[xlow, xhigh, zlow, zhigh])\n",
    "plt.colorbar()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d1027cdf",
   "metadata": {},
   "outputs": [],
   "source": [
    "# second version, this time using seaborn\n",
    "#\n",
    "import seaborn as sbn\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "fig, ax = plt.subplots()\n",
    "fig.set_size_inches(15, 9)\n",
    "\n",
    "# just sbn.heatmap(objective_values) would also work,\n",
    "# the other two options are just to match the colour\n",
    "# codes slightly better (I think) to the function values\n",
    "#\n",
    "sbn.heatmap(objective_values, robust=True, center=100)\n",
    "\n",
    "plt.xticks(fontsize=18, color=\"#322300\")\n",
    "plt.yticks(fontsize=18, color=\"#322300\")\n",
    "ax.set_xlabel(\"parameter x\", fontsize=24, color=\"#322300\")\n",
    "ax.set_ylabel(\"parameter z\", fontsize=24, color=\"#322300\")\n",
    "\n",
    "xsetoff, xstep = 16, 20\n",
    "xrange = range(xsetoff, len(x_param_values), xstep)\n",
    "zsetoff, zstep = 4, 20\n",
    "zrange = range(zsetoff, len(z_param_values), zstep)\n",
    "\n",
    "plt.xticks([i for i in xrange], [round(x_param_values[i], 2) for i in xrange])\n",
    "plt.yticks([i for i in zrange], [round(z_param_values[i], 2) for i in zrange])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8be429ff",
   "metadata": {},
   "outputs": [],
   "source": [
    "import scipy.optimize as opt\n",
    "\n",
    "initial_point = [3, 3]  # also try what happens with [1, 1] -> it will find the other local minimum\n",
    "\n",
    "# this local minimization scheme (Nelder-Mead) has the\n",
    "# advantage that it does not require us to provide the gradient, etc.,\n",
    "# it only requires three arguments:\n",
    "#\n",
    "# - the function, which needs to be given such that it takes its parameters in list form\n",
    "# - the initial value, in list form\n",
    "# - the threshold for termination (as \"xatol\" value as given below)\n",
    "#\n",
    "local_minimum = opt.minimize(objlistf, initial_point, method='nelder-mead', options={'xatol': 0.0001})\n",
    "\n",
    "final_point = local_minimum.x\n",
    "\n",
    "print(\"Parameter values at minimum: x = \", final_point[0], \", z = \", final_point[1], sep=\"\")\n",
    "print(\"Objective value at minimum: y =\", objlistf(final_point))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bfb1c2df",
   "metadata": {},
   "outputs": [],
   "source": [
    "# constant coefficients used in the function below\n",
    "#\n",
    "base_maintenance = 50000\n",
    "external_cost_coefficient = 0.9\n",
    "internal_cost_coefficient = 0.75\n",
    "lifetime_coefficient = 20\n",
    "market_saturation_onset = 1000000\n",
    "min_investment = 100000\n",
    "productivity_coefficient = 3\n",
    "oversaturation_coefficient = 1.5e-07\n",
    "\n",
    "# yearly deficit of some industrial operation, as a function of:\n",
    "#    - investment i\n",
    "#    - production volume p, overall market value of the products in an undersaturated market\n",
    "#    - depreciation period d (in years): time that the equipment is in use\n",
    "#\n",
    "# taking into account:\n",
    "#    - operating cost (repairs etc. increase if the equipment remains in use for more time)\n",
    "#    - depreciation cost (investment cost per year, i.e., divided by depreciation period)\n",
    "#    - production cost, including labour cost (increases if you need to rely on external providers)\n",
    "#    - income from sales (this is a negative contribution to cost)\n",
    "#\n",
    "def yearly_deficit(i, p, d, debug_output):\n",
    "    if(d < 1):\n",
    "        d = 1  # less than one year is not feasible\n",
    "    operating_cost = base_maintenance + lifetime_coefficient*(d**0.5)*(p**0.5)\n",
    "    \n",
    "    if(i < 0):\n",
    "        i = 0  # making a profit by \"negative investment\" would be nice, but does not work\n",
    "    annual_depreciation = i/d  # the investment is done once and used over d years\n",
    "    \n",
    "    internal_production_limit = max((i - min_investment)*productivity_coefficient, 0)\n",
    "    if (p < 0):\n",
    "        p = 0  # you cannot produce a negative amount of goods\n",
    "    if p > internal_production_limit:\n",
    "        production_cost = internal_cost_coefficient*internal_production_limit \\\n",
    "                          + external_cost_coefficient*(p - internal_production_limit)\n",
    "    else:\n",
    "        production_cost = internal_cost_coefficient*p\n",
    "    \n",
    "    income_per_unit = 1.0\n",
    "    if(p > market_saturation_onset):\n",
    "        income_per_unit -= oversaturation_coefficient * (p - market_saturation_onset)\n",
    "    sales_income = income_per_unit * p\n",
    "    \n",
    "    balance = operating_cost + annual_depreciation + production_cost - sales_income\n",
    "    \n",
    "    if debug_output:\n",
    "        print(\"total investment\\ti = GBP \", round(i, 2), \"\\nproduction volume\\tp = GBP \", \\\n",
    "              round(p, 2), \" per year\\ndepreciation period\\td = \", round(d, 2), \\\n",
    "              \" years\", sep=\"\", end=\"\\n\\n\")\n",
    "        print(\"operating cost:\\tGBP \", round(operating_cost, 2), \" per year\\ndepreciation:\\tGBP \", \\\n",
    "              round(annual_depreciation, 2), \" per year\\nprod.cost:\\tGBP \", round(production_cost, 2), \\\n",
    "              \" per year\\nsales contrib.:\\tGBP \", round(-sales_income, 2), \\\n",
    "              \" per year\\n==\\ntotal deficit:\\tGBP \", round(balance, 2), \\\n",
    "              \" per year\\n==\\n\", sep=\"\", end=\"\\n\")\n",
    "    \n",
    "    return balance\n",
    "\n",
    "# the same as above, with arguments passed in list form, so that the optimizer can access it\n",
    "#\n",
    "def yearly_deficit_list_wrapper(x):\n",
    "    return yearly_deficit(x[0], x[1], x[2], False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8c0749c0",
   "metadata": {},
   "outputs": [],
   "source": [
    "# example parameter choices that yield a deficit;\n",
    "# we would like to see a profit (negative deficit) instead,\n",
    "# for which the opt.minimize() function can help.\n",
    "#\n",
    "investment = min_investment\n",
    "volume = market_saturation_onset\n",
    "years_of_use = 7\n",
    "\n",
    "our_deficit = yearly_deficit(investment, volume, years_of_use, True)"
   ]
  }
 ],
 "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.9.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}