{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# setup needed Py packages:\n",
    "import numpy as np\n",
    "from scipy.optimize import differential_evolution as DE\n",
    "import win32com.client"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Instantiate the main OpenDSS Object:\n",
    "try:\n",
    "    dssObj = win32com.client.Dispatch(\"OpenDSSengine.DSS\")\n",
    "except:\n",
    "    print (\"Unable to start the OpenDSS Engine\")\n",
    "    raise SystemExit"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Instantiate the other OpenDSS Object:\n",
    "dssText = dssObj.Text  # Definiranje tekstualnog objekta\n",
    "dssCircuit = dssObj.ActiveCircuit  # Definiranje objekta mreze\n",
    "dssSolution = dssCircuit.Solution  # Definiranje objekta rjesenja\n",
    "dssElem = dssCircuit.ActiveCktElement  # Definiranje objekta elementa mreze"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Compail the OpenDSS file\n",
    "#dssText.Command = r\"Compile 'C:\\Users\\mimit1\\Documents\\OpenDSS\\Exercise3dss.dss'\"\n",
    "dssText.Command = r\"Compile 'C:\\Users\\Barukcic\\Documents\\OpenDSS\\Exercise3dss.dss'\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "('sour.1', 'sour.2', 'sour.3', 'b2.1', 'b2.2', 'b2.3', 'b3.1', 'b3.2', 'b3.3', 'b4.1', 'b4.2', 'b4.3', 'b5.1', 'b5.2', 'b5.3', 'b6.1', 'b6.2', 'b6.3', 'b7.1', 'b7.2', 'b7.3', 'b8.1', 'b8.2', 'b8.3', 'b9.1', 'b9.2', 'b9.3', 'b10.1', 'b10.2', 'b10.3', 'b11.1', 'b11.2', 'b11.3', 'b12.1', 'b12.2', 'b12.3', 'b13.1', 'b13.2', 'b13.3', 'b14.1', 'b14.2', 'b14.3', 'b15.1', 'b15.2', 'b15.3', 'b16.1', 'b16.2', 'b16.3', 'b17.1', 'b17.2', 'b17.3', 'b18.1', 'b18.2', 'b18.3', 'b19.1', 'b19.2', 'b19.3', 'b20.1', 'b20.2', 'b20.3', 'b21.1', 'b21.2', 'b21.3', 'b22.1', 'b22.2', 'b22.3', 'b23.1', 'b23.2', 'b23.3', 'b24.1', 'b24.2', 'b24.3', 'b25.1', 'b25.2', 'b25.3', 'b26.1', 'b26.2', 'b26.3', 'b27.1', 'b27.2', 'b27.3', 'b28.1', 'b28.2', 'b28.3', 'b29.1', 'b29.2', 'b29.3', 'b30.1', 'b30.2', 'b30.3', 'b31.1', 'b31.2', 'b31.3', 'b32.1', 'b32.2', 'b32.3', 'b33.1', 'b33.2', 'b33.3', 'b34.1', 'b34.2', 'b34.3')\n"
     ]
    }
   ],
   "source": [
    "# Get node names:\n",
    "NodeN = dssCircuit.AllNodeNames\n",
    "# and print them:\n",
    "print (NodeN)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "('sour', 'b2', 'b3', 'b4', 'b5', 'b6', 'b7', 'b8', 'b9', 'b10', 'b11', 'b12', 'b13', 'b14', 'b15', 'b16', 'b17', 'b18', 'b19', 'b20', 'b21', 'b22', 'b23', 'b24', 'b25', 'b26', 'b27', 'b28', 'b29', 'b30', 'b31', 'b32', 'b33', 'b34')\n"
     ]
    }
   ],
   "source": [
    "# Repeat for bus names:\n",
    "BusN = dssCircuit.AllBusNames\n",
    "print (BusN)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Nominal voltage of load L2 is:  11.0 kV \n",
      " Nominal active power of L2 is:  230.0 kW \n",
      " Nominal reactive power of L2 is:  142.5 kvar \n",
      " and complex power of L2 is:  (230+142.5j) kVA\n"
     ]
    }
   ],
   "source": [
    "# Get some properties of the load element:\n",
    "dssCircuit.Loads.Name = 'L2'\n",
    "VL2 = dssCircuit.Loads.kV\n",
    "PL2 = dssCircuit.Loads.kW\n",
    "QL2 = dssCircuit.Loads.kvar\n",
    "# and print them:\n",
    "print ('Nominal voltage of load L2 is: ', VL2, 'kV \\n',\n",
    "       'Nominal active power of L2 is: ', PL2, 'kW \\n',\n",
    "       'Nominal reactive power of L2 is: ', QL2, 'kvar \\n',\n",
    "       'and complex power of L2 is: ', complex(PL2, QL2), 'kVA')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Nominal loads are, PL:  [230.0, 0.0, 230.0, 230.0, 0.0, 0.0, 230.0, 230.0, 0.0, 230.0, 137.0, 72.0, 72.0, 72.0, 13.5, 230.0, 230.0, 230.0, 230.0, 230.0, 230.0, 230.0, 230.0, 230.0, 230.0, 137.0, 75.0, 75.0, 75.0, 57.0, 57.0, 57.0, 57.0] \n",
      " and QL:  [142.5, 0.0, 142.5, 142.5, 0.0, 0.0, 142.5, 142.5, 0.0, 142.5, 84.0, 45.0, 45.0, 45.0, 7.5, 142.5, 142.5, 142.5, 142.5, 142.5, 142.5, 142.5, 142.5, 142.5, 142.5, 85.0, 48.0, 48.0, 48.0, 34.5, 34.5, 34.5, 34.5]\n"
     ]
    }
   ],
   "source": [
    "# Get nominal loads:\n",
    "LNP = []\n",
    "LNQ = []\n",
    "load = dssCircuit.Loads\n",
    "Li = load.First \n",
    "nL=0\n",
    "while Li:     \n",
    "    LNP.append(load.kW)\n",
    "    LNQ.append(load.kVAR)\n",
    "    nL +=1\n",
    "    Li = load.Next\n",
    "\n",
    "# print nominal loads:\n",
    "print ('Nominal loads are, PL: ',LNP,\n",
    "      '\\n and QL: ', LNQ)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Let's make the objective function of the optimization problem of optimal allocation of DG in the distribution network in form of PY function:\n",
    "def FC(ind, *addtnP):\n",
    "    # decoding DG allocation from DE individual\n",
    "    LOCdg = int(round(ind[0]))  # decode DG location from DE individual\n",
    "    # print (type(LOCdg)) - control line\n",
    "    Pdg = ind[1]  # decode DG powers from DE individual\n",
    "    Qdg = ind[2]\n",
    "    \n",
    "    # change DG location and power in OpenDSS script:\n",
    "    dssText.Command = 'Edit generator.DG' + ' bus1=' + BusN[LOCdg] + ' kw=' + str(Pdg) + ' kvar=' + str(Qdg)\n",
    "\n",
    "    # execute simulation with changed DG location and powers for each part of load shape curve:\n",
    "    EL = [] # initilzation of energy loss list\n",
    "    VmaxV = [] # initilzation of maximum voltage list\n",
    "    VminV = [] # initilzation of minimum voltage list\n",
    "    for ij in range(len(PL)):\n",
    "        load = dssCircuit.Loads\n",
    "        Li = load.First \n",
    "        nL=0\n",
    "        while Li:     \n",
    "            load.kW = PL[ij] * LNP[nL]\n",
    "            load.kvar = PL[ij] * LNQ[nL]\n",
    "            nL +=1\n",
    "            Li = load.Next\n",
    "        dssSolution.Solve()\n",
    "        # get nodal voltages, maximum and minimum voltages in [p.u]\n",
    "        Vpu = dssCircuit.AllBusVmagPU\n",
    "        Vmax = max(Vpu)\n",
    "        Vmin = min(Vpu)\n",
    "        VmaxV.append(Vmax)\n",
    "        VminV.append(Vmin)\n",
    "        # get network loosses\n",
    "        NL = dssCircuit.Losses[0]/1000  # divide by 1000 to get numerical value in [kW]\n",
    "        # check inequality constraints and penalize OF if they not satisfied (soft penalisation):\n",
    "        if Vmax > Vu:\n",
    "            NL = 2*(1+Vmax-Vu)*NL # increased penalisation by factor 2\n",
    "        if Vmin < Vl:\n",
    "            NL = 2*(1+Vl-Vmin)*NL\n",
    "        EL.append(NL*TL[ij]) # get energy losses for ij-th load level in duration of ij-th period in [kWh]\n",
    "    # use network losses as objective function value\n",
    "    return EL, VminV, VmaxV\n",
    "\n",
    "def OF(ind, *addtnP):\n",
    "    FNC = FC(ind, *addtnP)\n",
    "    return sum(FNC[0]) # total daily energy losses"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [],
   "source": [
    "# set voltage limits in [p.u] (inequality constraints)\n",
    "Vl = 0.9\n",
    "Vu = 1.1\n",
    "\n",
    "# define load shape: load power and time duration in form of Python list:\n",
    "PL = [0.4, 0.8, 1.0]\n",
    "TL = [8, 6, 10]\n",
    "\n",
    "# and forward them as additional parameters of objective function (together with nominal loads)\n",
    "addtnP = (Vl, Vu, PL, TL, LNP, LNQ)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Active network losses are:  2557.5520849898053 [kWh]\n",
      "Active network losses, min and max voltages are:  ([166.08644709898854, 634.7292304204398, 1756.736407470377], [0.980343250929517, 0.9560131126610388, 0.9437487446524763], [0.9994445497306317, 0.9986372653413943, 0.9982302096374378]) [kW, p.u., p.u.]\n"
     ]
    }
   ],
   "source": [
    "# Control cell, check if objective function work well\n",
    "indp = [9.8, 500.0, 350.0]\n",
    "NLp = OF(indp, addtnP)\n",
    "NLpp = FC(indp, addtnP)\n",
    "\n",
    "# total network losses are:\n",
    "print ('Active network losses are: ', NLp, '[kWh]')\n",
    "print ('Active network losses, min and max voltages are: ', NLpp, '[kW, p.u., p.u.]')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[(0.0, 33.4), (0.0, 5000.0), (-2000.0, 2000.0)]\n"
     ]
    }
   ],
   "source": [
    "# define optimization part: (tip: see the notebook you made in exercise 1)\n",
    "# Power limits of DG\n",
    "PDGmin = 0.0\n",
    "PDGmax = 5000.0\n",
    "QDGmin = -2000.0\n",
    "QDGmax = 2000.0\n",
    "\n",
    "# define box constraints (decision variable ranges (limits))\n",
    "BXC = [(0.0, len(BusN)-0.6),\n",
    "      (PDGmin, PDGmax),\n",
    "      (QDGmin, QDGmax)] # 0 - 1000 kW\n",
    "print (BXC)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "##################################################### \n",
      "      fun: 964.5369190867824\n",
      " message: 'Optimization terminated successfully.'\n",
      "    nfev: 12000\n",
      "     nit: 19\n",
      " success: True\n",
      "       x: array([  19.92646986, 2217.46909805, 1358.77310326]) \n",
      " ##################################################### \n",
      "\n",
      "\n",
      " Problem solution is:  [  19.92646986 2217.46909805 1358.77310326] \n",
      " , and total network lossess are:  964.5369190867824 kW \n",
      " for optimal allocation:  \n",
      " DG location in bus:  b21 \n",
      " injecting active power of :  2217.4690980530677 kW \n",
      " and rective power of :  1358.7731032604463 kvar\n",
      "\n",
      " Minimal and maximal nodal voltages are:  [1.00014167817882, 0.9849445585628107, 0.9769826293088005]  and  [1.0163383827716714, 0.9993577518685043, 0.9989585187008927]\n"
     ]
    }
   ],
   "source": [
    "# performe optimization by using DE\n",
    "solution = DE(OF, BXC, addtnP, maxiter=500, popsize=200, mutation=(0.5, 1), recombination=0.7, disp=False, polish=False, tol=0.001)\n",
    "print ('##################################################### \\n', solution, '\\n ##################################################### \\n')\n",
    "print ('\\n','Problem solution is: ', solution.x, '\\n', ', and total network lossess are: ', solution.fun, 'kW',\n",
    "      '\\n for optimal allocation: ', \n",
    "       '\\n DG location in bus: ', BusN[int(round(solution.x[0]))],\n",
    "      '\\n injecting active power of : ', solution.x[1], 'kW',\n",
    "      '\\n and rective power of : ', solution.x[2], 'kvar',)\n",
    "print ('\\n Minimal and maximal nodal voltages are: ', FC(solution.x, addtnP)[1], ' and ', FC(solution.x, addtnP)[2])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      " Energy losses for different load levels are:  [218.0279801321201, 180.93571217174804, 565.5781826046511]\n",
      "\n",
      " and power losses for different load levels are:  [27.25349751651501, 30.15595202862467, 56.55781826046511]\n"
     ]
    }
   ],
   "source": [
    "# get losses for different load levels:\n",
    "print ('\\n Energy losses for different load levels are: ', FC(solution.x, addtnP)[0])\n",
    "print ('\\n and power losses for different load levels are: ', [FC(solution.x, addtnP)[0][iL]/TL[iL] for iL in range(len(PL))])"
   ]
  },
  {
   "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
}