Merge branch 'develop' into getem-wf-temp

nlopt/timer.c

@@ -43,7 +43,9 @@
 #  include <windows.h>    
 #endif
 
+// header needed to build on MacOS "Call to undeclared function 'gettimeofday'; ISO C99 and later do not support implicit function declarations" compile error
 int     gettimeofday(struct timeval * __restrict, void * __restrict);
+
 /* return time in seconds since some arbitrary point in the past */
 double nlopt_seconds(void)
 {

shared/lib_battery_dispatch.cpp

@@ -681,7 +681,7 @@ dispatch_automatic_t::dispatch_automatic_t(
 
     _mode = dispatch_mode;
     _weather_forecast_mode = weather_forecast_mode;
-    _safety_factor = 0.03;
+    _safety_factor = 0.0;
 
 	m_batteryPower->canClipCharge = can_clip_charge;
 	m_batteryPower->canSystemCharge = can_charge;

shared/lib_battery_dispatch.h

@@ -57,7 +57,7 @@ class dispatch_t
 public:
 
 	enum FOM_MODES { FOM_AUTOMATED_ECONOMIC, FOM_PV_SMOOTHING, FOM_CUSTOM_DISPATCH, FOM_MANUAL };
-	enum BTM_MODES { PEAK_SHAVING, MAINTAIN_TARGET, CUSTOM_DISPATCH, MANUAL, FORECAST };
+	enum BTM_MODES { PEAK_SHAVING, MAINTAIN_TARGET, CUSTOM_DISPATCH, MANUAL, RETAIL_RATE, SELF_CONSUMPTION };
 	enum METERING { BEHIND, FRONT };
     enum WEATHER_FORECAST_CHOICE { WF_LOOK_AHEAD, WF_LOOK_BEHIND, WF_CUSTOM };
     enum LOAD_FORECAST_CHOICE { LOAD_LOOK_AHEAD, LOAD_LOOK_BEHIND, LOAD_CUSTOM };
@@ -182,7 +182,7 @@ class dispatch_t
 
 	/**
 	The dispatch mode.
-	For behind-the-meter dispatch: 0 = PEAK_SHAVING, 1 = MAINTAIN_TARGET, 2 = CUSTOM, 3 = MANUAL, 4 = FORECAST
+	For behind-the-meter dispatch: 0 = PEAK_SHAVING, 1 = MAINTAIN_TARGET, 2 = CUSTOM, 3 = MANUAL, 4 = FORECAST, 5 = SELF_CONSUMPTION
 	For front-of-meter dispatch: 0 = FOM_AUTOMATED_ECONOMIC, 1 = FOM_PV_SMOOTHING, 2 = FOM_CUSTOM_DISPATCH, 3 = FOM_MANUAL
 	*/
 	int _mode;

shared/lib_battery_dispatch_automatic_btm.cpp

@@ -70,7 +70,8 @@ dispatch_automatic_behind_the_meter_t::dispatch_automatic_behind_the_meter_t(
     bool chargeOnlySystemExceedLoad,
     bool dischargeOnlyLoadExceedSystem,
     bool behindTheMeterDischargeToGrid,
-    double SOC_min_outage
+    double SOC_min_outage,
+    int load_forecast_mode
 	) : dispatch_automatic_t(Battery, dt_hour, SOC_min, SOC_max, current_choice, Ic_max, Id_max, Pc_max_kwdc, Pd_max_kwdc, Pc_max_kwac, Pd_max_kwac,
 		t_min, dispatch_mode, weather_forecast_mode, pv_dispatch, nyears, look_ahead_hours, dispatch_update_frequency_hours, can_charge, can_clip_charge, can_grid_charge, can_fuelcell_charge,
         battReplacementCostPerkWh, battCycleCostChoice, battCycleCost, battOMCost, interconnection_limit, chargeOnlySystemExceedLoad, dischargeOnlyLoadExceedSystem,
@@ -81,6 +82,8 @@ dispatch_automatic_behind_the_meter_t::dispatch_automatic_behind_the_meter_t(
 	_P_target_use.reserve(_num_steps);
 	_P_battery_use.reserve(_num_steps);
 
+    _load_forecast_mode = load_forecast_mode;
+
 	grid.reserve(_num_steps);
 	sorted_grid.reserve(_num_steps);
 
@@ -170,7 +173,7 @@ double dispatch_automatic_behind_the_meter_t::power_grid_target() { return _P_ta
 
 void dispatch_automatic_behind_the_meter_t::setup_rate_forecast()
 {
-    if (_mode == dispatch_t::FORECAST)
+    if (_mode == dispatch_t::RETAIL_RATE)
     {
 
         forecast_setup rate_setup(_steps_per_hour, _nyears);
@@ -193,7 +196,7 @@ void dispatch_automatic_behind_the_meter_t::update_dispatch(size_t year, size_t
     // [kWh] - the maximum energy that can be cycled
     double E_max = 0;
 
-    if (_mode == dispatch_t::FORECAST)
+    if (_mode == dispatch_t::RETAIL_RATE)
     {
         // Hourly rolling forecast horizon
         if ((hour_of_year != _hour_last_updated) || m_outage_manager->recover_from_outage)
@@ -264,6 +267,20 @@ void dispatch_automatic_behind_the_meter_t::initialize(size_t hour_of_year, size
 	m_batteryPower->powerBatteryTarget = 0;
     _day_index = 0;
 
+    // Give all algorithms real data for the current step (overwrite forecast if needed)
+    if (_load_forecast_mode != LOAD_LOOK_AHEAD) {
+        _P_load_ac[lifetimeIndex] = m_batteryPower->powerLoad;
+    }
+    // Lookahead forecasts may better account for losses than the code below, so don't run this if lookahead
+    if (_weather_forecast_mode != WF_LOOK_AHEAD) {
+        if (m_batteryPower->connectionMode == AC_CONNECTED) {
+            _P_pv_ac[lifetimeIndex] = m_batteryPower->powerSystem;
+        }
+        else {
+            _P_pv_ac[lifetimeIndex] = m_batteryPower->powerSystem * m_batteryPower->sharedInverter->efficiencyAC;
+        }
+    }
+
 	// clean up vectors
     size_t lifetimeMax = _P_pv_ac.size();
 	for (size_t ii = 0; ii != _num_steps && lifetimeIndex < lifetimeMax; ii++)
@@ -274,6 +291,7 @@ void dispatch_automatic_behind_the_meter_t::initialize(size_t hour_of_year, size
 		_P_battery_use.push_back(0.);
         lifetimeIndex++;
 	}
+
 }
 bool dispatch_automatic_behind_the_meter_t::check_new_month(size_t hour_of_year, size_t step)
 {
@@ -320,6 +338,7 @@ void dispatch_automatic_behind_the_meter_t::sort_grid(size_t idx, FILE *p, const
 
 	// compute grid net from pv and load (no battery)
 	size_t count = 0;
+
 	for (size_t hour = 0; hour != 24 && idx < _P_load_ac.size(); hour++)
 	{
 		for (size_t step = 0; step != _steps_per_hour; step++)
@@ -447,7 +466,7 @@ void dispatch_automatic_behind_the_meter_t::target_power(double E_useful, size_t
 				if (sorted_grid[ii].Grid() > P_target_min)
 					break;
 
-				E_charge += (P_target_min - sorted_grid[ii].Grid())*_dt_hour;
+				E_charge += (P_target_min - sorted_grid[ii].Grid()) * _dt_hour * m_batteryPower->singlePointEfficiencyACToDC;
 			}
 			E_charge_vec.push_back(E_charge);
 			if (debug)
@@ -864,7 +883,7 @@ void dispatch_automatic_behind_the_meter_t::costToCycle()
             m_cycleCost = 0.01 * capacityPercentDamagePerCycle * m_battReplacementCostPerKWH[curr_year] * _Battery->get_params().nominal_energy;
         }
         else {
-            // Should only apply to BattWatts. BattWatts doesn't have price signal dispatch, so this is fine.
+            // Should only apply to BattWatts. BattWatts doesn't have retal rate dispatch, so this is fine.
             m_cycleCost = 0.0;
         }
     }

shared/lib_battery_dispatch_automatic_btm.h

@@ -37,7 +37,7 @@ OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 #include "lib_utility_rate.h"
 
 /*
- * Data for price signal dispatch (FORECAST) to compare dispatch plans in the cost_based_target_power function
+ * Data for retail rate dispatch to compare dispatch plans in the cost_based_target_power function
  */
 struct dispatch_plan
 {
@@ -95,7 +95,8 @@ class dispatch_automatic_behind_the_meter_t : public dispatch_automatic_t
         bool chargeOnlySystemExceedLoad,
         bool dischargeOnlyLoadExceedSystem,
         bool behindTheMeterDischargeToGrid,
-        double SOC_min_outage
+        double SOC_min_outage,
+        int load_forecast_mode
 		);
 
 	~dispatch_automatic_behind_the_meter_t() override {};
@@ -153,7 +154,7 @@ class dispatch_automatic_behind_the_meter_t : public dispatch_automatic_t
     void target_power(double E_max, size_t idx, FILE* p = NULL, bool debug = false);
     void apply_target_power(size_t day_index);
 
-    /*! Functions used by price signal dispatch */
+    /*! Functions used by retail rate dispatch */
     double compute_costs(size_t idx, size_t year, size_t hour_of_year, FILE* p = NULL, bool debug = false); // Initial computation of no-dispatch costs, assigned hourly to grid points
     void cost_based_target_power(size_t idx, size_t year, size_t hour_of_year, double no_dispatch_cost, double E_max, FILE* p = NULL, const bool debug = false); // Optimizing loop, runs twelve possible dispatch scenarios
     void plan_dispatch_for_cost(dispatch_plan& plan, size_t idx, double E_max, double startingEnergy); // Generates each dispatch plan (input argument)
@@ -166,6 +167,9 @@ class dispatch_automatic_behind_the_meter_t : public dispatch_automatic_t
 	/*! Full time-series of loads [kW] */
 	double_vec _P_load_ac;
 
+    /*! Forecast mode for loads (look ahead, look behind, custom) */
+    int _load_forecast_mode;
+
 	/*! Full time-series of target power [kW] */
 	double_vec _P_target_input;
 

shared/lib_battery_powerflow.cpp

@@ -351,6 +351,11 @@ void BatteryPowerFlow::calculateACConnected()
 
     // Code simplification to remove redundancy for code that should use either critical load or actual load
     double calc_load_ac = (m_BatteryPower->isOutageStep ? P_crit_load_ac : P_load_ac);
+    double P_required_for_load = calc_load_ac;
+
+    if (ac_loss_percent_post_battery < 1) { // Account for possible divide by zero
+        P_required_for_load /= (1 - ac_loss_percent_post_battery);
+    }
 
     // charging and idle
     if (P_battery_ac <= 0)
@@ -365,11 +370,11 @@ void BatteryPowerFlow::calculateACConnected()
         if (m_BatteryPower->chargeOnlySystemExceedLoad) {
             P_pv_to_load_ac = P_pv_ac;
 
-            if (P_pv_to_load_ac > calc_load_ac) {
-                P_pv_to_load_ac = calc_load_ac;
+            if (P_pv_to_load_ac > P_required_for_load) {
+                P_pv_to_load_ac = P_required_for_load;
             }
             // Fuel cell goes to load next
-            P_fuelcell_to_load_ac = std::fmin(calc_load_ac - P_pv_to_load_ac, P_fuelcell_ac);
+            P_fuelcell_to_load_ac = std::fmin(P_required_for_load - P_pv_to_load_ac, P_fuelcell_ac);
         }
 
         // Excess PV can go to battery, if PV can cover charging losses
@@ -398,11 +403,11 @@ void BatteryPowerFlow::calculateACConnected()
                 P_pv_to_load_ac = P_pv_ac - P_pv_to_batt_ac;
             }
 
-            if (P_pv_to_load_ac > calc_load_ac) {
-                P_pv_to_load_ac = calc_load_ac;
+            if (P_pv_to_load_ac > P_required_for_load) {
+                P_pv_to_load_ac = P_required_for_load;
             }
             // Fuel cell goes to load next
-            P_fuelcell_to_load_ac = std::fmin(calc_load_ac - P_pv_to_load_ac, P_fuelcell_ac);
+            P_fuelcell_to_load_ac = std::fmin(P_required_for_load - P_pv_to_load_ac, P_fuelcell_ac);
         }
 
         // Fuelcell can also charge battery
@@ -439,6 +444,7 @@ void BatteryPowerFlow::calculateACConnected()
     // discharging, not idle
     else
     {
+
         // Test if battery is discharging erroneously
         if (!m_BatteryPower->canDischarge && P_battery_ac > 0) {
             P_batt_to_grid_ac = P_batt_to_load_ac = 0;
@@ -448,9 +454,9 @@ void BatteryPowerFlow::calculateACConnected()
             P_pv_to_load_ac = P_pv_ac;
 
             // Excess PV production, no other component meets load
-            if (P_pv_ac >= calc_load_ac)
+            if (P_pv_ac >= P_required_for_load)
             {
-                P_pv_to_load_ac = calc_load_ac;
+                P_pv_to_load_ac = P_required_for_load;
                 P_fuelcell_to_load_ac = 0;
                 P_batt_to_load_ac = 0;
 
@@ -459,14 +465,14 @@ void BatteryPowerFlow::calculateACConnected()
                 P_fuelcell_to_grid_ac = P_fuelcell_ac;
             }
             else {
-                P_fuelcell_to_load_ac = std::fmin(P_fuelcell_ac, calc_load_ac - P_pv_to_load_ac);
-                P_batt_to_load_ac = std::fmin(P_battery_ac - P_system_loss_ac, calc_load_ac - P_pv_to_load_ac - P_fuelcell_to_load_ac);
+                P_fuelcell_to_load_ac = std::fmin(P_fuelcell_ac, P_required_for_load - P_pv_to_load_ac);
+                P_batt_to_load_ac = std::fmin(P_battery_ac - P_system_loss_ac, P_required_for_load - P_pv_to_load_ac - P_fuelcell_to_load_ac);
             }
         }
         else {
-            P_batt_to_load_ac = std::fmin(P_battery_ac, calc_load_ac);
-            P_fuelcell_to_load_ac = std::fmin(P_fuelcell_ac, calc_load_ac - P_batt_to_load_ac);
-            P_pv_to_load_ac = std::fmin(std::fmax(0, calc_load_ac - P_fuelcell_to_load_ac - P_batt_to_load_ac), P_pv_ac);
+            P_batt_to_load_ac = std::fmin(P_battery_ac, P_required_for_load);
+            P_fuelcell_to_load_ac = std::fmin(P_fuelcell_ac, P_required_for_load - P_batt_to_load_ac);
+            P_pv_to_load_ac = std::fmin(std::fmax(0, P_required_for_load - P_fuelcell_to_load_ac - P_batt_to_load_ac), P_pv_ac);
             P_pv_to_grid_ac = std::fmax(0, P_pv_ac - P_pv_to_load_ac);
             P_fuelcell_to_grid_ac = std::fmax(0, P_fuelcell_ac - P_fuelcell_to_load_ac);
         }
@@ -483,11 +489,8 @@ void BatteryPowerFlow::calculateACConnected()
             P_fuelcell_to_grid_ac = 0;
         }
 
+        // Preliminary batt to grid for DC losses
         P_batt_to_grid_ac = P_battery_ac - P_system_loss_ac - P_batt_to_load_ac - P_batt_to_pv_inverter;
-        if (m_BatteryPower->isOutageStep && P_batt_to_grid_ac > tolerance) {
-            m_BatteryPower->powerBatteryDC = (P_battery_ac - P_batt_to_grid_ac) / m_BatteryPower->singlePointEfficiencyDCToAC;
-            return calculateACConnected();
-        }
 
         P_fuelcell_to_grid_ac = P_fuelcell_ac - P_fuelcell_to_load_ac;
         P_batt_to_system_loss = P_system_loss_ac;
@@ -503,15 +506,23 @@ void BatteryPowerFlow::calculateACConnected()
     // Compute total system output and grid power flow
     P_gen_ac = P_pv_ac + P_fuelcell_ac + P_inverter_draw_ac + P_battery_ac - P_system_loss_ac;
 
+    // Final batt to grid for outage accounting
+    if (P_battery_ac > 0)
+    {
+        P_batt_to_grid_ac = P_battery_ac * (1 - ac_loss_percent_post_battery) - P_system_loss_ac - P_batt_to_load_ac - P_batt_to_pv_inverter;
+        if (m_BatteryPower->isOutageStep && P_batt_to_grid_ac > tolerance) {
+            m_BatteryPower->powerBatteryDC = (P_battery_ac - P_batt_to_grid_ac) / m_BatteryPower->singlePointEfficiencyDCToAC;
+            return calculateACConnected();
+        }
+    }
+
     // Apply AC losses to powerflow - note that these are applied to gen later
     P_pv_to_batt_ac *= (1 - ac_loss_percent_post_battery);
     P_pv_to_load_ac *= (1 - ac_loss_percent_post_battery);
     P_pv_to_batt_ac *= (1 - ac_loss_percent_post_battery);
     P_pv_to_grid_ac *= (1 - ac_loss_percent_post_battery);
     P_grid_to_batt_ac *= (1 - ac_loss_percent_post_battery);
-    P_grid_to_load_ac *= (1 - ac_loss_percent_post_battery);
     P_batt_to_load_ac *= (1 - ac_loss_percent_post_battery);
-    P_batt_to_grid_ac *= (1 - ac_loss_percent_post_battery);
     P_fuelcell_to_batt_ac *= (1 - ac_loss_percent_post_battery);
     P_fuelcell_to_load_ac *= (1 - ac_loss_percent_post_battery);
     P_fuelcell_to_grid_ac *= (1 - ac_loss_percent_post_battery);

shared/lib_battery_voltage.cpp

@@ -123,6 +123,9 @@ void voltage_table_t::initialize() {
             double V0 = params->voltage_table[i - 1][1];
             slope = (V - V0) / (DOD - DOD0);
             intercept = V0 - (slope * DOD0);
+
+            if (fabs(slope) < 1e-7)
+                throw std::runtime_error("voltage_table_t error: Battery voltage matrix cannot have two identical voltages.");
         }
         slopes.emplace_back(slope);
         intercepts.emplace_back(intercept);
@@ -178,31 +181,35 @@ voltage_t *voltage_table_t::clone() {
     return new voltage_table_t(*this);
 }
 
-double voltage_table_t::calculate_voltage(double DOD) {
+double voltage_table_t::calculate_voltage(double DOD, double I) {
     DOD = fmax(0., DOD);
     DOD = fmin(DOD, 100.);
 
     size_t row = 0;
     while (row < params->voltage_table.size() && DOD > params->voltage_table[row][0]) {
         row++;
     }
+    //
+    if (DOD < tolerance || DOD > 100. - tolerance) {
+        I = 0.0; // At full or empty, current must go to zero
+    }
 
-    return fmax(slopes[row] * DOD + intercepts[row], 0);
+    return fmax(slopes[row] * DOD + intercepts[row], 0) - I * params->resistance;
 }
 
 void voltage_table_t::set_initial_SOC(double init_soc) {
-    state->cell_voltage = calculate_voltage(100. - init_soc);
+    state->cell_voltage = calculate_voltage(100. - init_soc, 0.0);
 }
 
 double voltage_table_t::calculate_voltage_for_current(double I, double q, double qmax, double) {
     double DOD = (q - I * params->dt_hr) / qmax * 100.;
-    return calculate_voltage(DOD) * params->num_cells_series;
+    return calculate_voltage(DOD, I / params->num_strings) * params->num_cells_series;
 }
 
 
-void voltage_table_t::updateVoltage(double q, double qmax, double, const double, double) {
+void voltage_table_t::updateVoltage(double q, double qmax, double I, const double, double) {
     double DOD = 100. * (1 - q / qmax);
-    state->cell_voltage = calculate_voltage(DOD);
+    state->cell_voltage = calculate_voltage(DOD, I / params->num_strings);
 }
 
 // helper fx to calculate depth of discharge from current and max capacities
@@ -212,7 +219,7 @@ double voltage_table_t::calculate_max_charge_w(double q, double qmax, double, do
     double current = (q - qmax) / params->dt_hr;
     if (max_current)
         *max_current = current;
-    return calculate_voltage(0.) * current * params->num_cells_series;
+    return calculate_voltage(0., current / params->num_strings) * current * params->num_cells_series;
 }
 
 double voltage_table_t::calculate_max_discharge_w(double q, double qmax, double, double *max_current) {
@@ -227,7 +234,7 @@ double voltage_table_t::calculate_max_discharge_w(double q, double qmax, double,
         dod = fmin(100, dod);
         dod = fmax(0, dod);
         double current = qmax * ((1. - DOD0 / 100.) - (1. - dod / 100.)) / params->dt_hr;
-        double p = calculate_voltage(dod) * current;
+        double p = calculate_voltage(dod, current / params->num_strings) * current;
         if (p > max_P) {
             max_P = p;
             max_I = current;
@@ -253,9 +260,9 @@ double voltage_table_t::calculate_current_for_target_w(double P_watts, double q,
 
     P_watts /= params->num_cells_series;
     P_watts *= params->dt_hr;
-    double multiplier = 1.;
+    int multiplier = 1;
     if (P_watts < 0)
-        multiplier = -1.;
+        multiplier = -1;
 
     size_t row = 0;
     while (row < params->voltage_table.size() && DOD > params->voltage_table[row][0]) {
@@ -266,12 +273,12 @@ double voltage_table_t::calculate_current_for_target_w(double P_watts, double q,
     double B = qmax / 100.;
 
     double DOD_new = 0.;
-    double incr = 0;
-    double DOD_best = DOD_best = multiplier == -1. ? 0 : 100;
+    int incr = 0;
+    double DOD_best = DOD_best = (multiplier == -1) ? 0 : 100;
     double P_best = 0;
-    while (incr + row < slopes.size() && incr + row >= 0) {
-        size_t i = row + (size_t) incr;
-        incr += 1 * multiplier;
+    while (((incr + row) < slopes.size()) && ((incr + row) >= 0)) {
+        size_t i = row + incr;
+        incr += multiplier;
 
         double a = B * slopes[i];
         double b = A * slopes[i] + B * intercepts[i];
@@ -288,7 +295,8 @@ double voltage_table_t::calculate_current_for_target_w(double P_watts, double q,
         auto DOD_upper = params->voltage_table[upper][0];
         auto DOD_lower = params->voltage_table[lower][0];
         if (DOD_new <= DOD_upper && DOD_new >= DOD_lower) {
-            double P = (q - (100. - DOD_new) * qmax/100) * (a * DOD_new + b);
+            current = qmax * ((1. - DOD / 100.) - (1. - DOD_new / 100.)) / params->dt_hr;
+            double P = current * (a * DOD_new + b - current / params->num_strings * params->resistance);
             if (std::abs(P) > std::abs(P_best)) {
                 P_best = P;
                 DOD_best = DOD_new;