# PVL_TestScript1

Example script to demonstrate the modeling of a PV system

This script demonstrates how to use the PV_Lib Toolbox to model the power output of a fixed-tilt photovoltaic array in Albuquerque, New Mexico, USA. The script uses irradiance and weather data that is read in from a Typical Meteorological Year 3 (TMY3) formatted file that contains one year of hourly-averaged data.

## Contents

- PV_LIB Functions used in this example
- Step 0: Define PV System
- 0.1 Define the PV module
- 0.2 Define the Inverter
- 0.3 Define the Array configuration
- 0.4 Define Additional Array Parameters
- Step 1: Read in Irradiance and Weather
- 1.1 Define Time and Irradiance Variabiles
- 1.2 Define the Site Location
- 1.3 Calculate Sun Position
- 1.4 Calculate Air Mass
- Step 2: Calculate Incident Radiation
- 2.1 Calculate Solar Angle of Incidence
- 2.2 Calculate Beam Radiation Component on Array
- 2.3 Calculate SkyDiffuse Radiation Component on Array
- 2.4 Calculate Ground Reflected Radiation Component on Array
- Step 3 Shading and Soiling
- Step 4: Calculate Cell Temperature
- Step 5: Calculate Module/Array IV Performance
- Step 6: DC and Mismatch Losses
- Step 7: DC to DC Max Power Point Tracking
- Step 8: DC to AC Conversion

## PV_LIB Functions used in this example

- `pvl_sapmmoduledb`
- `pvl_readtmy3`
- `pvl_maketimestruct`
- `pvl_makelocationstruct`
- `pvl_ephemeris`
- `pvl_absoluteairmass`
- `pvl_relativeairmass`
- `pvl_getaoi`
- `pvl_kingdiffuse`
- `pvl_grounddiffuse`
- `pvl_sapmcelltemp`
- `pvl_sapm`
- `pvl_snlinverter`

    clear all; %clear memory before running script
    close all; %close all files and figures

## Step 0: Define PV System

Before modeling can even begin, the PV system must be defined so that it can be properly represented by the model. This system we use for this example is outlined below:

- Location: Albuquerque, New Mexico USA
- Module: Candadian Solar CS5P-220M
- Inverter: PV Powered PVP2500
- # of Modules per String: 9
- # of Parallel strings: 2
- Array Tilt Angle: 35.04 degrees
- Array Azimuth Angle: 180 degrees*

*Azumuth angle convention is decimal degrees East of North. With North = 0 deg, East = 90 deg, South = 180 deg, and West = 270 deg.

## 0.1 Define the PV module

Define the PV module from the Sandia PV Module Database.

This step uses the PV_LIB function `pvl_sapmmoduledb` to extract a selected module from the excel version of the Sandia PV Module Database. The module selected is found on row 123 in that file.

    ModuleParameters = pvl_sapmmoduledb(123,'SandiaModuleDatabase_20120925.xlsx')

    ModuleParameters =

            name: 'Canadian Solar CS5P-220M'
         vintage: 2009
        material: 'c-Si'
            area: 1.7010
        AlphaIsc: 3.9700e-04
        AlphaImp: 1.8100e-04
            Isc0: 5.0911
            Imp0: 4.5463
            Voc0: 59.2608
            Vmp0: 48.3156
         BetaVoc: -0.2170
         BetaVmp: -0.2355
        mBetaVoc: 0
        mBetaVmp: 0
              Ns: 96
              Np: 1
            delT: 3
              fd: 1
               n: 1.4032
             Ix0: 4.9760
            Ixx0: 3.1880
          a_wind: -3.4064
          b_wind: -0.0842
               c: [1.0128 -0.0128 0.2793 -7.2446 0.9964 0.0036 1.1554 -0.1554]
               a: [-6.9304e-05 0.0017 -0.0158 0.0681 0.9284]
               b: [-1.3590e-09 2.1120e-07 -1.2460e-05 3.1030e-04 -0.0024 1]

## 0.2 Define the Inverter

Inverter database maintained by NREL System Advisor Model has been converted to a structure and saved in \Required Data folder. This .mat file can be generated from the latest SAM file using the function: `pvl_SAMLibraryReader_SNLInverters.m`.

    load('SandiaInverterDatabaseSAM2014.1.14.mat')
    % PV Powered PVP2500 inverter is #793 in the InverterNames cell array
    Inverter = SNLInverterDB(793)
    clear InverterNames SNLInverterDB

    Inverter =

               name: {'PV Powered: PVP2500 240V [CEC 2006]'}
                Vac: 240
               Pac0: 2500
               Pdc0: 2.6301e+03
               Vdc0: 219.2317
                Ps0: 41.3618
                 C0: -1.2650e-05
                 C1: 6.6747e-05
                 C2: 0.0017
                 C3: 6.1375e-04
                Pnt: 3.9000
             Vdcmax: 500
             Idcmax: 20
            MPPTLow: 140
             MPPTHi: 450
        LibraryType: {'SandiaInverter'}
        LibraryName: {'Sandia Inverters'}

## 0.3 Define the Array configuration

    Array.Tilt = 23; % Array tilt angle (deg)
    Array.Azimuth = 180; %Array azimuth (180 deg indicates array faces South)
    Array.Ms = 9; %Number of modules in series
    Array.Mp = 2; %Number of paralell strings

## 0.4 Define Additional Array Parameters

Because we will be using the module temperature model from the Sandia Photovoltaic Array Performance Model (SAPM), we will define the necessary parameters `a` and `b` here.

    Array.a = -3.56;
    Array.b = -0.075;

## Step 1: Read in Irradiance and Weather

Read in weather data from Albuquerque TMY3 file: `723650TY.csv`. This step uses the PV_LIB function `pvl_readtmy3` to read in data from this TMY3 formatted file.

    TMYData = pvl_readtmy3('723650TY.csv');

## 1.1 Define Time and Irradiance Variabiles

- Create time structure (Time) by using the PV_LIB function `pvl_maketimestruct`.
- Define variables for direct normal irradiance (DNI), diffuse horizontal irradiance (DHI), and global horizontal irradiance (GHI)

    TimeMatlab = TMYData.DateNumber;
    Time = pvl_maketimestruct(TimeMatlab, ones(size(TimeMatlab))*TMYData.SiteTimeZone);
    DNI = TMYData.DNI;
    DHI = TMYData.DHI;
    GHI = TMYData.GHI;

Let's examine irradiance for a sample day (August 2)

    figure
    tfilter = and(Time.month == 8,Time.day == 2);
    plot(Time.hour(tfilter),DNI(tfilter),'-*')
    hold all
    plot(Time.hour(tfilter),DHI(tfilter),'-o')
    plot(Time.hour(tfilter),GHI(tfilter),'-x')
    legend('DNI','DHI','GHI')
    xlabel('Hour of Day')
    ylabel('Irradiance (W/m^2)')
    title('Albuquerque TMY3 - Aug 2','FontSize',14)

## 1.2 Define the Site Location

Create a Location structure using the PV_LIB function `pvl_makelocationstruct`.

    Location = pvl_makelocationstruct(TMYData.SiteLatitude,TMYData.SiteLongitude,TMYData.SiteElevation) %Altitude is optional

    Location =

         latitude: 35.0400
        longitude: -106.6200
         altitude: 1619

## 1.3 Calculate Sun Position

Calculate sun position using PV_LIB function `pvl_ephemeris`. If Pressure and Temperature are specified they must be in units of Pa and deg C, respectively.

    PresPa = TMYData.Pressure*100; %Convert pressure from mbar to Pa
    [SunAz, SunEl, AppSunEl, SolarTime] = pvl_ephemeris(Time,Location,PresPa,TMYData.DryBulb);

Let's examine a plot of sun position for our site on August 2

    figure
    tfilter = and(Time.month == 8,Time.day == 2);
    plot(Time.hour(tfilter),90-AppSunEl(tfilter),'-s')
    hold all
    plot(Time.hour(tfilter),SunAz(tfilter),'-o')
    legend('Zenith angle','Azimuth Angle')
    xlabel('Hour of Day')
    ylabel('Angle (deg)')
    title('Albuquerque Sun Position - Aug 2','FontSize',14)

## 1.4 Calculate Air Mass

Air mass represents the optical path length of the sun through the Earth's atmosphere. It is primarily a function of the zenith angle but also varies with surface elevation and atmospheric pressure. It is calculated using the PV_LIB finctions `pvl_absoluteairmass` and `pvl_relativeairmass`.

    AMa = pvl_absoluteairmass(pvl_relativeairmass(90-AppSunEl),PresPa);

## Step 2: Calculate Incident Radiation

This step includes accounting for so

## 2.1 Calculate Solar Angle of Incidence

This step calculates the angle between direct beam from the sun to the array surface using the PV_LIB function `pvl_getaoi`.

    AOI = pvl_getaoi(Array.Tilt, Array.Azimuth, 90-AppSunEl, SunAz);

Let's examine a plot of sun angle of incidence on the array on August 2

    figure
    tfilter = and(Time.month == 8,Time.day == 2);
    plot(Time.hour(tfilter),AOI(tfilter),'-s')
    legend('AOI','Location','SE')
    xlabel('Hour of Day')
    ylabel('Angle (deg)')
    title('Albuquerque Angle of Incidence - Aug 2','FontSize',14)

## 2.2 Calculate Beam Radiation Component on Array

The beam component of the irradiance, `Eb` is simply:

    Eb = 0*AOI; %Initiallize variable
    Eb(AOI<90) = DNI(AOI<90).*cosd(AOI(AOI<90)); %Only calculate when sun is in view of the plane of array

## 2.3 Calculate SkyDiffuse Radiation Component on Array

There are many models that have been developed to estimate the amount of diffuse radiation that reaches a tlited plane from the sky `EdiffSky`. The PV_LIB Toolbox provides several of the more popular ones to choose from. They are listed below.

- `pvl_isotropicsky`
- `pvl_haydavies1980`
- `pvl_kingdiffuse`
- `pvl_reindl1990`

We will use `pvl_isotropicsky` for this example.

    EdiffSky = pvl_isotropicsky(Array.Tilt,DHI);

## 2.4 Calculate Ground Reflected Radiation Component on Array

If the array is tilted there will be a small amount of diffuse radiation reflecting off the ground surfaxce and hitting the plane of the array `EdiffGround`. We calculate this component using the PV_LIB function `pvl_grounddiffuse`. The ground reflectance or albedo is an input to this function.

    Albedo = 0.2;
    EdiffGround = pvl_grounddiffuse(Array.Tilt,GHI, Albedo);
    E = Eb + EdiffSky + EdiffGround; % Total incident irradiance (W/m^2)
    Ediff = EdiffSky + EdiffGround; % Total diffuse incident irradiance (W/m^2)

The total diffuse radiation reaching the array is calculated as the sum of the sky diffuse and the ground diffuse components. The total Plane of Array (POA) irradiance is the sum of `Eb` + `EdiffSky` + `EdiffGround`.

The following plot shows all of these components of August 2.

    figure
    tfilter = and(Time.month == 8,Time.day == 2);
    plot(Time.hour(tfilter),Eb(tfilter),'-s')
    hold all
    plot(Time.hour(tfilter),EdiffSky(tfilter),'-o')
    plot(Time.hour(tfilter),EdiffGround(tfilter),'-x')
    legend('Eb','EdiffSky','EdiffGround')
    xlabel('Hour of Day')
    ylabel('Irradiance (W/m^2')
    title('Albuquerque POA Irradiance Components - Aug 2','FontSize',14)

## Step 3 Shading and Soiling

Soiling reduces the amount of light reaching the array. For this example we will assume that soling levels are 2%

    SF=0.98;

## Step 4: Calculate Cell Temperature

Cell temperature is calculated using the PV_LIB function: `pvl_sapmcelltemp`, which is based on a method developed by King et al., 2004 (Eqs. 11 & 12).

    E0 = 1000; %Reference irradiance (1000 W/m^2)
    celltemp = pvl_sapmcelltemp(E, E0, Array.a, Array.b, TMYData.Wspd, TMYData.DryBulb, ModuleParameters.delT);

## Step 5: Calculate Module/Array IV Performance

For this example we use the Sandia PV Array Performance Model to represent the IV performance of the array at each time step. We can apply this model using the PV_LIB function: `pvl_sapm`. This function requires that we define the "effective irradiance" `Ee` as defined by King et al., 2004 in Eq. 7. `Ee` is esentially the POA irradiance corrected for spectral mismatch losses (`F1`) and angle of incidence losses (`F2`). The `F1` and `F2` functions are defined as 5th and 6th degree polynomials, respectively.

The output of this step is the IV characteristics of the array at each time step. The Sandia PV Array Performance Model used here only provides five points on the IV curve (Isc, Ix, Imp, Ixx, and Voc).

    F1 = max(0,polyval(ModuleParameters.a,AMa)); %Spectral loss function
    F2 = max(0,polyval(ModuleParameters.b,AOI)); % Angle of incidence loss function
    Ee = F1.*((Eb.*F2+ModuleParameters.fd.*Ediff)/E0)*SF; %Effective irradiance
    Ee(isnan(Ee))=0; % Set any NaNs to zero
    mSAPMResults = pvl_sapm(ModuleParameters, Ee, celltemp);
    aSAPMResults.Vmp = Array.Ms  *mSAPMResults.Vmp;
    aSAPMResults.Imp = Array.Mp  *mSAPMResults.Imp;
    aSAPMResults.Pmp = aSAPMResults.Vmp .* aSAPMResults.Imp;

The following plot shows the PV DC output predicted for the system on August 2.

    figure('OuterPosition',[100 100 600 800])
    tfilter = and(Time.month == 8,Time.day == 2);
    subplot(3,1,1)
    plot(Time.hour(tfilter),aSAPMResults.Pmp(tfilter),'-sr')
    legend('Pmp')
    %xlabel('Hour of Day')
    ylabel('DC Power (W)')
    title('Albuquerque DC Array Output - Aug 2','FontSize',14)
    subplot(3,1,2)
    plot(Time.hour(tfilter),aSAPMResults.Imp(tfilter),'-ob')
    legend('Imp')
    %xlabel('Hour of Day')
    ylabel('DC Current (A)')
    %title('Albuquerque Vmp and Imp - Aug 2','FontSize',12)
    subplot(3,1,3)
    plot(Time.hour(tfilter),aSAPMResults.Vmp(tfilter),'-xg')
    legend('Vmp')
    xlabel('Hour of Day')
    ylabel('DC Volatage (V)')
    %title('Albuquerque Vmp and Imp - Aug 2','FontSize',12)

## Step 6: DC and Mismatch Losses

In this step any losses between the PV array and the max power point tracker are included. In the present small system example, we will assume that these losses are neglegible. If you want to include such losses (e.g., wireing losses), you would have to module each DC circuit separately and adjust the current values.

## Step 7: DC to DC Max Power Point Tracking

In this step of the process the functionality of the maximum power point tracking (MPPT) is simulated. In most cases standard arrays, it is assumed that the inverter (with MPPT capability) is able to perfectly track the max power point (MPP) and thus there is no calculation required other than to calculate Pmp as was done in Step 5 above.

## Step 8: DC to AC Conversion

This step acconts for the power conversion of the inverter and results in the AC power of the PV system. In this example we use the `Sandia Performance Model for Grid-Connected Photovoltaic Inverters` (King et al. 2007), which is implemented in the PV_LIB function: `pvl_snlinverter`.

    ACPower = pvl_snlinverter(Inverter, mSAPMResults.Vmp*Array.Ms, mSAPMResults.Pmp*Array.Ms*Array.Mp);

The following plot compares the DC and AC power outputs predicted for the system on August 2.

    figure
    tfilter = and(Time.month == 8,Time.day == 2);
    plot(Time.hour(tfilter),aSAPMResults.Pmp(tfilter),'-sr')
    hold all
    plot(Time.hour(tfilter),ACPower(tfilter),'-ob')
    legend('DC Pmp', 'AC Power')
    %xlabel('Hour of Day')
    ylabel('Power (W)')
    title('Albuquerque DC and AC Power Output - Aug 2','FontSize',14)

Note that the inverter is undersized compared with the DC array and thus the inverter "clips" the AC power.

Copyright 2014 Sandia National Laboratories