proof of an implicit solver

diode-model-schematic.png

explicit-pv-model.ipynb

PV-system-implicit-solver.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Example of implicit approach to solving for current of a PV system given voltage\n",
    "\n",
    "This example is part of a series of notebooks that explore using an implicit approach to solving a PV system versus an explicit approach. Specifically this example finds the current given a specified system voltage.\n",
    "\n",
    "First set up the notebook to make plots inline, then import NumPy to vector math, and the solvers we're going to use from SciPy optimize package."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# imports\n",
    "%matplotlib inline\n",
    "\n",
    "import numpy as np\n",
    "from scipy.optimize import minimize_scalar, fsolve\n",
    "# NOTE: fsolve is the same as root() set to 'hybr'\n",
    "from matplotlib import pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# CONSTANTS\n",
    "IL = 7.0  # [A] photogenerated current\n",
    "I0 = 1.0e-6  # [A] dark current\n",
    "RSH = 20.0  # [ohms] shunt resistance\n",
    "RS = 0.001  # [ohms] series resistance\n",
    "GAMMA = 1.23  # diode ideality\n",
    "VT = 0.026  # [V] thermal voltage at 300[K]\n",
    "VBYPASS = -0.5  # bypass diode trigger voltage\n",
    "VSYS = 240.0  # [V] system voltage\n",
    "NSERIES_CELLS = 72  # number of series cells per module\n",
    "NMODS = 8  # number of modules per string\n",
    "NSTRINGS = 3  # number of strings in system\n",
    "NBYPASS = 3  # number of bypass diodes per module, assumes cells divided evenly"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Callback\n",
    "\n",
    "The SciPy solver takes a callback function which takes a vector, `x`, of solution guesses, and a tuple of extra arguments, `args`, to pass in the independent variables that define the conditions such as number of series cells per module, number of modules per string, number of strings, how many bypass diodes there are per module, and the single diode model parameters."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def isys_from_voltage(x, *args):\n",
    "    nseries_cells = args[0]  # number of series cells per module\n",
    "    nmods = args[1]  # number of modules per string\n",
    "    nstrings = args[2]  # number of strings\n",
    "    nbypass = args[3]  # number of bypass diodes per module, assumes evenly divided\n",
    "    # single diode parameters\n",
    "    il, i0, rsh, rs, gamma, vt, vbypass, vsys = args[4:]\n",
    "\n",
    "    # cells\n",
    "    idx0 = nseries_cells * nmods * nstrings\n",
    "    vcells = x[:idx0].reshape(-1, nstrings)  # reshape strings into columns\n",
    "\n",
    "    # submodules\n",
    "    idx1 = idx0 + nbypass * nmods * nstrings\n",
    "    vsubs = x[idx0:idx1].reshape(-1, nstrings)  # reshape strings into columns\n",
    "\n",
    "    # strings\n",
    "    idx2 = idx1 + nstrings\n",
    "    istrings = x[idx1:idx2]\n",
    "\n",
    "    # system\n",
    "    isys = x[-1]\n",
    "\n",
    "    # residuals\n",
    "    residuals = np.zeros(x.shape)\n",
    "\n",
    "    # cells constraints\n",
    "    vdiode = vcells + istrings*rs\n",
    "    # this is the \"single diode model\" of a solar cell\n",
    "    residuals[:idx0] = np.ravel(\n",
    "        il - i0*(np.exp(vdiode/gamma/vt) - 1.0) - vdiode/rsh - istrings\n",
    "    )\n",
    "\n",
    "    # bypass diode constraint\n",
    "    ncells_bypass = nseries_cells / nbypass\n",
    "    vsubmodule = np.zeros(vsubs.shape)\n",
    "    for string in range(nstrings):\n",
    "        for submodule in range(nbypass * nmods):\n",
    "            jdx = int(submodule*ncells_bypass)\n",
    "            vsubmodule[submodule, string] = np.sum(vcells[jdx:int(jdx+ncells_bypass), string])\n",
    "    vsubmodule[vsubmodule < vbypass] = vbypass\n",
    "    residuals[idx0:idx1] = np.ravel(vsubmodule - vsubs)\n",
    "\n",
    "    # string constraints\n",
    "    vstrings = [np.sum(vsubs[:, string]) for string in range(nstrings)]\n",
    "    residuals[idx1:idx2] = vsys - np.array(vstrings)\n",
    "\n",
    "    # system constraint\n",
    "    residuals[-1] = np.sum(istrings) - isys\n",
    "    \n",
    "    return residuals"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "({'nfev': 1820,\n",
       "  'fjac': array([[-9.99624138e-01,  8.01234430e-04,  3.05076907e-05, ...,\n",
       "           3.38372362e-15,  2.30579607e-12,  3.92898781e-17],\n",
       "         [-6.65622409e-05, -9.99613737e-01,  7.59709723e-05, ...,\n",
       "           3.38676626e-15,  3.97123134e-12,  4.56848262e-18],\n",
       "         [-6.08932842e-05, -8.13107856e-05, -9.99345635e-01, ...,\n",
       "           3.38782842e-15, -1.29754627e-12, -6.15498768e-16],\n",
       "         ...,\n",
       "         [ 1.20318649e-05,  1.21131317e-03,  6.15106427e-04, ...,\n",
       "          -1.57970479e-01, -1.45660662e-02, -6.34531601e-02],\n",
       "         [ 2.01752822e-06,  3.11726575e-06,  1.08245364e-03, ...,\n",
       "           7.10489987e-02, -1.72674442e-01, -3.75597748e-02],\n",
       "         [ 1.57046806e-03,  1.57071003e-03,  1.57501669e-03, ...,\n",
       "          -4.40901236e-02, -4.40884819e-02,  9.21943067e-01]]),\n",
       "  'r': array([ 2.32598960e+02, -1.99255704e-01, -1.65809651e-02, ...,\n",
       "         -7.43576114e-01,  3.77448978e-02, -9.21179634e-01]),\n",
       "  'qtf': array([-2.81372496e-10, -1.20942828e-09,  4.83815989e-09, ...,\n",
       "          9.01168166e-10,  1.52157082e-09, -1.78711456e-11]),\n",
       "  'fvec': array([-1.49120982e-09, -1.43849999e-09,  1.87828331e-09, ...,\n",
       "          0.00000000e+00,  2.84217094e-14,  0.00000000e+00])},\n",
       " 1,\n",
       " 'The solution converged.')"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# make an initial guess\n",
    "varb = 0.5  # [V] some voltage - use voc_est with Rsh=np.Inf & Rs=0\n",
    "iarb = 6.0  # [A] some current\n",
    "v0 = [varb] * (NSTRINGS*NMODS*NSERIES_CELLS)  # [V] cell voltages guesses\n",
    "vsub0 = [varb * NSERIES_CELLS/NBYPASS] * (NSTRINGS*NMODS*NBYPASS)  # [V] submodule voltages guesses\n",
    "istr0 = [iarb] * NSTRINGS  # [A] string currents guesses\n",
    "isys0 = sum(istr0)  # [A] initial system current guess\n",
    "X0 = np.array(v0 + vsub0 + istr0 + [isys0])  # inital guess\n",
    "\n",
    "# add some noise to the light current\n",
    "noise = np.random.randn(NMODS*NSERIES_CELLS, NSTRINGS) / 1000\n",
    "\n",
    "# solve it\n",
    "ARGS = (NSERIES_CELLS, NMODS, NSTRINGS, NBYPASS, IL+noise,\n",
    "        I0, RSH, RS, GAMMA, VT, VBYPASS, VSYS)\n",
    "full_output = fsolve(isys_from_voltage, x0=X0, args=ARGS, full_output=True)\n",
    "# if full output then results is a tuple\n",
    "results = full_output[0]\n",
    "full_output[1:]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[0.41670409, 0.41669829, 0.41659487],\n",
       "       [0.41668225, 0.41652674, 0.41669965],\n",
       "       [0.41664967, 0.41669136, 0.41672096],\n",
       "       ...,\n",
       "       [0.41668421, 0.41666032, 0.41674502],\n",
       "       [0.41676059, 0.41670815, 0.41664592],\n",
       "       [0.41671843, 0.41663004, 0.41672851]])"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "idx0 = NSTRINGS*NMODS*NSERIES_CELLS\n",
    "voltages = results[:idx0].reshape(-1, NSTRINGS)  # cell voltages\n",
    "voltages"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[10.00008429, 10.00004239, 10.00042427],\n",
       "       [ 9.99981484, 10.00012778,  9.99970785],\n",
       "       [ 9.99994136,  9.99952761,  9.99967285],\n",
       "       [ 9.99995088, 10.00000426,  9.99957088],\n",
       "       [10.00007806,  9.99988728,  9.99968472],\n",
       "       [10.00051033,  9.99985692,  9.99959002],\n",
       "       [ 9.99982495,  9.999643  ,  9.99994605],\n",
       "       [10.00016988, 10.00013181,  9.99999447],\n",
       "       [ 9.99966196, 10.0005762 , 10.00009929],\n",
       "       [10.0001879 ,  9.99999273, 10.00005471],\n",
       "       [ 9.999977  , 10.00002123, 10.00033012],\n",
       "       [10.00004141,  9.99995615, 10.00050979],\n",
       "       [10.00020996,  9.99947279, 10.0005103 ],\n",
       "       [ 9.99983156, 10.00009486, 10.00010905],\n",
       "       [10.0000343 , 10.0004712 , 10.00009353],\n",
       "       [10.0000557 ,  9.99985237, 10.00006917],\n",
       "       [ 9.99973493, 10.00027908, 10.0004533 ],\n",
       "       [ 9.99960626,  9.99998118, 10.00001254],\n",
       "       [ 9.99975987,  9.99996528, 10.00002193],\n",
       "       [10.00040464,  9.9996972 ,  9.99989118],\n",
       "       [10.0000792 ,  9.99978051,  9.99948485],\n",
       "       [ 9.99966007, 10.00021479,  9.99955789],\n",
       "       [10.00060179, 10.0004249 , 10.00021889],\n",
       "       [ 9.99977887,  9.99999847,  9.99999236]])"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "idx1 = idx0 + NSTRINGS*NMODS*NBYPASS\n",
    "submodules = results[idx0:idx1].reshape(-1, NSTRINGS)\n",
    "submodules"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([6.42207435, 6.42213142, 6.42219736])"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "strings = results[idx1:(idx1+NSTRINGS)]\n",
    "strings"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "19.26640313192219"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "system = results[-1]\n",
    "system"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "# optimization objective function\n",
    "def mpp(vsys, nseries_cells, nmods, nstrings, nbypass,\n",
    "        il, i0, rsh, rs, gamma, vt, vbypass, guess):\n",
    "    r = fsolve(isys_from_voltage, x0=guess,\n",
    "               args=(nseries_cells, nmods, nstrings, nbypass,\n",
    "                     il, i0, rsh, rs, gamma, vt, vbypass, vsys),\n",
    "               full_output=False)\n",
    "    psys = r[-1]*vsys\n",
    "    return -psys"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "290.3329375367371"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# estimate upper limit of system voltage\n",
    "VOC_EST = (GAMMA*VT * np.log(IL / I0 + 1.0))*NMODS*NSERIES_CELLS\n",
    "VOC_EST"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "     fun: -4625.825388207809\n",
       "    nfev: 16\n",
       "     nit: 12\n",
       " success: True\n",
       "       x: 238.24038460561962"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# find max power\n",
    "optimum = minimize_scalar(mpp, (0, VSYS, VOC_EST),\n",
    "                          args=(NSERIES_CELLS, NMODS, NSTRINGS, NBYPASS,\n",
    "                                IL+noise, I0, RSH, RS, GAMMA, VT, VBYPASS, X0))\n",
    "optimum"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "238.24038460561962"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Vmp = optimum.x\n",
    "Vmp"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "4625.825388207809"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Pmp = -optimum.fun\n",
    "Pmp"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "19.41662995493122"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Imp = Pmp / Vmp\n",
    "Imp"
   ]
  }
 ],
 "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
}

PV-system-MPP-implicit-solver.ipynb

pvmismatch-example.ipynb

solar-convex-optimum.png

Note: finding the max power point may be an non-convex optimization problem, therefore to guarantee convergence, find the index of
the max power point by calling numpy.argmax(power). If the spacing of points on the IV curve is too course, then use brentq in a convex trust region around the index of the max power point.