{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Ch-17, Generating Capacity Reliability Evaluation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## example 17.1 Page 421"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "No. of Units out\t\tCapacity Out MW\t\tCapacity Available MW\t\tProbability\n",
      "0 \t\t\t\t0 \t\t\t\t80 \t\t\t0.9409\n",
      "1 \t\t\t\t40 \t\t\t\t40 \t\t\t0.0582\n",
      "2 \t\t\t\t80 \t\t\t\t0 \t\t\t0.0009\n"
     ]
    }
   ],
   "source": [
    "from math import factorial\n",
    "#given\n",
    "n=2 #number of generating station\n",
    "f=0.03 #F.O.R\n",
    "a=1-f\n",
    "p=40 #generation station power \n",
    "g = range(0,3)\n",
    "Pg = range(0,3) \n",
    "def comb(m,r):\n",
    "    y = factorial(m)/(factorial(m-r)*factorial(r))\n",
    "    return y\n",
    "print \"No. of Units out\\t\\tCapacity Out MW\\t\\tCapacity Available MW\\t\\tProbability\"\n",
    "for i in g:\n",
    "    Pg[i] = comb(n,i)*f**i*a**(n-i)\n",
    "    \n",
    "    print i,'\\t\\t\\t\\t',p*i,'\\t\\t\\t\\t',p*(n-i),'\\t\\t\\t',Pg[i]\n",
    "    \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## example 17.2 Page: 422"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Part(a):\n",
      "\n",
      "number of units out 0 ,capacity out  0MW ,capacity available  30MW ,probability 0.97 \n",
      "\n",
      "number of units out 1 ,capacity out  30MW ,capacity available  0MW ,probability 0.03 \n",
      "Part(b):\n",
      "\n",
      "for exp 17 1 \n",
      "\n",
      "number of units out 0 ,capacity out  0MW ,capacity available 80MW ,probability 0.940900 \n",
      "\n",
      "number of units out 1 ,capacity out  40MW ,capacity available 40MW ,probability 0.058200 \n",
      "combination of 2 set of stations\n",
      "\n",
      "capacity out 0MW ,capacity available 110MW ,individual state probability 0.912673 ,cumulative probability 1.000000\n",
      "\n",
      "capacity out 40MW ,capacity available 70MW ,individual state probability 0.056454 ,cumulative probability 0.087327\n"
     ]
    }
   ],
   "source": [
    "from math import factorial\n",
    "#given\n",
    "n1=2 #number of generating station\n",
    "f1=0.03 #F.O.R\n",
    "a1=1-f1\n",
    "p1=40 #genetaion station power\n",
    "n2=1  #number of genreting station \n",
    "f2=0.03 #F.O.R for second set\n",
    "a2=1-f2\n",
    "p2=30 #generating station power in second set \n",
    "def comb(m,r):\n",
    "    y=factorial(m)/(factorial(m-r)*factorial(r))\n",
    "    return y\n",
    "pg2=range(0,n2+1)\n",
    "co2=range(0,n2+1)\n",
    "ca2=range(0,n2+1)\n",
    "print \"Part(a):\"\n",
    "for i in range(0,n2+1):\n",
    "    pg2[i]=comb(n2,i)*((f2)**i)*((a2)**(n2-i))\n",
    "    co2[i]=p2*i ;ca2[i]=p2*(n2-i)\n",
    "    print \"\\nnumber of units out %d ,capacity out  %dMW ,capacity available  %dMW ,probability %0.2f \"%(i,co2[i],ca2[i],pg2[i])\n",
    "\n",
    "pg1=range(0,n1)\n",
    "co1=range(0,n1)\n",
    "ca1=range(0,n1)\n",
    "print \"Part(b):\"\n",
    "print \"\\nfor exp 17 1 \"\n",
    "for i in range(0,n1):\n",
    "    pg1[i]=comb(n1,i)*((f1)**i)*((a1)**(n1-i))\n",
    "    co1[i]=p1*i ;ca1[i]=p1*(n1-i)\n",
    "    print \"\\nnumber of units out %d ,capacity out  %dMW ,capacity available %dMW ,probability %4f \"%(i,co1[i],ca1[i],pg1[i])\n",
    "\n",
    "print \"combination of 2 set of stations\"\n",
    "tp=1\n",
    "pocg=0\n",
    "for i in range(0,n1):\n",
    "    \n",
    "    for j in range(0,n2):\n",
    "        og=co1[i]+co2[j] #now total system capacity out\n",
    "        cg=ca1[i]+ca2[j]  #now total system capacity available\n",
    "        tp=tp-pocg\n",
    "        pocg=pg1[i]*pg2[j] #individual stste probability\n",
    "        print \"\\ncapacity out %dMW ,capacity available %dMW ,individual state probability %.6f ,cumulative probability %.6f\"%(og,cg,pocg,tp)\n",
    "    \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## example 17.3: Page 423"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "number of units out 0 ,capacity out 0MW ,capacity available 200MW ,probability 0.814506 \n",
      "\n",
      "number of units out 1 ,capacity out 50MW ,capacity available 150MW ,probability 0.171475 \n",
      "\n",
      "number of units out 2 ,capacity out 100MW ,capacity available 100MW ,probability 0.013538 \n",
      "\n",
      "number of units out 3 ,capacity out 150MW ,capacity available 50MW ,probability 0.000475 \n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAEZCAYAAACNebLAAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAG5BJREFUeJzt3X+8HHV97/HXOwkUEE2AaCAQDFoSgoIKEnJpaVa9WKsW\nvCmPBjGAgD76KG2I96oQqE0Oty3BC4qplse9lR9GeRDgqqTJvYBEymKESkATiIbkgIokUEIiQRBi\nCZxP/5g5YbPZc87uOTs7s7vv5+OxD2ZnZ2c+u2z2c2be851VRGBmZlZtVN4FmJlZMblBmJlZTW4Q\nZmZWkxuEmZnV5AZhZmY1uUGYmVlNbhDWtiT1SPpW3nUMRtLhkl6UpLxrMWuUG4RlQtITkj6Q8WYG\nHMQjqSSpL/1yflHSJkm3SHpvlgWlr/v9uwqMeDIi3hgZDDhS4kJJ6yT9Nn2Nt0p6Z7O3Zd3JDcKy\nEgzyBd4iT6Vfzm8EZgAbgFWVX+CNkDS6jsUCaNXewmLgQmAucAAwBVgGfKTRFUka09zSrBO4QVhL\nSfo9SV+R9FR6u1rS3ulj4yT9P0nPSnpO0gpJh1Y89whJ90p6QdJdwPh6txsRT0XEQuBa4Ivp+ian\nexm7/h1IKks6P53+pKT7JH1Z0jZgoaS3SfpXSdskbZV0o6Sx6fLfAg4HVqR7LZ+r3oakiZKWS/q1\npMckfapi2z3pHsCS9DX+VNLxA7yPRwIXAGdERDkidkbEjoi4KSL6X9+u11LxelZV3O+TdIGkXqBX\n0jWSrqzazr9I+u8VtX8n/f/zC0lz633/rT25QVir/Q0wHXhXepsOfCF9bBRwHcmX7OHADuBrFc+9\nCXgQOAj4O+AcGt9LuQ04TtK+AzxeveczHfg58BbgcpK9g38ADgGmAZOAHoCIOAt4EvhouudyVY31\n35wucwhwOnC5pPdVPP6nwFJgLLCc3V9/pQ8AmyLioUFeaz17caelr3Faut3Z/Q9IOgA4BViaNrgV\nwBpgYrr9z0j64BDrtzbmBmGtdibwPyNiW0RsAy4DzgKIiOci4raI+F1E/JbkC3kmJGEv8F7gb9O/\nlleRfGE1ejjn6fQ54+pdPiL+KSL60rp+HhF3pzVsA67ur3EokiYBJwEXR8QrEfEwyR7N2RWLrYqI\nO9PM4kaSJlrLQcAzdb6GwSyKiOcj4j+AHwIh6eT0sdOB+yPiGeAEYHxE/H1EvBoRv0xrP6MJNVhB\n+bijtdpE4FcV959M5yFpP5Iv3D8mOaYOsH96BtBEYHtE7Kh47q9I/oJvxKEkf1U/D0yoY/lNlXck\nTSA59v+HwBtJ/sh6rs5tTwSei4iXKuY9SdL4+m2pmH4Z2EfSqIjoq1rXr0n2QkZq1+uLiJB0M/Bx\nYBVJM/9m+vBbgYmStlc8dzTwgybUYAXlPQhrtaeByRX3DweeSqc/SxK0To+IsSR/mSu9/TtwQNpE\n+r2Vxg8x/Tfgx2mj6f+irlznwVXLV6//cuA14J1pjWex+7+jwep5GjhQ0v4V8w4HNtdZe6W7gcMG\nyihSLwFvqLhf/dpgz3qXAqdLeivJoafvpPOfBH4ZEQdU3N4UER8dRu3WJtwgLEt7S9qn4jaG5Avo\nC5LGSxoPLCA5lAKwP0nu8BtJBwIL+1cUEb8CHgIuk7SXpD8E6vpySk8HPVTSQuB84NJ0nVtJmtNZ\nkkZLOg94+xCr25/ki/eFNED/fNXjWwZaR0RsAu4HFqVh/bHAeRWvv24R8RhwDUk+MFNS/3t9hqSL\n08XWArMk7Svp90le+1DrXQtsIzl8dGdEvJA+tBp4UdJF6fpGS3qnMj5t2PLlBmFZup3kMEn/bQHw\n9yRf9I+kt4fSeQBfAfYl+YK6H7iD3f/CPRM4keSQzgJgySDbDpJDIi8CL5J8wb0DmBkR369Y7tMk\nX/LbgKOB+6rWUf0X9mXAccBvSDKQ71Qts4ikAW6X9D8q1tPv4yR7UE8D3wUWRMS/DrK9AfdIIuJC\nkhD7n4DtwOMkofPydJGrgVdImtYNJI2ocn0Drfsm4P3pf/u31UfSkN8N/ALYCvwz8KaB6rP2p6x+\nMEjS9STnYz8bEcdUzJ9Lcnrea8D/j4iL0/mXkPw19RpwYUTclUlhZmZWlyxD6huAr/J6yEV6Ot+p\nwLERsVPSm9P5R5OcXnc0SYj4fUlTagRzZmbWIpkdYkpPQ9xeNfsvSU6r25kuszWdfxqwND118AmS\nXeXpWdVmZmZDa3UGcSTwR5J+lI7y7A+4JrL7mRybSfYkzMwsJ60eBzEGOCAiZkg6AbgVeNsAy+Z9\nHR8zs67W6gaxmeTMDSLiwfRaMONJTjWsHPB0GK+fG7+LJDcNM7NhiIiGLyLZ6kNMy0hOn0PSFGDv\n9HIFy4Ez0nO5jyA5FLW61goiwrcm3RYuXJh7DZ108/vp97Oot+HKbA9C0lKSkbAHSdpEct769cD1\nktaRnJ99NkBErJd0K7AeeBW4IEbyqszMbMQyaxAR8fEBHjprgOUvJ7mMgZmZFYBHUnexUqmUdwkd\nxe9nc/n9zF9mI6mzIMlHnszMGiSJaIOQ2szM2oQbhJmZ1eQGYWZmNblBmJlZTW4QZmZWkxuEmZnV\n5AZhZmY1uUGYmVlNbhBmZlaTG4SZmdXkBmFmZjW1XYPwpZjMzFqj7RrEjBmwalXeVZiZdb62axAX\nXghz5sCsWdDbm3c1Zmadq+0axCc+ARs2wPTpcNJJMG8ebNuWd1VmZp2n7RoEwL77wvz58Oij8Npr\nMG0aXHkl/O53eVdmZtY52rJB9Hvzm+FrX4Mf/jC5TZsGN9/sINvMrBk66hflymX43Odg9Gi46io4\n+eTW1WZmVlTD/UW5jmoQAH19sHQpXHopHH88XHEFTJnSogLNzArIPzmaGjXKQbaZWTN0XIPo5yDb\nzGxkOrZB9HOQbWY2PB2XQQzFQbaZdRuH1A1wkG1m3cQhdQMcZJuZDa0rG0Q/B9lmZgPLrEFIul7S\nFknrajz2WUl9kg6smHeJpMckbZD0wazqqsVBtpnZnjLLICSdDPwW+GZEHFMxfxLwdWAqcHxEPCfp\naOAm4ATgUOD7wJSI6KtaZ1MyiKE4yDazTlK4DCIiVgHbazz0ZeCiqnmnAUsjYmdEPAE8DkzPqrah\nlEqwerUvLW5m3a2lGYSk04DNEfFI1UMTgc0V9zeT7EnkxkG2mXW7ljUISfsBlwILK2cP8pRCJAAO\nss2sW41p4bbeDkwGHpYEcBjwY0knAk8BkyqWPSydt4eenp5d06VSiVKplEmx1fqD7Llz4aKL4Jpr\nYNEimD0b1PCRPTOz7JTLZcrl8ojXk+lAOUmTgRWVIXXFY79kz5B6Oq+H1L9fnUi3KqSuh4NsM2sX\nhQupJS0F7gemSNok6dyqRXZ900fEeuBWYD1wB3BBYTrBABxkm1mn68pLbTTbjh2weHGyJ3HmmbBg\nAYwfn3dVZmaJwu1BdJPKILuvz0G2mXUGN4gm8ohsM+skPsSUIQfZZlYEvtx3QfnS4maWN2cQBVVr\nRPaFF3pEtpkVnxtEizjINrN24wbRYg6yzaxdOIPImYNsM8uaQ+o25iDbzLLkkLqNOcg2syJygygQ\nB9lmViRuEAXkINvMisAZRBtwkG1mI+GQusM5yDaz4XJI3eEcZJtZq7lBtBkH2WbWKm4Qbao6yD7q\nqOQQVJcegTOzDDiD6BAOss1sIA6pzUG2mdXkkNocZJtZU7lBdCAH2WbWDG4QHcxBtpmNhDOILuIg\n26w7OaS2ujjINus+DqmtLg6yzaxebhBdykG2mQ3FDaLLVQbZ993nINvMXucMwnbjINus8ziktqZx\nkG3WWQoXUku6XtIWSesq5l0p6VFJD0v6rqSxFY9dIukxSRskfTCrumxoDrLNDLLNIG4APlQ17y7g\nHRHxLqAXuARA0tHAbODo9DnXSHI+kjMH2WbdLbMv4YhYBWyvmrcyIvrSuw8Ah6XTpwFLI2JnRDwB\nPA5Mz6o2a4yDbLPulOdf6ecBt6fTE4HNFY9tBg5teUU2qKlTYdky+MY34EtfghkzYNWqvKsys6yM\nyWOjkv4GeCUibhpksZp/n/b09OyaLpVKlEqlptZmQyuVYPXqZC9izhwH2WZFUy6XKZfLI15Ppmcx\nSZoMrIiIYyrmfRL4NPCBiPhdOm8+QERckd6/E1gYEQ9Urc9nMRXMjh2weHFySuyZZ8KCBTB+fN5V\nmVmlwp3FVIukDwGfB07rbw6p5cAZkvaWdARwJLC6lbXZ8DjINutcWZ7muhS4H5gqaZOk84CvAvsD\nKyWtkXQNQESsB24F1gN3ABd4V6G9DBRk9/UN/VwzKyYPlLNM9I/IHjUqCbQ9ItssPx5JbYXjEdlm\nxdAWGYR1F4/INmtvbhCWOQfZZu3JDcJaxkG2WXtxBmG5cZBt1hoOqa0tOcg2y55DamtLDrLNissN\nwgrBQbZZ8bhBWKE4yDYrDmcQVmgOss1GziG1dSwH2WYj45DaOpaDbLN8uEFY23CQbdZabhDWdhxk\nm7WGMwhrew6yzQaXeQYh6WBJ3uOwwun/jex585LfyJ41C3p7867KrP3V9YUv6UDgl8Cp2ZZjNjwO\nss2ar949gk8AK4HzM6zFbMQcZJs1T70N4lzgr4BJkg7JsB6zpnCQbTZyQzYISe8FtkbEJuBbwCez\nLsqsWaZOhWXL4BvfSALsGTNg1aq8qzJrD0OexSTpfwP3RMQtkt4C3BsR01pS3Z61+CwmGzaPyLZu\nlclZTJLeAPwxcBtARDwLbJRUGk6RZnlykG3WmKEOMb0CzIiIVyrmnQ38JLuSzLLlINusPoM2iIjY\nGRFb+u9L+mhEvBARL2Rfmlm2HGSbDa6hkdSS1kTEezKsZ6jtO4OwzHhEtnUqX83VbIQ8Ittsd402\niL/IpAqzgnCQbfa6RhvEpzKpwqxgHGSbNd4gTqh3QUnXS9oiaV3FvAMlrZTUK+kuSeMqHrtE0mOS\nNkj6YIN1mWXCQbZ1s0YbxJahF9nlBuBDVfPmAysjYgpwd3ofSUcDs4Gj0+dc4yvHWpF4RLZ1o4a+\nhCOi+gt/sGVXAdurZp8KLEmnlwAfS6dPA5amp9U+ATwOTG+kNrNWcJBt3aTVf6VPqBhXsQWYkE5P\nBDZXLLcZOLSVhZnVy0G2dYsxeW04IkLSYIMaaj7W09Oza7pUKlEqlZpbmFmd+oPs88+Hnp4kyL7o\nIpg7F/bZJ+/qrJuVy2XK5fKI1zPgQDlJB0RE9SGixlYuTQZWRMQx6f0NQCkinkkvG35PRBwlaT5A\nRFyRLncnsDAiHqhanwfKWWFt3AgXXwxr18KiRTB7drK3YZa3LAbK9Up6VNLXJZ0rqRnXvVwOnJNO\nnwMsq5h/hqS9JR0BHAmsbsL2zFqmP8hessRBtnWGQS+1IWkqcFJ6+y/AW4B/A+6PiC8OumJpKTAT\nGE+SNywA/gW4FTgceAL484h4Pl3+UuA84FVgXkR8r8Y6vQdhbcGXFrciGe4eRN3XYpL0duAjwDzg\n0Iho+VFWNwhrNzt2wOLFcNVVcOaZsGABjB+fd1XWbZp+iEnSH0j6vKTvSnoQuBwYTfL71GOHX6pZ\n9/CIbGtng4XUfcAa4Grgtoh4qZWF1eI9CGt3DrItD00/xJSeZdSfPUwH9gJ+TJJB/FtE/GL45Q6P\nG4R1invvhc9+1pcWt9ZoRQaxH0mI/BngiIgY3ejGRsoNwjqJg2xrlSwyiLGS/kTS30m6m2R081nA\nCpLrJpnZCHhEthXdYIeYtpGe0grcBzwUES+3sLZaNXkPwjrW1q1w2WVwyy0ekW3NlfkhpiJwg7Bu\n4CDbmi2LkHoFyfWQaq00IuLURjc2Um4Q1k0cZFuzZNEgtpLkDkuB/msi9W8gIuLe4RQ6Em4Q1m0c\nZFszZHEtpkOAS4F3Al8BTgG2RkQ5j+Zg1o0cZFueBmwQEfFqRNwREWcDM0h+xOdeSX/dsurMDPCI\nbMvHoNGXpH0k/RlwI/BXwGLgtlYUZmZ78m9kWysNlkF8C3gHcDtwS0Ssa2VhtTiDMNudg2yrRxYh\ndR8w0PWXIiLe1OjGRsoNwmxPDrJtKE0PqSNiVES8cYBby5uDmdXmINuy4uE3Zh3CQbY1mxuEWYdx\nkG3N4kttmHU4B9nmazGZ2YAcZHe3LEZSm1mHcJBtw+EGYdZFHGRbI9wgzLqQg2yrhzMIM3OQ3eEc\nUpvZiDjI7lwOqc1sRBxkWzU3CDPbjYNs6+cGYWY1Ocg2ZxBmVhcH2e2rrTIISZdI+pmkdZJukvR7\nkg6UtFJSr6S7JI3LozYzq23mTFi9GubNgzlzYNYs6O3NuyrLUssbhKTJwKeB4yLiGGA0cAYwH1gZ\nEVOAu9P7ZlYgDrK7Sx57EC8AO4H9JI0B9gOeBk4FlqTLLAE+lkNtZlYHB9ndoeUNIiKeA74EPEnS\nGJ6PiJXAhIjYki62BZjQ6trMrDEOsjvbmFZvUNLbgc8Ak4HfAP9X0pzKZSIiJNVMo3t6enZNl0ol\nSqVSVqWaWZ2mToVly14Psq++2kF2nsrlMuVyecTraflZTJJmA6dExKfS+2cBM4D3A++LiGckHQLc\nExFHVT3XZzGZFZxHZBdPO53FtAGYIWlfSQL+K7AeWAGcky5zDrAsh9rMbIQcZHeOPDKIh4FvAg8B\nj6Sz/xm4AjhFUi/J3sQVra7NzJrHQXb780A5M2uJjRvh4oth7VpYtAhmz072Nix7vpqrmbUFj8hu\nPTcIM2sbDrJbq51CajPrcg6y24MbhJnlxkF2sblBmFnuPCK7mJxBmFnhOMhuLofUZtZRHGQ3j0Nq\nM+solUH2iSc6yM6DG4SZFdq++yYD7Bxkt54bhJm1BQfZrecMwszakoPs+jmkNrOu4yC7Pg6pzazr\nOMjOlhuEmbU9B9nZcIMws47hILu5nEGYWcdykJ1wSG1mVoODbIfUZmY1OcgePjcIM+sKDrIb5wZh\nZl3FQXb9nEGYWVfrhiDbIbWZ2TB1epDtkNrMbJgcZNfmBmFmlnKQvTs3CDOzKg6yE84gzMyG0O5B\ntkNqM7MMtXOQ7ZDazCxD3Rhk59IgJI2T9G1Jj0paL+lESQdKWimpV9JdksblUZuZ2WC6KcjOaw9i\nMXB7REwDjgU2APOBlRExBbg7vW9mVkjdEGS3PIOQNBZYExFvq5q/AZgZEVskHQyUI+KoqmWcQZhZ\nIRU5yG6nDOIIYKukGyT9RNLXJb0BmBARW9JltgATcqjNzGxYZs6E1ath3jyYMwdmzYLe3ryrGpk8\nGsQY4Djgmog4DniJqsNJ6W6CdxXMrK10WpA9JodtbgY2R8SD6f1vA5cAz0g6OCKekXQI8GytJ/f0\n9OyaLpVKlEqlbKs1M2tQf5B93nlw2WVJkH3RRTB3LuyzT/bbL5fLlMvlEa8nl3EQkn4AfCoieiX1\nAPulD/06Ir4oaT4wLiLmVz3PGYSZtZ2NG5OGsXYtLFoEs2cnexut0lYD5SS9C7gW2Bv4OXAuMBq4\nFTgceAL484h4vup5bhBm1rbyCrLbqkEMlxuEmbW7PEZkt9NZTGZmXaudgmw3CDOzHLTDiGw3CDOz\nHBV5RLYzCDOzAskiyHZIbWbWIZodZDukNjPrEEUJst0gzMwKKu8g2w3CzKzg8gqynUGYmbWZRoNs\nh9RmZl2kkSDbIbWZWRdpRZDtBmFm1sayDLLdIMzMOkAWQbYzCDOzDlQZZD/4oENqMzOr0B9kz5nj\nBmFmZjX4LCYzM2sqNwgzM6vJDcLMzGpygzAzs5rcIMzMrCY3CDMzq8kNwszManKDMDOzmtwgzMys\nJjcIMzOryQ3CzMxqcoMwM7OacmsQkkZLWiNpRXr/QEkrJfVKukvSuLxqMzOzfPcg5gHrgf7Ls84H\nVkbEFODu9L5lqFwu511CR/H72Vx+P/OXS4OQdBjwYeBaoP8StKcCS9LpJcDHciitq/gfYHP5/Wwu\nv5/5y2sP4mrg80Dlj+FNiIgt6fQWYELLqzIzs11a3iAkfRR4NiLW8Prew27SXwXyLwOZmeWo5b8o\nJ+ly4CzgVWAf4E3Ad4ETgFJEPCPpEOCeiDiq6rluGmZmw9B2PzkqaSbwuYj4U0n/C/h1RHxR0nxg\nXEQ4qDYzy0kRxkH0d6grgFMk9QLvT++bmVlOct2DMDOz4irCHsQeJH1I0gZJj0m6eIBl/jF9/GFJ\n72l1je1kqPdTUknSb9KBi2skfSGPOtuBpOslbZG0bpBl/Nms01Dvpz+b9ZM0SdI9kn4m6aeSLhxg\nufo/nxFRqBswGngcmAzsBawFplUt82Hg9nT6ROBHeddd1Fud72cJWJ53re1wA04G3gOsG+Bxfzab\n+376s1n/e3kw8O50en9g40i/O4u4BzEdeDwinoiIncDNwGlVy+waVBcRDwDjJHncRG31vJ8wwCnH\ntruIWAVsH2QRfzYbUMf7Cf5s1iUinomIten0b4FHgYlVizX0+SxigzgU2FRxf3M6b6hlDsu4rnZV\nz/sZwEnpLuftko5uWXWdx5/N5vJncxgkTSbZM3ug6qGGPp9jml1YE9Sbmlf/VeG0vbZ63pefAJMi\n4mVJfwIsA6ZkW1ZH82ezefzZbJCk/YFvA/PSPYk9Fqm6P+Dns4h7EE8BkyruTyLpcoMtc1g6z/Y0\n5PsZES9GxMvp9B3AXpIObF2JHcWfzSbyZ7MxkvYCvgPcGBHLaizS0OeziA3iIeBISZMl7Q3MBpZX\nLbMcOBtA0gzg+Xj9Ok62uyHfT0kTJCmdnk5y+vNzrS+1I/iz2UT+bNYvfZ+uA9ZHxFcGWKyhz2fh\nDjFFxKuS/hr4HskZONdFxKOS/iJ9/P9ExO2SPizpceAl4NwcSy60et5P4HTgLyW9CrwMnJFbwQUn\naSkwExgvaROwkOTsMH82h2Go9xN/NhvxB8Ac4BFJa9J5lwKHw/A+nx4oZ2ZmNRXxEJOZmRWAG4SZ\nmdXkBmFmZjW5QZiZWU1uEGZmVpMbhJmZ1eQGYV1J0kEVl5D+d0mb0+kXJX0t7/rMisDjIKzrSVoI\nvBgRX85p+/sBO9Or7ZoVhvcgzBL9l3MoSVqRTvdIWiLpB5KekDRL0lWSHpF0h6Qx6XLHSypLekjS\nnZIObnDbU4GNkq6UdFRzX5bZ8LlBmA3uCOB9JNfRvxFYGRHHAjuAj6QXR/sq8GcR8V7gBuAfGtlA\nRKwBjgU2ANdKWiXpk5Le0MTXYdawwl2LyaxAArgjIl6T9FNgVER8L31sHcmv9E0B3gF8P72m3Gjg\n6YY3lFyW+TrgOknT0unFwNiRvgiz4XKDMBvcKwAR0SepMiPoI/n3I+BnEXHSQCuQNBr4MUnDWQ6s\nIbkoHcD5EfGTdLnJwDkkF6RbCyxo5gsxa5QbhNnA6vmpy43AmyXNiIgfpYecjoyI9f0LRMRrwLur\nnrfrWv1pY7gWOAi4HjgpIob6GU6zzLlBmCWi4r+1pmHPX96KiNgp6XTgHyWNJfk3dTWwnvq9CsyP\niIcaL9ssOz7N1czMavJZTGZmVpMbhJmZ1eQGYWZmNblBmJlZTW4QZmZWkxuEmZnV5AZhZmY1uUGY\nmVlN/wkVxRcaayYErQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7fbaec2a1190>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "number of units out 0 ,capacity out 0MW ,capacity available 200MW ,probability 0.814506 ,tg in percentage 0.00 ,expected load 0.000000MW\n",
      "\n",
      "number of units out 1 ,capacity out 50MW ,capacity available 150MW ,probability 0.171475 ,tg in percentage 0.00 ,expected load 0.000000MW\n",
      "\n",
      "number of units out 2 ,capacity out 100MW ,capacity available 100MW ,probability 0.013538 ,tg in percentage 200.00 ,expected load 2.707500MW\n",
      "\n",
      "\n",
      "expected loss of load is 2.707500MW percent of time. assuming 365 days in a year, then expected loss of load is 9.882MW days per year\n"
     ]
    }
   ],
   "source": [
    "from __future__ import division\n",
    "from math import factorial\n",
    "from numpy import arange\n",
    "%matplotlib inline\n",
    "from matplotlib import pyplot as plt\n",
    "#given\n",
    "n=4 #number of generating station\n",
    "f=0.05 #F.O.R\n",
    "a=1-f\n",
    "p=50 #generation station power \n",
    "mp=150 #maximum alowable  power\n",
    "lf=50 #load factor in percentage\n",
    "def comb(m,r):\n",
    "    y=factorial(m)/(factorial(m-r)*factorial(r))\n",
    "    return y\n",
    "pg = range(0,n)\n",
    "co = range(0,n)\n",
    "ca = range(0,n)\n",
    "for i in range(0,n):\n",
    "    pg[i]=comb(n,i)*((f)**i)*((a)**(n-i))\n",
    "    co[i]=p*i ;ca[i]=p*(n-i)\n",
    "    print \"\\nnumber of units out %d ,capacity out %dMW ,capacity available %dMW ,probability %4f \"%(i,co[i],ca[i],pg[i])\n",
    "\n",
    "ld=arange(mp,0,-lf)\n",
    "#[m n]=size(ld)\n",
    "n =len(ld)\n",
    "plt.plot(ld)\n",
    "plt.title('Load Duration Curve')\n",
    "plt.xlabel('Time -->')\n",
    "plt.ylabel('MW -->')\n",
    "plt.show()\n",
    "tg = range(0,n+1)\n",
    "tg[n-2]=round(10000/(n-2))/100\n",
    "tg[n-1]=tg[n-2]*2\n",
    "tg[n]=100\n",
    "tg[1]=0 ;tg[0]=0 #maximum load limit\n",
    "el = range(0,n)\n",
    "for i in range(0,n):\n",
    "    el[i]=pg[i]*tg[i]\n",
    "    print \"\\nnumber of units out %d ,capacity out %dMW ,capacity available %dMW ,probability %4f ,tg in percentage %.2f ,expected load %.6fMW\"%(i,co[i],ca[i],pg[i],tg[i],el[i])\n",
    "\n",
    "lt=sum(el)\n",
    "print \"\\n\\nexpected loss of load is %.6fMW percent of time. assuming 365 days in a year, then expected loss of load is %.3fMW days per year\"%(lt,lt*365/100)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## example 17.4 Page 425"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "number of units out 0 ,capacity out 0MW ,capacity available 200MW ,probability 0.9223682 \n",
      "\n",
      "number of units out 1 ,capacity out 50MW ,capacity available 150MW ,probability 0.0752954 \n",
      "\n",
      "number of units out 2 ,capacity out 100MW ,capacity available 100MW ,probability 0.0023050 \n",
      "\n",
      "number of units out 3 ,capacity out 150MW ,capacity available 50MW ,probability 0.0000314 \n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAEZCAYAAACNebLAAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAHE1JREFUeJzt3X20HHWd5/H3RwKjBHkIWZE8YDIa0PDkEwEZkXZYWUad\nRI07RHcvQnTWs8wCa5xRiLPkuo4g6grMOOTsKOEgHq5CUEhcntWWrDy4KMRgJkMiRBLAQEjQkMEh\nkO/+UXWTTtP33u57ux66+/M6555UV1VXfbtyb/36U79fdSsiMDMzq/eKogswM7NycgNhZmYNuYEw\nM7OG3ECYmVlDbiDMzKwhNxBmZtaQGwjrWJL6JV1TdB3DkXSYpG2SVHQtZq1yA2GZkLRe0ikZ72bI\nm3gkVSTtTE/O2yRtkPRdSW/PsqD0df/prgIjHouIV0cGNxwpca6kVZKeS1/jdZKOave+rDe5gbCs\nBMOcwHPyeHpyfjVwArAGWFF7Am+FpL2aWC2AvNLC5cC5wDnAQcDhwI3A+1rdkKRx7S3NuoEbCMuV\npD+SdJmkx9OfSyXtky47UNIPJD0laYuk5ZIm1zx3uqSfSPq9pNuBic3uNyIej4hFwDeBS9LtTUtT\nxq6/A0lVSR9Pp8+U9FNJX5O0GVgk6Y8l/UjSZklPS/q2pAPS9a8BDgOWp6nlr+v3IWmSpGWSnpG0\nVtInavbdnyaAq9PX+JCktw1xHGcAZwPzIqIaETsi4vmIuDYiBl/frtdS83pW1DzeKelsSQ8DD0u6\nQtJX6vZzk6RP1dR+Q/r/84ikc5o9/taZ3EBY3j4HzAKOTX9mAX+bLnsFcCXJSfYw4Hng6zXPvRb4\nf8DBwBeAj9F6Svk+8FZJrxpieX3ymQX8GngNcBFJOvgicCjwJmAq0A8QEX3AY8D70+Ty1Qbb/066\nzqHAh4GLJL27ZvmfAwPAAcAy9nz9tU4BNkTE/cO81mZS3Jz0Nb4p3e/pgwskHQS8BxhIG7jlwAPA\npHT//13SqSNs3zqYGwjL20eB/xkRmyNiM/B5oA8gIrZExPcj4g8R8RzJCflkSDp7gbcD/yN9t7yC\n5ITV6uWcJ9LnHNjs+hHxjxGxM63r1xHxw7SGzcClgzWORNJU4ETgsxHxQkSsJEk0Z9SstiIibk37\nLL5N0og2cjDw2yZfw3AujohnI+LfgP8LhKST0mUfBu6OiN8CxwETI+LvIuLFiHg0rX1eG2qwkvJ1\nR8vbJOA3NY8fS+chaV+SE+5/ILmmDrBfOgJoErA1Ip6vee5vSN7Bt2IyybvqZ4FDmlh/Q+0DSYeQ\nXPt/J/BqkjdZW5rc9yRgS0Rsr5n3GEnDN2hTzfS/Aq+U9IqI2Fm3rWdIUshY7Xp9ERGSvgN8BFhB\n0ph/K138OmCSpK01z90LuKsNNVhJOUFY3p4AptU8Pgx4PJ3+NElH66yIOIDknbnSnyeBg9JGZNDr\naP0S0weBn6cNzeCJunabr61bv377FwEvAUelNfax59/RcPU8AUyQtF/NvMOAjU3WXuuHwJSh+ihS\n24HxNY/rXxu8vN4B4MOSXkdy6emGdP5jwKMRcVDNz/4R8f5R1G4dwg2EZWkfSa+s+RlHcgL6W0kT\nJU0ELiS5lAKwH0m/w+8kTQAWDW4oIn4D3A98XtLekt4JNHVySoeDTpa0CPg4sDDd5tMkjVOfpL0k\nzQdeP8Lm9iM58f4+7UD/m7rlm4baRkRsAO4GLk47648B5te8/qZFxFrgCpL+gZMlDR7reZI+m672\nIPAhSa+S9AaS1z7Sdh8ENpNcPro1In6fLvoZsE3SZ9Lt7SXpKGU8bNiK5QbCsnQzyWWSwZ8Lgb8j\nOdH/Mv25P50HcBnwKpIT1N3ALez5DvejwPEkl3QuBK4eZt9BcklkG7CN5AR3JHByRNxZs95fkpzk\nNwMzgZ/WbaP+HfbngbcCvyPpA7mhbp2LSRrArZIW1Gxn0EdIEtQTwPeACyPiR8Psb8hEEhHnknRi\n/yOwFVhH0um8LF3lUuAFkkbrKpKGqHZ7Q237WuBP038H97WTpEF+M/AI8DTwT8D+Q9VnnU9ZfWGQ\npCUk47Gfioij03n9wCdIfrkAFkbELemyC0jeTb0EnBsRt2dSmJmZNSXLBuIk4DngWzUNxCJgW0R8\nrW7dmSTvVo4j6US8Ezi8QcecmZnlJLNLTOkwxK0NFjUaljgHGEiHDq4nicqzsqrNzMxGVkQfxDmS\nVkq6UtLgWPRJ7DmSYyNJkjAzs4Lk3UAsBqaTdHQ9CfyvYdYt+nN8zMx6Wq43ykXEU4PTkr5JMgoE\nkqGGtTc8TWH32HhqnuNGw8xsFCKi5Q+RzDVBSKq98/ODwKp0ehkwLx3LPR2YQTIs8WWeey4499xg\n0qTgppuCiN78WbRoUeE1lOXHx8LHwsdi+J/RyixBSBoguRN2oqQNJDc9VSS9meTy0aPAJwEiYrWk\n64DVwIvA2THEqxo/Hi6/HObOhfnz4frrk8cTJmT1SszMelOWo5g+EhGTImKfiJgaEUsi4oyIOCYi\njo2ID0TEppr1L4qIN0TEGyPitpG2/653wcqVScNw9NGwbNlIzzAzs1Z09J3Ug2liYAAWLIC+PtjS\n7MemdbhKpVJ0CaXhY7Gbj8VuPhZjl9mNclmQNNSVJ7Zvh4ULYelSWLwYZs/OuTgzs5KSRIyik7pr\nGohBd92V9E284x3umzAzg9E3EB19iakR902YmbVH1yWIWk4TZmZOEA05TZiZjV5XJ4haThNm1quc\nIEbgNGFm1pqeSRC1nCbMrJc4QbTAacLMbGQ9mSBqOU2YWbdzghglpwkzs8Z6PkHUcpows27kBNEG\nThNmZrs5QQzBacLMuoUTRJs5TZhZr3OCaILThJl1MieIDDlNmFkvcoJokdOEmXUaJ4icOE2YWa9w\nghgDpwkz6wROEAVwmjCzbuYE0SZOE2ZWVk4QBXOaMLNu4wSRAacJMysTJ4gScZows27gBJExpwkz\nK5oTREk5TZhZp8qsgZC0RNImSasaLPu0pJ2SJtTMu0DSWklrJJ2aVV1FGD8+SQ8DA7BgAfT1wZYt\nRVdlZja8LBPEVcBp9TMlTQXeA/ymZt5M4HRgZvqcKyR1XbpxmjCzTpLZSTgiVgBbGyz6GvCZunlz\ngIGI2BER64F1wKysaiuS04SZdYpc36VLmgNsjIhf1i2aBGysebwRmJxbYQVwmjCzssutgZC0L7AQ\nWFQ7e5indNZwpVFwmjCzMhuX475eD0wDVkoCmAL8XNLxwOPA1Jp1p6TzXqa/v3/XdKVSoVKpZFJs\nngbTxMKFSZpYvBhmzy66KjPrVNVqlWq1OubtZHofhKRpwPKIOLrBskeBt0XElrST+lqSfofJwJ3A\nG+pveujE+yBa5fsmzKzdSncfhKQB4G7gcEkbJJ1Vt8quM31ErAauA1YDtwBnd31LMAT3TZhZWfhO\n6hJzmjCzdihdgrCxc5owsyI5QXQIpwkzGy0niC7nNGFmeXOC6EBOE2bWCieIHuI0YWZ5cILocE4T\nZjYSJ4ge5TRhZllxgugiThNm1ogThDlNmFlbOUF0KacJMxvkBGF7cJows7FygugBThNmvc0Jwobk\nNGFmo+EE0WOcJsx6jxOENcVpwsya5QTRw5wmzHqDE4S1zGnCzIbjBGGA04RZN3OCsDFxmjCzek4Q\n9jJOE2bdxQnC2sZpwszACcJG4DRh1vmcICwTThNmvcsJwprmNGHWmZwgLHNOE2a9xQnCRsVpwqxz\nOEFYrpwmzLqfE4SNmdOEWbk5QVhhnCbMulNmDYSkJZI2SVpVM+8LklZKelDSDyVNrVl2gaS1ktZI\nOjWruiwb48cn6WFgABYsgL4+2LKl6KrMbCyyTBBXAafVzftyRBwbEW8GbgQWAUiaCZwOzEyfc4Uk\np5sO5DRh1j0yOwlHxApga928bTUP9wM2p9NzgIGI2BER64F1wKysarNsOU2YdYfc36VL+qKkx4Az\ngYvT2ZOAjTWrbQQm51yatZnThFlnG5f3DiPic8DnJJ0PXAacNdSqjWb29/fvmq5UKlQqlTZXaO00\nmCbmzk1GOl1/vUc6mWWtWq1SrVbHvJ1Mh7lKmgYsj4ijGyw7DLg5Io5KGwsi4kvpsluBRRFxX91z\nPMy1g23fDgsXwtKlsHgxzJ5ddEVmvaEjhrlKmlHzcA7wQDq9DJgnaR9J04EZwM/yrM2y574Js86S\n5TDXAeBu4AhJGyTNBy6WtErSg0AF+DRARKwGrgNWA7cAZzsqdC/3TZh1Bt9JbYXyXdhm2euIS0xm\n9ZwmzMrLCcJKw2nCLBtOENbxnCbMysUJwkrJacKsfZwgrKs4TZgVzwnCSs9pwmxsnCCsazlNmBXD\nCcI6itOEWeucIKwnOE2Y5ccJwjqW04RZc5wgrOc4TZhlywnCuoLThNnQnCCspzlNmLWfE4R1HacJ\nsz05QZilnCbM2qPpBCHptcBTEbEz25KGrcEJwlriNGGWcYKQNAF4FPC3CFtHcZowG72mEoSkc4D3\npOv/eeZVDV2HE4SNmtOE9aqs+yDOAv4KmCrp0FZ3YlYGThNmrRmxgZD0duDpiNgAXAOcmXVRZlkZ\nPz5JDwMDsGAB9PXBli1FV2VWTs0kiE8AS9Lpa4AzsivHLB9OE2YjG7YPQtJ44CHgiIh4IZ13I3BZ\nRFRzqXDPetwHYW3nvgnrdln1QbwAnDDYOKTOAH7R6o7Mysppwqyxlu6klvT+iPhBhvWMtH8nCMuU\n04R1o7zupP5Cqzsw6yROE2a7tZogHoiIt2RYz0j7d4Kw3DhNWLfIK0F8stUdmHUqpwnrda0miH+K\niP+SYT0j7d8JwgrhNGGdLK8EcVyzK0paImmTpFU1874i6Z8lrZT0PUkH1Cy7QNJaSWskndpiXWaZ\ncpqwXtRqA7GphXWvAk6rm3c7cGREHAs8DFwAIGkmcDowM33OFZL8UeRWKr4L23pNSyfhiKg/4Q+3\n7gpga928O2o+Lvw+YEo6PQcYiIgdEbEeWAfMaqU2s7w4TVivKPJd+nzg5nR6ErCxZtlGYHLuFZk1\nyWnCesG4InYq6XPACxFx7TCrNeyN7u/v3zVdqVSoVCptrc2sFYNpYuHCJE0sXgyz/a0pVrBqtUq1\nWh3zdoYcxSTpoIjY2nBhsxuXpgHLI+LomnlnAn8JnBIRf0jnnQ8QEV9KH98KLIqI++q251FMVloe\n6WRllcUopofTEUffkHSWpMPHUB8Akk4D/gaYM9g4pJYB8yTtI2k6MAP42Vj3Z5Yn901Ytxnp01yP\nAE5Mf94BvAa4B7g7Ii4ZdsPSAHAyMJFk9NMiklFL+wCDV2vviYiz0/UXkvRLvAicFxG3NdimE4R1\nBKcJK5PRJoimb5ST9HrgfcB5wOSIeGWrOxsrNxDWSbZvT/omli5134QVq+0NhKQ/YXdymAo8AtxL\nkiAeiIh/G325o+MGwjqR04QVLYs+iBXAPOB7QCUiTo+ISyPi3iIaB7NO5b4J61TDJYhD2Z0gZgF7\nAz8nSRD3RMQjeRVZU5MThHU0pwkrQtsTREQ8GRE3RMRfR8S7gFOANcDngbWjL9WsdzlNWCcZLkEc\nwO4RTCcCbyFpGO4GfhoRS/MqsqYmJwjrGk4Tlpcs+iB+DZwN/CtJapgSEcdHxKeKaBzMuo3ThJVd\nS98HUTQnCOtWThOWpSyGuS4n+TykRhuNiMh9VLcbCOtmvm/CspJFA/E0yaeqDpB8NDfsbiwiIn4y\nmkLHwg2E9QKnCWu3LPogDgUWAkcBlwHvAZ6OiGoRjYNZr3DfhJVFU30Qkv4I+AjwVaA/Ir6edWFD\n1OEEYT3FacLaIZPvpJb0SklzgW8DfwVcDnx/dCWaWaucJqxIw/VBXAMcSfKtb9+NiFV5FtaIE4T1\nMqcJG60sOql3AtuHeF5ExP6t7mys3EBYr/NIJxuNzD/uuwzcQJglnCasFZn0QZhZOblvwvLgBGHW\n4ZwmbCROEGY9ymnCsuIEYdZFnCasEScIM3OasLZygjDrUk4TNsgJwsz24DRhY+UEYdYDnCZ6mxOE\nmQ3JacJGwwnCrMc4TfQeJwgza4rThDXLCcKshzlN9AYnCDNrmdOEDSezBkLSEkmbJK2qmfcfJf1K\n0kuS3lq3/gWS1kpaI+nUrOoysz2NH5+kh4EBWLAA+vpgy5aiq7IyyDJBXAWcVjdvFfBB4K7amZJm\nAqcDM9PnXCHJ6cYsR04TVi+zk3BErAC21s1bExEPN1h9DjAQETsiYj2wDpiVVW1m1pjThNUqy7v0\nScDGmscbgckF1WLW85wmDGBc0QUMo+Fwpf7+/l3TlUqFSqWSUzlmvWUwTcydm4x0uv56j3TqFNVq\nlWq1OubtZDrMVdI0YHlEHF03/8fApyPiF+nj8wEi4kvp41uBRRFxX93zPMzVrAD+LuzO1onDXGuL\nXQbMk7SPpOnADOBnxZRlZvXcN9GbshzmOgDcDRwhaYOk+ZI+IGkDcALwfyTdAhARq4HrgNXALcDZ\njgpm5eO+id7iO6nNbFR8F3bn6MRLTGbWwZwmup8ThJmNmdNEuTlBmFlhnCa6kxOEmbWV00T5OEGY\nWSk4TXQPJwgzy4zTRDk4QZhZ6ThNdDYnCDPLhdNEcZwgzKzUnCY6jxOEmeXOaSJfThBm1jGcJjqD\nE4SZFcppIntOEGbWkZwmyssJwsxKw2kiG04QZtbxnCbKxQnCzErJaaJ9nCDMrKs4TRTPCcLMSs9p\nYmycIMysazlNFMMJwsw6itNE65wgzKwnOE3kxwnCzDqW00RznCDMrOc4TWTLCcLMuoLTxNCcIMys\npzlNtJ8ThJl1HaeJPTlBmJmlnCbawwnCzLqa00QJE4SkJZI2SVpVM2+CpDskPSzpdkkH1iy7QNJa\nSWsknZpVXWbWW5wmRi+zBCHpJOA54FsRcXQ678vA5oj4sqTPAgdFxPmSZgLXAscBk4E7gcMjYmfd\nNp0gzGzUejVNlC5BRMQKYGvd7NnA1en01cAH0uk5wEBE7IiI9cA6YFZWtZlZb3KaaE3endSHRMSm\ndHoTcEg6PQnYWLPeRpIkYWbWVuPHJ+lhYAAWLIC+PtiypeiqyqmwUUzptaLhrhf5WpKZZcZpYmTj\nct7fJkmvjYjfSjoUeCqd/zgwtWa9Kem8l+nv7981XalUqFQq2VRqZl1vME3MnZv0TVx/fXf0TVSr\nVarV6pi3k+kwV0nTgOV1ndTPRMQlks4HDqzrpJ7F7k7qN9T3SLuT2syysn07LFwIS5fC4sUwe3bR\nFbXPaDupsxzFNACcDEwk6W+4ELgJuA44DFgP/EVEPJuuvxCYD7wInBcRtzXYphsIM8tUN450Kl0D\nkQU3EGaWh25LE24gzMzarFvSROnugzAz63S9PtLJCcLMrAmdnCacIMzMMtSLacIJwsysRZ2WJpwg\nzMxy0itpwgnCzGwMOiFNOEGYmRWgm9OEE4SZWZuUNU04QZiZFazb0oQThJlZBsqUJpwgzMxKpBvS\nhBOEmVnGik4TThBmZiXVqWnCCcLMLEdFpAknCDOzDtBJacIJwsysIHmlCScIM7MOU/Y04QRhZlYC\nWaYJJwgzsw5WxjThBGFmVjLtThNOEGZmXaIsacIJwsysxNqRJpwgzMy6UJFpwgnCzKxDjDZNOEGY\nmXW5vNOEE4SZWQdqJU04QZiZ9ZA80kQhCULSecAnAAHfiIjLJU0Avgu8DlgP/EVEPFv3PCcIM7M6\nI6WJjkkQko4iaRyOA44F3i/p9cD5wB0RcTjww/SxDaFarRZdQmn4WOzmY7FbLx2LrNJEEZeY3gjc\nFxF/iIiXgJ8Ac4HZwNXpOlcDHyigto7RS7/8I/Gx2M3HYrdeOxbjxyfpYWAAFiyAvj7YsmVs2yyi\ngXgIOEnSBEn7Au8FpgCHRMSmdJ1NwCEF1GZm1tHamSbGta+s5kTEGkmXALcD24EHgZfq1glJ7mww\nMxuFwTQxd27SNzFahQ9zlfRFYCNwHlCJiN9KOhT4cUS8sW5dNxpmZqMwmk7q3BMEgKTXRMRTkg4D\nPgScAEwHPgZckv57Y/3zRvMCzcxsdIoa5noXcDCwA/hURPw4HeZ6HXAYQwxzNTOz/BR+icnMzMqp\nlHdSSzpN0hpJayV9doh1/j5dvlLSW/KuMS8jHQtJ/yk9Br+U9FNJxxRRZx6a+b1I1ztO0ouSPpRn\nfXlq8m+kIukBSQ9JquZcYm6a+BuZKOlWSQ+mx+LMAsrMnKQlkjZJWjXMOq2dNyOiVD/AXsA6YBqw\nN8kopzfVrfNe4OZ0+njg3qLrLvBYvAM4IJ0+rZePRc16PwJ+AMwtuu4Cfy8OBH4FTEkfTyy67gKP\nRT9w8eBxAJ4BxhVdewbH4iTgLcCqIZa3fN4sY4KYBayLiPURsQP4DjCnbp1dN9VFxH3AgZK68b6J\nEY9FRNwTEb9LH95Hck9JN2rm9wLgHGAp8HSexeWsmWPxUeCGiNgIEBGbc64xL80ciyeB/dPp/YFn\nIuLFHGvMRUSsALYOs0rL580yNhCTgQ01jzem80ZapxtPjM0ci1ofB27OtKLijHgsJE0mOTksTmd1\nawdbM78XM4AJkn4s6X5JfblVl69mjsU3gCMlPQGsJBlS34taPm8WMsx1BM3+UdcPee3Gk0HTr0nS\nu4H5wJ9kV06hmjkWlwHnR0RIEi//HekWzRyLvYG3AqcA+wL3SLo3ItZmWln+mjkWC4EHI6KSfu7b\nHZKOjYhtGddWRi2dN8vYQDwOTK15PJWkpRtunSnpvG7TzLEg7Zj+BnBaRAwXMTtZM8fibcB3kraB\nicCfSdoREQV95XtmmjkWG4DNEfE88Hw6tPxYoNsaiGaOxYnAFwEi4teSHgWOAO7PpcLyaPm8WcZL\nTPcDMyRNk7QPcDpQ/we+DDgDQNIJwLOx+3OcusmIxyK92fB7wH+OiHUF1JiXEY9FRPxxREyPiOkk\n/RD/tQsbB2jub+Qm4J2S9ko/8+x4YHXOdeahmWOxBvj3AOk19yOAR3KtshxaPm+WLkFExIuS/htw\nG8kIhSsj4p8lfTJd/r8j4mZJ75W0juTznM4qsOTMNHMsgAuBg4DF6TvnHRExq6ias9LksegJTf6N\nrJF0K/BLYCfJ9650XQPR5O/FRcBVklaSvCn+TESM8XNOy0fSAHAyMFHSBmARyaXGUZ83faOcmZk1\nVMZLTGZmVgJuIMzMrCE3EGZm1pAbCDMza8gNhJmZNeQGwszMGnIDYT1J0sHpR2E/IOlJSRvT6W2S\nvl50fWZl4PsgrOdJWgRsi4ivFbT/fUlucNxRxP7NhuIEYZYQ7PqSneXpdL+kqyXdJWm9pA9J+mr6\n5Uy3SBqXrvc2SdX0U1NvlfTaFvd9BPAvkr4i6Y3tfVlmo+cGwmx404F3k3yW/reBOyLiGOB54H2S\n9gb+geTLid4OXEX6wXDNiogHgGNIPjPom5JWSDpT0vg2vg6zlpXus5jMSiSAWyLiJUkPAa+IiNvS\nZatIvsXscOBI4M70s7D2Ap5oeUcRzwFXAldKelM6fTlwwFhfhNlouYEwG94LABGxU1JtH8FOkr8f\nAb+KiBOH2oCkvYCfkzQ4y4AHSD5IDeDjEfGLdL1pwMeAeSRfnXlhO1+IWavcQJgNrZkvHPoX4N9J\nOiEi7k0vOc2o/eTUiHgJeHPd827ctZOkYfgmcDCwBDixi7/XwzqIGwizRNT822gaXv7tWxEROyR9\nGPh7SQeQ/E1dSmvfvfAiyTfh9doX2FjJeZirmZk15FFMZmbWkBsIMzNryA2EmZk15AbCzMwacgNh\nZmYNuYEwM7OG3ECYmVlDbiDMzKyh/w8YEJUAVpFFJgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7fbaec2a1310>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "\n",
      "number of units out -1 ,capacity out 0MW ,capacity available 200MW ,probability 0.922368 ,tg in percentage 0.00 ,expected load 0.000000MW\n",
      "\n",
      "number of units out 0 ,capacity out 50MW ,capacity available 150MW ,probability 0.075295 ,tg in percentage 0.00 ,expected load 0.000000MW\n",
      "\n",
      "number of units out 1 ,capacity out 100MW ,capacity available 100MW ,probability 0.002305 ,tg in percentage 41.67 ,expected load 0.096040MW\n",
      "\n",
      "number of units out 2 ,capacity out 150MW ,capacity available 50MW ,probability 0.000031 ,tg in percentage 83.33 ,expected load 0.002613MW\n",
      "\n",
      "\n",
      "expected loss of load is 0.098653MW percent of time. assuming 365 days in a year, then expected loss of load is 0.360MW days per year ,some times the loss of load is also expressed as reciprocal of this figure and then the units are years per day this result is 2.7771MW years per day.\n"
     ]
    }
   ],
   "source": [
    "from math import factorial\n",
    "from numpy import arange\n",
    "%matplotlib inline\n",
    "from matplotlib import pyplot as plt\n",
    "\n",
    "#given\n",
    "n=4 #number of generating station\n",
    "f=0.02 #F.O.R\n",
    "a=1-f\n",
    "p=50 #generation station power \n",
    "mp=150 #maximum alowable  power\n",
    "minp=30 #minimum power \n",
    "lf=60 #load factor in percentage\n",
    "def comb(m,r):\n",
    "    y=factorial(m)/(factorial(m-r)*factorial(r))\n",
    "    return y\n",
    "\n",
    "pg = range(0,n)\n",
    "co = range(0,n)\n",
    "ca = range(0,n)\n",
    "for i in range(0,n):\n",
    "    pg[i]=comb(n,i)*((f)**i)*((a)**(n-i))\n",
    "    co[i]=p*i ;ca[i]=p*(n-i)\n",
    "    print \"\\nnumber of units out %d ,capacity out %dMW ,capacity available %dMW ,probability %.7f \"%(i,co[i],ca[i],pg[i])\n",
    "\n",
    "\n",
    "ld=arange(mp,minp,-lf)\n",
    "#[m n1]=size(ld)\n",
    "n1=len(ld)\n",
    "#[mm m]=max(co)\n",
    "m=max(co)\n",
    "\n",
    "plt.plot(ld)\n",
    "plt.title('Load Duration Curve')\n",
    "plt.xlabel('Time -->')\n",
    "plt.ylabel('MW -->')\n",
    "plt.show()\n",
    "\n",
    "tg = range(0,n)\n",
    "tg[0]=0\n",
    "for i in range(1,n):\n",
    "    \n",
    "    tg[i]=(mp-ca[i])*100/(2*lf) #percentage time \n",
    "\n",
    "print \"\"\n",
    "el = range(0,n)\n",
    "for i in range(0,n):\n",
    "    el[i]=pg[i]*tg[i]\n",
    "    print \"\\nnumber of units out %d ,capacity out %dMW ,capacity available %dMW ,probability %4f ,tg in percentage %.2f ,expected load %.6fMW\"%(i-1,co[i],ca[i],pg[i],tg[i],el[i])\n",
    "\n",
    "lt=sum(el)\n",
    "print \"\\n\\nexpected loss of load is %.6fMW percent of time. assuming 365 days in a year, then expected loss of load is %.3fMW days per year ,some times the loss of load is also expressed as reciprocal of this figure and then the units are years per day this result is %.4fMW years per day.\"%(lt,lt*365/100,100/(lt*365))"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}