{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# initialize necessary PY packages:\n",
    "import numpy as np\n",
    "from scipy.optimize import differential_evolution as DE"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# code DE individual\n",
    "# ind = [P1, P2, P3, P4]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# define objective function:\n",
    "def OF(ind, *par):\n",
    "    # parameters of generation costs and loads\n",
    "    a1, b1, c1 = par[:3]\n",
    "    a2, b2, c2 = par[3:3+3]\n",
    "    a3, b3, c3 = par[6:6+3]\n",
    "    a4, b4, c4 = par[9:9+3]\n",
    "    load = par[-1]\n",
    "    # decode of DE individual\n",
    "    P1, P2, P3, P4 = ind[0], ind[1], ind[2], ind[3]\n",
    "    # calculate production costs\n",
    "    C1 = a1*P1**2 + b1*P1 + c1\n",
    "    C2 = a2*P2**2 + b2*P2 + c2\n",
    "    C3 = a3*P3**2 + b3*P3 + c3\n",
    "    C4 = a4*P4**2 + b4*P4 + c4\n",
    "    # calculate OF value\n",
    "    OFv = C1 + C2 + C3 + C4# objective function value\n",
    "    # calculate equality constraint\n",
    "    EQC = (P1 + P2 + P3 + P4) - load\n",
    "    # penalize OF due to dissatisfaction of the equality constraint\n",
    "    OFp = OFv + 1e1*abs(EQC)\n",
    "    # hard penalization\n",
    "    #OFp = OFv\n",
    "    #if abs(EQC) > 0.1*par[-1]: # abs(EQC) > 1e-2, EQC != 0\n",
    "    #    OFp = 1e3*OFv\n",
    "    return OFp # measure of DE idividual quality"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# define objective function:\n",
    "def OF1(ind, *par):\n",
    "    # parameters of generation costs and loads\n",
    "    a1, b1, c1 = par[:3]\n",
    "    a2, b2, c2 = par[3:3+3]\n",
    "    a3, b3, c3 = par[6:6+3]\n",
    "    a4, b4, c4 = par[9:9+3]\n",
    "    load = par[-1]\n",
    "    # decode of DE individual\n",
    "    P1, P2, P3, P4 = ind[0], ind[1], ind[2], ind[3]\n",
    "    # calculate production costs\n",
    "    C1 = a1*P1**2 + b1*P1 + c1\n",
    "    C2 = a2*P2**2 + b2*P2 + c2\n",
    "    C3 = a3*P3**2 + b3*P3 + c3\n",
    "    C4 = a4*P4**2 + b4*P4 + c4\n",
    "    # calculate OF value\n",
    "    OFv = C1 + C2 + C3 + C4# objective function value\n",
    "    # calculate equality constraint\n",
    "    EQC = (P1 + P2 + P3 + P4) - load\n",
    "    # penalize OF due to dissatisfaction of the equality constraint\n",
    "    # OFp = OFv + 1e3*abs(EQC)\n",
    "    # hard penalization\n",
    "    OFp = OFv\n",
    "    if abs(EQC) > 0.1*par[-1]: # abs(EQC) > 1e-2, EQC != 0\n",
    "        OFp = 1e3*OFv\n",
    "    return OFp # measure of DE idividual quality"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[(0.0, 1000.0), (0.0, 1000.0), (0.0, 1000.0), (0.0, 1000.0)]"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# define box constraints (decision variable ranges (limits))\n",
    "BXC = 4*[(0.0, 1000.0)] # 0 - 1000 kW\n",
    "print (BXC)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "# define additional parameters:\n",
    "par = [3e-4, 2e-2, 7,\n",
    "      2e-4, 1.5e-2, 5,\n",
    "      8e-4, 2.1e-2, 4,\n",
    "      1e-4, 1.8e-2, 7,\n",
    "      2500]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "     fun: 402.06868577754096\n",
      " message: 'Optimization terminated successfully.'\n",
      "    nfev: 35600\n",
      "     nit: 88\n",
      " success: True\n",
      "       x: array([517.67316208, 787.52615128, 194.88088592, 999.91884713])\n",
      "\n",
      " Problem solution is:  [517.67316208 787.52615128 194.88088592 999.91884713] \n",
      " , and total production cost is:  402.06868577754096 EUR/h\n"
     ]
    }
   ],
   "source": [
    "# performe optimization by using DE\n",
    "solution = DE(OF, BXC, par, maxiter=500, popsize=100, mutation=(0.5, 1), recombination=0.7, disp=False, polish=False, tol=0.001)\n",
    "print (solution)\n",
    "print ('\\n','Problem solution is: ', solution.x, '\\n', ', and total production cost is: ', solution.fun, 'EUR/h')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-0.0009535962085465144"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# check equality constraint (EQC)\n",
    "print (sum(solution.x) - par[-1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "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.6.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}