From a0588d49f6cf35187ab0455ba9c762ae3a1d1525 Mon Sep 17 00:00:00 2001 From: joallen6 Date: Thu, 16 Jun 2016 09:54:23 -0700 Subject: [PATCH] Add files via upload --- PVLIB/pvl_absoluteairmass.m | 8 +- PVLIB/pvl_calcparams_desoto.m | 115 ++++++++--------------- PVLIB/pvl_grounddiffuse.m | 6 +- PVLIB/pvl_isotropicsky.m | 4 +- PVLIB/pvl_perez.m | 36 ++++---- PVLIB/pvl_sapmcelltemp.m | 38 +++----- PVLIB/pvl_singleaxis.m | 16 ++-- PVLIB/pvl_singlediode.m | 167 ++++++++++++++++++++++++++++------ PVLIB/pvl_snlinverter.m | 15 +-- 9 files changed, 230 insertions(+), 175 deletions(-) diff --git a/PVLIB/pvl_absoluteairmass.m b/PVLIB/pvl_absoluteairmass.m index e3ea43c..f858af9 100755 --- a/PVLIB/pvl_absoluteairmass.m +++ b/PVLIB/pvl_absoluteairmass.m @@ -1,6 +1,6 @@ function AMa = pvl_absoluteairmass(AMrelative,pressure) % PVL_ABSOLUTEAIRMASS Determine absolute (pressure corrected) airmass from relative airmass and pressure -% +% % Syntax % AMa = pvl_absoluteairmass(AMrelative, pressure) % @@ -29,11 +29,11 @@ % clear sky solar irradiance models using theoretical and measured data," % Solar Energy, vol. 51, pp. 121-138, 1993. % -% See also PVL_RELATIVEAIRMASS +% See also PVL_RELATIVEAIRMASS p=inputParser; -p.addRequired('AMrelative', @(x) all(isnumeric(x) | isnan(x))); -p.addRequired('pressure', @(x) all(isnumeric(x) & x>=0)); +p.addRequired('AMrelative', @(x) isnumeric(x)); +p.addRequired('pressure', @(x) isnumeric(x) & all(x>=0 | isnan(x))); p.parse(AMrelative,pressure); AMa = AMrelative.*pressure/101325; \ No newline at end of file diff --git a/PVLIB/pvl_calcparams_desoto.m b/PVLIB/pvl_calcparams_desoto.m index 16b5d7b..e849fff 100755 --- a/PVLIB/pvl_calcparams_desoto.m +++ b/PVLIB/pvl_calcparams_desoto.m @@ -1,46 +1,44 @@ function [IL, I0, Rs, Rsh, nNsVth] = pvl_calcparams_desoto(S,Tcell,alpha_isc,ModuleParameters, EgRef, dEgdT, varargin) -% PVL_CALCPARAMS_DESOTO Applies temperature and irradiance corrections to reference parameters per [1] +% PVL_CALCPARAMS_DESOTO Calculates five parameters to compute IV curves using the De Soto model. % % Syntax % [IL, I0, Rs, Rsh, nNsVth] = pvl_calcparams_desoto(S, Tcell, alpha_isc, ModuleParameters, EgRef, dEgdt) -% [IL, I0, Rs, Rsh, nNsVth] = pvl_calcparams_desoto(S, Tcell, alpha_isc, ModuleParameters, EgRef, dEgdt, M) -% [IL, I0, Rs, Rsh, nNsVth] = pvl_calcparams_desoto(S, Tcell, alpha_isc, ModuleParameters, EgRef, dEgdt, M, Sref) -% [IL, I0, Rs, Rsh, nNsVth] = pvl_calcparams_desoto(S, Tcell, alpha_isc, ModuleParameters, EgRef, dEgdt, M, Sref, Tref) +% [IL, I0, Rs, Rsh, nNsVth] = pvl_calcparams_desoto(S, Tcell, alpha_isc, ModuleParameters, EgRef, dEgdt) +% [IL, I0, Rs, Rsh, nNsVth] = pvl_calcparams_desoto(S, Tcell, alpha_isc, ModuleParameters, EgRef, dEgdt, Sref) +% [IL, I0, Rs, Rsh, nNsVth] = pvl_calcparams_desoto(S, Tcell, alpha_isc, ModuleParameters, EgRef, dEgdt, Sref, Tref) % [IL, I0, Rs, Rsh, nNsVth] = pvl_calcparams_desoto(S, Tcell, alpha_isc, ModuleParameters, EgRef, dEgdt, 'Sref', Sref, 'Tref', Tref) % % Description -% Applies the temperature and irradiance corrections to the IL, I0, -% Rs, Rsh, and a parameters at reference conditions (IL_ref, I0_ref, -% etc.) according to the De Soto et. al description given in [1]. The -% results of this correction procedure may be used in a single diode -% model to determine IV curves at irradiance = S, cell temperature = -% Tcell. +% Applies the temperature and irradiance corrections to calculate the +% five parameters for an IV curve according to the De Soto model [1]. +% The results of this procedure may be used to determine the IV curve +% at effective irradiance S and cell temperature Tcell. % % Input Parameters: -% S - The irradiance (in W/m^2) absorbed by the module. S must be >= 0. -% May be a vector of irradiances, but must be the same size as all -% other input vectors. Due to a division by S in the script, any S +% S - The effective irradiance (in W/m^2) absorbed by the module. S must be >= 0. +% May be a vector the same size as Tcell. Due to a division by S in the script, any S % value equal to 0 will be set to 1E-10. % Tcell - The average cell temperature of cells within a module in C. -% Tcell must be >= -273.15. May be a vector of cell temperatures, but -% must be the same size as all other input vectors. +% Tcell must be >= -273.15. May be a vector of the same size as S. % alpha_isc - The short-circuit current temperature coefficient of the -% module in units of 1/C. +% module in units of A/C (or A/K). % ModuleParameters - a struct with parameters describing PV module -% performance at reference conditions according to DeSoto's paper. A +% performance at reference conditions according to DeSoto's paper. % Parameters may be generated or found by lookup. For ease of use, a library -% of parameters may be found within the System Advisor Model (SAM) [2]. +% of modules A SAM library of modules is provided with PV_LIB (|\Required Data\CECModuleDatabaseSAM2014.1.14.mat|), +% which may be read by the SAM library reader +% function . may be found within the System Advisor Model (SAM) [2]. % The SAM library has been provided as a .mat file, or may % be read and converted into a .mat file by the SAM library reader % functions. The ModuleParameters struct must contain (at least) the -% following 5 fields: +% following fields: % ModuleParameters.a_ref - modified diode ideality factor parameter at -% reference conditions (units of eV), a_ref can be calculated from the +% reference conditions, a_ref can be calculated from the % usual diode ideality factor (n), number of cells in series (Ns), % and cell temperature (Tcell) per equation (2) in [1]. % ModuleParameters.IL_ref - Light-generated current (or photocurrent) % in amperes at reference conditions. This value is referred to -% as Iph in some literature. +% as Iphi in some literature. % ModuleParameters.I0_ref - diode reverse saturation current in amperes, % under reference conditions. % ModuleParameters.Rsh_ref - shunt resistance under reference conditions (ohms) @@ -51,14 +49,6 @@ % May be either a scalar value (e.g. -0.0002677 as in [1]) or a % vector of dEgdT values corresponding to each input condition (this % may be useful if dEgdT is a function of temperature). -% M - An optional airmass modifier, if omitted, M is given a value of 1, -% which assumes absolute (pressure corrected) airmass = 1.5. In this -% code, M is equal to M/Mref as described in [1] (i.e. Mref is assumed -% to be 1). Source [1] suggests that an appropriate value for M -% as a function absolute airmass (AMa) may be: -% M = polyval([-0.000126, 0.002816, -0.024459, 0.086257, 0.918093], AMa) -% M may be a vector, but must be of the same size as all other input -% vectors. % Sref - Optional reference irradiance in W/m^2. If omitted, a value of % 1000 is used. % Tref - Optional reference cell temperature in C. If omitted, a value of @@ -119,44 +109,37 @@ % Silicon (Si): % EgRef = 1.121 % dEgdT = -0.0002677 -% M = polyval([-0.000126 0.002816 -0.024459 0.086257 0.918093], AMa) % Source = Reference 1 % Cadmium Telluride (CdTe): % EgRef = 1.475 % dEgdT = -0.0003 -% M = polyval([-2.46E-5 9.607E-4 -0.0134 0.0716 0.9196], AMa) % Source = Reference 4 % Copper Indium diSelenide (CIS): % EgRef = 1.010 % dEgdT = -0.00011 -% M = polyval([-3.74E-5 0.00125 -0.01462 0.0718 0.9210], AMa) % Source = Reference 4 % Copper Indium Gallium diSelenide (CIGS): % EgRef = 1.15 % dEgdT = ???? -% M = polyval([-9.07E-5 0.0022 -0.0202 0.0652 0.9417], AMa) % Source = Wikipedia % Gallium Arsenide (GaAs): % EgRef = 1.424 % dEgdT = -0.000433 -% M = unknown % Source = Reference 4 p = inputParser; -p.addRequired('S',@(x) all(x>=0) & isnumeric(x) & isvector(x) ); -p.addRequired('Tcell',@(x) all(x>=-273.15) & isnumeric(x) & isvector(x) ); -p.addRequired('alpha_isc', @(x) (isnumeric(x) & isvector(x))); -p.addRequired('ModuleParameters', @(x) (isstruct(x))); -p.addRequired('EgRef', @(x) (isnumeric(x) & isvector(x) & x>0)); -p.addRequired('dEgdT', @(x) (isnumeric(x) & isvector(x))); -p.addOptional('M', 1, @(x) (isnumeric(x) & isvector(x))); -p.addOptional('Sref',1000, @(x) (all(x>0)& isnumeric(x) & isvector(x))); -p.addOptional('Tref',25, @(x) (all(x>-273.15) & isnumeric(x) & isvector(x))); +p.addRequired('S',@(x) isnumeric(x) && isvector(x) && all(x>=0 | isnan(x))); +p.addRequired('Tcell',@(x) isnumeric(x) && isvector(x) && all(x>=-273.15 | isnan(x))); +p.addRequired('alpha_isc', @(x) isnumeric(x) && isvector(x)); +p.addRequired('ModuleParameters', @(x) isstruct(x)); +p.addRequired('EgRef', @(x) isnumeric(x) && isvector(x) && all(x>0 | isnan(x))); +p.addRequired('dEgdT', @(x) isnumeric(x) && isvector(x)); +p.addOptional('Sref',1000, @(x) isnumeric(x) && isvector(x) && (all(x>0 | isnan(x)))); +p.addOptional('Tref',25, @(x) isnumeric(x) && isvector(x) && (all(x>-273.15 | isnan(x)))); p.parse(S, Tcell, alpha_isc, ModuleParameters, EgRef, dEgdT, varargin{:}); S = p.Results.S(:); -M = p.Results.M(:); -M = max(M,0); % Ensure that the airmass modifier is never negative + Sref = p.Results.Sref(:); Tref = p.Results.Tref(:); Tcell = p.Results.Tcell(:); @@ -165,57 +148,31 @@ dEgdT = p.Results.dEgdT(:); EgRef = p.Results.EgRef(:); -%Reference parameter should be taken and input from CEC database or similar a_ref=ModuleParameters.a_ref(:); IL_ref=ModuleParameters.IL_ref(:); I0_ref=ModuleParameters.I0_ref(:); Rsh_ref=ModuleParameters.Rsh_ref(:); -Rs_ref=ModuleParameters.Rs_ref(:); %R_s=R_s_ref -% -VectorSizes = [numel(S), numel(M), numel(Tcell), numel(dEgdT),... - numel(Sref), numel(Tref), numel(EgRef), numel(a_ref), numel(IL_ref),... - numel(I0_ref), numel(Rsh_ref), numel(Rs_ref), numel(alpha_isc)]; +Rs_ref=ModuleParameters.Rs_ref(:); + +VectorSizes = [numel(S), numel(Tcell)]; MaxVectorSize = max(VectorSizes); if not(all((VectorSizes==MaxVectorSize) | (VectorSizes==1))) - error(['Input vectors S, Tcell, alpha_isc, EgRef, dEgdT, M, Sref (if used), '... - 'Tref (if used), and all used components of ModuleParameters must '... + error(['Input vectors S and Tcell must '... 'either be scalars or vectors of the same length.']); end %k=1.3806488e-23; %Boltzman's constant in units of J/K k = 8.617332478e-5; % Boltzmann constant in units of eV/K - Tref_K=Tref+273.15; % Reference cell temperature in Kelvin Tcell_K=Tcell+273.15; % cell temperature in Kelvin - - - -% We divide by S later, so make S a very small number for any S == 0. -S(S==0) = 1E-10; - -% -% % If any input variable is not a scalar, then make any scalar input values -% % into a column vector of the correct size. -% if MaxVectorSize >1 -% if VectorSizes(1) < MaxVectorSize -% S = S.*ones(MaxVectorSize , 1); -% end -% if VectorSizes(2) < MaxVectorSize -% M = M.*ones(MaxVectorSize, 1); -% end -% if VectorSizes(3) < MaxVectorSize -% Tcell_K = Tcell_K.*ones(MaxVectorSize, 1); -% end -% end - -%These parameters (a, I_L, I_o, M, and R_sh) need to be vectors of equal -%length to the number of conditions (number of S,Tcell pairs). E_g=EgRef.*(1+dEgdT.*(Tcell_K-Tref_K)); % Equation 10 in [1] nNsVth=a_ref.*(Tcell_K./Tref_K); % Equation 8 in [1] -IL=S./Sref.*M.*(IL_ref+alpha_isc.*(Tcell_K-Tref_K)); -% Equation 9 in [1] +IL=S./Sref.*(IL_ref+alpha_isc.*(Tcell_K-Tref_K)); +IL(S <= 0) = 0; % If there is no light then no current is produced I0=I0_ref.*((Tcell_K./Tref_K).^3).*exp((EgRef./(k.*Tref_K))-(E_g./(k.*Tcell_K))); +I0(IL==0) = 0; % If there is no light-generated current, there is no reverse saturation current Rsh=Rsh_ref.*(Sref./S); % Equation 12 in [1] +Rsh(S <= 0) = inf; % Rsh is undefined if there is no current Rs = Rs_ref; diff --git a/PVLIB/pvl_grounddiffuse.m b/PVLIB/pvl_grounddiffuse.m index d702699..c827b3c 100755 --- a/PVLIB/pvl_grounddiffuse.m +++ b/PVLIB/pvl_grounddiffuse.m @@ -38,9 +38,9 @@ p = inputParser; -p.addRequired('GHI',@(x) all(isnumeric(x) & isvector(x) & x>=0)); -p.addRequired('Albedo', @(x) all(isnumeric(x) & isvector(x) & x>=0 & x<=1)); -p.addRequired('SurfTilt', @(x) all(isvector(x) & isnumeric(x))); +p.addRequired('GHI',@(x) isnumeric(x) && isvector(x) && all(x>=0 | isnan(x))); +p.addRequired('Albedo', @(x) isnumeric(x) && isvector(x) && all((x>=0 & x<=1) | isnan(x))); +p.addRequired('SurfTilt', @(x) isnumeric(x) && isvector(x)); p.parse(GHI,Albedo,SurfTilt); GHI = GHI(:); diff --git a/PVLIB/pvl_isotropicsky.m b/PVLIB/pvl_isotropicsky.m index 6c6e6af..e43d1c7 100755 --- a/PVLIB/pvl_isotropicsky.m +++ b/PVLIB/pvl_isotropicsky.m @@ -44,8 +44,8 @@ % p = inputParser; -p.addRequired('SurfTilt', @(x) all(isnumeric(x) & x<=180 & x>=0 & isvector(x))); -p.addRequired('DHI', @(x) all(isnumeric(x) & isvector(x) & x>=0)); +p.addRequired('SurfTilt', @(x) isnumeric(x) && isvector(x) && all((x>=0 & x<=180) | isnan(x))); +p.addRequired('DHI', @(x) isnumeric(x) && isvector(x) && all(x>=0 | isnan(x))); p.parse(SurfTilt, DHI); SkyDiffuse = DHI *(1+ cosd(SurfTilt)) * 0.5; diff --git a/PVLIB/pvl_perez.m b/PVLIB/pvl_perez.m index 28b97ca..0c8ae12 100755 --- a/PVLIB/pvl_perez.m +++ b/PVLIB/pvl_perez.m @@ -1,17 +1,16 @@ function SkyDiffuse = pvl_perez(SurfTilt, SurfAz, DHI, DNI, HExtra, SunZen, SunAz, AM, varargin) -% PVL_PEREZ Determine diffuse irradiance from the sky on a tilted surface using one of the Perez models +% PVL_PEREZ Determine diffuse irradiance from the sky on a tilted surface using the Perez model % % Syntax % SkyDiffuse = pvl_perez(SurfTilt, SurfAz, DHI, DNI, HExtra, SunZen, SunAz, AM) % SkyDiffuse = pvl_perez(SurfTilt, SurfAz, DHI, DNI, HExtra, SunZen, SunAz, AM, model) % % Description -% Perez models determine the diffuse irradiance from the sky (ground -% reflected irradiance is not included in this algorithm) on a tilted +% The Perez model [3] determines the sky diffuse irradiance on a tilted % surface using the surface tilt angle, surface azimuth angle, diffuse % horizontal irradiance, direct normal irradiance, extraterrestrial % irradiance, sun zenith angle, sun azimuth angle, and relative (not -% pressure-corrected) airmass. Optionally a selector may be used to use +% pressure-corrected) airmass. An optional selector may be used to specify % any of Perez's model coefficient sets. % % Inputs: @@ -64,17 +63,18 @@ % the input vector(s). % % References -% [1] Loutzenhiser P.G. et. al. "Empirical validation of models to compute -% solar irradiance on inclined surfaces for building energy simulation" -% 2007, Solar Energy vol. 81. pp. 254-267 +% [1] Loutzenhiser P.G. et. al., 2007. Empirical validation of models to compute +% solar irradiance on inclined surfaces for building energy simulation, +% Solar Energy vol. 81. pp. 254-267. % [2] Perez, R., Seals, R., Ineichen, P., Stewart, R., Menicucci, D., 1987. A new % simplified version of the Perez diffuse irradiance model for tilted % surfaces. Solar Energy 39 (3), 221–232. % [3] Perez, R., Ineichen, P., Seals, R., Michalsky, J., Stewart, R., 1990. % Modeling daylight availability and irradiance components from direct % and global irradiance. Solar Energy 44 (5), 271–289. -% [4] Perez, R. et. al 1988. "The Development and Verification of the -% Perez Diffuse Radiation Model". SAND88-7030 +% [4] Perez, R. et. al 1988. The Development and Verification of the +% Perez Diffuse Radiation Model,.SAND88-7030, Sandia National +% Laboratories. % % See also PVL_EPHEMERIS PVL_EXTRARADIATION PVL_ISOTROPICSKY % PVL_HAYDAVIES1980 PVL_REINDL1990 PVL_KLUCHER1979 PVL_KINGDIFFUSE @@ -84,15 +84,15 @@ % % -p = inputParser; -p.addRequired('SurfTilt', @(x) all(isnumeric(x) & x<=180 & x>=0 & isvector(x))); -p.addRequired('SurfAz', @(x) all(isnumeric(x) & x<=360 & x>=0 & isvector(x))); -p.addRequired('DHI', @(x) all(isnumeric(x) & isvector(x) & x>=0)); -p.addRequired('DNI', @(x) all(isnumeric(x) & isvector(x) & x>=0)); -p.addRequired('HExtra', @(x) all(isnumeric(x) & isvector(x) & x>=0)); -p.addRequired('SunZen', @(x) all(isnumeric(x) & x<=180 & x>=0 & isvector(x))); -p.addRequired('SunAz', @(x) all(isnumeric(x) & x<=360 & x>=0 & isvector(x))); -p.addRequired('AM', @(x) (all(((isnumeric(x) & x>=0) | isnan(x))) & isvector(x))); +p=inputParser; +p.addRequired('SurfTilt', @(x) isnumeric(x) && isvector(x) && all((x>=0 & x<=180) | isnan(x))); +p.addRequired('SurfAz', @(x) isnumeric(x) && isvector(x) && all((x>=0 & x<=360) | isnan(x))); +p.addRequired('DHI', @(x) isnumeric(x) && isvector(x) && all(x>=0 | isnan(x))); +p.addRequired('DNI', @(x) isnumeric(x) && isvector(x) && all(x>=0 | isnan(x))); +p.addRequired('HExtra', @(x) isnumeric(x) && isvector(x) && all(x>=0 | isnan(x))); +p.addRequired('SunZen', @(x) isnumeric(x) && isvector(x) && all((x>=0 & x<=180) | isnan(x))); +p.addRequired('SunAz', @(x) isnumeric(x) && isvector(x) && all((x>=0 & x<=360) | isnan(x))); +p.addRequired('AM', @(x) isnumeric(x) && isvector(x) && all(x>=0 | isnan(x))); p.addOptional('model', '1990', @(x) ischar(x)); p.parse(SurfTilt, SurfAz, DHI, DNI, HExtra, SunZen, SunAz, AM, varargin{:}); diff --git a/PVLIB/pvl_sapmcelltemp.m b/PVLIB/pvl_sapmcelltemp.m index 06edafc..154aafb 100755 --- a/PVLIB/pvl_sapmcelltemp.m +++ b/PVLIB/pvl_sapmcelltemp.m @@ -1,4 +1,4 @@ -function [Tcell Tmodule] = pvl_sapmcelltemp(E, E0, a, b, windspeed, Tamb, deltaT) +function [Tcell, Tmodule] = pvl_sapmcelltemp(E, E0, a, b, windspeed, Tamb, deltaT) % PVL_SAPMCELLTEMP Estimate cell temperature from irradiance, windspeed, ambient temperature, and module parameters (SAPM) % % Syntax @@ -41,33 +41,21 @@ % See also PVL_SAPM % p = inputParser; -p.addRequired('E', @(x) all(isnumeric(x) & (x >= 0) & isvector(x))); -p.addRequired('E0', @(x) all(isnumeric(x) & isscalar(x) & (x >= 0))); -p.addRequired('a', @(x) all(isscalar(x) & isnumeric(x))); -p.addRequired('b', @(x) all(isscalar(x) & isnumeric(x))); -p.addRequired('windspeed', @(x) all(isvector(x) & isnumeric(x))); -p.addRequired('Tamb', @(x) all(isnumeric(x) & (x >= -273.15) & isvector(x))); -p.addRequired('deltaT', @(x) all(isnumeric(x) & (x >= 0) & isscalar(x))); +p.addRequired('E', @(x) isnumeric(x) && isvector(x) && all(x>=0 | isnan(x))); +p.addRequired('E0', @(x) isnumeric(x) && isscalar(x) && all(x>=0 | isnan(x))); +p.addRequired('a', @(x) isnumeric(x) && isscalar(x)); +p.addRequired('b', @(x) isnumeric(x) && isscalar(x)); +p.addRequired('windspeed', @(x) isnumeric(x) && isvector(x)); +p.addRequired('Tamb', @(x) isnumeric(x) && isvector(x) && all(x >= -273.15 | isnan(x))); +p.addRequired('deltaT', @(x) isnumeric(x) && isscalar(x) && all(x>=0 | isnan(x))); p.parse(E, E0, a, b, windspeed, Tamb, deltaT) -E = E(:); - -if isscalar(p.Results.windspeed) - windspeed = p.Results.windspeed*ones(size(E)); -else - windspeed = p.Results.windspeed(:); -end - -if isscalar(p.Results.Tamb) - Tamb = p.Results.Tamb*ones(size(E)); -else - Tamb = p.Results.Tamb(:); +if ~(isscalar(windspeed) || numel(windspeed) == numel(E)) + error('Input windspeed must be scalar or vectors of same length as E.'); end - -if (numel(E) ~= numel(windspeed)) || (numel(E) ~= numel(Tamb)) || (numel(E) ~= numel(windspeed)) - error(['Error in pvl_kingcelltemp. Inputs E, Tamb, and windspeed must '... - 'be scalars or vectors of same length.']); +if ~(isscalar(Tamb) || numel(Tamb) == numel(E)) + error('Input Tamb must be scalar or vectors of same length as E.'); end Tmodule = E.*(exp(a+b.*windspeed))+Tamb; -Tcell = Tmodule + E./E0.*deltaT; \ No newline at end of file +Tcell = Tmodule + E./E0.*deltaT; diff --git a/PVLIB/pvl_singleaxis.m b/PVLIB/pvl_singleaxis.m index 035ffe6..bcb723d 100755 --- a/PVLIB/pvl_singleaxis.m +++ b/PVLIB/pvl_singleaxis.m @@ -80,14 +80,14 @@ % % p = inputParser; -p.addRequired('SunZen', @(x) all(isnumeric(x) & x>=0 & x<=180) & isvector(x)); -p.addRequired('SunAz', @(x) all(isnumeric(x) & x>=0 & x<=360) & isvector(x)); -p.addRequired('Latitude', @(x) (isnumeric(x) & isscalar(x))); -p.addRequired('AxisTilt', @(x) (isnumeric(x) & isscalar(x) & x>=0 & x<=90)); -p.addRequired('AxisAzimuth', @(x) (isnumeric(x) & isscalar(x) & x>=0 & x<=360)); -p.addRequired('MaxAngle', @(x) (isnumeric(x) & isscalar(x) & x<=180 & x>=0)); -p.addOptional('BackTrack', 0, @(x) (isscalar(x))); -p.addOptional('GCR',1/3.5, @(x) (isnumeric(x) & isscalar(x) & x<=1)); +p.addRequired('SunZen', @(x) isnumeric(x) && isvector(x) && all((x>=0 & x<=180) | isnan(x))); +p.addRequired('SunAz', @(x) isnumeric(x) && isvector(x) && all((x>=0 & x<=360) | isnan(x))); +p.addRequired('Latitude', @(x) isnumeric(x) && isscalar(x)); +p.addRequired('AxisTilt', @(x) isnumeric(x) && isscalar(x) && x>=0 && x<=90); +p.addRequired('AxisAzimuth', @(x) isnumeric(x) && isscalar(x) && x>=0 && x<=360); +p.addRequired('MaxAngle', @(x) isnumeric(x) && isscalar(x) && x>=0 && x<=180); +p.addOptional('BackTrack', 0, @(x) isscalar(x)); +p.addOptional('GCR',1/3.5, @(x) isnumeric(x) & isscalar(x) & x<=1); p.parse(SunZen, SunAz, Latitude, AxisTilt, AxisAzimuth, MaxAngle, varargin{:}) if Latitude<0 diff --git a/PVLIB/pvl_singlediode.m b/PVLIB/pvl_singlediode.m index b41f1ce..f5970ac 100755 --- a/PVLIB/pvl_singlediode.m +++ b/PVLIB/pvl_singlediode.m @@ -91,15 +91,15 @@ % Vth is the thermal voltage of the cell in VOLTS % n is the usual diode ideality factor, assumed to be linear -% nVth is n*Vth +% nNsVth is n*Ns*Vth p = inputParser; -p.addRequired('IL',@(x) all(x>=0) & isnumeric(x) & isvector(x) ); -p.addRequired('I0', @(x) all(x>=0) & isnumeric(x) & isvector(x)); -p.addRequired('Rs', @(x) all(x>=0) & isnumeric(x) & isvector(x)); -p.addRequired('Rsh', @(x) all(x>=0) & isnumeric(x) & isvector(x)); -p.addRequired('nNsVth',@(x) all(x>=0) & isnumeric(x) & isvector(x) ); -p.addOptional('NumPoints', 0, @(x) isfinite(x) & isscalar(x)); +p.addRequired('IL',@(x) isnumeric(x) && isvector(x) && all(x>=0 | isnan(x))); +p.addRequired('I0', @(x) isnumeric(x) && isvector(x) && all(x>=0 | isnan(x))); +p.addRequired('Rs', @(x) isnumeric(x) && isvector(x) && all(x>=0 | isnan(x))); +p.addRequired('Rsh', @(x) isnumeric(x) && isvector(x) && all(x>=0 | isnan(x))); +p.addRequired('nNsVth',@(x) isnumeric(x) && isvector(x) && all(x>=0 | isnan(x))); +p.addOptional('NumPoints', 0, @(x) isfinite(x) && isscalar(x)); p.parse(IL, I0, Rs, Rsh, nNsVth, varargin{:}); % Make all inputs into column vectors @@ -129,30 +129,34 @@ nNsVth = nNsVth.*ones(MaxVectorSize , 1); end +Imax = zeros(MaxVectorSize, 1); +Pmp = zeros(MaxVectorSize, 1); +Vmax = zeros(MaxVectorSize, 1); +Ix = zeros(MaxVectorSize, 1); +Ixx = zeros(MaxVectorSize, 1); +Voc = zeros(MaxVectorSize, 1); +Isc = zeros(MaxVectorSize, 1); + +u = IL>0; + % Take care of any pesky non-integers by rounding up NumPoints = ceil(NumPoints); +% Find Isc using Lambert W +Isc(u) = I_from_V(Rsh(u), Rs(u), nNsVth(u), 0, I0(u), IL(u)); -% Set the tolerance of the fminbnd function to be 1e-8 instead of the -% default 1e-4 -defaultoptions = optimset('fminbnd'); -options = optimset(defaultoptions, 'TolX', 1e-8); +% Find Voc using Lambert W +Voc(u) = V_from_I(Rsh(u), Rs(u), nNsVth(u), 0, I0(u), IL(u)); -% Find Isc using Lambert W -Isc = I_from_V(Rsh, Rs, nNsVth, 0, I0, IL); +% Calculate I, V and P at the maximum power point -% Find Voc using Lambert W -Voc=V_from_I(Rsh, Rs, nNsVth, 0, I0, IL); +[Imax(u), Vmax(u), Pmp(u)]=calc_Pmp_bisect(IL(u),I0(u),nNsVth(u),Rs(u),Rsh(u)); -% Invert the Power-Current curve. Find the current where the inverted power -% is minimized. This is Imax. -[Imax negPmp ~] = pvl_fminbnd_vec(@(x) V_from_I(Rsh, Rs, nNsVth, x, I0, IL).*x.*-1, 0, Isc, options); -Vmax = -1*negPmp./Imax; % Find Ix and Ixx using Lambert W -Ix = I_from_V(Rsh, Rs, nNsVth, 0.5*Voc, I0, IL); -Ixx = I_from_V(Rsh, Rs, nNsVth, 0.5*(Voc+Vmax), I0, IL); +Ix(u) = I_from_V(Rsh(u), Rs(u), nNsVth(u), 0.5*Voc(u), I0(u), IL(u)); +Ixx(u) = I_from_V(Rsh(u), Rs(u), nNsVth(u), 0.5*(Voc(u)+Vmax(u)), I0(u), IL(u)); % If the user says they want a curve of with number of points equal to @@ -164,9 +168,12 @@ % curve. Thus the nth (V,I) point of curve m would be found as follows: % (Result.V(m,n),Result.I(m,n)). if NumPoints >= 2 + Result.I = zeros(MaxVectorSize, NumPoints); + Result.V = zeros(MaxVectorSize, NumPoints); s = ones(1,NumPoints); % shaping vector to shape the column vector parameters into 2-D matrices Result.V = (Voc)*(0:1/(NumPoints-1):1); - Result.I = I_from_V(Rsh*s, Rs*s, nNsVth*s, Result.V, I0*s, IL*s); + Result.I(u,:) = I_from_V(Rsh(u)*s, Rs(u)*s, nNsVth(u)*s, Result.V(u,:), I0(u)*s, IL(u)*s); + Result.I(:,end) = 0; % Make sure that I at Voc (the last point) is always 0 end % Wrap up the results into the Result struct @@ -175,26 +182,23 @@ Result.Imp = Imax; Result.Ix = Ix; Result.Ixx = Ixx; -Result.Pmp = -negPmp; +Result.Pmp = Pmp; Result.Isc = Isc; end - - - function [V]=V_from_I(Rsh, Rs, nNsVth, I, Io, Iphi) % calculates V from I per Eq 3 Jain and Kapoor 2004 -% uses Lambert W implemented in wapr_vec.m +% uses Lambert W implemented in pvl_lambertw.m % Rsh, nVth, I, Io, Iphi can all be vectors % Rs can be a vector, but should be a scalar % Generate the argument of the LambertW function argW = (Io.*Rsh./nNsVth) .* ... exp(Rsh.*(-I+Iphi+Io)./(nNsVth)); -inputterm = wapr_vec(argW); % Get the LambertW output +inputterm = pvl_lambertw(argW); % Get the LambertW output f = isnan(inputterm); % If argW is too big, the LambertW result will be NaN and we have to go to logspace % If it is necessary to go to logspace @@ -223,16 +227,119 @@ function [I]=I_from_V(Rsh, Rs, nNsVth, V, Io, Iphi) % calculates I from V per Eq 2 Jain and Kapoor 2004 -% uses Lambert W implemented in wapr_vec.m +% uses Lambert W implemented in pvl_lambertw.m % Rsh, nVth, V, Io, Iphi can all be vectors % Rs can be a vector, but should be a scalar argW = Rs.*Io.*Rsh.*... exp(Rsh.*(Rs.*(Iphi+Io)+V)./(nNsVth.*(Rs+Rsh)))./... (nNsVth.*(Rs + Rsh)); -inputterm = wapr_vec(argW); +inputterm = pvl_lambertw(argW); % Eqn. 4 in Jain and Kapoor, 2004 I = -V./(Rs + Rsh) - (nNsVth./Rs) .* inputterm + ... Rsh.*(Iphi + Io)./(Rs + Rsh); +end + +function [W] = calc_phi_exact(Imp, IL, Io, a, Rsh) + +% calculates W(phi) where phi is the argument of the +% Lambert W function in V = V(I) at I=Imp ([2], Eq. 3). Formula for +% phi is given in code below as argw. + +% phi +argw = Rsh.*Io./a .*exp(Rsh.*(IL + Io - Imp)./a); + +% Screen out any negative values for argw +u = argw>0; +W(~u)=NaN; + +tmp = pvl_lambertw(argw(u)); + +ff = isnan(tmp); + +% take care of any numerical overflow by evaluating log(W(phi)) +if any(ff) + logargW = log(Rsh(u)) + log(Io(u)) - log(a(u)) + Rsh(u).*(IL(u) + Io(u) - Imp(u))./a(u); + % Three iterations of Newton-Raphson method to solve w+log(w)=logargW. + % The initial guess is w=logargW. Where direct evaluation (above) results + % in NaN from overflow, 3 iterations of Newton's method gives + % approximately 8 digits of precision. + x = logargW; + for i=1:5 + x = x.*((1-log(x)+logargW)./(1+x)); + end; + tmp(ff) = x(ff); +end; + +W(u) = tmp; + +end + +function Imp=calc_Imp_bisect(Iph,Io,a,Rs,Rsh) + +% calculates the value of Imp (current at maximum power point) for an IV +% curve with parameters Iph, Io, a, Rs, Rsh. Imp is found as the value of +% I for which g(I)=dP/dV (I) = 0. + +% Set up lower and upper bounds on Imp +A = 0*Iph; +B = Iph+Io; + +% Detect when lower and upper bounds are not consistent with finding the +% zero of dP/dV + +gA=g(A,Iph,Io,a,Rs,Rsh); +gB=g(B,Iph,Io,a,Rs,Rsh); + +if any(gA.*gB>0) + % where gA*gB>0, then there is a problem with the IV curve parameters. + % In the event where gA and gB have the same sign, alert the user with + % a warning and replace erroneous cases with NaN + errorvalues = gA .* gB > 0; + warning(['Warning: pvl_singlediode has found at least one case where' ... + ' the single diode parameters are such that dP/dV may not have' ... + ' a zero. A NaN value has been reported for all such cases.']) + A(errorvalues) = NaN; % This will set Imp values where gA*gB>0 to NaN +end + +% midpoint is initial guess for Imp +p = (A+B)./2; +err = g(p,Iph,Io,a,Rs,Rsh); % value of dP/dV at initial guess p + +while max(abs(B-A))>1e-6 % set precision of estimate of Imp to 1e-6 (A) + gA=g(A,Iph,Io,a,Rs,Rsh); % value of dP/dV at left endpoint + u=(gA.*err)<0; + B(u)=p(u); + A(~u)=p(~u); + p=(A+B)/2; + err = g(p,Iph,Io,a,Rs,Rsh); +end; +Imp = p; + +end + +function y=g(I,Iph,Io,a,Rs,Rsh) + +% calculates dP/dV exactly, using p=I*V=I*V(I), where V=V(I) uses the +% Lambert's W function W(phi) ([2], Eq. 3). + +[z]=calc_phi_exact(I,Iph,Io,a,Rsh); % calculate W(phi) +z = z(:); + +% calculate dP/dV +y = (Iph+Io-2*I).*Rsh - 2*I.*Rs - a.*z +I.*Rsh.*z./(1+z); + +end + +function [Imp, Vmp, Pmp]=calc_Pmp_bisect(Iph,Io,a,Rs,Rsh) + +% Returns Imp, Vmp, Pmp for the IV curve described by input parameters. +% Vectorized. + +Imp=calc_Imp_bisect(Iph,Io,a,Rs,Rsh); % find Imp +[z]=calc_phi_exact(Imp,Iph,Io,a,Rsh); % Calculate W(phi) at Imp, where W is Lambert's W function and phi is its argument ([2], Eq. 3) +z = z(:); +Vmp=(Iph+Io-Imp).*Rsh - Imp.*Rs - a.*z; % Compute V from Imp and W(phi) +Pmp = Vmp.*Imp; end \ No newline at end of file diff --git a/PVLIB/pvl_snlinverter.m b/PVLIB/pvl_snlinverter.m index e5f6e8c..dec5ba0 100755 --- a/PVLIB/pvl_snlinverter.m +++ b/PVLIB/pvl_snlinverter.m @@ -16,9 +16,12 @@ % Inverter - A struct defining the inverter to be used, giving the % inverter performance parameters according to the Sandia % Grid-Connected Photovoltaic Inverter Model (SAND 2007-5036) [1]. A set of -% inverter performance parameters are provided with PV_LIB, or may be -% generated from a System Advisor Model (SAM) [2] library using the SAM -% library reader functions. Required struct components are: +% inverter performance parameters are provided with PV_LIB at +% /Required Data/, +% or may be generated from a System Advisor Model (SAM) [2] library +% using the SAM library reader functions. +% +% Required struct components are: % Inverter.Pac0 - AC-power output from inverter based on input power and voltage, (W) % Inverter.Pdc0 - DC-power input to inverter, typically assumed to be equal to the PV array maximum % power, (W) @@ -61,12 +64,12 @@ % [2] System Advisor Model web page. https://sam.nrel.gov. % % See also -% PVL_SAPM PVL_SAMLIBRARYREADER_SNLINVERTERS +% PVL_SAPM PVL_SAMLIBRARYREADER_SNLINVERTERS PVL_ADRINVERTER p = inputParser; p.addRequired('Inverter',@(x) isstruct(x)) -p.addRequired('Vdc', @(x) all(isnumeric(x) & x>=0 & isvector(x))); -p.addRequired('Pdc', @(x) all(isnumeric(x) & x>=0 & isvector(x))); +p.addRequired('Vdc', @(x) isnumeric(x) && isvector(x) && all(x>=0 | isnan(x))); +p.addRequired('Pdc', @(x) isnumeric(x) && isvector(x) && all(x>=0 | isnan(x))); p.parse(Inverter, Vdc, Pdc); Pac0 = p.Results.Inverter.Pac0;