Merge branch 'develop' into hybrids

shared/lib_irradproc.cpp

@@ -1692,6 +1692,38 @@ perez(double, double dn, double df, double alb, double inc, double tilt, double
     }
 }
 
+void ineichen(double clearsky_results[3], double apparent_zenith, double absolute_airmass, double linke_turbidity = 1.0, double altitude = 0.0, double dni_extra = 1364.0, bool perez_enhancement = false) {
+    double cos_zenith = Max(cosd(apparent_zenith), 0);
+    double tl = linke_turbidity;
+
+    double fh1 = exp(-altitude / 8000.0);
+    double fh2 = exp(-altitude / 1250.0);
+    double cg1 = 5.09e-5 * altitude + 0.868;
+    double cg2 = 3.92e-5 * altitude + 0.0387;
+
+    double ghi = exp(-cg2 * absolute_airmass * (fh1 + fh2 * (tl - 1)));
+    if (perez_enhancement) ghi *= exp(0.01 * pow(absolute_airmass, 1.8));
+
+    ghi = cg1 * dni_extra * cos_zenith * tl / tl * Max(ghi, 0);
+
+    double b = 0.664 + 0.163 / fh1;
+
+    double bnci = b * exp(-0.09 * absolute_airmass * (tl - 1));
+    bnci = dni_extra * Max(bnci, 0);
+
+    double bnci_2 = ((1 - (0.1 - 0.2 * exp(-tl)) / (0.1 + 0.882 / fh1)) / cos_zenith);
+    bnci_2 = ghi * Min(Max(bnci_2, 0), 1e20);
+
+    double dni = Min(bnci, bnci_2);
+
+    double dhi = ghi - dni * cos_zenith;
+    clearsky_results[0] = ghi;
+    clearsky_results[1] = dni;
+    clearsky_results[2] = dhi;
+    return;
+}
+
+
 void irrad::setup() {
     year = month = day = hour = -999;
     minute = delt = latitudeDegrees = longitudeDegrees = timezone = -999;
@@ -1719,6 +1751,7 @@ void irrad::setup() {
     poaRearDirectDiffuse = 0.;
     poaRearRowReflections = 0.;
     poaRearSelfShaded = 0.;
+
 }
 
 irrad::irrad() {
@@ -1733,7 +1766,7 @@ irrad::irrad(weather_record wf, weather_header hdr,
              double groundCoverageRatioIn, double slopeTiltIn, double slopeAzmIn, std::vector<double> monthlyTiltDegrees,
              std::vector<double> userSpecifiedAlbedo,
              poaDecompReq *poaAllIn,
-             bool useSpatialAlbedos, const util::matrix_t<double>* userSpecifiedSpatialAlbedos) :
+             bool useSpatialAlbedos, const util::matrix_t<double>* userSpecifiedSpatialAlbedos, bool enableSubhourlyClipping) :
         skyModel(skyModelIn), radiationMode(radiationModeIn), trackingMode(trackModeIn),
         enableBacktrack(backtrackingEnabled), forceToStow(forceToStowIn),
         delt(dtHour), tiltDegrees(tiltDegreesIn), surfaceAzimuthDegrees(azimuthDegreesIn),
@@ -1760,6 +1793,8 @@ irrad::irrad(weather_record wf, weather_header hdr,
     set_optional(hdr.elev, wf.pres, wf.tdry);
     set_sky_model(skyModel, albedo, albedoSpatial);
 
+    set_subhourly_clipping(enableSubhourlyClipping);
+
     if (radiationMode == irrad::DN_DF) set_beam_diffuse(wf.dn, wf.df);
     else if (radiationMode == irrad::DN_GH) set_global_beam(wf.gh, wf.dn);
     else if (radiationMode == irrad::GH_DF) set_global_diffuse(wf.gh, wf.df);
@@ -1852,6 +1887,17 @@ void irrad::get_poa(double *beam, double *skydiff, double *gnddiff,
     if (horizon != 0) *horizon = diffuseIrradianceFront[2];
 }
 
+void irrad::get_poa_clearsky(double* beam, double* skydiff, double* gnddiff,
+    double* isotrop, double* circum, double* horizon) {
+    if (beam != 0) *beam = planeOfArrayIrradianceFrontCS[0];
+    if (skydiff != 0) *skydiff = planeOfArrayIrradianceFrontCS[1];
+    if (gnddiff != 0) *gnddiff = planeOfArrayIrradianceFrontCS[2];
+    if (isotrop != 0) *isotrop = diffuseIrradianceFrontCS[0];
+    if (circum != 0) *circum = diffuseIrradianceFrontCS[1];
+    if (horizon != 0) *horizon = diffuseIrradianceFrontCS[2];
+}
+
+
 double irrad::get_poa_rear() {
     return planeOfArrayIrradianceRearAverage;
 }
@@ -1937,6 +1983,11 @@ void irrad::set_optional(double elev, double pres, double t_amb) //defaults of 0
         this->tamb = t_amb;
 }
 
+void irrad::set_subhourly_clipping(bool enable)
+{
+    if (enable) this->enableSubhourlyClipping = true;
+}
+
 void irrad::set_sky_model(int sm, double alb, const std::vector<double> &albSpatial) {
     this->skyModel = sm;
     this->albedo = alb;
@@ -2107,9 +2158,17 @@ int irrad::calc() {
         timeStepSunPosition[2] = 0;
     }
 
+    //clearsky
+    if (enableSubhourlyClipping) {
+        ineichen(clearskyIrradiance, RTOD * sunAnglesRadians[1], 1.5, 1.0, elevation);
+    }
+
 
     planeOfArrayIrradianceFront[0] = planeOfArrayIrradianceFront[1] = planeOfArrayIrradianceFront[2] = 0;
+    planeOfArrayIrradianceFrontCS[0] = planeOfArrayIrradianceFrontCS[1] = planeOfArrayIrradianceFrontCS[2] = 0;
     diffuseIrradianceFront[0] = diffuseIrradianceFront[1] = diffuseIrradianceFront[2] = 0;
+    diffuseIrradianceFrontCS[0] = diffuseIrradianceFrontCS[1] = diffuseIrradianceFrontCS[2] = 0;
+
     surfaceAnglesRadians[0] = surfaceAnglesRadians[1] = surfaceAnglesRadians[2] = surfaceAnglesRadians[3] = surfaceAnglesRadians[4] = 0;
 
     // do irradiance calculations if sun is up
@@ -2174,13 +2233,39 @@ int irrad::calc() {
                           diffuseIrradianceFront);
                     break;
             }
+
+            if (enableSubhourlyClipping) {
+                switch (skyModel) {
+                    case 0:
+                        isotropic(hextra, clearskyIrradiance[1], clearskyIrradiance[2], albedo,
+                            surfaceAnglesRadians[0], surfaceAnglesRadians[1], sunAnglesRadians[1],
+                            planeOfArrayIrradianceFrontCS, diffuseIrradianceFrontCS);
+                        break;
+                    case 1:
+                        hdkr(hextra, clearskyIrradiance[1], clearskyIrradiance[2], albedo, surfaceAnglesRadians[0],
+                            surfaceAnglesRadians[1], sunAnglesRadians[1], planeOfArrayIrradianceFrontCS,
+                            diffuseIrradianceFrontCS);
+                        break;
+                    default:
+                        perez(hextra, clearskyIrradiance[1], clearskyIrradiance[2], albedo, surfaceAnglesRadians[0],
+                            surfaceAnglesRadians[1], sunAnglesRadians[1], planeOfArrayIrradianceFrontCS,
+                            diffuseIrradianceFrontCS);
+                        break;
+                }
+            }
         }
         else { // Sev 2015/09/11 - perform a POA decomp.
             int errorcode = poaDecomp(weatherFilePOA, surfaceAnglesRadians, sunAnglesRadians, albedo, poaAll,
                                       directNormal, diffuseHorizontal, globalHorizontal, planeOfArrayIrradianceFront,
                                       diffuseIrradianceFront);
+            if (enableSubhourlyClipping) {
+                int errorcode_cs = poaDecomp(weatherFilePOA, surfaceAnglesRadians, sunAnglesRadians, albedo, poaAll,
+                    clearskyIrradiance[1], clearskyIrradiance[2], clearskyIrradiance[0], planeOfArrayIrradianceFrontCS,
+                    diffuseIrradianceFrontCS);
+            }
             calculatedDirectNormal = directNormal;
             calculatedDiffuseHorizontal = diffuseHorizontal;
+
             return errorcode; //this will return 0 if successful, otherwise 40, 41, or 42 if calculated decomposed dni, dhi, or ghi are negative
         }
     }

shared/lib_irradproc.h

@@ -759,7 +759,26 @@ void incidence(int mode, double tilt, double sazm, double rlim, double zen, doub
 * \param[out] diffc[2] horizon brightening
 */
 void perez(double hextra, double dn, double df, double alb, double inc, double tilt, double zen, double poa[3], double diffc[3] /* can be NULL */);
-
+/**
+* Ineichen model for calculating clearsky irradiance components
+* solar radiation + ground reflected radiation for a tilted surface and returns the total plane-of-array irradiance(poa),
+* see also isotropic(), hdkr().
+*
+* Function does not check all input for valid entries; consequently, this should be
+* done before calling the function.  (Reference: Perez et al, Solar Energy Vol. 44, No.5, pp.271-289,1990.)
+*
+* \param[in] apparent_zenith apparent solar zenith angle in degrees
+* \param[in] absolute_airmass pressure corrected airmass
+* \param[in] linke_turbidity Linke turbidity
+* \param[in] altitude site altitude in meters
+* \param[in] dni_extra extraterrestrial direct normal irradiance (W/m2)
+* \param[in] perez_enhancement is the perez enhcancement factor applied (y/n). Setting to true may produce spurious results for times when the sun is near the horizon and the airmass is high. 
+* \param[out] clear_sky_resuls calculated clearsky  irradiances (W/m2)
+* \param[out] clearsky_results[0] clear sky ghi
+* \param[out] clearsky_results[1] clear sky dni
+* \param[out] clearsky_results[2] clear sky dhi
+*/
+void ineichen(double clearsky_results[3], double apparent_zenith, double absolute_airmass, double linke_turbidity, double altitude, double dni_extra, bool perez_enhancement);
 /**
 * Isotropic sky model for diffuse irradiance on a tilted surface, see also perez(), hdkr().
 *
@@ -779,6 +798,7 @@ void perez(double hextra, double dn, double df, double alb, double inc, double t
 * \param[out] diffc[1] circumsolar diffuse
 * \param[out] diffc[2] horizon brightening
 */
+
 void isotropic(double hextra, double dn, double df, double alb, double inc, double tilt, double zen, double poa[3], double diffc[3] /* can be NULL */);
 
 /**
@@ -968,6 +988,9 @@ class irrad
     int year, month, day, hour;
     double minute, delt;
 
+    //Enable subhourly clipping correction
+    bool enableSubhourlyClipping;
+
     // Subarray properties
     double tiltDegrees;				///< Surface tilt of subarray in degrees
     double surfaceAzimuthDegrees;	///< Surface azimuth of subarray in degrees
@@ -1000,11 +1023,14 @@ class irrad
     double sunAnglesRadians[9];				///< Sun angles in radians calculated from solarpos()	
     double surfaceAnglesRadians[5];			///< Surface angles in radians calculated from incidence()
     double planeOfArrayIrradianceFront[3];	///< Front-side plane-of-array irradiance for beam, sky diffuse, ground diffuse (W/m2)
+    double planeOfArrayIrradianceFrontCS[3]; ///< Front-side plane-of-array clearsky irradiance for beam, sky diffuse, ground diffuse (W/m2)
     double planeOfArrayIrradianceRear[3];	///< Rear-side plane-of-array irradiance for beam, sky diffuse, ground diffuse (W/m2)
     double diffuseIrradianceFront[3];		///< Front-side diffuse irradiance for isotropic, circumsolar, and horizon (W/m2)
+    double diffuseIrradianceFrontCS[3];		///< Front-side diffuse clearsky irradiance for isotropic, circumsolar, and horizon (W/m2)
     double diffuseIrradianceRear[3];		///< Rear-side diffuse irradiance for isotropic, circumsolar, and horizon (W/m2)
     int timeStepSunPosition[3];				///< [0] effective hour of day used for sun position, [1] effective minute of hour used for sun position, [2] is sun up?  (0=no, 1=midday, 2=sunup, 3=sundown)
     double planeOfArrayIrradianceRearAverage; ///< Average rear side plane-of-array irradiance (W/m2)
+    double clearskyIrradiance[3];           /// [0] clearsky GHI, [1] clearsky DNI, [2] clearsky GHI from Ineichen model (W/m2);
     std::vector<double> planeOfArrayIrradianceRearSpatial;  ///< Spatial rear side plane-of-array irradiance (W/m2), where index 0 is at row bottom
     std::vector<double> groundIrradianceSpatial;            ///< Spatial irradiance incident on the ground in between rows, where index 0 is towards front of array
 
@@ -1029,7 +1055,7 @@ class irrad
         double dtHour, double tiltDegrees, double azimuthDegrees, double trackerRotationLimitDegrees, double stowAngleDegreesIn,
         double groundCoverageRatio, double slopeTilt, double slopeAzm, std::vector<double> monthlyTiltDegrees, std::vector<double> userSpecifiedAlbedo,
         poaDecompReq* poaAllIn,
-        bool useSpatialAlbedos = false, const util::matrix_t<double>* userSpecifiedSpatialAlbedos = nullptr);
+        bool useSpatialAlbedos = false, const util::matrix_t<double>* userSpecifiedSpatialAlbedos = nullptr, bool enableSubhourlyClipping = false);
 
     /// Construct the irrad class with an Irradiance_IO() object and Subarray_IO() object
     irrad();
@@ -1049,6 +1075,9 @@ class irrad
     // Set optional parameters for solarpos_spa calculation
     void set_optional(double elev = 0, double pres = 1013.25, double t_amb = 15);
 
+    //Set whether to use subhourly clipping model
+    void set_subhourly_clipping(bool enable = false);
+
     /// Set the sky model for the irradiance processor, using \link Irradiance_IO::SKYMODEL 
     void set_sky_model(int skymodel, double albedo, const std::vector<double> &albedoSpatial = std::vector<double>());
 
@@ -1105,6 +1134,10 @@ class irrad
     void get_poa(double* beam, double* skydiff, double* gnddiff,
         double* isotrop, double* circum, double* horizon);
 
+    /// Return the front-side plane-of-array irradiance and diffuse components of irradiation
+    void get_poa_clearsky(double* beam, double* skydiff, double* gnddiff,
+        double* isotrop, double* circum, double* horizon);
+
     /// Return the rear-side average total plane-of-array irradiance
     double get_poa_rear();
 

shared/lib_pv_io_manager.cpp

@@ -678,6 +678,12 @@ PVSystem_IO::PVSystem_IO(compute_module* cm, std::string cmName, Simulation_IO*
     }
 
     numberOfInverters = cm->as_integer("inverter_count");
+
+    dcNameplate = cm->as_double("system_capacity");
+    //numberOfInvertersClipping = cm->as_integer("num_inverter_subhourly_clipping");
+    numberOfInvertersClipping = dcNameplate / (Inverter->ratedACOutput / 1000);
+
+
     ratedACOutput = Inverter->ratedACOutput * numberOfInverters;
     acDerate = 1 - cm->as_double("acwiring_loss") / 100;
     acLossPercent = (1 - acDerate) * 100;
@@ -689,7 +695,7 @@ PVSystem_IO::PVSystem_IO(compute_module* cm, std::string cmName, Simulation_IO*
     enableSnowModel = cm->as_boolean("en_snow_model");
 
     // The shared inverter of the PV array and a tightly-coupled DC connected battery
-    std::unique_ptr<SharedInverter> tmpSharedInverter(new SharedInverter(Inverter->inverterType, numberOfInverters, &Inverter->sandiaInverter, &Inverter->partloadInverter, &Inverter->ondInverter));
+    std::unique_ptr<SharedInverter> tmpSharedInverter(new SharedInverter(Inverter->inverterType, numberOfInverters, &Inverter->sandiaInverter, &Inverter->partloadInverter, &Inverter->ondInverter, numberOfInvertersClipping));
     m_sharedInverter = std::move(tmpSharedInverter);
 
     // Register shared inverter with inverter_IO
@@ -859,6 +865,10 @@ void PVSystem_IO::AllocateOutputs(compute_module* cm)
             p_derateSelfShading.push_back(cm->allocate(prefix + "ss_derate", numberOfWeatherFileRecords));
             p_derateSelfShadingDiffuse.push_back(cm->allocate(prefix + "ss_diffuse_derate", numberOfWeatherFileRecords));
             p_derateSelfShadingReflected.push_back(cm->allocate(prefix + "ss_reflected_derate", numberOfWeatherFileRecords));
+            p_DNIIndex.push_back(cm->allocate(prefix + "dni_index", numberOfWeatherFileRecords));
+            p_poaBeamFrontCS.push_back(cm->allocate(prefix + "poa_beam_front_cs", numberOfWeatherFileRecords));
+            p_poaDiffuseFrontCS.push_back(cm->allocate(prefix + "poa_diffuse_front_cs", numberOfWeatherFileRecords));
+            p_poaGroundFrontCS.push_back(cm->allocate(prefix + "poa_ground_front_cs", numberOfWeatherFileRecords));
 
             if (enableSnowModel) {
                 p_snowLoss.push_back(cm->allocate(prefix + "snow_loss", numberOfWeatherFileRecords));
@@ -921,13 +931,18 @@ void PVSystem_IO::AllocateOutputs(compute_module* cm)
 
     p_inverterACOutputPreLoss = cm->allocate("ac_gross", numberOfWeatherFileRecords);
     p_acWiringLoss = cm->allocate("ac_wiring_loss", numberOfWeatherFileRecords);
+    p_ClippingPotential = cm->allocate("clipping_potential", numberOfWeatherFileRecords);
     p_transmissionLoss = cm->allocate("ac_transmission_loss", numberOfWeatherFileRecords);
     p_acPerfAdjLoss = cm->allocate("ac_perf_adj_loss", numberOfWeatherFileRecords);
     p_acLifetimeLoss = cm->allocate("ac_lifetime_loss", numberOfWeatherFileRecords);
     p_dcLifetimeLoss = cm->allocate("dc_lifetime_loss", numberOfWeatherFileRecords);
     p_systemDCPower = cm->allocate("dc_net", numberOfLifetimeRecords);
     p_systemACPower = cm->allocate("gen", numberOfLifetimeRecords);
 
+    p_systemDCPowerCS = cm->allocate("dc_net_clearsky", numberOfLifetimeRecords);
+    p_subhourlyClippingLoss = cm->allocate("subhourly_clipping_loss", numberOfLifetimeRecords);
+    p_subhourlyClippingLossFactor = cm->allocate("subhourly_clipping_loss_factor", numberOfLifetimeRecords);
+
     if (Simulation->useLifetimeOutput)
     {
         p_dcDegradationFactor = cm->allocate("dc_degrade_factor", numberOfYears);

shared/lib_pv_io_manager.h

@@ -283,6 +283,8 @@ struct PVSystem_IO
 
 	size_t numberOfSubarrays;
 	size_t numberOfInverters;
+    size_t numberOfInvertersClipping;
+    double dcNameplate;
 
 	Irradiance_IO * Irradiance;
 	Simulation_IO * Simulation;
@@ -346,6 +348,11 @@ struct PVSystem_IO
     std::vector<ssc_number_t *> p_poaRearSpatial;
     std::vector<ssc_number_t *> p_groundRear;
 
+    std::vector<ssc_number_t*> p_poaBeamFrontCS;
+    std::vector<ssc_number_t*> p_poaDiffuseFrontCS;
+    std::vector<ssc_number_t*> p_poaGroundFrontCS;
+    std::vector<ssc_number_t*> p_DNIIndex;
+
 	// MPPT level outputs
 	std::vector<ssc_number_t *> p_mpptVoltage; /// An output vector containing input DC voltage in V to each mppt input
 	std::vector<ssc_number_t *> p_dcPowerNetPerMppt; /// An output vector containing Net DC Power in W for each mppt input
@@ -409,7 +416,12 @@ struct PVSystem_IO
     ssc_number_t *p_dcLifetimeLoss; // kWdc
 
 	ssc_number_t *p_systemDCPower; // kWdc
+    ssc_number_t* p_systemDCPowerCS; // kWdc
 	ssc_number_t *p_systemACPower; // kWac
+
+    ssc_number_t *p_subhourlyClippingLoss;
+    ssc_number_t* p_subhourlyClippingLossFactor;
+    ssc_number_t* p_ClippingPotential;
 };
 
 /**
@@ -493,6 +505,9 @@ struct Subarray_IO
 		double poaBeamFront;	/// POA due to beam irradiance on the front of the subarray [W/m2]
 		double poaDiffuseFront; /// POA due to diffuse irradiance on the front of the subarray [W/m2]
 		double poaGroundFront;  /// POA due to ground reflection on the front of the subarray [W/m2]
+        double poaBeamFrontCS;	/// POA due to clearsky beam irradiance on the front of the subarray [W/m2]
+        double poaDiffuseFrontCS; /// POA due to clearsky diffuse irradiance on the front of the subarray [W/m2]
+        double poaGroundFrontCS;  /// POA due to clearsky ground reflection on the front of the subarray [W/m2]
 		double poaRear;			/// POA total irradiance on the back of the subarray if bifacial modules [W/m2]
 		double poaTotal;		/// POA total of front and rear side of array [W/m2]
 		bool sunUp;				/// Flag indicating whether the sun is up or not
@@ -507,6 +522,7 @@ struct Subarray_IO
 
 	//calculated- subarray power
 	double dcPowerSubarray; /// DC power for this subarray [W]
+    double dcPowerSubarrayCS;
 
 };
 
@@ -567,7 +583,7 @@ struct Module_IO
 	double temperatureCellCelcius; /// The weighted moving average  cell temperature of the module [C]
 	double temperatureCellCelciusSS; /// The SS average cell temperature of the module [C]
 	double angleOfIncidenceModifier; /// The angle of incidence modifier on the total poa front-side irradiance [0-1]
-
+    double dcPowerWCS;
 };