bishop88_init_voc.py

# -*- coding: utf-8 -*-
"""
Created on Thu Apr 25 14:39:40 2024

@author: cliff
"""
import numpy as np
import pandas as pd
import pvlib
from pvlib.singlediode import bishop88, estimate_voc, _prepare_newton_inputs
from scipy.optimize import brentq, newton

# parameter sets from smmeredith 
# https://github.com/pvlib/pvlib-python/issues/2013#issuecomment-2073155300
params = {
  "0": {
    "alpha_sc": 0.001,
    "gamma_ref": 1.4084485593678782,
    "mu_gamma": 0.0008751660156250001,
    "I_L_ref": 2.510492627609538,
    "I_o_ref": 3.955007235783788e-10,
    "R_sh_ref": 4900,
    "R_sh_0": 8800,
    "R_s": 6.76,
    "cells_in_series": 264,
    "R_sh_exp": 7.0,
    "EgRef": 1.5
  },
  "1": {
    "alpha_sc": 0.001016,
    "gamma_ref": 1.1941564983432582,
    "mu_gamma": 0.00035303955078124936,
    "I_L_ref": 2.552079086883624,
    "I_o_ref": 4.231194870537287e-12,
    "R_sh_ref": 7100,
    "R_sh_0": 7100,
    "R_s": 5.164,
    "cells_in_series": 264,
    "R_sh_exp": 5.5,
    "EgRef": 1.5
  },
  "2": {
    "alpha_sc": 0.00102,
    "gamma_ref": 1.169653055295,
    "mu_gamma": 0.00031853759765625,
    "I_L_ref": 2.5523401529493297,
    "I_o_ref": 2.3970457796738275e-12,
    "R_sh_ref": 6400,
    "R_sh_0": 6400,
    "R_s": 5.277,
    "cells_in_series": 264,
    "R_sh_exp": 5.5,
    "EgRef": 1.5
  },
  "3": {
    "alpha_sc": 0.00102,
    "gamma_ref": 1.1717140423666776,
    "mu_gamma": 0.0003388037109374999,
    "I_L_ref": 2.5604199291711414,
    "I_o_ref": 2.4215577240761377e-12,
    "R_sh_ref": 6900,
    "R_sh_0": 6900,
    "R_s": 5.068,
    "cells_in_series": 264,
    "R_sh_exp": 5.5,
    "EgRef": 1.5
  },
  "4": {
    "alpha_sc": 0.00106,
    "gamma_ref": 1.369626412489172,
    "mu_gamma": 0.001527337646484375,
    "I_L_ref": 2.662739113120348,
    "I_o_ref": 2.1030694268456338e-10,
    "R_sh_ref": 9100,
    "R_sh_0": 13000,
    "R_s": 5.402,
    "cells_in_series": 271,
    "R_sh_exp": 5.5,
    "EgRef": 1.5
  },
  "5": {
    "alpha_sc": 0.001036,
    "gamma_ref": 1.1894167823312154,
    "mu_gamma": 0.0002741040039062494,
    "I_L_ref": 2.595323159794087,
    "I_o_ref": 2.512992175286693e-12,
    "R_sh_ref": 13200,
    "R_sh_0": 13200,
    "R_s": 3.878,
    "cells_in_series": 264,
    "R_sh_exp": 5.5,
    "EgRef": 1.5
  },
  "6": {
    "alpha_sc": 0.001036,
    "gamma_ref": 1.2354899524789542,
    "mu_gamma": 0.00045881187318457043,
    "I_L_ref": 2.591148415191749,
    "I_o_ref": 7.213105050899882e-12,
    "R_sh_ref": 14100,
    "R_sh_0": 14100,
    "R_s": 4.017,
    "cells_in_series": 264,
    "R_sh_exp": 4.5,
    "EgRef": 1.5
  },
  "7": {
    "alpha_sc": 0.0010400000000000001,
    "gamma_ref": 1.2107613754021975,
    "mu_gamma": 0.00034734741210937497,
    "I_L_ref": 2.60238977108264,
    "I_o_ref": 3.740548935305017e-12,
    "R_sh_ref": 15000,
    "R_sh_0": 15000,
    "R_s": 3.888,
    "cells_in_series": 264,
    "R_sh_exp": 5.5,
    "EgRef": 1.5
  },
  "8": {
    "alpha_sc": 0.001044,
    "gamma_ref": 1.1975609455461838,
    "mu_gamma": 0.0002969702148437495,
    "I_L_ref": 2.6123968333667733,
    "I_o_ref": 2.6124981962119807e-12,
    "R_sh_ref": 14000,
    "R_sh_0": 14500,
    "R_s": 3.641,
    "cells_in_series": 264,
    "R_sh_exp": 5.5,
    "EgRef": 1.5
  },
  "9": {
    "alpha_sc": 0.00107,
    "gamma_ref": 1.2003103709150689,
    "mu_gamma": 0.00094519775390625,
    "I_L_ref": 2.671958847107402,
    "I_o_ref": 4.7590638988688474e-12,
    "R_sh_ref": 13700,
    "R_sh_0": 13700,
    "R_s": 4.161,
    "cells_in_series": 271,
    "R_sh_exp": 5.5,
    "EgRef": 1.5
  },
  "10": {
    "alpha_sc": 0.0012052,
    "gamma_ref": 1.4514249591272739,
    "mu_gamma": 0.00024388427734374946,
    "I_L_ref": 3.042521297278528,
    "I_o_ref": 4.911461636489296e-10,
    "R_sh_ref": 11500,
    "R_sh_0": 12000,
    "R_s": 3.885,
    "cells_in_series": 268,
    "R_sh_exp": 5.5,
    "EgRef": 1.5
  },
  "11": {
    "alpha_sc": 0.0012079999999999999,
    "gamma_ref": 1.4347976874113653,
    "mu_gamma": 0.00019899658203124953,
    "I_L_ref": 3.048129655613767,
    "I_o_ref": 3.5982220311717015e-10,
    "R_sh_ref": 12000,
    "R_sh_0": 12500,
    "R_s": 3.755,
    "cells_in_series": 268,
    "R_sh_exp": 5.5,
    "EgRef": 1.5
  },
  "12": {
    "alpha_sc": 0.0012112,
    "gamma_ref": 1.3407235630860397,
    "mu_gamma": 0.0006101342773437489,
    "I_L_ref": 3.0376398708690666,
    "I_o_ref": 7.524792039031038e-11,
    "R_sh_ref": 5500,
    "R_sh_0": 6000,
    "R_s": 5.263,
    "cells_in_series": 268,
    "R_sh_exp": 5.5,
    "EgRef": 1.5
  },
  "13": {
    "alpha_sc": 0.001214,
    "gamma_ref": 1.4415843664096415,
    "mu_gamma": 0.0008313085937499984,
    "I_L_ref": 3.0393075212182685,
    "I_o_ref": 4.0183201130931835e-10,
    "R_sh_ref": 9000,
    "R_sh_0": 9500,
    "R_s": 4.818,
    "cells_in_series": 268,
    "R_sh_exp": 5.5,
    "EgRef": 1.5
  },
  "14": {
    "alpha_sc": 0.0012168,
    "gamma_ref": 1.3630431939833358,
    "mu_gamma": 0.0006750585937499987,
    "I_L_ref": 3.046848984680672,
    "I_o_ref": 1.0236166994787753e-10,
    "R_sh_ref": 7000,
    "R_sh_0": 7000,
    "R_s": 4.957,
    "cells_in_series": 268,
    "R_sh_exp": 5.5,
    "EgRef": 1.5
  },
  "15": {
    "alpha_sc": 0.00122,
    "gamma_ref": 1.3539877202291934,
    "mu_gamma": 0.0006272265624999989,
    "I_L_ref": 3.054640150849507,
    "I_o_ref": 8.24804073383895e-11,
    "R_sh_ref": 7500,
    "R_sh_0": 7500,
    "R_s": 4.781,
    "cells_in_series": 268,
    "R_sh_exp": 5.5,
    "EgRef": 1.5
  },
  "16": {
    "alpha_sc": 0.0012228,
    "gamma_ref": 1.3544935976803756,
    "mu_gamma": 0.0005923730468749989,
    "I_L_ref": 3.0613273003417176,
    "I_o_ref": 7.883312366503028e-11,
    "R_sh_ref": 8500,
    "R_sh_0": 8500,
    "R_s": 4.543,
    "cells_in_series": 268,
    "R_sh_exp": 5.5,
    "EgRef": 1.5
  },
  "17": {
    "alpha_sc": 0.0012256,
    "gamma_ref": 1.2916241612804187,
    "mu_gamma": 0.00047308959960937403,
    "I_L_ref": 3.068717040806731,
    "I_o_ref": 2.2691248021217617e-11,
    "R_sh_ref": 7000,
    "R_sh_0": 7000,
    "R_s": 4.602,
    "cells_in_series": 268,
    "R_sh_exp": 5.5,
    "EgRef": 1.5
  }
}
    

def bishop88_mpp(photocurrent, saturation_current, resistance_series,
                 resistance_shunt, nNsVth, d2mutau=0, NsVbi=np.Inf,
                 breakdown_factor=0., breakdown_voltage=-5.5,
                 breakdown_exp=3.28, method='newton', method_kwargs=None):
    """
    Find max power point.

    Parameters
    ----------
    photocurrent : numeric
        photogenerated current (Iph or IL) [A]
    saturation_current : numeric
        diode dark or saturation current (Io or Isat) [A]
    resistance_series : numeric
        series resistance (Rs) in [Ohm]
    resistance_shunt : numeric
        shunt resistance (Rsh) [Ohm]
    nNsVth : numeric
        product of diode ideality factor (n), number of series cells (Ns), and
        thermal voltage (Vth = k_b * T / q_e) in volts [V]
    d2mutau : numeric, default 0
        PVsyst parameter for cadmium-telluride (CdTe) and amorphous-silicon
        (a-Si) modules that accounts for recombination current in the
        intrinsic layer. The value is the ratio of intrinsic layer thickness
        squared :math:`d^2` to the diffusion length of charge carriers
        :math:`\\mu \\tau`. [V]
    NsVbi : numeric, default np.inf
        PVsyst parameter for cadmium-telluride (CdTe) and amorphous-silicon
        (a-Si) modules that is the product of the PV module number of series
        cells ``Ns`` and the builtin voltage ``Vbi`` of the intrinsic layer.
        [V].
    breakdown_factor : numeric, default 0
        fraction of ohmic current involved in avalanche breakdown :math:`a`.
        Default of 0 excludes the reverse bias term from the model. [unitless]
    breakdown_voltage : numeric, default -5.5
        reverse breakdown voltage of the photovoltaic junction :math:`V_{br}`
        [V]
    breakdown_exp : numeric, default 3.28
        avalanche breakdown exponent :math:`m` [unitless]
    method : str, default 'newton'
       Either ``'newton'`` or ``'brentq'``. ''method'' must be ``'newton'``
       if ``breakdown_factor`` is not 0.
    method_kwargs : dict, optional
        Keyword arguments passed to root finder method. See
        :py:func:`scipy:scipy.optimize.brentq` and
        :py:func:`scipy:scipy.optimize.newton` parameters.
        ``'full_output': True`` is allowed, and ``optimizer_output`` would be
        returned. See examples section.

    Returns
    -------
    tuple
        max power current ``i_mp`` [A], max power voltage ``v_mp`` [V], and
        max power ``p_mp`` [W]
    optimizer_output : tuple, optional, if specified in ``method_kwargs``
        see root finder documentation for selected method.
        Found root is diode voltage in [1]_.

    Examples
    --------
    Using the following arguments that may come from any
    `calcparams_.*` function in :py:mod:`pvlib.pvsystem`:

    >>> args = {'photocurrent': 1., 'saturation_current': 9e-10, 'nNsVth': 4.,
    ...         'resistance_series': 4., 'resistance_shunt': 5000.0}

    Use default values:

    >>> i_mp, v_mp, p_mp = bishop88_mpp(**args)

    Specify tolerances and maximum number of iterations:

    >>> i_mp, v_mp, p_mp = bishop88_mpp(**args, method='newton',
    ...     method_kwargs={'tol': 1e-3, 'rtol': 1e-3, 'maxiter': 20})

    Retrieve full output from the root finder:

    >>> (i_mp, v_mp, p_mp), method_output = bishop88_mpp(**args,
    ...     method='newton', method_kwargs={'full_output': True})
    """
    # collect args
    args = (photocurrent, saturation_current,
            resistance_series, resistance_shunt, nNsVth, d2mutau, NsVbi,
            breakdown_factor, breakdown_voltage, breakdown_exp)
    method = method.lower()

    # method_kwargs create dict if not provided
    # this pattern avoids bugs with Mutable Default Parameters
    if not method_kwargs:
        method_kwargs = {}

    # first bound the search using voc
    voc_est = estimate_voc(photocurrent, saturation_current, nNsVth)

    def fmpp(x, *a):
        return bishop88(x, *a, gradients=True)[6]

    if method == 'brentq':
        # break out arguments for numpy.vectorize to handle broadcasting
        xp = np.where(voc_est < NsVbi, voc_est, 0.99*NsVbi)
        vec_fun = np.vectorize(
            lambda voc, iph, isat, rs, rsh, gamma, d2mutau, NsVbi, vbr_a, vbr,
            vbr_exp: brentq(fmpp, 0.0, voc,
                            args=(iph, isat, rs, rsh, gamma, d2mutau, NsVbi,
                                  vbr_a, vbr, vbr_exp),
                            **method_kwargs)
        )
        vd = vec_fun(xp, *args)
    elif method == 'newton':
        # make sure all args are numpy arrays if max size > 1
        # if voc_est is an array, then make a copy to use for initial guess, v0
        xp = np.where(voc_est < NsVbi, voc_est, 0.99*NsVbi)
        x0, args, method_kwargs = \
            _prepare_newton_inputs(xp, args, method_kwargs)
        vd = newton(func=fmpp, x0=x0,
                    fprime=lambda x, *a: bishop88(x, *a, gradients=True)[7],
                    args=args, **method_kwargs)
    else:
        raise NotImplementedError("Method '%s' isn't implemented" % method)

    # When 'full_output' parameter is specified, returned 'vd' is a tuple with
    # many elements, where the root is the first one. So we use it to output
    # the bishop88 result and return
    # tuple(tuple with bishop88 solution, tuple with method results)
    if method_kwargs.get('full_output') is True:
        return (bishop88(vd[0], *args), vd)
    else:
        return bishop88(vd, *args)


def check_solution(i, v, il, io, rs, rsh, a, NsVbi=np.inf, d2mutau=0.):
    vd = v + rs * i
    return il - io*np.expm1(vd/a) - vd/rsh - il*d2mutau/(NsVbi - vd) - i


#%% set up weather
irrad = np.arange(20, 1100, 20)
tc = np.arange(-25, 74, 1)
weather = np.array(np.meshgrid(irrad, tc)).T.reshape(-1, 2)


#%% Newton method. Have to handle newton and brentq separately because brentq
# doesn't fail gracefully

method = 'newton'

print(method+'\n')
print('==================================')


keys = params.keys()
# Loop over parameter sets
for k in keys:
    p = params[k]
    sde_params = pvlib.pvsystem.calcparams_pvsyst(
        weather[:, 0], weather[:, 1], **p)
    
    # line up quantities
    df = pd.DataFrame(data=np.array(sde_params).T,
                      columns=['IL', 'Io', 'Rs', 'Rsh', 'a'])        
    # Use PVsyst defaults for Vbi and d2mutau
    # https://www.pvsyst.com/help/index.html?thinfilm_recombinationloss.htm
    # assume single junction, built-in voltage of 0.9V per cell
    # assume d2mutau of 1.4 V
    NsVbi = params[k]['cells_in_series'] * 0.9
    d2mutau = 1.4
    
    voc_est = pvlib.singlediode.estimate_voc(sde_params[0], sde_params[1],
                                             sde_params[-1])

    print('Before fix')
    #with np.errstate(over='ignore', divide='ignore', invalid='ignore'):
    result = pvlib.singlediode.bishop88_mpp(
            *sde_params, d2mutau=d2mutau, NsVbi=NsVbi, method=method)

    imp, vmp, pmp = result
    err = check_solution(imp, vmp, df['IL'], df['Io'], df['Rs'],
                         df['Rsh'], df['a'], NsVbi, d2mutau)
    bad_results = np.isnan(pmp) | (pmp < 0) | (err.abs() > 0.00001)
    print('  Failed MPP: ' + str(bad_results.sum()))
    
    print('After fix')
    result2 = bishop88_mpp(*sde_params, d2mutau=d2mutau, NsVbi=NsVbi)
    imp2, vmp2, pmp2 = result2
    
    err2 = check_solution(imp2, vmp2, df['IL'], df['Io'], df['Rs'],
                          df['Rsh'], df['a'], NsVbi, d2mutau)
    bad_results2 = np.isnan(pmp2) | (pmp2 < 0) | (err2.abs() > 0.00001)
    print('  Failed MPP: ' + str(bad_results2.sum()))
    print('------------------------------------\n')


#%% Brentq method, take a while to run since we have to loop over weather
# to handle failures

method = 'brentq'

print(method+'\n')
print('==================================')
for i, k in enumerate(keys):
    p = params[k]
    sde_params = pvlib.pvsystem.calcparams_pvsyst(
        weather[:, 0], weather[:, 1], **p)
    
    # line up quantities
    df = pd.DataFrame(data=np.array(sde_params).T,
                      columns=['IL', 'Io', 'Rs', 'Rsh', 'a'])        
    # Use PVsyst defaults for Vbi and d2mutau
    # https://www.pvsyst.com/help/index.html?thinfilm_recombinationloss.htm
    # assume single junction, built-in voltage of 0.9V per cell
    # assume d2mutau of 1.4 V
    NsVbi = params[k]['cells_in_series'] * 0.9
    d2mutau = 1.4
    
    voc_est = pvlib.singlediode.estimate_voc(sde_params[0], sde_params[1],
                                             sde_params[-1])
    failed_pmp = 0
    print('Before fix')
    #with np.errstate(over='ignore', divide='ignore', invalid='ignore'):
    for j in df.index:
        try:
            result = pvlib.singlediode.bishop88_mpp(
                *df.loc[j,:], d2mutau=d2mutau, NsVbi=NsVbi, method=method)
            imp, vmp, pmp = result
            if np.isnan(pmp) | (pmp<0):
                failed_pmp += 1
            else:
                err = check_solution(imp, vmp, df.loc[j,'IL'], df.loc[j,'Io'],
                                     df.loc[j,'Rs'], df.loc[j,'Rsh'],
                                     df.loc[j,'a'], NsVbi, d2mutau)
                if np.abs(err)>0.0001:
                    failed_pmp += 1
        except ValueError:
            failed_pmp += 1
            continue
    print('  Failed MPP: ' + str(failed_pmp))

    
    print('After fix')
    result2 = bishop88_mpp(*sde_params, d2mutau=d2mutau, NsVbi=NsVbi)
    imp2, vmp2, pmp2 = result2
    
    err2 = check_solution(imp2, vmp2, df['IL'], df['Io'], df['Rs'],
                          df['Rsh'], df['a'], NsVbi, d2mutau)
    bad_results2 = np.isnan(pmp2) | (pmp2 < 0) | (err2.abs() > 0.00001)
    print('  Failed MPP: ' + str(bad_results2.sum()))
    print('------------------------------------\n')