From 799f9237818594566353b95cce73f2bd5f340471 Mon Sep 17 00:00:00 2001
From: Margaret McCall <majemc@gmail.com>
Date: Mon, 18 May 2020 13:31:19 -0700
Subject: [PATCH 01/14] first edits to new opti

---
 v2gsim/post_simulation/netload_opti_plus.py | 578 ++++++++++++++++++++
 1 file changed, 578 insertions(+)
 create mode 100644 v2gsim/post_simulation/netload_opti_plus.py

diff --git a/v2gsim/post_simulation/netload_opti_plus.py b/v2gsim/post_simulation/netload_opti_plus.py
new file mode 100644
index 0000000..9888530
--- /dev/null
+++ b/v2gsim/post_simulation/netload_opti_plus.py
@@ -0,0 +1,578 @@
+# -*- coding: utf-8 -*-
+from __future__ import division
+from pyomo.opt import SolverFactory
+from pyomo.environ import *
+import time
+import pandas
+import numpy
+import matplotlib.pyplot as plt
+import seaborn as sns
+import v2gsim.model
+import v2gsim.result
+import datetime
+
+
+class CentralOptimization(object):
+    """Creates an object to perform optimization.
+    The object contains some general parameters for the optimization
+    """
+    def __init__(self, project, optimization_timestep, date_from,
+                 date_to, minimum_SOC=0.1, maximum_SOC=0.95):
+        # All the variables are at the project timestep except for the model variables
+        # optimization_timestep is in minutes
+        self.optimization_timestep = optimization_timestep
+        # Set minimum_SOC
+        self.minimum_SOC = minimum_SOC
+        self.maximum_SOC = maximum_SOC
+        # Set date boundaries, should be same as the one used during the simulation
+        self.date_from = date_from
+        self.date_to = date_to
+        self.SOC_index_from = int((date_from - project.date).total_seconds() / project.timestep)
+        self.SOC_index_to = int((date_to - project.date).total_seconds() / project.timestep)
+
+    def solve(self, project, net_load, real_number_of_vehicle, SOC_margin=0.02,
+              SOC_offset=0.0, peak_shaving='peak_shaving', penalization=5, beta=None, plot=False):
+        """Launch the optimization and the post_processing fucntion. Results
+        and assumptions are appended to a data frame.
+
+        Args:
+            project (Project): project
+            net_load (pandas.DataFrame): data frame with date index and a 'net_load' column in [W]
+            real_number_of_vehicle (int): number of vehicle expected on the net load,
+                False if number is the same as in project
+            SOC_margin (float): SOC margin that can be used by the optimization at the end of the day [0, 1]
+            SOC_offset (float): energy offset [0, 1]
+            peak_shaving (bolean): if True ramping constraints are not taking in account within the objective else it is.
+        """
+        # Reset model
+        self.times = []
+        self.vehicles = []
+        self.d = {}
+        self.pmax = {}
+        self.pmin = {}
+        self.emin = {}
+        self.emax = {}
+        self.efinal = {}
+
+        # Set the variables for the optimization
+        new_net_load = self.initialize_net_load(net_load, real_number_of_vehicle, project)
+        self.initialize_model(project, new_net_load, SOC_margin, SOC_offset)
+
+        # Run the optimization
+        timer = time.time()
+        opti_model, result = self.process(self.times, self.vehicles, self.d, self.pmax,
+                                          self.pmin, self.emin, self.emax,
+                                          self.efinal, peak_shaving, penalization)
+        timer2 = time.time()
+        print('The optimization duration was ' + str((timer2 - timer) / 60) + ' minutes')
+        print('')
+
+        # Post process results
+        return self.post_process(project, net_load, opti_model, result, plot)
+
+    def initialize_net_load(self, net_load, real_number_of_vehicle, project):
+        """Make sure that the net load has the right size and scale the net
+        load for the optimization scale.
+
+        Args:
+            net_load (pandas.DataFrame): data frame with date index and a 'net_load' column in [W]
+            net_load_pmax (int): maximum power on the scaled net load
+        """
+        # Make sure we are not touching the initial data
+        new_net_load = net_load.copy()
+
+        # Resample the net load
+        new_net_load = new_net_load.resample(str(self.optimization_timestep) + 'T').first()
+
+        # Check against the actual lenght it should have
+        diff = (len(new_net_load) -
+                int((self.date_to - self.date_from).total_seconds() / (60 * self.optimization_timestep)))
+        if diff > 0:
+            # We should trim the net load with diff elements (normaly 1 because of slicing inclusion)
+            new_net_load.drop(new_net_load.tail(diff).index, inplace=True)
+        elif diff < 0:
+            print('The net load does not contain enough data points')
+
+        if real_number_of_vehicle:
+            # Set scaling factor
+            scaling_factor = len(project.vehicles) / real_number_of_vehicle
+
+            # Scale the temp net load
+            new_net_load['netload'] *= scaling_factor
+
+        return new_net_load
+
+    def check_energy_constraints_feasible(self, vehicle, SOC_init, SOC_final, SOC_offset, verbose=False):
+        """Make sure that SOC final can be reached from SOC init under uncontrolled
+        charging (best case scenario). Print details when conditions are not met.
+
+        Args:
+            vehicle (Vehicle): vehicle
+            SOC_init (float): state of charge at the begining of the optimization [0, 1]
+            SOC_final (float): state of charge at the end of the optimization [0, 1]
+            SOC_offset (float): energy offset [0, 1]
+
+        Return:
+            (Boolean)
+        """
+        def print_status(SOC_final, SOC_init, diff):
+                print('Set battery gain: ' + str((SOC_final - SOC_init) * 100) + '%')
+                print('Simulation battery gain: ' + str(diff * 100))
+
+        # Check if below minimum SOC at any time
+        if (min(vehicle.SOC[self.SOC_index_from:self.SOC_index_to]) - SOC_offset) <= self.minimum_SOC:
+            if verbose:
+                print('Vehicle: ' + str(vehicle.id) + ' has a minimum SOC of ' +
+                      str(min(vehicle.SOC[self.SOC_index_from:self.SOC_index_to]) * 100) + '%')
+            return False
+
+        # Check SOC difference between date_from and date_to ?
+        # Diff represent the minimum loss or the maximum gain
+        diff = vehicle.SOC[self.SOC_index_to] - vehicle.SOC[self.SOC_index_from]
+        # Simulation show a battery gain
+        if diff > 0:
+            # Gain should be greater than the one we set up
+            if SOC_final - SOC_init < diff:
+                # Good to go
+                return True
+            else:
+                # Set final SOC to be under diff
+                print_status(SOC_final, SOC_init, diff)
+                return False
+        # Simulation show battery loss
+        else:
+            # energy balance should be negative
+            if SOC_final - SOC_init > 0:
+                print_status(SOC_final, SOC_init, diff)
+                return False
+            # Loss should be smaller than the one we set (less negative)
+            if diff < SOC_init - SOC_final:
+                return True
+            else:
+                # Set final SOC to include at least the lost
+                print_status(SOC_final, SOC_init, diff)
+                return False
+
+    def initialize_time_index(self, net_load):
+        """Replace date index by time ids
+
+        Args:
+            net_load (pandas.DataFrame): data frame with date index and a 'net_load' column in [W]
+
+        Retrun:
+            net_load as a dictionary with time ids, time ids
+        """
+        temp_index = pandas.DataFrame(range(0, len(net_load)), columns=['index'])
+        # Set temp_index
+        temp_net_load = net_load.copy()
+        temp_net_load = temp_net_load.set_index(temp_index['index'])
+        # Return a dictionary
+        return temp_net_load.to_dict()['netload'], temp_index.index.values.tolist()
+
+    def get_initial_SOC(self, vehicle, SOC_offset, SOC_init=None):
+        """Get the initial SOC with which people start the optimization
+        """
+        if SOC_init is not None:
+            return SOC_init
+        else:
+            return vehicle.SOC[self.SOC_index_from] - SOC_offset
+
+    def get_final_SOC(self, vehicle, SOC_margin, SOC_offset, SOC_end=None):
+        """Get final SOC that vehicle must reached at the end of the optimization
+        """
+        if SOC_end is not None:
+            return SOC_end
+        else:
+            return vehicle.SOC[self.SOC_index_to] - SOC_offset - SOC_margin
+
+    def initialize_model(self, project, net_load, SOC_margin, SOC_offset):
+        """Select the vehicles that were plugged at controlled chargers and create
+        the optimization variables (see inputs of optimization)
+
+        Args:
+            project (Project): project
+
+        Return:
+            times, vehicles, d, pmax, pmin, emin, emax, efinal
+        """
+        # Create a dict with the net load and get time index in a data frame
+        self.d, self.times = self.initialize_time_index(net_load)
+        vehicle_to_optimize = 0
+        unfeasible_vehicle = 0
+
+        for vehicle in project.vehicles:
+            if vehicle.result is not None:
+                # Get SOC init and SOC end
+                SOC_init = self.get_initial_SOC(vehicle, SOC_offset)
+                SOC_final = self.get_final_SOC(vehicle, SOC_margin, SOC_offset)
+
+                # Find out if vehicle itinerary is feasible
+                if not self.check_energy_constraints_feasible(vehicle, SOC_init, SOC_final, SOC_offset):
+                    # Reset vehicle result to None
+                    vehicle.result = None
+                    unfeasible_vehicle += 1
+                    continue
+
+                # Add vehicle id to a list
+                self.vehicles.append(vehicle.id)
+                vehicle_to_optimize += 1
+
+                # Resample vehicle result
+                temp_vehicle_result = vehicle.result.resample(str(self.optimization_timestep) + 'T').first()
+
+                # Set time_vehicle_index
+                temp_vehicle_result = temp_vehicle_result.set_index(pandas.DataFrame(
+                    index=[(time, vehicle.id) for time in self.times]).index)
+
+                # Push pmax and pmin with vehicle and time key
+                self.pmin.update(temp_vehicle_result.to_dict()['p_min'])
+                self.pmax.update(temp_vehicle_result.to_dict()['p_max'])
+
+                # Push emin and emax with vehicle and time key
+                # Units! if project.timestep in seconds, self.timestep in minutes and battery in Wh
+                # Units! Wproject.timestep --> Wself.timestep * (project.timestep / (60 * self.timestep))
+                # Units! Wtimestep --> Wh * (60 / self.timestep)
+                temp_vehicle_result['emin'] = (temp_vehicle_result.energy * (project.timestep / (60 * self.optimization_timestep)) +
+                                               (self.minimum_SOC - SOC_init) * vehicle.car_model.battery_capacity *
+                                               (60 / self.optimization_timestep))
+                self.emin.update(temp_vehicle_result.to_dict()['emin'])
+                temp_vehicle_result['emax'] = (temp_vehicle_result.energy * (project.timestep / (60 * self.optimization_timestep)) + 10000 +
+                                               (self.maximum_SOC - SOC_init) * vehicle.car_model.battery_capacity *
+                                               (60 / self.optimization_timestep))
+                self.emax.update(temp_vehicle_result.to_dict()['emax'])
+
+                # Push efinal with vehicle key
+                self.efinal.update({vehicle.id: (temp_vehicle_result.tail(1).energy.values[0] * (project.timestep / (60 * self.optimization_timestep)) +
+                                                 (SOC_final - SOC_init) * vehicle.car_model.battery_capacity *
+                                                 (60 / self.optimization_timestep))})
+
+        print('There is ' + str(vehicle_to_optimize) + ' vehicle participating in the optimization (' +
+              str(vehicle_to_optimize * 100 / len(project.vehicles)) + '%)')
+        print('There is ' + str(unfeasible_vehicle) + ' unfeasible vehicle.')
+        print('')
+
+    def process(self, times, vehicles, d, price, pmax, pmin, emin, emax,
+                efinal, peak_shaving, penalization, solver="gurobi"):
+        """The process function creates the pyomo model and solve it.
+        Minimize sum( net_load(t) + sum(power_demand(t, v)))**2
+        subject to:
+        pmin(t, v) <= power_demand(t, v) <= pmax(t, v)
+        emin(t, v) <= sum(power_demand(t, v)) <= emax(t, v)
+        sum(power_demand(t, v)) >= efinal(v)
+        rampmin(t) <= net_load_ramp(t) + power_demand_ramp(t, v) <= rampmax(t)
+
+        Args:
+            times (list): timestep list
+            vehicles (list): unique list of vehicle ids
+            d (dict): time - net load at t
+            price (dict): time - power price at t ###
+            pmax (dict): (time, id) - power maximum at t for v
+            pmin (dict): (time, id) - power minimum at t for v
+            emin (dict): (time, id) - energy minimum at t for v
+            emax (dict): (time, id) - energy maximum at t for v
+            efinal (dict): id - final SOC
+            solver (string): name of the solver to use (default is gurobi)
+
+        Return:
+            model (ConcreteModel), result
+        """
+
+        # Select gurobi solver
+        with SolverFactory(solver) as opt:
+            # Solver option see Gurobi website
+            # opt.options['Method'] = 1
+
+            # Creation of a Concrete Model
+            model = ConcreteModel()
+
+            # ###### Set
+            model.t = Set(initialize=times, doc='Time', ordered=True)
+            last_t = model.t.last()
+            model.v = Set(initialize=vehicles, doc='Vehicles')
+
+            # ###### Parameters
+            # Net load
+            model.d = Param(model.t, initialize=d, doc='Net load')
+
+            # Price
+            model.price = Param(model.t, initialize=price, doc='Prices')
+
+            # Power
+            model.p_max = Param(model.t, model.v, initialize=pmax, doc='P max')
+            model.p_min = Param(model.t, model.v, initialize=pmin, doc='P min')
+
+            # Energy
+            model.e_min = Param(model.t, model.v, initialize=emin, doc='E min')
+            model.e_max = Param(model.t, model.v, initialize=emax, doc='E max')
+
+            model.e_final = Param(model.v, initialize=efinal, doc='final energy balance')
+            # model.beta = Param(initialize=beta, doc='beta')
+
+            # ###### Variable
+            model.u = Var(model.t, model.v, domain=Integers, doc='Power used')
+
+            # ###### Rules
+            def maximum_power_rule(model, t, v):
+                return model.u[t, v] <= model.p_max[t, v]
+            model.power_max_rule = Constraint(model.t, model.v, rule=maximum_power_rule, doc='P max rule')
+
+            def minimum_power_rule(model, t, v):
+                return model.u[t, v] >= model.p_min[t, v]
+            model.power_min_rule = Constraint(model.t, model.v, rule=minimum_power_rule, doc='P min rule')
+
+            def minimum_energy_rule(model, t, v):
+                return sum(model.u[i, v] for i in range(0, t + 1)) >= model.e_min[t, v]
+            model.minimum_energy_rule = Constraint(model.t, model.v, rule=minimum_energy_rule, doc='E min rule')
+
+            def maximum_energy_rule(model, t, v):
+                return sum(model.u[i, v] for i in range(0, t + 1)) <= model.e_max[t, v]
+            model.maximum_energy_rule = Constraint(model.t, model.v, rule=maximum_energy_rule, doc='E max rule')
+
+            def final_energy_balance(model, v):
+                return sum(model.u[i, v] for i in model.t) >= model.e_final[v]
+            model.final_energy_rule = Constraint(model.v, rule=final_energy_balance, doc='E final rule')
+
+            # Set the objective to be either peak shaving or ramp mitigation
+            if peak_shaving == 'peak_shaving':
+                def objective_rule(model):
+                    return sum([(model.d[t] + sum([model.u[t, v] for v in model.v]))**2 for t in model.t])
+                model.objective = Objective(rule=objective_rule, sense=minimize, doc='Define objective function')
+
+            elif peak_shaving == 'penalized_peak_shaving':
+                def objective_rule(model):
+                    return (sum([(model.d[t] + sum([model.u[t, v] for v in model.v]))**2 for t in model.t]) +
+                            penalization * sum([sum([model.u[t, v]**2 for v in model.v]) for t in model.t]))
+                model.objective = Objective(rule=objective_rule, sense=minimize, doc='Define objective function')
+
+            elif peak_shaving == 'ramp_mitigation':
+                def objective_rule(model):
+                    return sum([(model.d[t + 1] - model.d[t] + sum([model.u[t + 1, v] - model.u[t, v] for v in model.v]))**2 for t in model.t if t != last_t])
+                model.objective = Objective(rule=objective_rule, sense=minimize, doc='Define objective function')
+
+            ### MM additions
+            # Hybrid
+            # Definition: optimizes for both ramp mitigation and peak shaving
+            # No new constraints needed, but need to define lambda
+            elif peak_shaving == 'hybrid': 
+                def objective_rule(model):
+                    lamb = .5 #weights peak shaving and ramp mitigation goals equally
+                    return lamb * sum([(model.d[t] + sum([model.u[t, v] for v in model.v]))**2 for t in model.t]) + \
+                        (1-lamb) * sum([(model.d[t + 1] - model.d[t] + sum([model.u[t + 1, v] - model.u[t, v] for v in model.v]))**2 for t in model.t if t != last_t])
+                model.objective = Objective(rule=objective_rule, sense=minimize, doc='Define objective function')
+
+            # Emergency
+            # Definition: Emergency V2G; mostly V1G but each vehicle can get called for V2G a certain number of times/year
+            # Clearly a new constraint, but do I need to create a new vehicle attribute to track this...?
+            # Can be used either for peak shaving or ramp mitigation; for now just implementing peak shaving
+            elif peak_shaving == 'emergency':
+                def objective_rule(model):
+                    return sum([(model.d[t] + sum([model.u[t, v] for v in model.v]))**2 for t in model.t])
+                model.objective = Objective(rule=objective_rule, sense=minimize, doc='Define objective function')
+
+                def max_call_rule(model, v):
+                    max_calls = 2 #max V2G calls in a week...need to calculate based on timesteps fed in
+                    return sum(model.u[i, v] for i in model.t < 0) <= max_calls #hoping the first gives # of times power is discharging
+                model.max_call_rule = Constraint(model.v, rule=max_call_rule, doc='Max call rule')
+
+            # Cost-based
+            # Definition: minimizes charging cost based on given cost stack
+            # Should probably only have visibility 24/48 hours into the future, ideally
+            # V1G vs V2G will make a big difference here--incentivizing selling on-peak
+            # Need to feed in p (prices at t)
+            elif peak_shaving == 'cost':
+                def objective_rule(model):
+                    return sum((model.u[t, v] * model.price[t] for v in model.v) for t in model.t)
+                model.objective = Objective(rule=objective_rule, sense=minimize, doc='Define objective function')
+
+            
+            
+            results = opt.solve(model)
+            # results.write()
+
+        return model, results
+
+    def plot_result(self, model):
+        """Create a plot showing the power constraints, the energy constraints and the ramp
+        constraints as well as the final net load.
+        """
+        # Set the graph style
+        sns.set_style("whitegrid")
+        sns.despine()
+
+        result = pandas.DataFrame()
+
+        # Get the result
+        df = pandas.DataFrame(index=['power'], data=model.u.get_values()).transpose().groupby(level=0).sum()
+        # Ramp of the result
+        mylist = [0]
+        mylist.extend(list(numpy.diff(df['power'].values)))
+        df['ramppower'] = mylist
+        result = pandas.concat([result, df], axis=1)
+        # cum sum of the result
+        df = pandas.DataFrame(
+            pandas.DataFrame(
+                index=['anything'],
+                data=model.u.get_values()).transpose().unstack().cumsum().sum(axis=1), columns=['powercum']) * (5 / 60)
+        result = pandas.concat([result, df], axis=1)
+
+        # Get pmax and pmin units of power [W]
+        df = pandas.DataFrame(index=['pmax'], data=model.p_max.extract_values()).transpose().groupby(level=0).sum()
+        result = pandas.concat([result, df], axis=1)
+        df = pandas.DataFrame(index=['pmin'], data=model.p_min.extract_values()).transpose().groupby(level=0).sum()
+        result = pandas.concat([result, df], axis=1)
+
+        # Get emin and emax units of energy [Wh]
+        df = pandas.DataFrame(index=['emax'], data=model.e_max.extract_values()).transpose().groupby(level=0).sum() * (5 / 60)
+        result = pandas.concat([result, df], axis=1)
+        df = pandas.DataFrame(index=['emin'], data=model.e_min.extract_values()).transpose().groupby(level=0).sum() * (5 / 60)
+        result = pandas.concat([result, df], axis=1)
+
+        # Get the minimum final energy quantity
+        e_final = pandas.DataFrame(index=['efinal'], data=model.e_final.extract_values()).transpose().sum() * (5 / 60)
+
+        # Get the actual ramp
+        mylist = [0]
+        df = pandas.DataFrame(index=['net_load'], data=model.d.extract_values()).transpose()
+        mylist.extend(list(numpy.diff(df['net_load'].values)))
+        df['rampnet_load'] = mylist
+        result = pandas.concat([result, df], axis=1)
+
+        # Plot power constraints
+        plt.subplot(411)
+        plt.plot(result.index.values, result.pmax.values, label='pmax')
+        plt.plot(result.index.values, result.power.values, label='power')
+        plt.plot(result.index.values, result.pmin.values, label='pmin')
+        plt.legend(loc=0)
+
+        # Plot energy constraints
+        plt.subplot(412)
+        plt.plot(result.index.values, result.emax.values, label='emax')
+        plt.plot(result.index.values, result.powercum.values, label='powercum')
+        plt.plot(result.index.values, result.emin.values, label='emin')
+        plt.plot(result.index[-1], e_final, '*', markersize=15, label='efinal')
+        plt.legend(loc=0)
+
+        # Plot the result ramp
+        plt.subplot(413)
+        plt.plot(result.index.values, result.ramppower.values + result.rampnet_load.values, label='result ramp')
+        plt.plot(result.index.values, result.rampnet_load.values, label='net_load ramp')
+        plt.legend(loc=0)
+
+        # Plot the power demand results
+        plt.subplot(414)
+        plt.plot(result.index.values, result.net_load.values, label='net_load')
+        plt.plot(result.index.values, result.net_load.values + result.power.values, label='net_load + vehicle')
+        plt.legend(loc=0)
+        plt.show()
+
+    def post_process(self, project, netload, model, result, plot):
+        """Recompute SOC profiles and compute new total power demand
+
+        Args:
+            project (Project): project
+
+        Note: Should check that 'vehicle before' and 'after' contain the same number of vehicles
+        """
+        if plot:
+            self.plot_result(model)
+
+        temp = pandas.DataFrame()
+        first = True
+        for vehicle in project.vehicles:
+            if vehicle.result is not None:
+                if first:
+                    temp['vehicle_before'] = vehicle.result['power_demand']
+                    first = False
+                else:
+                    temp['vehicle_before'] += vehicle.result['power_demand']
+
+        temp2 = pandas.DataFrame(index=['vehicle_after'], data=model.u.get_values()).transpose().groupby(level=0).sum()
+        i = pandas.date_range(start=self.date_from, end=self.date_to,
+                              freq=str(self.optimization_timestep) + 'T', closed='left')
+        temp2 = temp2.set_index(i)
+        temp2 = temp2.resample(str(project.timestep) + 'S')
+        temp2 = temp2.fillna(method='ffill').fillna(method='bfill')
+
+        final_result = pandas.DataFrame()
+        final_result = pandas.concat([temp['vehicle_before'], temp2['vehicle_after']], axis=1)
+        final_result = final_result.fillna(method='ffill').fillna(method='bfill')
+
+        return final_result
+
+
+def save_vehicle_state_for_optimization(vehicle, timestep, date_from,
+                                        date_to, activity=None, power_demand=None,
+                                        SOC=None, detail=None, nb_interval=None, init=False,
+                                        run=False, post=False):
+    """Save results for individual vehicles. Power demand is positive when charging
+    negative when driving. Energy consumption is positive when driving and negative
+    when charging. Charging station that offer after simulation processing should
+    have activity.charging_station.post_simulation True.
+    """
+    if run:
+        if vehicle.result is not None:
+            activity_index1, activity_index2, location_index1, location_index2, save = v2gsim.result._map_index(
+                activity.start, activity.end, date_from, date_to, len(power_demand),
+                len(vehicle.result['power_demand']), timestep)
+            # Time frame are matching
+            if save:
+                # If driving pmin and pmax are equal to 0 since we are not plugged
+                if isinstance(activity, v2gsim.model.Driving):
+                    vehicle.result['p_max'][location_index1:location_index2] -= (
+                        [0.0] * (activity_index2 - activity_index1))
+                    vehicle.result['p_min'][location_index1:location_index2] -= (
+                        [0.0] * (activity_index2 - activity_index1))
+                    # Energy consumed is directly the power demand (sum later)
+                    vehicle.result['energy'][location_index1:location_index2] += (
+                        power_demand[activity_index1:activity_index2])
+                    # Power demand on the grid is 0 since we are driving
+                    vehicle.result['power_demand'][location_index1:location_index2] -= (
+                        [0.0] * (activity_index2 - activity_index1))
+
+                # If parked pmin and pmax are not necessary the same
+                if isinstance(activity, v2gsim.model.Parked):
+                    # Save the positive power demand of this specific vehicle
+                    vehicle.result['power_demand'][location_index1:location_index2] += (
+                        power_demand[activity_index1:activity_index2])
+                    if activity.charging_station.post_simulation:
+                        # Find if vehicle or infra is limiting
+                        pmax = min(activity.charging_station.maximum_power,
+                                   vehicle.car_model.maximum_power)
+                        pmin = max(activity.charging_station.minimum_power,
+                                   vehicle.car_model.minimum_power)
+                        vehicle.result['p_max'][location_index1:location_index2] += (
+                            [pmax] * (activity_index2 - activity_index1))
+                        vehicle.result['p_min'][location_index1:location_index2] += (
+                            [pmin] * (activity_index2 - activity_index1))
+                        # Energy consumed is 0 the optimization will decide
+                        vehicle.result['energy'][location_index1:location_index2] -= (
+                            [0.0] * (activity_index2 - activity_index1))
+                    else:
+                        vehicle.result['p_max'][location_index1:location_index2] += (
+                            power_demand[activity_index1:activity_index2])
+                        vehicle.result['p_min'][location_index1:location_index2] += (
+                            power_demand[activity_index1:activity_index2])
+                        # Energy is 0.0 because it's already accounted in power_demand
+                        vehicle.result['energy'][location_index1:location_index2] -= (
+                            [0.0] * (activity_index2 - activity_index1))
+
+    elif init:
+        vehicle.SOC = [vehicle.SOC[0]]
+        vehicle.result = None
+        for activity in vehicle.activities:
+            if isinstance(activity, v2gsim.model.Parked):
+                if activity.charging_station.post_simulation:
+                    # Initiate a dictionary of numpy array to hold result (faster than DataFrame)
+                    vehicle.result = {'power_demand': numpy.array([0.0] * int((date_to - date_from).total_seconds() / timestep)),
+                                      'p_max': numpy.array([0.0] * int((date_to - date_from).total_seconds() / timestep)),
+                                      'p_min': numpy.array([0.0] * int((date_to - date_from).total_seconds() / timestep)),
+                                      'energy': numpy.array([0.0] * int((date_to - date_from).total_seconds() / timestep))}
+                    # Leave the init function
+                    return
+    elif post:
+        if vehicle.result is not None:
+            # Convert location result back into pandas DataFrame (faster that way)
+            i = pandas.date_range(start=date_from, end=date_to,
+                                  freq=str(timestep) + 's', closed='left')
+            vehicle.result = pandas.DataFrame(index=i, data=vehicle.result)
+            vehicle.result['energy'] = vehicle.result['energy'].cumsum()

From 1d95f9fb787dd8266b78a5e1aaed01eaba4bedc3 Mon Sep 17 00:00:00 2001
From: Margaret McCall <majemc@gmail.com>
Date: Thu, 21 May 2020 10:49:50 -0700
Subject: [PATCH 02/14] first opti updates

---
 v2gsim/post_simulation/netload_opti_plus.py | 20 ++++++++++++--------
 1 file changed, 12 insertions(+), 8 deletions(-)

diff --git a/v2gsim/post_simulation/netload_opti_plus.py b/v2gsim/post_simulation/netload_opti_plus.py
index 9888530..562c97e 100644
--- a/v2gsim/post_simulation/netload_opti_plus.py
+++ b/v2gsim/post_simulation/netload_opti_plus.py
@@ -69,6 +69,7 @@ def solve(self, project, net_load, real_number_of_vehicle, SOC_margin=0.02,
 
         # Post process results
         return self.post_process(project, net_load, opti_model, result, plot)
+        #return opti_model, result
 
     def initialize_net_load(self, net_load, real_number_of_vehicle, project):
         """Make sure that the net load has the right size and scale the net
@@ -251,7 +252,7 @@ def initialize_model(self, project, net_load, SOC_margin, SOC_offset):
         print('There is ' + str(unfeasible_vehicle) + ' unfeasible vehicle.')
         print('')
 
-    def process(self, times, vehicles, d, price, pmax, pmin, emin, emax,
+    def process(self, times, vehicles, d, pmax, pmin, emin, emax,
                 efinal, peak_shaving, penalization, solver="gurobi"):
         """The process function creates the pyomo model and solve it.
         Minimize sum( net_load(t) + sum(power_demand(t, v)))**2
@@ -290,12 +291,16 @@ def process(self, times, vehicles, d, price, pmax, pmin, emin, emax,
             last_t = model.t.last()
             model.v = Set(initialize=vehicles, doc='Vehicles')
 
+
             # ###### Parameters
             # Net load
             model.d = Param(model.t, initialize=d, doc='Net load')
 
             # Price
-            model.price = Param(model.t, initialize=price, doc='Prices')
+            #model.price = Param(model.t, initialize=price, doc='Prices')
+
+            # Lambda for hybrid model
+            model.lamb = Param(initialize = .5)
 
             # Power
             model.p_max = Param(model.t, model.v, initialize=pmax, doc='P max')
@@ -308,9 +313,11 @@ def process(self, times, vehicles, d, price, pmax, pmin, emin, emax,
             model.e_final = Param(model.v, initialize=efinal, doc='final energy balance')
             # model.beta = Param(initialize=beta, doc='beta')
 
+
             # ###### Variable
             model.u = Var(model.t, model.v, domain=Integers, doc='Power used')
 
+
             # ###### Rules
             def maximum_power_rule(model, t, v):
                 return model.u[t, v] <= model.p_max[t, v]
@@ -355,9 +362,8 @@ def objective_rule(model):
             # No new constraints needed, but need to define lambda
             elif peak_shaving == 'hybrid': 
                 def objective_rule(model):
-                    lamb = .5 #weights peak shaving and ramp mitigation goals equally
-                    return lamb * sum([(model.d[t] + sum([model.u[t, v] for v in model.v]))**2 for t in model.t]) + \
-                        (1-lamb) * sum([(model.d[t + 1] - model.d[t] + sum([model.u[t + 1, v] - model.u[t, v] for v in model.v]))**2 for t in model.t if t != last_t])
+                    return model.lamb * sum([(model.d[t] + sum([model.u[t, v] for v in model.v]))**2 for t in model.t]) + \
+                        (1-model.lamb) * sum([(model.d[t + 1] - model.d[t] + sum([model.u[t + 1, v] - model.u[t, v] for v in model.v]))**2 for t in model.t if t != last_t])
                 model.objective = Objective(rule=objective_rule, sense=minimize, doc='Define objective function')
 
             # Emergency
@@ -370,7 +376,7 @@ def objective_rule(model):
                 model.objective = Objective(rule=objective_rule, sense=minimize, doc='Define objective function')
 
                 def max_call_rule(model, v):
-                    max_calls = 2 #max V2G calls in a week...need to calculate based on timesteps fed in
+                    max_calls = 2 #max V2G calls per vehicle in a week...need to calculate based on timesteps fed in
                     return sum(model.u[i, v] for i in model.t < 0) <= max_calls #hoping the first gives # of times power is discharging
                 model.max_call_rule = Constraint(model.v, rule=max_call_rule, doc='Max call rule')
 
@@ -383,8 +389,6 @@ def max_call_rule(model, v):
                 def objective_rule(model):
                     return sum((model.u[t, v] * model.price[t] for v in model.v) for t in model.t)
                 model.objective = Objective(rule=objective_rule, sense=minimize, doc='Define objective function')
-
-            
             
             results = opt.solve(model)
             # results.write()

From 2a079a95c5cbf0f14f457f77647ce123b833a3b0 Mon Sep 17 00:00:00 2001
From: Margaret McCall <majemc@gmail.com>
Date: Thu, 21 May 2020 12:12:48 -0700
Subject: [PATCH 03/14] hybrid optimization working

---
 v2gsim/post_simulation/netload_opti_plus.py | 11 +++++++----
 1 file changed, 7 insertions(+), 4 deletions(-)

diff --git a/v2gsim/post_simulation/netload_opti_plus.py b/v2gsim/post_simulation/netload_opti_plus.py
index 562c97e..973dc9c 100644
--- a/v2gsim/post_simulation/netload_opti_plus.py
+++ b/v2gsim/post_simulation/netload_opti_plus.py
@@ -31,7 +31,7 @@ def __init__(self, project, optimization_timestep, date_from,
         self.SOC_index_to = int((date_to - project.date).total_seconds() / project.timestep)
 
     def solve(self, project, net_load, real_number_of_vehicle, SOC_margin=0.02,
-              SOC_offset=0.0, peak_shaving='peak_shaving', penalization=5, beta=None, plot=False):
+              SOC_offset=0.0, peak_shaving='peak_shaving', penalization=5, beta=None, plot=False, peak_scalar = .00137):
         """Launch the optimization and the post_processing fucntion. Results
         and assumptions are appended to a data frame.
 
@@ -62,7 +62,7 @@ def solve(self, project, net_load, real_number_of_vehicle, SOC_margin=0.02,
         timer = time.time()
         opti_model, result = self.process(self.times, self.vehicles, self.d, self.pmax,
                                           self.pmin, self.emin, self.emax,
-                                          self.efinal, peak_shaving, penalization)
+                                          self.efinal, peak_shaving, penalization, peak_scalar)
         timer2 = time.time()
         print('The optimization duration was ' + str((timer2 - timer) / 60) + ' minutes')
         print('')
@@ -253,7 +253,7 @@ def initialize_model(self, project, net_load, SOC_margin, SOC_offset):
         print('')
 
     def process(self, times, vehicles, d, pmax, pmin, emin, emax,
-                efinal, peak_shaving, penalization, solver="gurobi"):
+                efinal, peak_shaving, penalization, peak_scalar, solver="gurobi"):
         """The process function creates the pyomo model and solve it.
         Minimize sum( net_load(t) + sum(power_demand(t, v)))**2
         subject to:
@@ -302,6 +302,9 @@ def process(self, times, vehicles, d, pmax, pmin, emin, emax,
             # Lambda for hybrid model
             model.lamb = Param(initialize = .5)
 
+            # Scaling factor for peak-shaving sub-objective for hybrid model
+            model.peak_scalar = Param(initialize = peak_scalar)
+
             # Power
             model.p_max = Param(model.t, model.v, initialize=pmax, doc='P max')
             model.p_min = Param(model.t, model.v, initialize=pmin, doc='P min')
@@ -362,7 +365,7 @@ def objective_rule(model):
             # No new constraints needed, but need to define lambda
             elif peak_shaving == 'hybrid': 
                 def objective_rule(model):
-                    return model.lamb * sum([(model.d[t] + sum([model.u[t, v] for v in model.v]))**2 for t in model.t]) + \
+                    return model.lamb * model.peak_scalar * sum([(model.d[t] + sum([model.u[t, v] for v in model.v]))**2 for t in model.t]) + \
                         (1-model.lamb) * sum([(model.d[t + 1] - model.d[t] + sum([model.u[t + 1, v] - model.u[t, v] for v in model.v]))**2 for t in model.t if t != last_t])
                 model.objective = Objective(rule=objective_rule, sense=minimize, doc='Define objective function')
 

From b6d015a9fd9d57b144c13ad7373f471b5dca01ff Mon Sep 17 00:00:00 2001
From: Margaret McCall <majemc@gmail.com>
Date: Thu, 21 May 2020 12:58:03 -0700
Subject: [PATCH 04/14] beginning emergency optimization

---
 v2gsim/post_simulation/netload_opti_plus.py | 33 ++++++++++++++-------
 1 file changed, 23 insertions(+), 10 deletions(-)

diff --git a/v2gsim/post_simulation/netload_opti_plus.py b/v2gsim/post_simulation/netload_opti_plus.py
index 973dc9c..00df9f1 100644
--- a/v2gsim/post_simulation/netload_opti_plus.py
+++ b/v2gsim/post_simulation/netload_opti_plus.py
@@ -299,12 +299,6 @@ def process(self, times, vehicles, d, pmax, pmin, emin, emax,
             # Price
             #model.price = Param(model.t, initialize=price, doc='Prices')
 
-            # Lambda for hybrid model
-            model.lamb = Param(initialize = .5)
-
-            # Scaling factor for peak-shaving sub-objective for hybrid model
-            model.peak_scalar = Param(initialize = peak_scalar)
-
             # Power
             model.p_max = Param(model.t, model.v, initialize=pmax, doc='P max')
             model.p_min = Param(model.t, model.v, initialize=pmin, doc='P min')
@@ -316,9 +310,21 @@ def process(self, times, vehicles, d, pmax, pmin, emin, emax,
             model.e_final = Param(model.v, initialize=efinal, doc='final energy balance')
             # model.beta = Param(initialize=beta, doc='beta')
 
+            #HYBRID
+            # Lambda for hybrid model
+            model.lamb = Param(initialize = .5)
+            # Scaling factor for peak-shaving sub-objective for hybrid model
+            model.peak_scalar = Param(initialize = peak_scalar)
+
+            #EMERGENCY
+            # Maximum number of times a vehicle can be called for V2G in emergency scenario
+            model.max_calls = Param(model.v, initialize = 2) #trying to initialize max_calls for each vehicle being 2
+            model.zero = Param(model.v, initialize = 0) #trying to get vector to assess if u is negative
+
 
             # ###### Variable
             model.u = Var(model.t, model.v, domain=Integers, doc='Power used')
+            model.num_export = Var(model.t, model.v, domain=Binary, doc='Power exported at timestep?')
 
 
             # ###### Rules
@@ -362,7 +368,6 @@ def objective_rule(model):
             ### MM additions
             # Hybrid
             # Definition: optimizes for both ramp mitigation and peak shaving
-            # No new constraints needed, but need to define lambda
             elif peak_shaving == 'hybrid': 
                 def objective_rule(model):
                     return model.lamb * model.peak_scalar * sum([(model.d[t] + sum([model.u[t, v] for v in model.v]))**2 for t in model.t]) + \
@@ -371,16 +376,24 @@ def objective_rule(model):
 
             # Emergency
             # Definition: Emergency V2G; mostly V1G but each vehicle can get called for V2G a certain number of times/year
-            # Clearly a new constraint, but do I need to create a new vehicle attribute to track this...?
             # Can be used either for peak shaving or ramp mitigation; for now just implementing peak shaving
             elif peak_shaving == 'emergency':
                 def objective_rule(model):
                     return sum([(model.d[t] + sum([model.u[t, v] for v in model.v]))**2 for t in model.t])
                 model.objective = Objective(rule=objective_rule, sense=minimize, doc='Define objective function')
 
+                """def count_exports_rule(model, t, v):
+                    for t in model.t:
+                        for v in model.v:
+                            if value(model.u[t, v]) < 0:
+                                model.num_export[t, v] == 1  
+                            elif value(model.u[u, v]) >= 0:
+                                model.num_export[t, v] == 0
+                    return model.num_export
+                model.count_exports_rule = Constraint(model.t, model.v, rule=count_exports_rule, doc='Count exports rule')"""
+
                 def max_call_rule(model, v):
-                    max_calls = 2 #max V2G calls per vehicle in a week...need to calculate based on timesteps fed in
-                    return sum(model.u[i, v] for i in model.t < 0) <= max_calls #hoping the first gives # of times power is discharging
+                    return sum(value(model.u[t, v]) < 0 for t in model.t) <= model.max_calls[v]
                 model.max_call_rule = Constraint(model.v, rule=max_call_rule, doc='Max call rule')
 
             # Cost-based

From 751792954f927ffee4ce64f766cbf80b692fab9c Mon Sep 17 00:00:00 2001
From: Margaret McCall <majemc@gmail.com>
Date: Tue, 26 May 2020 12:48:05 -0700
Subject: [PATCH 05/14] added cost optimization, modified hybrid

---
 v2gsim/post_simulation/netload_opti_plus.py | 71 +++++++++++++++------
 1 file changed, 53 insertions(+), 18 deletions(-)

diff --git a/v2gsim/post_simulation/netload_opti_plus.py b/v2gsim/post_simulation/netload_opti_plus.py
index 00df9f1..165ccb4 100644
--- a/v2gsim/post_simulation/netload_opti_plus.py
+++ b/v2gsim/post_simulation/netload_opti_plus.py
@@ -31,7 +31,7 @@ def __init__(self, project, optimization_timestep, date_from,
         self.SOC_index_to = int((date_to - project.date).total_seconds() / project.timestep)
 
     def solve(self, project, net_load, real_number_of_vehicle, SOC_margin=0.02,
-              SOC_offset=0.0, peak_shaving='peak_shaving', penalization=5, beta=None, plot=False, peak_scalar = .00137):
+              SOC_offset=0.0, peak_shaving='peak_shaving', penalization=5, beta=None, plot=False, peak_scalar = .01, peak_subtractor = 50000, price=pandas.DataFrame()): #price can't be array, is dict; can't have non-default arg follow default
         """Launch the optimization and the post_processing fucntion. Results
         and assumptions are appended to a data frame.
 
@@ -42,7 +42,7 @@ def solve(self, project, net_load, real_number_of_vehicle, SOC_margin=0.02,
                 False if number is the same as in project
             SOC_margin (float): SOC margin that can be used by the optimization at the end of the day [0, 1]
             SOC_offset (float): energy offset [0, 1]
-            peak_shaving (bolean): if True ramping constraints are not taking in account within the objective else it is.
+            peak_shaving (boolean): if True ramping constraints are not taking in account within the objective else it is.
         """
         # Reset model
         self.times = []
@@ -57,12 +57,15 @@ def solve(self, project, net_load, real_number_of_vehicle, SOC_margin=0.02,
         # Set the variables for the optimization
         new_net_load = self.initialize_net_load(net_load, real_number_of_vehicle, project)
         self.initialize_model(project, new_net_load, SOC_margin, SOC_offset)
+        
+        if peak_shaving=='cost':
+            price = self.initialize_price(price, net_load)
 
         # Run the optimization
         timer = time.time()
         opti_model, result = self.process(self.times, self.vehicles, self.d, self.pmax,
                                           self.pmin, self.emin, self.emax,
-                                          self.efinal, peak_shaving, penalization, peak_scalar)
+                                          self.efinal, peak_shaving, penalization, peak_scalar, peak_subtractor, price)
         timer2 = time.time()
         print('The optimization duration was ' + str((timer2 - timer) / 60) + ' minutes')
         print('')
@@ -103,6 +106,30 @@ def initialize_net_load(self, net_load, real_number_of_vehicle, project):
 
         return new_net_load
 
+    def initialize_price(self, price, net_load): ##>>>>
+        """Resample power price data to be on same frequency as net load data, and returns dictionary with timesteps instead of timestamps.
+        Args:
+            price (pandas.DataFrame): data frame with date index and a 'price' column (in $/MWh)
+            net_load (pandas.DataFrame): data frame with date index and a 'net_load' column in [W]
+        
+        Returns dictionary with 
+            """
+        # Check that price and net_load dataframes have same length
+        assert len(price) == len(net_load)
+        # Make sure we are not touching the initial data
+        new_price = price.copy()
+        # Resample the net load
+        new_price = new_price.resample(str(self.optimization_timestep) + 'T').first()
+
+
+        # Convert to timestep index instead of timestamps
+        temp_index = pandas.DataFrame(range(0, len(new_price)), columns=['index'])
+        # Set temp_index
+        temp_new_price = new_price.copy()
+        temp_new_price= temp_new_price.set_index(temp_index['index'])
+        # Return a dictionary
+        return temp_new_price.to_dict()['price']
+
     def check_energy_constraints_feasible(self, vehicle, SOC_init, SOC_final, SOC_offset, verbose=False):
         """Make sure that SOC final can be reached from SOC init under uncontrolled
         charging (best case scenario). Print details when conditions are not met.
@@ -160,7 +187,7 @@ def initialize_time_index(self, net_load):
         Args:
             net_load (pandas.DataFrame): data frame with date index and a 'net_load' column in [W]
 
-        Retrun:
+        Return:
             net_load as a dictionary with time ids, time ids
         """
         temp_index = pandas.DataFrame(range(0, len(net_load)), columns=['index'])
@@ -253,7 +280,8 @@ def initialize_model(self, project, net_load, SOC_margin, SOC_offset):
         print('')
 
     def process(self, times, vehicles, d, pmax, pmin, emin, emax,
-                efinal, peak_shaving, penalization, peak_scalar, solver="gurobi"):
+                efinal, peak_shaving, penalization, peak_scalar, peak_subtractor, price, solver="gurobi"):
+
         """The process function creates the pyomo model and solve it.
         Minimize sum( net_load(t) + sum(power_demand(t, v)))**2
         subject to:
@@ -296,9 +324,6 @@ def process(self, times, vehicles, d, pmax, pmin, emin, emax,
             # Net load
             model.d = Param(model.t, initialize=d, doc='Net load')
 
-            # Price
-            #model.price = Param(model.t, initialize=price, doc='Prices')
-
             # Power
             model.p_max = Param(model.t, model.v, initialize=pmax, doc='P max')
             model.p_min = Param(model.t, model.v, initialize=pmin, doc='P min')
@@ -311,15 +336,25 @@ def process(self, times, vehicles, d, pmax, pmin, emin, emax,
             # model.beta = Param(initialize=beta, doc='beta')
 
             #HYBRID
-            # Lambda for hybrid model
-            model.lamb = Param(initialize = .5)
-            # Scaling factor for peak-shaving sub-objective for hybrid model
-            model.peak_scalar = Param(initialize = peak_scalar)
+            if peak_shaving == 'hybrid':
+                # Lambda for hybrid model
+                model.lamb = Param(initialize = .5)
+                # Scaling factor for peak-shaving sub-objective for hybrid model
+                model.peak_scalar = Param(initialize = peak_scalar)
+                # Normalization factor for peak-shaving sub-objective
+                model.peak_subtractor = Param(initialize = peak_subtractor)
 
             #EMERGENCY
-            # Maximum number of times a vehicle can be called for V2G in emergency scenario
-            model.max_calls = Param(model.v, initialize = 2) #trying to initialize max_calls for each vehicle being 2
-            model.zero = Param(model.v, initialize = 0) #trying to get vector to assess if u is negative
+            if peak_shaving == 'emergency':
+                # Maximum number of times a vehicle can be called for V2G in emergency scenario
+                model.max_calls = Param(model.v, initialize = 2) #trying to initialize max_calls for each vehicle being 2
+                model.zero = Param(model.v, initialize = 0) #trying to get vector to assess if u is negative
+
+            #COST
+            if peak_shaving == 'cost':
+                # Price
+                model.price = Param(model.t, initialize=price, doc='Prices') #energy price at each timestep
+
 
 
             # ###### Variable
@@ -370,7 +405,7 @@ def objective_rule(model):
             # Definition: optimizes for both ramp mitigation and peak shaving
             elif peak_shaving == 'hybrid': 
                 def objective_rule(model):
-                    return model.lamb * model.peak_scalar * sum([(model.d[t] + sum([model.u[t, v] for v in model.v]))**2 for t in model.t]) + \
+                    return model.lamb * model.peak_scalar * (sum([(model.d[t] + sum([model.u[t, v] for v in model.v]))**2 for t in model.t]) - model.peak_subtractor) + \
                         (1-model.lamb) * sum([(model.d[t + 1] - model.d[t] + sum([model.u[t + 1, v] - model.u[t, v] for v in model.v]))**2 for t in model.t if t != last_t])
                 model.objective = Objective(rule=objective_rule, sense=minimize, doc='Define objective function')
 
@@ -400,10 +435,10 @@ def max_call_rule(model, v):
             # Definition: minimizes charging cost based on given cost stack
             # Should probably only have visibility 24/48 hours into the future, ideally
             # V1G vs V2G will make a big difference here--incentivizing selling on-peak
-            # Need to feed in p (prices at t)
             elif peak_shaving == 'cost':
                 def objective_rule(model):
-                    return sum((model.u[t, v] * model.price[t] for v in model.v) for t in model.t)
+                    return sum([sum([model.u[t, v] * model.price[t] for t in model.t]) for v in model.v])
+                    #return sum((model.u[t, v] * model.price[t] for v in model.v) for t in model.t)
                 model.objective = Objective(rule=objective_rule, sense=minimize, doc='Define objective function')
             
             results = opt.solve(model)

From 8460a2e15ee910ffb30dddbebc080ad5357e1892 Mon Sep 17 00:00:00 2001
From: Margaret McCall <majemc@gmail.com>
Date: Thu, 28 May 2020 13:42:01 -0700
Subject: [PATCH 06/14] added cost optimization, modified hybrid

---
 v2gsim/post_simulation/netload_opti_plus.py | 34 ++++++++-------------
 1 file changed, 13 insertions(+), 21 deletions(-)

diff --git a/v2gsim/post_simulation/netload_opti_plus.py b/v2gsim/post_simulation/netload_opti_plus.py
index 165ccb4..148f618 100644
--- a/v2gsim/post_simulation/netload_opti_plus.py
+++ b/v2gsim/post_simulation/netload_opti_plus.py
@@ -31,7 +31,7 @@ def __init__(self, project, optimization_timestep, date_from,
         self.SOC_index_to = int((date_to - project.date).total_seconds() / project.timestep)
 
     def solve(self, project, net_load, real_number_of_vehicle, SOC_margin=0.02,
-              SOC_offset=0.0, peak_shaving='peak_shaving', penalization=5, beta=None, plot=False, peak_scalar = .01, peak_subtractor = 50000, price=pandas.DataFrame()): #price can't be array, is dict; can't have non-default arg follow default
+              SOC_offset=0.0, peak_shaving='peak_shaving', penalization=5, beta=None, plot=False, peak_scalar = .01, peak_subtractor = 50000, price=pandas.DataFrame(), calls=2): #price can't be array, is dict; can't have non-default arg follow default
         """Launch the optimization and the post_processing fucntion. Results
         and assumptions are appended to a data frame.
 
@@ -53,10 +53,11 @@ def solve(self, project, net_load, real_number_of_vehicle, SOC_margin=0.02,
         self.emin = {}
         self.emax = {}
         self.efinal = {}
+        self.max_calls = {}
 
         # Set the variables for the optimization
         new_net_load = self.initialize_net_load(net_load, real_number_of_vehicle, project)
-        self.initialize_model(project, new_net_load, SOC_margin, SOC_offset)
+        self.initialize_model(project, new_net_load, SOC_margin, SOC_offset, calls)
         
         if peak_shaving=='cost':
             price = self.initialize_price(price, net_load)
@@ -65,14 +66,14 @@ def solve(self, project, net_load, real_number_of_vehicle, SOC_margin=0.02,
         timer = time.time()
         opti_model, result = self.process(self.times, self.vehicles, self.d, self.pmax,
                                           self.pmin, self.emin, self.emax,
-                                          self.efinal, peak_shaving, penalization, peak_scalar, peak_subtractor, price)
+                                          self.efinal, peak_shaving, penalization, peak_scalar, peak_subtractor, price, calls)
         timer2 = time.time()
         print('The optimization duration was ' + str((timer2 - timer) / 60) + ' minutes')
         print('')
 
         # Post process results
-        return self.post_process(project, net_load, opti_model, result, plot)
-        #return opti_model, result
+        #return self.post_process(project, net_load, opti_model, result, plot)
+        return opti_model, result
 
     def initialize_net_load(self, net_load, real_number_of_vehicle, project):
         """Make sure that the net load has the right size and scale the net
@@ -213,7 +214,7 @@ def get_final_SOC(self, vehicle, SOC_margin, SOC_offset, SOC_end=None):
         else:
             return vehicle.SOC[self.SOC_index_to] - SOC_offset - SOC_margin
 
-    def initialize_model(self, project, net_load, SOC_margin, SOC_offset):
+    def initialize_model(self, project, net_load, SOC_margin, SOC_offset, calls):
         """Select the vehicles that were plugged at controlled chargers and create
         the optimization variables (see inputs of optimization)
 
@@ -274,13 +275,16 @@ def initialize_model(self, project, net_load, SOC_margin, SOC_offset):
                                                  (SOC_final - SOC_init) * vehicle.car_model.battery_capacity *
                                                  (60 / self.optimization_timestep))})
 
+                # Create max_calls for each vehicle for emergency optimization ######################
+                self.max_calls.update({vehicle.id: calls})
+
         print('There is ' + str(vehicle_to_optimize) + ' vehicle participating in the optimization (' +
               str(vehicle_to_optimize * 100 / len(project.vehicles)) + '%)')
         print('There is ' + str(unfeasible_vehicle) + ' unfeasible vehicle.')
         print('')
 
     def process(self, times, vehicles, d, pmax, pmin, emin, emax,
-                efinal, peak_shaving, penalization, peak_scalar, peak_subtractor, price, solver="gurobi"):
+                efinal, peak_shaving, penalization, peak_scalar, peak_subtractor, price, max_calls, solver="gurobi"):
 
         """The process function creates the pyomo model and solve it.
         Minimize sum( net_load(t) + sum(power_demand(t, v)))**2
@@ -333,7 +337,6 @@ def process(self, times, vehicles, d, pmax, pmin, emin, emax,
             model.e_max = Param(model.t, model.v, initialize=emax, doc='E max')
 
             model.e_final = Param(model.v, initialize=efinal, doc='final energy balance')
-            # model.beta = Param(initialize=beta, doc='beta')
 
             #HYBRID
             if peak_shaving == 'hybrid':
@@ -347,8 +350,7 @@ def process(self, times, vehicles, d, pmax, pmin, emin, emax,
             #EMERGENCY
             if peak_shaving == 'emergency':
                 # Maximum number of times a vehicle can be called for V2G in emergency scenario
-                model.max_calls = Param(model.v, initialize = 2) #trying to initialize max_calls for each vehicle being 2
-                model.zero = Param(model.v, initialize = 0) #trying to get vector to assess if u is negative
+                model.max_calls = Param(model.v, initialize = max_calls) #trying to initialize max_calls for each vehicle being 2...need to make look like efinal
 
             #COST
             if peak_shaving == 'cost':
@@ -417,18 +419,8 @@ def objective_rule(model):
                     return sum([(model.d[t] + sum([model.u[t, v] for v in model.v]))**2 for t in model.t])
                 model.objective = Objective(rule=objective_rule, sense=minimize, doc='Define objective function')
 
-                """def count_exports_rule(model, t, v):
-                    for t in model.t:
-                        for v in model.v:
-                            if value(model.u[t, v]) < 0:
-                                model.num_export[t, v] == 1  
-                            elif value(model.u[u, v]) >= 0:
-                                model.num_export[t, v] == 0
-                    return model.num_export
-                model.count_exports_rule = Constraint(model.t, model.v, rule=count_exports_rule, doc='Count exports rule')"""
-
                 def max_call_rule(model, v):
-                    return sum(value(model.u[t, v]) < 0 for t in model.t) <= model.max_calls[v]
+                    return len([x for x in [model.u[t, v] for t in model.t] if x < 0]) <= model.max_calls[v]
                 model.max_call_rule = Constraint(model.v, rule=max_call_rule, doc='Max call rule')
 
             # Cost-based

From 9a2be31f99515afb353195647b09123710589def Mon Sep 17 00:00:00 2001
From: Margaret McCall <majemc@gmail.com>
Date: Mon, 15 Jun 2020 11:29:44 -0700
Subject: [PATCH 07/14] placeholder before working on hybrid and emergency

---
 v2gsim/post_simulation/netload_opti_plus.py | 94 ++++++++++++++++++---
 1 file changed, 81 insertions(+), 13 deletions(-)

diff --git a/v2gsim/post_simulation/netload_opti_plus.py b/v2gsim/post_simulation/netload_opti_plus.py
index 148f618..a489a3e 100644
--- a/v2gsim/post_simulation/netload_opti_plus.py
+++ b/v2gsim/post_simulation/netload_opti_plus.py
@@ -1,4 +1,5 @@
 # -*- coding: utf-8 -*-
+
 from __future__ import division
 from pyomo.opt import SolverFactory
 from pyomo.environ import *
@@ -72,8 +73,8 @@ def solve(self, project, net_load, real_number_of_vehicle, SOC_margin=0.02,
         print('')
 
         # Post process results
-        #return self.post_process(project, net_load, opti_model, result, plot)
-        return opti_model, result
+        return self.post_process(project, net_load, opti_model, result, plot)
+        #return opti_model, result
 
     def initialize_net_load(self, net_load, real_number_of_vehicle, project):
         """Make sure that the net load has the right size and scale the net
@@ -308,9 +309,11 @@ def process(self, times, vehicles, d, pmax, pmin, emin, emax,
 
         Return:
             model (ConcreteModel), result
-        """
+            """
+        
 
         # Select gurobi solver
+         
         with SolverFactory(solver) as opt:
             # Solver option see Gurobi website
             # opt.options['Method'] = 1
@@ -347,21 +350,24 @@ def process(self, times, vehicles, d, pmax, pmin, emin, emax,
                 # Normalization factor for peak-shaving sub-objective
                 model.peak_subtractor = Param(initialize = peak_subtractor)
 
-            #EMERGENCY
-            if peak_shaving == 'emergency':
-                # Maximum number of times a vehicle can be called for V2G in emergency scenario
-                model.max_calls = Param(model.v, initialize = max_calls) #trying to initialize max_calls for each vehicle being 2...need to make look like efinal
-
             #COST
             if peak_shaving == 'cost':
                 # Price
                 model.price = Param(model.t, initialize=price, doc='Prices') #energy price at each timestep
 
+            #EMERGENCY
+            if peak_shaving == 'emergency':
+                # Maximum Wh a vehicle can export over course of optimization (battery capacity * given fraction)
+                model.max_calls = Param(model.v, initialize = max_calls)
 
 
             # ###### Variable
+            
+
+            #model.b = Var(model.t, model.v, domain=Binary, doc='Power exported at timestep?') #should it be within=Binary?
+            #model.u_neg = Var(model.t, model.v, domain=Integers, doc='Power used')
+            #model.u_pos = Var(model.t, model.v, domain=Integers, doc='Power used') 
             model.u = Var(model.t, model.v, domain=Integers, doc='Power used')
-            model.num_export = Var(model.t, model.v, domain=Binary, doc='Power exported at timestep?')
 
 
             # ###### Rules
@@ -372,6 +378,27 @@ def maximum_power_rule(model, t, v):
             def minimum_power_rule(model, t, v):
                 return model.u[t, v] >= model.p_min[t, v]
             model.power_min_rule = Constraint(model.t, model.v, rule=minimum_power_rule, doc='P min rule')
+            # you'd have four rules instead: 1 for charging (>0, <pmax), 1 for discharging (<0, >pmin)
+            # u is just the sum of those two things. basically the same thing, but now 3 decision vars
+
+            #do we need to worry about big neg + big +? constrain one to be 0 if other is not 0?
+            # if i do whole problem w u_neg and u_post instead of u, takes care of the issue
+            # where u_neg is super neg & vv
+            # bc when calculating the amount of energy in the battery, there's an efficiency term
+            # assigned to the in and the out...bc if both were values, there'd be more losses
+            # 
+
+            #def maximum_power_rule(model, t, v):
+            #    return model.u_pos[t, v] <= model.p_max[t, v]
+            #model.power_max_rule = Constraint(model.t, model.v, rule=maximum_power_rule, doc='P max rule')
+
+            #def positive_power_rule(model, t, v):
+             #   return model.u_pos[t, v] >= 0
+            #rederfine, and do for 2 negative ones -- model.power_max_rule = Constraint(model.t, model.v, rule=maximum_power_rule, doc='P max rule')
+
+            #def u_sum(model, t, v):
+             #   return model.u = model.u_neg + model.u_pos #...
+
 
             def minimum_energy_rule(model, t, v):
                 return sum(model.u[i, v] for i in range(0, t + 1)) >= model.e_min[t, v]
@@ -385,6 +412,23 @@ def final_energy_balance(model, v):
                 return sum(model.u[i, v] for i in model.t) >= model.e_final[v]
             model.final_energy_rule = Constraint(model.v, rule=final_energy_balance, doc='E final rule')
 
+            #def not_two_values:
+             #   return minimum( maximum(u_neg, 0), minimum(u_pos, 0))==0 # too hard...can you even put min/max in constraints? no?
+                # have to somehow penalize them with an efficiency thing to prevent simultaneous charging/discharging
+                # need variable that's gonna be relative to the losses...add a little bit of each charge/discharge
+                # to objective function (an additive term...just like adding a constant...doesn't have to be on same scale, but if scale too small, could be ~epsilon). 3 orders of magnitude
+                # def wanna look at u_neg/u_pos afterwards to see if one is 0 and other not*****
+                # check out google pyomo forum; respnsive
+                # if optimization works, it works
+                # could do fake losses terms where you multiply u_neg*0.9 &
+
+            #def power_export_rule(model, t, v):
+             #   return model.b[t, v] for model.u[t, v] < 0 == 1 ###this is the tricky part ## b is min of 0 or u. if else ok in constr
+                ## if u is +...if t is something, but NOT variables
+                # constraint: model.b <0
+                # constraint: model.b > u
+            
+
             # Set the objective to be either peak shaving or ramp mitigation
             if peak_shaving == 'peak_shaving':
                 def objective_rule(model):
@@ -417,11 +461,34 @@ def objective_rule(model):
             elif peak_shaving == 'emergency':
                 def objective_rule(model):
                     return sum([(model.d[t] + sum([model.u[t, v] for v in model.v]))**2 for t in model.t])
+                    # minimize obj function that has model.b term in it
+                    # u should be cut into two parts--charging and discharging
                 model.objective = Objective(rule=objective_rule, sense=minimize, doc='Define objective function')
 
-                def max_call_rule(model, v):
-                    return len([x for x in [model.u[t, v] for t in model.t] if x < 0]) <= model.max_calls[v]
-                model.max_call_rule = Constraint(model.v, rule=max_call_rule, doc='Max call rule')
+                
+                ### Put constraint on negative part!!
+                # define another variable: model.c (counting) 
+
+                #model.c = Var(...)
+
+              #  rule count(): # try to force count to be 0 unless it can be below 0...model.c > model.u...so when u is positive
+                ### make it follow u, and set clear limits (can't be below 0), and ensure that optimization will
+                # try to make it 0 most of the time, except when it has the ability to not make it 0, only when
+                # you want to count one of the session. only time optimization can minimize this term (ie, be something /= 0)
+                # is when u is negative
+
+#                def negative_power(model, v):
+             #       return sum([model.u_neg for ])
+
+                # Could use mod to constrain sessions of discharging...as soon as you have X timesteps
+                # consecutive that you're discharging...harder
+                # if you want to count how many times u_neg is below zero, easy--if < 0, +1 this var
+                # this would include counting variable in objective function...
+                
+                def max_calls_rule(model, v):
+                    return sum([model.b[t, v] for t in model.t]) <= model.max_calls[v] 
+                    #return len([x for x in [model.u[t, v] for t in model.t] if x < 0]) <= model.max_calls[v] 
+                model.max_calls_rule = Constraint(model.v, rule=max_calls_rule, doc='Max output rule')
 
             # Cost-based
             # Definition: minimizes charging cost based on given cost stack
@@ -433,7 +500,8 @@ def objective_rule(model):
                     #return sum((model.u[t, v] * model.price[t] for v in model.v) for t in model.t)
                 model.objective = Objective(rule=objective_rule, sense=minimize, doc='Define objective function')
             
-            results = opt.solve(model)
+            solver_parameters = "ResultFile=model.ilp"
+            results = opt.solve(model, options_string=solver_parameters, symbolic_solver_labels=True)
             # results.write()
 
         return model, results

From b991c04efe2673d13f781e0e6d4606f9271eca39 Mon Sep 17 00:00:00 2001
From: Margaret McCall <majemc@gmail.com>
Date: Mon, 15 Jun 2020 11:57:14 -0700
Subject: [PATCH 08/14] hybrid isolated, works

---
 v2gsim/post_simulation/netload_opti_hybrid.py | 693 ++++++++++++++++++
 1 file changed, 693 insertions(+)
 create mode 100644 v2gsim/post_simulation/netload_opti_hybrid.py

diff --git a/v2gsim/post_simulation/netload_opti_hybrid.py b/v2gsim/post_simulation/netload_opti_hybrid.py
new file mode 100644
index 0000000..6fd5a52
--- /dev/null
+++ b/v2gsim/post_simulation/netload_opti_hybrid.py
@@ -0,0 +1,693 @@
+# -*- coding: utf-8 -*-
+
+from __future__ import division
+from pyomo.opt import SolverFactory
+from pyomo.environ import *
+import time
+import pandas
+import numpy
+import matplotlib.pyplot as plt
+import seaborn as sns
+import v2gsim.model
+import v2gsim.result
+import datetime
+
+
+class CentralOptimization(object):
+    """Creates an object to perform optimization.
+    The object contains some general parameters for the optimization
+    """
+    def __init__(self, project, optimization_timestep, date_from,
+                 date_to, minimum_SOC=0.1, maximum_SOC=0.95):
+        # All the variables are at the project timestep except for the model variables
+        # optimization_timestep is in minutes
+        self.optimization_timestep = optimization_timestep
+        # Set minimum_SOC
+        self.minimum_SOC = minimum_SOC
+        self.maximum_SOC = maximum_SOC
+        # Set date boundaries, should be same as the one used during the simulation
+        self.date_from = date_from
+        self.date_to = date_to
+        self.SOC_index_from = int((date_from - project.date).total_seconds() / project.timestep)
+        self.SOC_index_to = int((date_to - project.date).total_seconds() / project.timestep)
+
+    def solve(self, project, net_load, real_number_of_vehicle, SOC_margin=0.02,
+              SOC_offset=0.0, peak_shaving='peak_shaving', penalization=5, beta=None, plot=False, peak_scalar = .01, peak_subtractor = 50000, price=pandas.DataFrame(), calls=2): #price can't be array, is dict; can't have non-default arg follow default
+        """Launch the optimization and the post_processing fucntion. Results
+        and assumptions are appended to a data frame.
+
+        Args:
+            project (Project): project
+            net_load (pandas.DataFrame): data frame with date index and a 'net_load' column in [W]
+            real_number_of_vehicle (int): number of vehicle expected on the net load,
+                False if number is the same as in project
+            SOC_margin (float): SOC margin that can be used by the optimization at the end of the day [0, 1]
+            SOC_offset (float): energy offset [0, 1]
+            peak_shaving (boolean): if True ramping constraints are not taking in account within the objective else it is.
+        """
+        # Reset model
+        self.times = []
+        self.vehicles = []
+        self.d = {}
+        self.pmax = {}
+        self.pmin = {}
+        self.emin = {}
+        self.emax = {}
+        self.efinal = {}
+        self.max_calls = {}
+
+        # Set the variables for the optimization
+        new_net_load = self.initialize_net_load(net_load, real_number_of_vehicle, project)
+        self.initialize_model(project, new_net_load, SOC_margin, SOC_offset, calls)
+        
+        if peak_shaving=='cost':
+            price = self.initialize_price(price, net_load)
+
+        # Run the optimization
+        timer = time.time()
+        opti_model, result = self.process(self.times, self.vehicles, self.d, self.pmax,
+                                          self.pmin, self.emin, self.emax,
+                                          self.efinal, peak_shaving, penalization, peak_scalar, peak_subtractor, price, calls)
+        timer2 = time.time()
+        print('The optimization duration was ' + str((timer2 - timer) / 60) + ' minutes')
+        print('')
+
+        # Post process results
+        return self.post_process(project, net_load, opti_model, result, plot)
+        #return opti_model, result
+
+    def initialize_net_load(self, net_load, real_number_of_vehicle, project):
+        """Make sure that the net load has the right size and scale the net
+        load for the optimization scale.
+
+        Args:
+            net_load (pandas.DataFrame): data frame with date index and a 'net_load' column in [W]
+            net_load_pmax (int): maximum power on the scaled net load
+        """
+        # Make sure we are not touching the initial data
+        new_net_load = net_load.copy()
+
+        # Resample the net load
+        new_net_load = new_net_load.resample(str(self.optimization_timestep) + 'T').first()
+
+        # Check against the actual lenght it should have
+        diff = (len(new_net_load) -
+                int((self.date_to - self.date_from).total_seconds() / (60 * self.optimization_timestep)))
+        if diff > 0:
+            # We should trim the net load with diff elements (normaly 1 because of slicing inclusion)
+            new_net_load.drop(new_net_load.tail(diff).index, inplace=True)
+        elif diff < 0:
+            print('The net load does not contain enough data points')
+
+        if real_number_of_vehicle:
+            # Set scaling factor
+            scaling_factor = len(project.vehicles) / real_number_of_vehicle
+
+            # Scale the temp net load
+            new_net_load['netload'] *= scaling_factor
+
+        return new_net_load
+
+    def initialize_price(self, price, net_load): ##>>>>
+        """Resample power price data to be on same frequency as net load data, and returns dictionary with timesteps instead of timestamps.
+        Args:
+            price (pandas.DataFrame): data frame with date index and a 'price' column (in $/MWh)
+            net_load (pandas.DataFrame): data frame with date index and a 'net_load' column in [W]
+        
+        Returns dictionary with 
+            """
+        # Check that price and net_load dataframes have same length
+        assert len(price) == len(net_load)
+        # Make sure we are not touching the initial data
+        new_price = price.copy()
+        # Resample the net load
+        new_price = new_price.resample(str(self.optimization_timestep) + 'T').first()
+
+
+        # Convert to timestep index instead of timestamps
+        temp_index = pandas.DataFrame(range(0, len(new_price)), columns=['index'])
+        # Set temp_index
+        temp_new_price = new_price.copy()
+        temp_new_price= temp_new_price.set_index(temp_index['index'])
+        # Return a dictionary
+        return temp_new_price.to_dict()['price']
+
+    def check_energy_constraints_feasible(self, vehicle, SOC_init, SOC_final, SOC_offset, verbose=False):
+        """Make sure that SOC final can be reached from SOC init under uncontrolled
+        charging (best case scenario). Print details when conditions are not met.
+
+        Args:
+            vehicle (Vehicle): vehicle
+            SOC_init (float): state of charge at the begining of the optimization [0, 1]
+            SOC_final (float): state of charge at the end of the optimization [0, 1]
+            SOC_offset (float): energy offset [0, 1]
+
+        Return:
+            (Boolean)
+        """
+        def print_status(SOC_final, SOC_init, diff):
+                print('Set battery gain: ' + str((SOC_final - SOC_init) * 100) + '%')
+                print('Simulation battery gain: ' + str(diff * 100))
+
+        # Check if below minimum SOC at any time
+        if (min(vehicle.SOC[self.SOC_index_from:self.SOC_index_to]) - SOC_offset) <= self.minimum_SOC:
+            if verbose:
+                print('Vehicle: ' + str(vehicle.id) + ' has a minimum SOC of ' +
+                      str(min(vehicle.SOC[self.SOC_index_from:self.SOC_index_to]) * 100) + '%')
+            return False
+
+        # Check SOC difference between date_from and date_to ?
+        # Diff represent the minimum loss or the maximum gain
+        diff = vehicle.SOC[self.SOC_index_to] - vehicle.SOC[self.SOC_index_from]
+        # Simulation show a battery gain
+        if diff > 0:
+            # Gain should be greater than the one we set up
+            if SOC_final - SOC_init < diff:
+                # Good to go
+                return True
+            else:
+                # Set final SOC to be under diff
+                print_status(SOC_final, SOC_init, diff)
+                return False
+        # Simulation show battery loss
+        else:
+            # energy balance should be negative
+            if SOC_final - SOC_init > 0:
+                print_status(SOC_final, SOC_init, diff)
+                return False
+            # Loss should be smaller than the one we set (less negative)
+            if diff < SOC_init - SOC_final:
+                return True
+            else:
+                # Set final SOC to include at least the lost
+                print_status(SOC_final, SOC_init, diff)
+                return False
+
+    def initialize_time_index(self, net_load):
+        """Replace date index by time ids
+
+        Args:
+            net_load (pandas.DataFrame): data frame with date index and a 'net_load' column in [W]
+
+        Return:
+            net_load as a dictionary with time ids, time ids
+        """
+        temp_index = pandas.DataFrame(range(0, len(net_load)), columns=['index'])
+        # Set temp_index
+        temp_net_load = net_load.copy()
+        temp_net_load = temp_net_load.set_index(temp_index['index'])
+        # Return a dictionary
+        return temp_net_load.to_dict()['netload'], temp_index.index.values.tolist()
+
+    def get_initial_SOC(self, vehicle, SOC_offset, SOC_init=None):
+        """Get the initial SOC with which people start the optimization
+        """
+        if SOC_init is not None:
+            return SOC_init
+        else:
+            return vehicle.SOC[self.SOC_index_from] - SOC_offset
+
+    def get_final_SOC(self, vehicle, SOC_margin, SOC_offset, SOC_end=None):
+        """Get final SOC that vehicle must reached at the end of the optimization
+        """
+        if SOC_end is not None:
+            return SOC_end
+        else:
+            return vehicle.SOC[self.SOC_index_to] - SOC_offset - SOC_margin
+
+    def initialize_model(self, project, net_load, SOC_margin, SOC_offset, calls):
+        """Select the vehicles that were plugged at controlled chargers and create
+        the optimization variables (see inputs of optimization)
+
+        Args:
+            project (Project): project
+
+        Return:
+            times, vehicles, d, pmax, pmin, emin, emax, efinal
+        """
+        # Create a dict with the net load and get time index in a data frame
+        self.d, self.times = self.initialize_time_index(net_load)
+        vehicle_to_optimize = 0
+        unfeasible_vehicle = 0
+
+        for vehicle in project.vehicles:
+            if vehicle.result is not None:
+                # Get SOC init and SOC end
+                SOC_init = self.get_initial_SOC(vehicle, SOC_offset)
+                SOC_final = self.get_final_SOC(vehicle, SOC_margin, SOC_offset)
+
+                # Find out if vehicle itinerary is feasible
+                if not self.check_energy_constraints_feasible(vehicle, SOC_init, SOC_final, SOC_offset):
+                    # Reset vehicle result to None
+                    vehicle.result = None
+                    unfeasible_vehicle += 1
+                    continue
+
+                # Add vehicle id to a list
+                self.vehicles.append(vehicle.id)
+                vehicle_to_optimize += 1
+
+                # Resample vehicle result
+                temp_vehicle_result = vehicle.result.resample(str(self.optimization_timestep) + 'T').first()
+
+                # Set time_vehicle_index
+                temp_vehicle_result = temp_vehicle_result.set_index(pandas.DataFrame(
+                    index=[(time, vehicle.id) for time in self.times]).index)
+
+                # Push pmax and pmin with vehicle and time key
+                self.pmin.update(temp_vehicle_result.to_dict()['p_min'])
+                self.pmax.update(temp_vehicle_result.to_dict()['p_max'])
+
+                # Push emin and emax with vehicle and time key
+                # Units! if project.timestep in seconds, self.timestep in minutes and battery in Wh
+                # Units! Wproject.timestep --> Wself.timestep * (project.timestep / (60 * self.timestep))
+                # Units! Wtimestep --> Wh * (60 / self.timestep)
+                temp_vehicle_result['emin'] = (temp_vehicle_result.energy * (project.timestep / (60 * self.optimization_timestep)) +
+                                               (self.minimum_SOC - SOC_init) * vehicle.car_model.battery_capacity *
+                                               (60 / self.optimization_timestep))
+                self.emin.update(temp_vehicle_result.to_dict()['emin'])
+                temp_vehicle_result['emax'] = (temp_vehicle_result.energy * (project.timestep / (60 * self.optimization_timestep)) + 10000 +
+                                               (self.maximum_SOC - SOC_init) * vehicle.car_model.battery_capacity *
+                                               (60 / self.optimization_timestep))
+                self.emax.update(temp_vehicle_result.to_dict()['emax'])
+
+                # Push efinal with vehicle key
+                self.efinal.update({vehicle.id: (temp_vehicle_result.tail(1).energy.values[0] * (project.timestep / (60 * self.optimization_timestep)) +
+                                                 (SOC_final - SOC_init) * vehicle.car_model.battery_capacity *
+                                                 (60 / self.optimization_timestep))})
+
+                # Create max_calls for each vehicle for emergency optimization ######################
+                self.max_calls.update({vehicle.id: calls})
+
+        print('There is ' + str(vehicle_to_optimize) + ' vehicle participating in the optimization (' +
+              str(vehicle_to_optimize * 100 / len(project.vehicles)) + '%)')
+        print('There is ' + str(unfeasible_vehicle) + ' unfeasible vehicle.')
+        print('')
+
+    def process(self, times, vehicles, d, pmax, pmin, emin, emax,
+                efinal, peak_shaving, penalization, peak_scalar, peak_subtractor, price, max_calls, solver="gurobi"):
+
+        """The process function creates the pyomo model and solve it.
+        Minimize sum( net_load(t) + sum(power_demand(t, v)))**2
+        subject to:
+        pmin(t, v) <= power_demand(t, v) <= pmax(t, v)
+        emin(t, v) <= sum(power_demand(t, v)) <= emax(t, v)
+        sum(power_demand(t, v)) >= efinal(v)
+        rampmin(t) <= net_load_ramp(t) + power_demand_ramp(t, v) <= rampmax(t)
+
+        Args:
+            times (list): timestep list
+            vehicles (list): unique list of vehicle ids
+            d (dict): time - net load at t
+            price (dict): time - power price at t ###
+            pmax (dict): (time, id) - power maximum at t for v
+            pmin (dict): (time, id) - power minimum at t for v
+            emin (dict): (time, id) - energy minimum at t for v
+            emax (dict): (time, id) - energy maximum at t for v
+            efinal (dict): id - final SOC
+            solver (string): name of the solver to use (default is gurobi)
+
+        Return:
+            model (ConcreteModel), result
+            """
+        
+
+        # Select gurobi solver
+         
+        with SolverFactory(solver) as opt:
+            # Solver option see Gurobi website
+            # opt.options['Method'] = 1
+
+            # Creation of a Concrete Model
+            model = ConcreteModel()
+
+            # ###### Set
+            model.t = Set(initialize=times, doc='Time', ordered=True)
+            last_t = model.t.last()
+            model.v = Set(initialize=vehicles, doc='Vehicles')
+
+
+            # ###### Parameters
+            # Net load
+            model.d = Param(model.t, initialize=d, doc='Net load')
+
+            # Power
+            model.p_max = Param(model.t, model.v, initialize=pmax, doc='P max')
+            model.p_min = Param(model.t, model.v, initialize=pmin, doc='P min')
+
+            # Energy
+            model.e_min = Param(model.t, model.v, initialize=emin, doc='E min')
+            model.e_max = Param(model.t, model.v, initialize=emax, doc='E max')
+
+            model.e_final = Param(model.v, initialize=efinal, doc='final energy balance')
+
+            #HYBRID
+            if peak_shaving == 'hybrid':
+                # Lambda for hybrid model
+                model.lamb = Param(initialize = .5)
+                # Scaling factor for peak-shaving sub-objective for hybrid model
+                model.peak_scalar = Param(initialize = peak_scalar)
+                # Normalization factor for peak-shaving sub-objective
+                model.peak_subtractor = Param(initialize = peak_subtractor)
+
+            #COST
+            if peak_shaving == 'cost':
+                # Price
+                model.price = Param(model.t, initialize=price, doc='Prices') #energy price at each timestep
+
+            #EMERGENCY
+            if peak_shaving == 'emergency':
+                # Maximum Wh a vehicle can export over course of optimization (battery capacity * given fraction)
+                model.max_calls = Param(model.v, initialize = max_calls)
+
+
+            # ###### Variable
+            
+
+            #model.b = Var(model.t, model.v, domain=Binary, doc='Power exported at timestep?') #should it be within=Binary?
+            #model.u_neg = Var(model.t, model.v, domain=Integers, doc='Power used')
+            #model.u_pos = Var(model.t, model.v, domain=Integers, doc='Power used') 
+            model.u = Var(model.t, model.v, domain=Integers, doc='Power used')
+
+
+            # ###### Rules
+            def maximum_power_rule(model, t, v):
+                return model.u[t, v] <= model.p_max[t, v]
+            model.power_max_rule = Constraint(model.t, model.v, rule=maximum_power_rule, doc='P max rule')
+
+            def minimum_power_rule(model, t, v):
+                return model.u[t, v] >= model.p_min[t, v]
+            model.power_min_rule = Constraint(model.t, model.v, rule=minimum_power_rule, doc='P min rule')
+            # you'd have four rules instead: 1 for charging (>0, <pmax), 1 for discharging (<0, >pmin)
+            # u is just the sum of those two things. basically the same thing, but now 3 decision vars
+
+            #do we need to worry about big neg + big +? constrain one to be 0 if other is not 0?
+            # if i do whole problem w u_neg and u_post instead of u, takes care of the issue
+            # where u_neg is super neg & vv
+            # bc when calculating the amount of energy in the battery, there's an efficiency term
+            # assigned to the in and the out...bc if both were values, there'd be more losses
+            # 
+
+            #def maximum_power_rule(model, t, v):
+            #    return model.u_pos[t, v] <= model.p_max[t, v]
+            #model.power_max_rule = Constraint(model.t, model.v, rule=maximum_power_rule, doc='P max rule')
+
+            #def positive_power_rule(model, t, v):
+             #   return model.u_pos[t, v] >= 0
+            #rederfine, and do for 2 negative ones -- model.power_max_rule = Constraint(model.t, model.v, rule=maximum_power_rule, doc='P max rule')
+
+            #def u_sum(model, t, v):
+             #   return model.u = model.u_neg + model.u_pos #...
+
+
+            def minimum_energy_rule(model, t, v):
+                return sum(model.u[i, v] for i in range(0, t + 1)) >= model.e_min[t, v]
+            model.minimum_energy_rule = Constraint(model.t, model.v, rule=minimum_energy_rule, doc='E min rule')
+
+            def maximum_energy_rule(model, t, v):
+                return sum(model.u[i, v] for i in range(0, t + 1)) <= model.e_max[t, v]
+            model.maximum_energy_rule = Constraint(model.t, model.v, rule=maximum_energy_rule, doc='E max rule')
+
+            def final_energy_balance(model, v):
+                return sum(model.u[i, v] for i in model.t) >= model.e_final[v]
+            model.final_energy_rule = Constraint(model.v, rule=final_energy_balance, doc='E final rule')
+
+            #def not_two_values:
+             #   return minimum( maximum(u_neg, 0), minimum(u_pos, 0))==0 # too hard...can you even put min/max in constraints? no?
+                # have to somehow penalize them with an efficiency thing to prevent simultaneous charging/discharging
+                # need variable that's gonna be relative to the losses...add a little bit of each charge/discharge
+                # to objective function (an additive term...just like adding a constant...doesn't have to be on same scale, but if scale too small, could be ~epsilon). 3 orders of magnitude
+                # def wanna look at u_neg/u_pos afterwards to see if one is 0 and other not*****
+                # check out google pyomo forum; respnsive
+                # if optimization works, it works
+                # could do fake losses terms where you multiply u_neg*0.9 &
+
+            #def power_export_rule(model, t, v):
+             #   return model.b[t, v] for model.u[t, v] < 0 == 1 ###this is the tricky part ## b is min of 0 or u. if else ok in constr
+                ## if u is +...if t is something, but NOT variables
+                # constraint: model.b <0
+                # constraint: model.b > u
+            
+
+            # Set the objective to be either peak shaving or ramp mitigation
+            if peak_shaving == 'peak_shaving':
+                def objective_rule(model):
+                    return sum([(model.d[t] + sum([model.u[t, v] for v in model.v]))**2 for t in model.t])
+                model.objective = Objective(rule=objective_rule, sense=minimize, doc='Define objective function')
+
+            elif peak_shaving == 'penalized_peak_shaving':
+                def objective_rule(model):
+                    return (sum([(model.d[t] + sum([model.u[t, v] for v in model.v]))**2 for t in model.t]) +
+                            penalization * sum([sum([model.u[t, v]**2 for v in model.v]) for t in model.t]))
+                model.objective = Objective(rule=objective_rule, sense=minimize, doc='Define objective function')
+
+            elif peak_shaving == 'ramp_mitigation':
+                def objective_rule(model):
+                    return sum([(model.d[t + 1] - model.d[t] + sum([model.u[t + 1, v] - model.u[t, v] for v in model.v]))**2 for t in model.t if t != last_t])
+                model.objective = Objective(rule=objective_rule, sense=minimize, doc='Define objective function')
+
+            ### MM additions
+            # Hybrid
+            # Definition: optimizes for both ramp mitigation and peak shaving
+            elif peak_shaving == 'hybrid': 
+                def objective_rule(model):
+                    return model.lamb * model.peak_scalar * (sum([(model.d[t] - model.peak_subtractor + sum([model.u[t, v] for v in model.v]))**2 for t in model.t]) ) + \
+                        (1-model.lamb) * sum([(model.d[t + 1] - model.d[t] + sum([model.u[t + 1, v] - model.u[t, v] for v in model.v]))**2 for t in model.t if t != last_t])
+                model.objective = Objective(rule=objective_rule, sense=minimize, doc='Define objective function')
+
+            # Emergency
+            # Definition: Emergency V2G; mostly V1G but each vehicle can get called for V2G a certain number of times/year
+            # Can be used either for peak shaving or ramp mitigation; for now just implementing peak shaving
+            elif peak_shaving == 'emergency':
+                def objective_rule(model):
+                    return sum([(model.d[t] + sum([model.u[t, v] for v in model.v]))**2 for t in model.t])
+                    # minimize obj function that has model.b term in it
+                    # u should be cut into two parts--charging and discharging
+                model.objective = Objective(rule=objective_rule, sense=minimize, doc='Define objective function')
+
+                
+                ### Put constraint on negative part!!
+                # define another variable: model.c (counting) 
+
+                #model.c = Var(...)
+
+              #  rule count(): # try to force count to be 0 unless it can be below 0...model.c > model.u...so when u is positive
+                ### make it follow u, and set clear limits (can't be below 0), and ensure that optimization will
+                # try to make it 0 most of the time, except when it has the ability to not make it 0, only when
+                # you want to count one of the session. only time optimization can minimize this term (ie, be something /= 0)
+                # is when u is negative
+
+#                def negative_power(model, v):
+             #       return sum([model.u_neg for ])
+
+                # Could use mod to constrain sessions of discharging...as soon as you have X timesteps
+                # consecutive that you're discharging...harder
+                # if you want to count how many times u_neg is below zero, easy--if < 0, +1 this var
+                # this would include counting variable in objective function...
+                
+                def max_calls_rule(model, v):
+                    return sum([model.b[t, v] for t in model.t]) <= model.max_calls[v] 
+                    #return len([x for x in [model.u[t, v] for t in model.t] if x < 0]) <= model.max_calls[v] 
+                model.max_calls_rule = Constraint(model.v, rule=max_calls_rule, doc='Max output rule')
+
+            # Cost-based
+            # Definition: minimizes charging cost based on given cost stack
+            # Should probably only have visibility 24/48 hours into the future, ideally
+            # V1G vs V2G will make a big difference here--incentivizing selling on-peak
+            elif peak_shaving == 'cost':
+                def objective_rule(model):
+                    return sum([sum([model.u[t, v] * model.price[t] for t in model.t]) for v in model.v])
+                    #return sum((model.u[t, v] * model.price[t] for v in model.v) for t in model.t)
+                model.objective = Objective(rule=objective_rule, sense=minimize, doc='Define objective function')
+            
+            solver_parameters = "ResultFile=model.ilp"
+            results = opt.solve(model, options_string=solver_parameters, symbolic_solver_labels=True)
+            # results.write()
+
+        return model, results
+
+    def plot_result(self, model):
+        """Create a plot showing the power constraints, the energy constraints and the ramp
+        constraints as well as the final net load.
+        """
+        # Set the graph style
+        sns.set_style("whitegrid")
+        sns.despine()
+
+        result = pandas.DataFrame()
+
+        # Get the result
+        df = pandas.DataFrame(index=['power'], data=model.u.get_values()).transpose().groupby(level=0).sum()
+        # Ramp of the result
+        mylist = [0]
+        mylist.extend(list(numpy.diff(df['power'].values)))
+        df['ramppower'] = mylist
+        result = pandas.concat([result, df], axis=1)
+        # cum sum of the result
+        df = pandas.DataFrame(
+            pandas.DataFrame(
+                index=['anything'],
+                data=model.u.get_values()).transpose().unstack().cumsum().sum(axis=1), columns=['powercum']) * (5 / 60)
+        result = pandas.concat([result, df], axis=1)
+
+        # Get pmax and pmin units of power [W]
+        df = pandas.DataFrame(index=['pmax'], data=model.p_max.extract_values()).transpose().groupby(level=0).sum()
+        result = pandas.concat([result, df], axis=1)
+        df = pandas.DataFrame(index=['pmin'], data=model.p_min.extract_values()).transpose().groupby(level=0).sum()
+        result = pandas.concat([result, df], axis=1)
+
+        # Get emin and emax units of energy [Wh]
+        df = pandas.DataFrame(index=['emax'], data=model.e_max.extract_values()).transpose().groupby(level=0).sum() * (5 / 60)
+        result = pandas.concat([result, df], axis=1)
+        df = pandas.DataFrame(index=['emin'], data=model.e_min.extract_values()).transpose().groupby(level=0).sum() * (5 / 60)
+        result = pandas.concat([result, df], axis=1)
+
+        # Get the minimum final energy quantity
+        e_final = pandas.DataFrame(index=['efinal'], data=model.e_final.extract_values()).transpose().sum() * (5 / 60)
+
+        # Get the actual ramp
+        mylist = [0]
+        df = pandas.DataFrame(index=['net_load'], data=model.d.extract_values()).transpose()
+        mylist.extend(list(numpy.diff(df['net_load'].values)))
+        df['rampnet_load'] = mylist
+        result = pandas.concat([result, df], axis=1)
+
+        # Plot power constraints
+        plt.subplot(411)
+        plt.plot(result.index.values, result.pmax.values, label='pmax')
+        plt.plot(result.index.values, result.power.values, label='power')
+        plt.plot(result.index.values, result.pmin.values, label='pmin')
+        plt.legend(loc=0)
+
+        # Plot energy constraints
+        plt.subplot(412)
+        plt.plot(result.index.values, result.emax.values, label='emax')
+        plt.plot(result.index.values, result.powercum.values, label='powercum')
+        plt.plot(result.index.values, result.emin.values, label='emin')
+        plt.plot(result.index[-1], e_final, '*', markersize=15, label='efinal')
+        plt.legend(loc=0)
+
+        # Plot the result ramp
+        plt.subplot(413)
+        plt.plot(result.index.values, result.ramppower.values + result.rampnet_load.values, label='result ramp')
+        plt.plot(result.index.values, result.rampnet_load.values, label='net_load ramp')
+        plt.legend(loc=0)
+
+        # Plot the power demand results
+        plt.subplot(414)
+        plt.plot(result.index.values, result.net_load.values, label='net_load')
+        plt.plot(result.index.values, result.net_load.values + result.power.values, label='net_load + vehicle')
+        plt.legend(loc=0)
+        plt.show()
+
+    def post_process(self, project, netload, model, result, plot):
+        """Recompute SOC profiles and compute new total power demand
+
+        Args:
+            project (Project): project
+
+        Note: Should check that 'vehicle before' and 'after' contain the same number of vehicles
+        """
+        if plot:
+            self.plot_result(model)
+
+        temp = pandas.DataFrame()
+        first = True
+        for vehicle in project.vehicles:
+            if vehicle.result is not None:
+                if first:
+                    temp['vehicle_before'] = vehicle.result['power_demand']
+                    first = False
+                else:
+                    temp['vehicle_before'] += vehicle.result['power_demand']
+
+        temp2 = pandas.DataFrame(index=['vehicle_after'], data=model.u.get_values()).transpose().groupby(level=0).sum()
+        i = pandas.date_range(start=self.date_from, end=self.date_to,
+                              freq=str(self.optimization_timestep) + 'T', closed='left')
+        temp2 = temp2.set_index(i)
+        temp2 = temp2.resample(str(project.timestep) + 'S')
+        temp2 = temp2.fillna(method='ffill').fillna(method='bfill')
+
+        final_result = pandas.DataFrame()
+        final_result = pandas.concat([temp['vehicle_before'], temp2['vehicle_after']], axis=1)
+        final_result = final_result.fillna(method='ffill').fillna(method='bfill')
+
+        return final_result
+
+
+def save_vehicle_state_for_optimization(vehicle, timestep, date_from,
+                                        date_to, activity=None, power_demand=None,
+                                        SOC=None, detail=None, nb_interval=None, init=False,
+                                        run=False, post=False):
+    """Save results for individual vehicles. Power demand is positive when charging
+    negative when driving. Energy consumption is positive when driving and negative
+    when charging. Charging station that offer after simulation processing should
+    have activity.charging_station.post_simulation True.
+    """
+    if run:
+        if vehicle.result is not None:
+            activity_index1, activity_index2, location_index1, location_index2, save = v2gsim.result._map_index(
+                activity.start, activity.end, date_from, date_to, len(power_demand),
+                len(vehicle.result['power_demand']), timestep)
+            # Time frame are matching
+            if save:
+                # If driving pmin and pmax are equal to 0 since we are not plugged
+                if isinstance(activity, v2gsim.model.Driving):
+                    vehicle.result['p_max'][location_index1:location_index2] -= (
+                        [0.0] * (activity_index2 - activity_index1))
+                    vehicle.result['p_min'][location_index1:location_index2] -= (
+                        [0.0] * (activity_index2 - activity_index1))
+                    # Energy consumed is directly the power demand (sum later)
+                    vehicle.result['energy'][location_index1:location_index2] += (
+                        power_demand[activity_index1:activity_index2])
+                    # Power demand on the grid is 0 since we are driving
+                    vehicle.result['power_demand'][location_index1:location_index2] -= (
+                        [0.0] * (activity_index2 - activity_index1))
+
+                # If parked pmin and pmax are not necessary the same
+                if isinstance(activity, v2gsim.model.Parked):
+                    # Save the positive power demand of this specific vehicle
+                    vehicle.result['power_demand'][location_index1:location_index2] += (
+                        power_demand[activity_index1:activity_index2])
+                    if activity.charging_station.post_simulation:
+                        # Find if vehicle or infra is limiting
+                        pmax = min(activity.charging_station.maximum_power,
+                                   vehicle.car_model.maximum_power)
+                        pmin = max(activity.charging_station.minimum_power,
+                                   vehicle.car_model.minimum_power)
+                        vehicle.result['p_max'][location_index1:location_index2] += (
+                            [pmax] * (activity_index2 - activity_index1))
+                        vehicle.result['p_min'][location_index1:location_index2] += (
+                            [pmin] * (activity_index2 - activity_index1))
+                        # Energy consumed is 0 the optimization will decide
+                        vehicle.result['energy'][location_index1:location_index2] -= (
+                            [0.0] * (activity_index2 - activity_index1))
+                    else:
+                        vehicle.result['p_max'][location_index1:location_index2] += (
+                            power_demand[activity_index1:activity_index2])
+                        vehicle.result['p_min'][location_index1:location_index2] += (
+                            power_demand[activity_index1:activity_index2])
+                        # Energy is 0.0 because it's already accounted in power_demand
+                        vehicle.result['energy'][location_index1:location_index2] -= (
+                            [0.0] * (activity_index2 - activity_index1))
+
+    elif init:
+        vehicle.SOC = [vehicle.SOC[0]]
+        vehicle.result = None
+        for activity in vehicle.activities:
+            if isinstance(activity, v2gsim.model.Parked):
+                if activity.charging_station.post_simulation:
+                    # Initiate a dictionary of numpy array to hold result (faster than DataFrame)
+                    vehicle.result = {'power_demand': numpy.array([0.0] * int((date_to - date_from).total_seconds() / timestep)),
+                                      'p_max': numpy.array([0.0] * int((date_to - date_from).total_seconds() / timestep)),
+                                      'p_min': numpy.array([0.0] * int((date_to - date_from).total_seconds() / timestep)),
+                                      'energy': numpy.array([0.0] * int((date_to - date_from).total_seconds() / timestep))}
+                    # Leave the init function
+                    return
+    elif post:
+        if vehicle.result is not None:
+            # Convert location result back into pandas DataFrame (faster that way)
+            i = pandas.date_range(start=date_from, end=date_to,
+                                  freq=str(timestep) + 's', closed='left')
+            vehicle.result = pandas.DataFrame(index=i, data=vehicle.result)
+            vehicle.result['energy'] = vehicle.result['energy'].cumsum()

From e003242bd4d3cc38c7571db0d2bccb876d8fefb8 Mon Sep 17 00:00:00 2001
From: Margaret McCall <majemc@gmail.com>
Date: Mon, 15 Jun 2020 15:12:10 -0700
Subject: [PATCH 09/14] working hybrid model

---
 v2gsim/post_simulation/netload_opti_hybrid.py | 55 ++++++++++++++++---
 1 file changed, 47 insertions(+), 8 deletions(-)

diff --git a/v2gsim/post_simulation/netload_opti_hybrid.py b/v2gsim/post_simulation/netload_opti_hybrid.py
index 6fd5a52..185625f 100644
--- a/v2gsim/post_simulation/netload_opti_hybrid.py
+++ b/v2gsim/post_simulation/netload_opti_hybrid.py
@@ -31,8 +31,8 @@ def __init__(self, project, optimization_timestep, date_from,
         self.SOC_index_from = int((date_from - project.date).total_seconds() / project.timestep)
         self.SOC_index_to = int((date_to - project.date).total_seconds() / project.timestep)
 
-    def solve(self, project, net_load, real_number_of_vehicle, SOC_margin=0.02,
-              SOC_offset=0.0, peak_shaving='peak_shaving', penalization=5, beta=None, plot=False, peak_scalar = .01, peak_subtractor = 50000, price=pandas.DataFrame(), calls=2): #price can't be array, is dict; can't have non-default arg follow default
+    def solve(self, project, net_load, unoptimized_demand, real_number_of_vehicle, SOC_margin=0.02,
+              SOC_offset=0.0, peak_shaving='peak_shaving', penalization=5, beta=None, plot=False, peak_scalar=.01, peak_subtractor=30000, price=pandas.DataFrame(), calls=2): #price can't be array, is dict; can't have non-default arg follow default
         """Launch the optimization and the post_processing fucntion. Results
         and assumptions are appended to a data frame.
 
@@ -57,12 +57,17 @@ def solve(self, project, net_load, real_number_of_vehicle, SOC_margin=0.02,
         self.max_calls = {}
 
         # Set the variables for the optimization
-        new_net_load = self.initialize_net_load(net_load, real_number_of_vehicle, project)
+        new_net_load, new_unoptimized_demand = self.initialize_net_load(net_load, real_number_of_vehicle, project, unoptimized_demand)
         self.initialize_model(project, new_net_load, SOC_margin, SOC_offset, calls)
         
         if peak_shaving=='cost':
             price = self.initialize_price(price, net_load)
 
+        elif peak_shaving=='hybrid':
+            peak_scalar, peak_subtractor = self.get_peak_scalar(new_net_load, new_unoptimized_demand)
+            print(peak_scalar)
+            print(peak_subtractor)
+
         # Run the optimization
         timer = time.time()
         opti_model, result = self.process(self.times, self.vehicles, self.d, self.pmax,
@@ -76,7 +81,37 @@ def solve(self, project, net_load, real_number_of_vehicle, SOC_margin=0.02,
         return self.post_process(project, net_load, opti_model, result, plot)
         #return opti_model, result
 
-    def initialize_net_load(self, net_load, real_number_of_vehicle, project):
+    def get_peak_scalar(self, net_load_updated, unoptimized_demand):
+        """
+        Args: 
+            net_load_updated: series of net load without vehicle loads
+            unoptimized_demand: series of unoptimized vehicle loads, scaled appropriately
+        Returns:
+            peak_scalar (multiplier for peak-minimization argument in objective fxn)
+            peak_subtractor (to be subtracted from peak-min arg in obj fxn, to center around 0)
+        
+        """
+        #Combining net load and unoptimized demand
+        new_load = pandas.DataFrame()
+        new_load['load'] = numpy.array(net_load_updated['netload']) + numpy.array(unoptimized_demand)
+        
+        # Normalizing peak load (subtracting mean)
+        peak_subtractor = numpy.mean(new_load['load'])
+        norm_load = new_load['load'] - peak_subtractor
+
+        # Scaling factor for peak shaving: sum of squared values at each timestep
+        peak_scale = sum(x**2 for x in norm_load)
+
+        # Scaling factor for ramp mitigation: sum of squared ramp rate between all timesteps
+        ramp_rate = new_load.subtract(new_load.shift(-1)).dropna()
+        ramp_scale = sum(x**2 for x in ramp_rate['load'])
+
+        # Getting comparative scaling factor for peak-shaving side of hybrid optimization objective
+        peak_scalar = ramp_scale / peak_scale
+        
+        return peak_scalar, peak_subtractor
+    
+    def initialize_net_load(self, net_load, real_number_of_vehicle, project, unoptimized_demand):
         """Make sure that the net load has the right size and scale the net
         load for the optimization scale.
 
@@ -87,8 +122,12 @@ def initialize_net_load(self, net_load, real_number_of_vehicle, project):
         # Make sure we are not touching the initial data
         new_net_load = net_load.copy()
 
-        # Resample the net load
+        # Resample the net load & unoptimized demand
         new_net_load = new_net_load.resample(str(self.optimization_timestep) + 'T').first()
+        new_unoptimized_demand = unoptimized_demand.resample(str(self.optimization_timestep) + 'T').first()
+
+        # Check that new net load and new unoptimized demand are same length
+        assert len(new_net_load) == len(new_unoptimized_demand)
 
         # Check against the actual lenght it should have
         diff = (len(new_net_load) -
@@ -106,7 +145,7 @@ def initialize_net_load(self, net_load, real_number_of_vehicle, project):
             # Scale the temp net load
             new_net_load['netload'] *= scaling_factor
 
-        return new_net_load
+        return new_net_load, new_unoptimized_demand
 
     def initialize_price(self, price, net_load): ##>>>>
         """Resample power price data to be on same frequency as net load data, and returns dictionary with timesteps instead of timestamps.
@@ -343,8 +382,8 @@ def process(self, times, vehicles, d, pmax, pmin, emin, emax,
 
             #HYBRID
             if peak_shaving == 'hybrid':
-                # Lambda for hybrid model
-                model.lamb = Param(initialize = .5)
+                # Lambda for hybrid model--1 is all weight on peak shaving
+                model.lamb = Param(initialize = .9)
                 # Scaling factor for peak-shaving sub-objective for hybrid model
                 model.peak_scalar = Param(initialize = peak_scalar)
                 # Normalization factor for peak-shaving sub-objective

From 4789ad64a1cbf8106f7541699db46e8bb6a385f9 Mon Sep 17 00:00:00 2001
From: Margaret McCall <majemc@gmail.com>
Date: Mon, 15 Jun 2020 17:19:24 -0700
Subject: [PATCH 10/14] emergency model in progress

---
 .../post_simulation/netload_opti_emergency.py | 689 ++++++++++++++++++
 1 file changed, 689 insertions(+)
 create mode 100644 v2gsim/post_simulation/netload_opti_emergency.py

diff --git a/v2gsim/post_simulation/netload_opti_emergency.py b/v2gsim/post_simulation/netload_opti_emergency.py
new file mode 100644
index 0000000..e5cb03e
--- /dev/null
+++ b/v2gsim/post_simulation/netload_opti_emergency.py
@@ -0,0 +1,689 @@
+# -*- coding: utf-8 -*-
+
+from __future__ import division
+from pyomo.opt import SolverFactory
+from pyomo.environ import *
+import time
+import pandas
+import numpy
+import matplotlib.pyplot as plt
+import seaborn as sns
+import v2gsim.model
+import v2gsim.result
+import datetime
+
+
+class CentralOptimization(object):
+    """Creates an object to perform optimization.
+    The object contains some general parameters for the optimization
+    """
+    def __init__(self, project, optimization_timestep, date_from,
+                 date_to, minimum_SOC=0.1, maximum_SOC=0.95):
+        # All the variables are at the project timestep except for the model variables
+        # optimization_timestep is in minutes
+        self.optimization_timestep = optimization_timestep
+        # Set minimum_SOC
+        self.minimum_SOC = minimum_SOC
+        self.maximum_SOC = maximum_SOC
+        # Set date boundaries, should be same as the one used during the simulation
+        self.date_from = date_from
+        self.date_to = date_to
+        self.SOC_index_from = int((date_from - project.date).total_seconds() / project.timestep)
+        self.SOC_index_to = int((date_to - project.date).total_seconds() / project.timestep)
+
+    def solve(self, project, net_load, real_number_of_vehicle, SOC_margin=0.02,
+              SOC_offset=0.0, peak_shaving='peak_shaving', penalization=5, beta=None, plot=False, peak_scalar=.01, peak_subtractor=30000, 
+              unoptimized_demand = pandas.DataFrame(), price=pandas.DataFrame(), calls=-3000): #price can't be array, is dict; can't have non-default arg follow default
+        """Launch the optimization and the post_processing fucntion. Results
+        and assumptions are appended to a data frame.
+
+        Args:
+            project (Project): project
+            net_load (pandas.DataFrame): data frame with date index and a 'net_load' column in [W]
+            real_number_of_vehicle (int): number of vehicle expected on the net load,
+                False if number is the same as in project
+            SOC_margin (float): SOC margin that can be used by the optimization at the end of the day [0, 1]
+            SOC_offset (float): energy offset [0, 1]
+            peak_shaving (boolean): if True ramping constraints are not taking in account within the objective else it is.
+        """
+        # Reset model
+        self.times = []
+        self.vehicles = []
+        self.d = {}
+        self.pmax = {}
+        self.pmin = {}
+        self.emin = {}
+        self.emax = {}
+        self.efinal = {}
+        self.max_calls = {}
+
+        # Set the variables for the optimization
+        new_net_load, new_unoptimized_demand = self.initialize_net_load(net_load, real_number_of_vehicle, project, unoptimized_demand)
+        self.initialize_model(project, new_net_load, SOC_margin, SOC_offset, calls)
+        
+        if peak_shaving=='cost':
+            price = self.initialize_price(price, net_load)
+
+        elif peak_shaving=='hybrid':
+            peak_scalar, peak_subtractor = self.get_peak_scalar(new_net_load, new_unoptimized_demand)
+            print(peak_scalar)
+            print(peak_subtractor)
+
+        # Run the optimization
+        timer = time.time()
+        opti_model, result = self.process(self.times, self.vehicles, self.d, self.pmax,
+                                          self.pmin, self.emin, self.emax,
+                                          self.efinal, peak_shaving, penalization, peak_scalar, peak_subtractor, price, calls)
+        timer2 = time.time()
+        print('The optimization duration was ' + str((timer2 - timer) / 60) + ' minutes')
+        print('')
+
+        # Post process results
+        return self.post_process(project, net_load, opti_model, result, plot)
+        #return opti_model, result
+
+    def get_peak_scalar(self, net_load_updated, unoptimized_demand):
+        """
+        Args: 
+            net_load_updated: series of net load without vehicle loads
+            unoptimized_demand: series of unoptimized vehicle loads, scaled appropriately
+        Returns:
+            peak_scalar (multiplier for peak-minimization argument in objective fxn)
+            peak_subtractor (to be subtracted from peak-min arg in obj fxn, to center around 0)
+        
+        """
+        #Combining net load and unoptimized demand
+        new_load = pandas.DataFrame()
+        new_load['load'] = numpy.array(net_load_updated['netload']) + numpy.array(unoptimized_demand)
+        
+        # Normalizing peak load (subtracting mean)
+        peak_subtractor = numpy.mean(new_load['load'])
+        norm_load = new_load['load'] - peak_subtractor
+
+        # Scaling factor for peak shaving: sum of squared values at each timestep
+        peak_scale = sum(x**2 for x in norm_load)
+
+        # Scaling factor for ramp mitigation: sum of squared ramp rate between all timesteps
+        ramp_rate = new_load.subtract(new_load.shift(-1)).dropna()
+        ramp_scale = sum(x**2 for x in ramp_rate['load'])
+
+        # Getting comparative scaling factor for peak-shaving side of hybrid optimization objective
+        peak_scalar = ramp_scale / peak_scale
+        
+        return peak_scalar, peak_subtractor
+    
+    def initialize_net_load(self, net_load, real_number_of_vehicle, project, unoptimized_demand):
+        """Make sure that the net load has the right size and scale the net
+        load for the optimization scale.
+
+        Args:
+            net_load (pandas.DataFrame): data frame with date index and a 'net_load' column in [W]
+            net_load_pmax (int): maximum power on the scaled net load
+        """
+        # Make sure we are not touching the initial data
+        new_net_load = net_load.copy()
+
+        # Resample the net load & unoptimized demand
+        new_net_load = new_net_load.resample(str(self.optimization_timestep) + 'T').first()
+        new_unoptimized_demand = unoptimized_demand.resample(str(self.optimization_timestep) + 'T').first()
+
+        # Check that new net load and new unoptimized demand are same length
+        assert len(new_net_load) == len(new_unoptimized_demand)
+
+        # Check against the actual lenght it should have
+        diff = (len(new_net_load) -
+                int((self.date_to - self.date_from).total_seconds() / (60 * self.optimization_timestep)))
+        if diff > 0:
+            # We should trim the net load with diff elements (normaly 1 because of slicing inclusion)
+            new_net_load.drop(new_net_load.tail(diff).index, inplace=True)
+        elif diff < 0:
+            print('The net load does not contain enough data points')
+
+        if real_number_of_vehicle:
+            # Set scaling factor
+            scaling_factor = len(project.vehicles) / real_number_of_vehicle
+
+            # Scale the temp net load
+            new_net_load['netload'] *= scaling_factor
+
+        return new_net_load, new_unoptimized_demand
+
+    def initialize_price(self, price, net_load): ##>>>>
+        """Resample power price data to be on same frequency as net load data, and returns dictionary with timesteps instead of timestamps.
+        Args:
+            price (pandas.DataFrame): data frame with date index and a 'price' column (in $/MWh)
+            net_load (pandas.DataFrame): data frame with date index and a 'net_load' column in [W]
+        
+        Returns dictionary with 
+            """
+        # Check that price and net_load dataframes have same length
+        assert len(price) == len(net_load)
+        # Make sure we are not touching the initial data
+        new_price = price.copy()
+        # Resample the net load
+        new_price = new_price.resample(str(self.optimization_timestep) + 'T').first()
+
+
+        # Convert to timestep index instead of timestamps
+        temp_index = pandas.DataFrame(range(0, len(new_price)), columns=['index'])
+        # Set temp_index
+        temp_new_price = new_price.copy()
+        temp_new_price= temp_new_price.set_index(temp_index['index'])
+        # Return a dictionary
+        return temp_new_price.to_dict()['price']
+
+    def check_energy_constraints_feasible(self, vehicle, SOC_init, SOC_final, SOC_offset, verbose=False):
+        """Make sure that SOC final can be reached from SOC init under uncontrolled
+        charging (best case scenario). Print details when conditions are not met.
+
+        Args:
+            vehicle (Vehicle): vehicle
+            SOC_init (float): state of charge at the begining of the optimization [0, 1]
+            SOC_final (float): state of charge at the end of the optimization [0, 1]
+            SOC_offset (float): energy offset [0, 1]
+
+        Return:
+            (Boolean)
+        """
+        def print_status(SOC_final, SOC_init, diff):
+                print('Set battery gain: ' + str((SOC_final - SOC_init) * 100) + '%')
+                print('Simulation battery gain: ' + str(diff * 100))
+
+        # Check if below minimum SOC at any time
+        if (min(vehicle.SOC[self.SOC_index_from:self.SOC_index_to]) - SOC_offset) <= self.minimum_SOC:
+            if verbose:
+                print('Vehicle: ' + str(vehicle.id) + ' has a minimum SOC of ' +
+                      str(min(vehicle.SOC[self.SOC_index_from:self.SOC_index_to]) * 100) + '%')
+            return False
+
+        # Check SOC difference between date_from and date_to ?
+        # Diff represent the minimum loss or the maximum gain
+        diff = vehicle.SOC[self.SOC_index_to] - vehicle.SOC[self.SOC_index_from]
+        # Simulation show a battery gain
+        if diff > 0:
+            # Gain should be greater than the one we set up
+            if SOC_final - SOC_init < diff:
+                # Good to go
+                return True
+            else:
+                # Set final SOC to be under diff
+                print_status(SOC_final, SOC_init, diff)
+                return False
+        # Simulation show battery loss
+        else:
+            # energy balance should be negative
+            if SOC_final - SOC_init > 0:
+                print_status(SOC_final, SOC_init, diff)
+                return False
+            # Loss should be smaller than the one we set (less negative)
+            if diff < SOC_init - SOC_final:
+                return True
+            else:
+                # Set final SOC to include at least the lost
+                print_status(SOC_final, SOC_init, diff)
+                return False
+
+    def initialize_time_index(self, net_load):
+        """Replace date index by time ids
+
+        Args:
+            net_load (pandas.DataFrame): data frame with date index and a 'net_load' column in [W]
+
+        Return:
+            net_load as a dictionary with time ids, time ids
+        """
+        temp_index = pandas.DataFrame(range(0, len(net_load)), columns=['index'])
+        # Set temp_index
+        temp_net_load = net_load.copy()
+        temp_net_load = temp_net_load.set_index(temp_index['index'])
+        # Return a dictionary
+        return temp_net_load.to_dict()['netload'], temp_index.index.values.tolist()
+
+    def get_initial_SOC(self, vehicle, SOC_offset, SOC_init=None):
+        """Get the initial SOC with which people start the optimization
+        """
+        if SOC_init is not None:
+            return SOC_init
+        else:
+            return vehicle.SOC[self.SOC_index_from] - SOC_offset
+
+    def get_final_SOC(self, vehicle, SOC_margin, SOC_offset, SOC_end=None):
+        """Get final SOC that vehicle must reached at the end of the optimization
+        """
+        if SOC_end is not None:
+            return SOC_end
+        else:
+            return vehicle.SOC[self.SOC_index_to] - SOC_offset - SOC_margin
+
+    def initialize_model(self, project, net_load, SOC_margin, SOC_offset, calls):
+        """Select the vehicles that were plugged at controlled chargers and create
+        the optimization variables (see inputs of optimization)
+
+        Args:
+            project (Project): project
+
+        Return:
+            times, vehicles, d, pmax, pmin, emin, emax, efinal
+        """
+        # Create a dict with the net load and get time index in a data frame
+        self.d, self.times = self.initialize_time_index(net_load)
+        vehicle_to_optimize = 0
+        unfeasible_vehicle = 0
+
+        for vehicle in project.vehicles:
+            if vehicle.result is not None:
+                # Get SOC init and SOC end
+                SOC_init = self.get_initial_SOC(vehicle, SOC_offset)
+                SOC_final = self.get_final_SOC(vehicle, SOC_margin, SOC_offset)
+
+                # Find out if vehicle itinerary is feasible
+                if not self.check_energy_constraints_feasible(vehicle, SOC_init, SOC_final, SOC_offset):
+                    # Reset vehicle result to None
+                    vehicle.result = None
+                    unfeasible_vehicle += 1
+                    continue
+
+                # Add vehicle id to a list
+                self.vehicles.append(vehicle.id)
+                vehicle_to_optimize += 1
+
+                # Resample vehicle result
+                temp_vehicle_result = vehicle.result.resample(str(self.optimization_timestep) + 'T').first()
+
+                # Set time_vehicle_index
+                temp_vehicle_result = temp_vehicle_result.set_index(pandas.DataFrame(
+                    index=[(time, vehicle.id) for time in self.times]).index)
+
+                # Push pmax and pmin with vehicle and time key
+                self.pmin.update(temp_vehicle_result.to_dict()['p_min'])
+                self.pmax.update(temp_vehicle_result.to_dict()['p_max'])
+
+                # Push emin and emax with vehicle and time key
+                # Units! if project.timestep in seconds, self.timestep in minutes and battery in Wh
+                # Units! Wproject.timestep --> Wself.timestep * (project.timestep / (60 * self.timestep))
+                # Units! Wtimestep --> Wh * (60 / self.timestep)
+                temp_vehicle_result['emin'] = (temp_vehicle_result.energy * (project.timestep / (60 * self.optimization_timestep)) +
+                                               (self.minimum_SOC - SOC_init) * vehicle.car_model.battery_capacity *
+                                               (60 / self.optimization_timestep))
+                self.emin.update(temp_vehicle_result.to_dict()['emin'])
+                temp_vehicle_result['emax'] = (temp_vehicle_result.energy * (project.timestep / (60 * self.optimization_timestep)) + 10000 +
+                                               (self.maximum_SOC - SOC_init) * vehicle.car_model.battery_capacity *
+                                               (60 / self.optimization_timestep))
+                self.emax.update(temp_vehicle_result.to_dict()['emax'])
+
+                # Push efinal with vehicle key
+                self.efinal.update({vehicle.id: (temp_vehicle_result.tail(1).energy.values[0] * (project.timestep / (60 * self.optimization_timestep)) +
+                                                 (SOC_final - SOC_init) * vehicle.car_model.battery_capacity *
+                                                 (60 / self.optimization_timestep))})
+
+                # Create max_calls for each vehicle for emergency optimization ######################
+                self.max_calls.update({vehicle.id: calls})
+
+        print('There is ' + str(vehicle_to_optimize) + ' vehicle participating in the optimization (' +
+              str(vehicle_to_optimize * 100 / len(project.vehicles)) + '%)')
+        print('There is ' + str(unfeasible_vehicle) + ' unfeasible vehicle.')
+        print('')
+
+    def process(self, times, vehicles, d, pmax, pmin, emin, emax,
+                efinal, peak_shaving, penalization, peak_scalar, peak_subtractor, price, max_calls, solver="gurobi"):
+
+        """The process function creates the pyomo model and solve it.
+        Minimize sum( net_load(t) + sum(power_demand(t, v)))**2
+        subject to:
+        pmin(t, v) <= power_demand(t, v) <= pmax(t, v)
+        emin(t, v) <= sum(power_demand(t, v)) <= emax(t, v)
+        sum(power_demand(t, v)) >= efinal(v)
+        rampmin(t) <= net_load_ramp(t) + power_demand_ramp(t, v) <= rampmax(t)
+
+        Args:
+            times (list): timestep list
+            vehicles (list): unique list of vehicle ids
+            d (dict): time - net load at t
+            price (dict): time - power price at t ###
+            pmax (dict): (time, id) - power maximum at t for v
+            pmin (dict): (time, id) - power minimum at t for v
+            emin (dict): (time, id) - energy minimum at t for v
+            emax (dict): (time, id) - energy maximum at t for v
+            efinal (dict): id - final SOC
+            solver (string): name of the solver to use (default is gurobi)
+
+        Return:
+            model (ConcreteModel), result
+            """
+        
+
+        # Select gurobi solver
+         
+        with SolverFactory(solver) as opt:
+            # Solver option see Gurobi website
+            # opt.options['Method'] = 1
+
+            # Creation of a Concrete Model
+            model = ConcreteModel()
+
+            # ###### Set
+            model.t = Set(initialize=times, doc='Time', ordered=True)
+            last_t = model.t.last()
+            model.v = Set(initialize=vehicles, doc='Vehicles')
+
+
+            # ###### Parameters
+            # Net load
+            model.d = Param(model.t, initialize=d, doc='Net load')
+
+            # Power
+            model.p_max = Param(model.t, model.v, initialize=pmax, doc='P max')
+            model.p_min = Param(model.t, model.v, initialize=pmin, doc='P min')
+
+            # Energy
+            model.e_min = Param(model.t, model.v, initialize=emin, doc='E min')
+            model.e_max = Param(model.t, model.v, initialize=emax, doc='E max')
+
+            model.e_final = Param(model.v, initialize=efinal, doc='final energy balance')
+
+            #HYBRID
+            if peak_shaving == 'hybrid':
+                # Lambda for hybrid model--1 is all weight on peak shaving
+                model.lamb = Param(initialize = .5)
+                # Scaling factor for peak-shaving sub-objective for hybrid model
+                model.peak_scalar = Param(initialize = peak_scalar)
+                # Normalization factor for peak-shaving sub-objective
+                model.peak_subtractor = Param(initialize = peak_subtractor)
+
+            #COST
+            if peak_shaving == 'cost':
+                # Price
+                model.price = Param(model.t, initialize=price, doc='Prices') #energy price at each timestep
+
+            #EMERGENCY
+            if peak_shaving == 'emergency':
+                # Maximum Wh a vehicle can export over course of optimization (battery capacity * given fraction)
+                model.max_calls = Param(model.v, initialize = max_calls)
+
+
+            # ###### Variable
+            model.u_neg = Var(model.t, model.v, domain=Integers, doc='Power exported')
+            model.u_pos = Var(model.t, model.v, domain=Integers, doc='Power imported') 
+            model.u = Var(model.t, model.v, domain=Integers, doc='Power used')
+
+
+            # ###### Rules
+
+            # Power max/min
+            def max_charging_rule(model, t, v):
+                return model.u_pos[t, v] <= model.p_max[t, v]
+            model.max_charging_rule = Constraint(model.t, model.v, rule=max_charging_rule, doc='u_pos <= p_max')
+
+            def min_charging_rule(model, t, v):
+                return model.u_pos[t, v] >= 0
+            model.min_charging_rule = Constraint(model.t, model.v, rule=min_charging_rule, doc='u_pos >= 0')
+
+            def max_discharging_rule(model, t, v):
+                return model.u_neg[t, v] >= model.p_min[t, v]
+            model.max_discharging_rule = Constraint(model.t, model.v, rule=max_discharging_rule, doc='u_neg >= p_min')
+
+            def min_discharging_rule(model, t, v):
+                return model.u_neg[t, v] <= 0
+            model.min_discharging_rule = Constraint(model.t, model.v, rule=min_discharging_rule, doc='u_neg <= 0')
+
+            def total_power_rule(model, t, v):
+                return model.u[t, v] == model.u_neg[t, v] + model.u_pos[t, v]
+            model.total_power_rule = Constraint(model.t, model.v, rule=total_power_rule, doc='u = u_neg + u_pos')
+
+            # Energy max/min
+            def minimum_energy_rule(model, t, v):
+                return sum(model.u[i, v] for i in range(0, t + 1)) >= model.e_min[t, v]
+            model.minimum_energy_rule = Constraint(model.t, model.v, rule=minimum_energy_rule, doc='E min rule')
+
+            def maximum_energy_rule(model, t, v):
+                return sum(model.u[i, v] for i in range(0, t + 1)) <= model.e_max[t, v]
+            model.maximum_energy_rule = Constraint(model.t, model.v, rule=maximum_energy_rule, doc='E max rule')
+
+            def final_energy_balance(model, v):
+                return sum(model.u[i, v] for i in model.t) >= model.e_final[v]
+            model.final_energy_rule = Constraint(model.v, rule=final_energy_balance, doc='E final rule')
+            
+            
+
+            # Set the objective to be either peak shaving or ramp mitigation
+            if peak_shaving == 'peak_shaving':
+                def objective_rule(model):
+                    return sum([(model.d[t] + sum([model.u[t, v] for v in model.v]))**2 for t in model.t])
+                model.objective = Objective(rule=objective_rule, sense=minimize, doc='Define objective function')
+
+            elif peak_shaving == 'penalized_peak_shaving':
+                def objective_rule(model):
+                    return (sum([(model.d[t] + sum([model.u[t, v] for v in model.v]))**2 for t in model.t]) +
+                            penalization * sum([sum([model.u[t, v]**2 for v in model.v]) for t in model.t]))
+                model.objective = Objective(rule=objective_rule, sense=minimize, doc='Define objective function')
+
+            elif peak_shaving == 'ramp_mitigation':
+                def objective_rule(model):
+                    return sum([(model.d[t + 1] - model.d[t] + sum([model.u[t + 1, v] - model.u[t, v] for v in model.v]))**2 for t in model.t if t != last_t])
+                model.objective = Objective(rule=objective_rule, sense=minimize, doc='Define objective function')
+
+            ### MM additions
+            # Hybrid
+            # Definition: optimizes for both ramp mitigation and peak shaving
+            elif peak_shaving == 'hybrid': 
+                def objective_rule(model):
+                    return model.lamb * model.peak_scalar * (sum([(model.d[t] - model.peak_subtractor + sum([model.u[t, v] for v in model.v]))**2 for t in model.t]) ) + \
+                        (1-model.lamb) * sum([(model.d[t + 1] - model.d[t] + sum([model.u[t + 1, v] - model.u[t, v] for v in model.v]))**2 for t in model.t if t != last_t])
+                model.objective = Objective(rule=objective_rule, sense=minimize, doc='Define objective function')
+
+            # Emergency
+            # Definition: Emergency V2G; mostly V1G but each vehicle can get called for V2G a certain number of times/year
+            # Can be used either for peak shaving or ramp mitigation; for now just implementing peak shaving
+            elif peak_shaving == 'emergency':
+                # Peak shaving objective
+                def objective_rule(model):
+                    return sum([(model.d[t] + sum([model.u[t, v] for v in model.v]))**2 for t in model.t])
+                model.objective = Objective(rule=objective_rule, sense=minimize, doc='Define objective function')
+
+                # V2G (export) limits
+                def max_export_rule(model, v): 
+                    return sum(model.u_neg[i, v] for i in model.t) >= model.max_calls[v]
+                model.max_export_rule = Constraint(model.v, rule=max_export_rule, doc='Max export rule')
+
+
+            # Cost-based
+            # Definition: minimizes charging cost based on given cost stack
+            # Should probably only have visibility 24/48 hours into the future, ideally
+            # V1G vs V2G will make a big difference here--incentivizing selling on-peak
+            elif peak_shaving == 'cost':
+                def objective_rule(model):
+                    return sum([sum([model.u[t, v] * model.price[t] for t in model.t]) for v in model.v])
+                    #return sum((model.u[t, v] * model.price[t] for v in model.v) for t in model.t)
+                model.objective = Objective(rule=objective_rule, sense=minimize, doc='Define objective function')
+            
+            solver_parameters = "ResultFile=model.ilp"
+            results = opt.solve(model, options_string=solver_parameters, symbolic_solver_labels=True)
+            # results.write()
+
+        return model, results
+
+    def plot_result(self, model):
+        """Create a plot showing the power constraints, the energy constraints and the ramp
+        constraints as well as the final net load.
+        """
+        # Set the graph style
+        sns.set_style("whitegrid")
+        sns.despine()
+
+        result = pandas.DataFrame()
+
+        # Get the result
+        df = pandas.DataFrame(index=['power'], data=model.u.get_values()).transpose().groupby(level=0).sum()
+        # Ramp of the result
+        mylist = [0]
+        mylist.extend(list(numpy.diff(df['power'].values)))
+        df['ramppower'] = mylist
+        result = pandas.concat([result, df], axis=1)
+        # cum sum of the result
+        df = pandas.DataFrame(
+            pandas.DataFrame(
+                index=['anything'],
+                data=model.u.get_values()).transpose().unstack().cumsum().sum(axis=1), columns=['powercum']) * (5 / 60)
+        result = pandas.concat([result, df], axis=1)
+
+        # Get pmax and pmin units of power [W]
+        df = pandas.DataFrame(index=['pmax'], data=model.p_max.extract_values()).transpose().groupby(level=0).sum()
+        result = pandas.concat([result, df], axis=1)
+        df = pandas.DataFrame(index=['pmin'], data=model.p_min.extract_values()).transpose().groupby(level=0).sum()
+        result = pandas.concat([result, df], axis=1)
+
+        # Get emin and emax units of energy [Wh]
+        df = pandas.DataFrame(index=['emax'], data=model.e_max.extract_values()).transpose().groupby(level=0).sum() * (5 / 60)
+        result = pandas.concat([result, df], axis=1)
+        df = pandas.DataFrame(index=['emin'], data=model.e_min.extract_values()).transpose().groupby(level=0).sum() * (5 / 60)
+        result = pandas.concat([result, df], axis=1)
+
+        # Get the minimum final energy quantity
+        e_final = pandas.DataFrame(index=['efinal'], data=model.e_final.extract_values()).transpose().sum() * (5 / 60)
+
+        # Get the actual ramp
+        mylist = [0]
+        df = pandas.DataFrame(index=['net_load'], data=model.d.extract_values()).transpose()
+        mylist.extend(list(numpy.diff(df['net_load'].values)))
+        df['rampnet_load'] = mylist
+        result = pandas.concat([result, df], axis=1)
+
+        # Plot power constraints
+        plt.subplot(411)
+        plt.plot(result.index.values, result.pmax.values, label='pmax')
+        plt.plot(result.index.values, result.power.values, label='power')
+        plt.plot(result.index.values, result.pmin.values, label='pmin')
+        plt.legend(loc=0)
+
+        # Plot energy constraints
+        plt.subplot(412)
+        plt.plot(result.index.values, result.emax.values, label='emax')
+        plt.plot(result.index.values, result.powercum.values, label='powercum')
+        plt.plot(result.index.values, result.emin.values, label='emin')
+        plt.plot(result.index[-1], e_final, '*', markersize=15, label='efinal')
+        plt.legend(loc=0)
+
+        # Plot the result ramp
+        plt.subplot(413)
+        plt.plot(result.index.values, result.ramppower.values + result.rampnet_load.values, label='result ramp')
+        plt.plot(result.index.values, result.rampnet_load.values, label='net_load ramp')
+        plt.legend(loc=0)
+
+        # Plot the power demand results
+        plt.subplot(414)
+        plt.plot(result.index.values, result.net_load.values, label='net_load')
+        plt.plot(result.index.values, result.net_load.values + result.power.values, label='net_load + vehicle')
+        plt.legend(loc=0)
+        plt.show()
+
+    def post_process(self, project, netload, model, result, plot):
+        """Recompute SOC profiles and compute new total power demand
+
+        Args:
+            project (Project): project
+
+        Note: Should check that 'vehicle before' and 'after' contain the same number of vehicles
+        """
+        if plot:
+            self.plot_result(model)
+
+        temp = pandas.DataFrame()
+        first = True
+        for vehicle in project.vehicles:
+            if vehicle.result is not None:
+                if first:
+                    temp['vehicle_before'] = vehicle.result['power_demand']
+                    first = False
+                else:
+                    temp['vehicle_before'] += vehicle.result['power_demand']
+
+        temp2 = pandas.DataFrame(index=['vehicle_after'], data=model.u.get_values()).transpose().groupby(level=0).sum()
+        i = pandas.date_range(start=self.date_from, end=self.date_to,
+                              freq=str(self.optimization_timestep) + 'T', closed='left')
+        temp2 = temp2.set_index(i)
+        temp2 = temp2.resample(str(project.timestep) + 'S')
+        temp2 = temp2.fillna(method='ffill').fillna(method='bfill')
+
+        final_result = pandas.DataFrame()
+        final_result = pandas.concat([temp['vehicle_before'], temp2['vehicle_after']], axis=1)
+        final_result = final_result.fillna(method='ffill').fillna(method='bfill')
+
+        return final_result
+
+
+def save_vehicle_state_for_optimization(vehicle, timestep, date_from,
+                                        date_to, activity=None, power_demand=None,
+                                        SOC=None, detail=None, nb_interval=None, init=False,
+                                        run=False, post=False):
+    """Save results for individual vehicles. Power demand is positive when charging
+    negative when driving. Energy consumption is positive when driving and negative
+    when charging. Charging station that offer after simulation processing should
+    have activity.charging_station.post_simulation True.
+    """
+    if run:
+        if vehicle.result is not None:
+            activity_index1, activity_index2, location_index1, location_index2, save = v2gsim.result._map_index(
+                activity.start, activity.end, date_from, date_to, len(power_demand),
+                len(vehicle.result['power_demand']), timestep)
+            # Time frame are matching
+            if save:
+                # If driving pmin and pmax are equal to 0 since we are not plugged
+                if isinstance(activity, v2gsim.model.Driving):
+                    vehicle.result['p_max'][location_index1:location_index2] -= (
+                        [0.0] * (activity_index2 - activity_index1))
+                    vehicle.result['p_min'][location_index1:location_index2] -= (
+                        [0.0] * (activity_index2 - activity_index1))
+                    # Energy consumed is directly the power demand (sum later)
+                    vehicle.result['energy'][location_index1:location_index2] += (
+                        power_demand[activity_index1:activity_index2])
+                    # Power demand on the grid is 0 since we are driving
+                    vehicle.result['power_demand'][location_index1:location_index2] -= (
+                        [0.0] * (activity_index2 - activity_index1))
+
+                # If parked pmin and pmax are not necessary the same
+                if isinstance(activity, v2gsim.model.Parked):
+                    # Save the positive power demand of this specific vehicle
+                    vehicle.result['power_demand'][location_index1:location_index2] += (
+                        power_demand[activity_index1:activity_index2])
+                    if activity.charging_station.post_simulation:
+                        # Find if vehicle or infra is limiting
+                        pmax = min(activity.charging_station.maximum_power,
+                                   vehicle.car_model.maximum_power)
+                        pmin = max(activity.charging_station.minimum_power,
+                                   vehicle.car_model.minimum_power)
+                        vehicle.result['p_max'][location_index1:location_index2] += (
+                            [pmax] * (activity_index2 - activity_index1))
+                        vehicle.result['p_min'][location_index1:location_index2] += (
+                            [pmin] * (activity_index2 - activity_index1))
+                        # Energy consumed is 0 the optimization will decide
+                        vehicle.result['energy'][location_index1:location_index2] -= (
+                            [0.0] * (activity_index2 - activity_index1))
+                    else:
+                        vehicle.result['p_max'][location_index1:location_index2] += (
+                            power_demand[activity_index1:activity_index2])
+                        vehicle.result['p_min'][location_index1:location_index2] += (
+                            power_demand[activity_index1:activity_index2])
+                        # Energy is 0.0 because it's already accounted in power_demand
+                        vehicle.result['energy'][location_index1:location_index2] -= (
+                            [0.0] * (activity_index2 - activity_index1))
+
+    elif init:
+        vehicle.SOC = [vehicle.SOC[0]]
+        vehicle.result = None
+        for activity in vehicle.activities:
+            if isinstance(activity, v2gsim.model.Parked):
+                if activity.charging_station.post_simulation:
+                    # Initiate a dictionary of numpy array to hold result (faster than DataFrame)
+                    vehicle.result = {'power_demand': numpy.array([0.0] * int((date_to - date_from).total_seconds() / timestep)),
+                                      'p_max': numpy.array([0.0] * int((date_to - date_from).total_seconds() / timestep)),
+                                      'p_min': numpy.array([0.0] * int((date_to - date_from).total_seconds() / timestep)),
+                                      'energy': numpy.array([0.0] * int((date_to - date_from).total_seconds() / timestep))}
+                    # Leave the init function
+                    return
+    elif post:
+        if vehicle.result is not None:
+            # Convert location result back into pandas DataFrame (faster that way)
+            i = pandas.date_range(start=date_from, end=date_to,
+                                  freq=str(timestep) + 's', closed='left')
+            vehicle.result = pandas.DataFrame(index=i, data=vehicle.result)
+            vehicle.result['energy'] = vehicle.result['energy'].cumsum()

From e45d5275ce92fef45afd66ea6d6ed41f751274cf Mon Sep 17 00:00:00 2001
From: Margaret McCall <majemc@gmail.com>
Date: Mon, 22 Jun 2020 10:07:59 -0700
Subject: [PATCH 11/14] commmitting before porting emergency mods into main
 file

---
 v2gsim/__init__.py                            |  5 +-
 .../post_simulation/netload_opti_emergency.py | 20 +++----
 v2gsim/post_simulation/netload_opti_hybrid.py |  2 +-
 v2gsim/post_simulation/netload_opti_plus.py   | 55 ++++++++++++++++---
 .../post_simulation/netload_optimization.py   | 28 +++++++++-
 5 files changed, 88 insertions(+), 22 deletions(-)

diff --git a/v2gsim/__init__.py b/v2gsim/__init__.py
index f73c81e..259590d 100644
--- a/v2gsim/__init__.py
+++ b/v2gsim/__init__.py
@@ -10,6 +10,9 @@
 import v2gsim.model
 import v2gsim.tool
 import v2gsim.post_simulation.netload_optimization
+import v2gsim.post_simulation.netload_opti_plus
+import v2gsim.post_simulation.netload_opti_hybrid
+import v2gsim.post_simulation.netload_opti_emergency
 import v2gsim.post_simulation.result
 import v2gsim.charging.uncontrolled
 import v2gsim.charging.controlled
@@ -25,4 +28,4 @@
            'v2gsim.charging.uncontrolled', 'v2gsim.charging.controlled', 'v2gsim.charging.station',
            'v2gsim.driving.basic_powertrain', 'v2gsim.driving.drivecycle.generator',
            'v2gsim.driving.detailed.power_train', 'v2gsim.driving.detailed.init_model',
-           'v2gsim.post_simulation.result', 'v2gsim.battery_degradation.BatteryDegradation']
+           'v2gsim.post_simulation.result', 'v2gsim.battery_degradation.BatteryDegradation']
\ No newline at end of file
diff --git a/v2gsim/post_simulation/netload_opti_emergency.py b/v2gsim/post_simulation/netload_opti_emergency.py
index e5cb03e..9890d66 100644
--- a/v2gsim/post_simulation/netload_opti_emergency.py
+++ b/v2gsim/post_simulation/netload_opti_emergency.py
@@ -33,7 +33,7 @@ def __init__(self, project, optimization_timestep, date_from,
 
     def solve(self, project, net_load, real_number_of_vehicle, SOC_margin=0.02,
               SOC_offset=0.0, peak_shaving='peak_shaving', penalization=5, beta=None, plot=False, peak_scalar=.01, peak_subtractor=30000, 
-              unoptimized_demand = pandas.DataFrame(), price=pandas.DataFrame(), calls=-3000): #price can't be array, is dict; can't have non-default arg follow default
+              unoptimized_demand = pandas.DataFrame(), price=pandas.DataFrame(), max_export=-3000): #price can't be array, is dict; can't have non-default arg follow default
         """Launch the optimization and the post_processing fucntion. Results
         and assumptions are appended to a data frame.
 
@@ -55,11 +55,11 @@ def solve(self, project, net_load, real_number_of_vehicle, SOC_margin=0.02,
         self.emin = {}
         self.emax = {}
         self.efinal = {}
-        self.max_calls = {}
+        self.export_limit = {}
 
         # Set the variables for the optimization
         new_net_load, new_unoptimized_demand = self.initialize_net_load(net_load, real_number_of_vehicle, project, unoptimized_demand)
-        self.initialize_model(project, new_net_load, SOC_margin, SOC_offset, calls)
+        self.initialize_model(project, new_net_load, SOC_margin, SOC_offset, max_export)
         
         if peak_shaving=='cost':
             price = self.initialize_price(price, net_load)
@@ -73,7 +73,7 @@ def solve(self, project, net_load, real_number_of_vehicle, SOC_margin=0.02,
         timer = time.time()
         opti_model, result = self.process(self.times, self.vehicles, self.d, self.pmax,
                                           self.pmin, self.emin, self.emax,
-                                          self.efinal, peak_shaving, penalization, peak_scalar, peak_subtractor, price, calls)
+                                          self.efinal, peak_shaving, penalization, peak_scalar, peak_subtractor, price, max_export)
         timer2 = time.time()
         print('The optimization duration was ' + str((timer2 - timer) / 60) + ' minutes')
         print('')
@@ -255,7 +255,7 @@ def get_final_SOC(self, vehicle, SOC_margin, SOC_offset, SOC_end=None):
         else:
             return vehicle.SOC[self.SOC_index_to] - SOC_offset - SOC_margin
 
-    def initialize_model(self, project, net_load, SOC_margin, SOC_offset, calls):
+    def initialize_model(self, project, net_load, SOC_margin, SOC_offset, max_export):
         """Select the vehicles that were plugged at controlled chargers and create
         the optimization variables (see inputs of optimization)
 
@@ -316,8 +316,8 @@ def initialize_model(self, project, net_load, SOC_margin, SOC_offset, calls):
                                                  (SOC_final - SOC_init) * vehicle.car_model.battery_capacity *
                                                  (60 / self.optimization_timestep))})
 
-                # Create max_calls for each vehicle for emergency optimization ######################
-                self.max_calls.update({vehicle.id: calls})
+                # Create export_limit for each vehicle for emergency optimization ######################
+                self.export_limit.update({vehicle.id: max_export})
 
         print('There is ' + str(vehicle_to_optimize) + ' vehicle participating in the optimization (' +
               str(vehicle_to_optimize * 100 / len(project.vehicles)) + '%)')
@@ -325,7 +325,7 @@ def initialize_model(self, project, net_load, SOC_margin, SOC_offset, calls):
         print('')
 
     def process(self, times, vehicles, d, pmax, pmin, emin, emax,
-                efinal, peak_shaving, penalization, peak_scalar, peak_subtractor, price, max_calls, solver="gurobi"):
+                efinal, peak_shaving, penalization, peak_scalar, peak_subtractor, price, export_limit, solver="gurobi"):
 
         """The process function creates the pyomo model and solve it.
         Minimize sum( net_load(t) + sum(power_demand(t, v)))**2
@@ -398,7 +398,7 @@ def process(self, times, vehicles, d, pmax, pmin, emin, emax,
             #EMERGENCY
             if peak_shaving == 'emergency':
                 # Maximum Wh a vehicle can export over course of optimization (battery capacity * given fraction)
-                model.max_calls = Param(model.v, initialize = max_calls)
+                model.export_limit = Param(model.v, initialize = export_limit)
 
 
             # ###### Variable
@@ -482,7 +482,7 @@ def objective_rule(model):
 
                 # V2G (export) limits
                 def max_export_rule(model, v): 
-                    return sum(model.u_neg[i, v] for i in model.t) >= model.max_calls[v]
+                    return sum(model.u_neg[i, v] for i in model.t) >= model.export_limit[v]
                 model.max_export_rule = Constraint(model.v, rule=max_export_rule, doc='Max export rule')
 
 
diff --git a/v2gsim/post_simulation/netload_opti_hybrid.py b/v2gsim/post_simulation/netload_opti_hybrid.py
index 185625f..1d342b4 100644
--- a/v2gsim/post_simulation/netload_opti_hybrid.py
+++ b/v2gsim/post_simulation/netload_opti_hybrid.py
@@ -383,7 +383,7 @@ def process(self, times, vehicles, d, pmax, pmin, emin, emax,
             #HYBRID
             if peak_shaving == 'hybrid':
                 # Lambda for hybrid model--1 is all weight on peak shaving
-                model.lamb = Param(initialize = .9)
+                model.lamb = Param(initialize = .5)
                 # Scaling factor for peak-shaving sub-objective for hybrid model
                 model.peak_scalar = Param(initialize = peak_scalar)
                 # Normalization factor for peak-shaving sub-objective
diff --git a/v2gsim/post_simulation/netload_opti_plus.py b/v2gsim/post_simulation/netload_opti_plus.py
index a489a3e..1d342b4 100644
--- a/v2gsim/post_simulation/netload_opti_plus.py
+++ b/v2gsim/post_simulation/netload_opti_plus.py
@@ -31,8 +31,8 @@ def __init__(self, project, optimization_timestep, date_from,
         self.SOC_index_from = int((date_from - project.date).total_seconds() / project.timestep)
         self.SOC_index_to = int((date_to - project.date).total_seconds() / project.timestep)
 
-    def solve(self, project, net_load, real_number_of_vehicle, SOC_margin=0.02,
-              SOC_offset=0.0, peak_shaving='peak_shaving', penalization=5, beta=None, plot=False, peak_scalar = .01, peak_subtractor = 50000, price=pandas.DataFrame(), calls=2): #price can't be array, is dict; can't have non-default arg follow default
+    def solve(self, project, net_load, unoptimized_demand, real_number_of_vehicle, SOC_margin=0.02,
+              SOC_offset=0.0, peak_shaving='peak_shaving', penalization=5, beta=None, plot=False, peak_scalar=.01, peak_subtractor=30000, price=pandas.DataFrame(), calls=2): #price can't be array, is dict; can't have non-default arg follow default
         """Launch the optimization and the post_processing fucntion. Results
         and assumptions are appended to a data frame.
 
@@ -57,12 +57,17 @@ def solve(self, project, net_load, real_number_of_vehicle, SOC_margin=0.02,
         self.max_calls = {}
 
         # Set the variables for the optimization
-        new_net_load = self.initialize_net_load(net_load, real_number_of_vehicle, project)
+        new_net_load, new_unoptimized_demand = self.initialize_net_load(net_load, real_number_of_vehicle, project, unoptimized_demand)
         self.initialize_model(project, new_net_load, SOC_margin, SOC_offset, calls)
         
         if peak_shaving=='cost':
             price = self.initialize_price(price, net_load)
 
+        elif peak_shaving=='hybrid':
+            peak_scalar, peak_subtractor = self.get_peak_scalar(new_net_load, new_unoptimized_demand)
+            print(peak_scalar)
+            print(peak_subtractor)
+
         # Run the optimization
         timer = time.time()
         opti_model, result = self.process(self.times, self.vehicles, self.d, self.pmax,
@@ -76,7 +81,37 @@ def solve(self, project, net_load, real_number_of_vehicle, SOC_margin=0.02,
         return self.post_process(project, net_load, opti_model, result, plot)
         #return opti_model, result
 
-    def initialize_net_load(self, net_load, real_number_of_vehicle, project):
+    def get_peak_scalar(self, net_load_updated, unoptimized_demand):
+        """
+        Args: 
+            net_load_updated: series of net load without vehicle loads
+            unoptimized_demand: series of unoptimized vehicle loads, scaled appropriately
+        Returns:
+            peak_scalar (multiplier for peak-minimization argument in objective fxn)
+            peak_subtractor (to be subtracted from peak-min arg in obj fxn, to center around 0)
+        
+        """
+        #Combining net load and unoptimized demand
+        new_load = pandas.DataFrame()
+        new_load['load'] = numpy.array(net_load_updated['netload']) + numpy.array(unoptimized_demand)
+        
+        # Normalizing peak load (subtracting mean)
+        peak_subtractor = numpy.mean(new_load['load'])
+        norm_load = new_load['load'] - peak_subtractor
+
+        # Scaling factor for peak shaving: sum of squared values at each timestep
+        peak_scale = sum(x**2 for x in norm_load)
+
+        # Scaling factor for ramp mitigation: sum of squared ramp rate between all timesteps
+        ramp_rate = new_load.subtract(new_load.shift(-1)).dropna()
+        ramp_scale = sum(x**2 for x in ramp_rate['load'])
+
+        # Getting comparative scaling factor for peak-shaving side of hybrid optimization objective
+        peak_scalar = ramp_scale / peak_scale
+        
+        return peak_scalar, peak_subtractor
+    
+    def initialize_net_load(self, net_load, real_number_of_vehicle, project, unoptimized_demand):
         """Make sure that the net load has the right size and scale the net
         load for the optimization scale.
 
@@ -87,8 +122,12 @@ def initialize_net_load(self, net_load, real_number_of_vehicle, project):
         # Make sure we are not touching the initial data
         new_net_load = net_load.copy()
 
-        # Resample the net load
+        # Resample the net load & unoptimized demand
         new_net_load = new_net_load.resample(str(self.optimization_timestep) + 'T').first()
+        new_unoptimized_demand = unoptimized_demand.resample(str(self.optimization_timestep) + 'T').first()
+
+        # Check that new net load and new unoptimized demand are same length
+        assert len(new_net_load) == len(new_unoptimized_demand)
 
         # Check against the actual lenght it should have
         diff = (len(new_net_load) -
@@ -106,7 +145,7 @@ def initialize_net_load(self, net_load, real_number_of_vehicle, project):
             # Scale the temp net load
             new_net_load['netload'] *= scaling_factor
 
-        return new_net_load
+        return new_net_load, new_unoptimized_demand
 
     def initialize_price(self, price, net_load): ##>>>>
         """Resample power price data to be on same frequency as net load data, and returns dictionary with timesteps instead of timestamps.
@@ -343,7 +382,7 @@ def process(self, times, vehicles, d, pmax, pmin, emin, emax,
 
             #HYBRID
             if peak_shaving == 'hybrid':
-                # Lambda for hybrid model
+                # Lambda for hybrid model--1 is all weight on peak shaving
                 model.lamb = Param(initialize = .5)
                 # Scaling factor for peak-shaving sub-objective for hybrid model
                 model.peak_scalar = Param(initialize = peak_scalar)
@@ -451,7 +490,7 @@ def objective_rule(model):
             # Definition: optimizes for both ramp mitigation and peak shaving
             elif peak_shaving == 'hybrid': 
                 def objective_rule(model):
-                    return model.lamb * model.peak_scalar * (sum([(model.d[t] + sum([model.u[t, v] for v in model.v]))**2 for t in model.t]) - model.peak_subtractor) + \
+                    return model.lamb * model.peak_scalar * (sum([(model.d[t] - model.peak_subtractor + sum([model.u[t, v] for v in model.v]))**2 for t in model.t]) ) + \
                         (1-model.lamb) * sum([(model.d[t + 1] - model.d[t] + sum([model.u[t + 1, v] - model.u[t, v] for v in model.v]))**2 for t in model.t if t != last_t])
                 model.objective = Objective(rule=objective_rule, sense=minimize, doc='Define objective function')
 
diff --git a/v2gsim/post_simulation/netload_optimization.py b/v2gsim/post_simulation/netload_optimization.py
index 47e5424..804532b 100644
--- a/v2gsim/post_simulation/netload_optimization.py
+++ b/v2gsim/post_simulation/netload_optimization.py
@@ -227,6 +227,7 @@ def initialize_model(self, project, net_load, SOC_margin, SOC_offset):
                 # Push pmax and pmin with vehicle and time key
                 self.pmin.update(temp_vehicle_result.to_dict()['p_min'])
                 self.pmax.update(temp_vehicle_result.to_dict()['p_max'])
+                
 
                 # Push emin and emax with vehicle and time key
                 # Units! if project.timestep in seconds, self.timestep in minutes and battery in Wh
@@ -345,7 +346,8 @@ def objective_rule(model):
                     return sum([(model.d[t + 1] - model.d[t] + sum([model.u[t + 1, v] - model.u[t, v] for v in model.v]))**2 for t in model.t if t != last_t])
                 model.objective = Objective(rule=objective_rule, sense=minimize, doc='Define objective function')
 
-            results = opt.solve(model)
+            solver_parameters = "ResultFile=model2.ilp" ##new
+            results = opt.solve(model, options_string=solver_parameters, symbolic_solver_labels=True) ##new
             # results.write()
 
         return model, results
@@ -468,6 +470,7 @@ def save_vehicle_state_for_optimization(vehicle, timestep, date_from,
     when charging. Charging station that offer after simulation processing should
     have activity.charging_station.post_simulation True.
     """
+
     if run:
         if vehicle.result is not None:
             activity_index1, activity_index2, location_index1, location_index2, save = v2gsim.result._map_index(
@@ -510,11 +513,32 @@ def save_vehicle_state_for_optimization(vehicle, timestep, date_from,
                         vehicle.result['p_max'][location_index1:location_index2] += (
                             power_demand[activity_index1:activity_index2])
                         vehicle.result['p_min'][location_index1:location_index2] += (
-                            power_demand[activity_index1:activity_index2])
+                            power_demand[activity_index1:activity_index2]) #if net load scaled too small...careful w units
+                            #bc kW vs kWh, but charging every 10m...5 kW, 5 kWh battery...no should be same
+                            #things are running as it should; some numerical issue...
+                            #final energy balance not something you set, just same as beginning
+                            #maybe really obvious somehow that one of constraints hasn't been scaled...e_final is issue??
+                            #JC has worked with one giant car before
+                            #interesting to plot optimization boundaries...
+                            #max energy rule is potentially happening; final energy balance...more likely #1, since final e is just trying to 
+                            #be above certain threshold
+                            #did i scale chargers wrong compared to cars??
+                            #weird bc it should stop the charging when it reached maximum
+                            #check if below minimum SOC at any time...dont check for above SOC, bc charging function not supposed
+                            #to charge above...high rates might make this weird
+                            #woudl be intersting to see if we could check if we dont go above max SOC at any time
+                            #if all constraints are set properly, shouldnt matter if i'm doing one large car
+                            #look at vehicles, look at SOC, see if you find anything weird--SOC is going above 1
+                            #charging function uncontrolled--charging uncontrolled . poi -- if SOC < maximum
+                            #ah--if SOC is below maximum, you charge; else, you don't charge. Have a chance that SOC is slightly below maximum,
+                            # but as soonas you add giant value one more time, you're above maximum
+                            #maybe make cars really small??? let him know how it goes, if looking at charging fxn isnt answer
                         # Energy is 0.0 because it's already accounted in power_demand
                         vehicle.result['energy'][location_index1:location_index2] -= (
                             [0.0] * (activity_index2 - activity_index1))
 
+
+
     elif init:
         vehicle.SOC = [vehicle.SOC[0]]
         vehicle.result = None

From 6a192cbedf60ecbb648b2d5c733044e87834f3c4 Mon Sep 17 00:00:00 2001
From: Margaret McCall <majemc@gmail.com>
Date: Thu, 13 Aug 2020 12:34:23 -0500
Subject: [PATCH 12/14] working opti plus model

---
 v2gsim/post_simulation/netload_opti_plus.py | 151 +++++++++-----------
 1 file changed, 70 insertions(+), 81 deletions(-)

diff --git a/v2gsim/post_simulation/netload_opti_plus.py b/v2gsim/post_simulation/netload_opti_plus.py
index 1d342b4..01d5e53 100644
--- a/v2gsim/post_simulation/netload_opti_plus.py
+++ b/v2gsim/post_simulation/netload_opti_plus.py
@@ -31,8 +31,9 @@ def __init__(self, project, optimization_timestep, date_from,
         self.SOC_index_from = int((date_from - project.date).total_seconds() / project.timestep)
         self.SOC_index_to = int((date_to - project.date).total_seconds() / project.timestep)
 
-    def solve(self, project, net_load, unoptimized_demand, real_number_of_vehicle, SOC_margin=0.02,
-              SOC_offset=0.0, peak_shaving='peak_shaving', penalization=5, beta=None, plot=False, peak_scalar=.01, peak_subtractor=30000, price=pandas.DataFrame(), calls=2): #price can't be array, is dict; can't have non-default arg follow default
+    def solve(self, project, net_load, unoptimized_demand, real_number_of_vehicle, SOC_margin=0.02, 
+              SOC_offset=0.0, peak_shaving='peak_shaving', penalization=5, beta=None, plot=False, peak_scalar=.01, peak_subtractor=30000, 
+              price=pandas.DataFrame(), max_export=-3000): #price can't be array, is dict; can't have non-default arg follow default
         """Launch the optimization and the post_processing fucntion. Results
         and assumptions are appended to a data frame.
 
@@ -54,11 +55,12 @@ def solve(self, project, net_load, unoptimized_demand, real_number_of_vehicle, S
         self.emin = {}
         self.emax = {}
         self.efinal = {}
-        self.max_calls = {}
+        self.export_limit = {}
+
 
         # Set the variables for the optimization
         new_net_load, new_unoptimized_demand = self.initialize_net_load(net_load, real_number_of_vehicle, project, unoptimized_demand)
-        self.initialize_model(project, new_net_load, SOC_margin, SOC_offset, calls)
+        self.initialize_model(project, new_net_load, SOC_margin, SOC_offset, max_export) 
         
         if peak_shaving=='cost':
             price = self.initialize_price(price, net_load)
@@ -72,13 +74,13 @@ def solve(self, project, net_load, unoptimized_demand, real_number_of_vehicle, S
         timer = time.time()
         opti_model, result = self.process(self.times, self.vehicles, self.d, self.pmax,
                                           self.pmin, self.emin, self.emax,
-                                          self.efinal, peak_shaving, penalization, peak_scalar, peak_subtractor, price, calls)
+                                          self.efinal, peak_shaving, penalization, peak_scalar, peak_subtractor, price, max_export)
         timer2 = time.time()
         print('The optimization duration was ' + str((timer2 - timer) / 60) + ' minutes')
         print('')
 
         # Post process results
-        return self.post_process(project, net_load, opti_model, result, plot)
+        return self.post_process(project, net_load, opti_model, result, plot), opti_model ##ADDED OPTI_MODEL
         #return opti_model, result
 
     def get_peak_scalar(self, net_load_updated, unoptimized_demand):
@@ -94,6 +96,8 @@ def get_peak_scalar(self, net_load_updated, unoptimized_demand):
         #Combining net load and unoptimized demand
         new_load = pandas.DataFrame()
         new_load['load'] = numpy.array(net_load_updated['netload']) + numpy.array(unoptimized_demand)
+        plt.plot(numpy.array(net_load_updated['netload']))
+        plt.plot(numpy.array(unoptimized_demand))
         
         # Normalizing peak load (subtracting mean)
         peak_subtractor = numpy.mean(new_load['load'])
@@ -188,6 +192,11 @@ def print_status(SOC_final, SOC_init, diff):
                 print('Set battery gain: ' + str((SOC_final - SOC_init) * 100) + '%')
                 print('Simulation battery gain: ' + str(diff * 100))
 
+        print('CHECKING SOC for vehicle '+str(vehicle.id))
+        print('SOC_init: '+str(SOC_init))
+        print('SOC_final: '+str(SOC_final))
+        
+
         # Check if below minimum SOC at any time
         if (min(vehicle.SOC[self.SOC_index_from:self.SOC_index_to]) - SOC_offset) <= self.minimum_SOC:
             if verbose:
@@ -198,6 +207,7 @@ def print_status(SOC_final, SOC_init, diff):
         # Check SOC difference between date_from and date_to ?
         # Diff represent the minimum loss or the maximum gain
         diff = vehicle.SOC[self.SOC_index_to] - vehicle.SOC[self.SOC_index_from]
+        print('Diff between opti SOC diff and unopti SOC diff: '+str(SOC_final - SOC_init - diff))
         # Simulation show a battery gain
         if diff > 0:
             # Gain should be greater than the one we set up
@@ -254,7 +264,7 @@ def get_final_SOC(self, vehicle, SOC_margin, SOC_offset, SOC_end=None):
         else:
             return vehicle.SOC[self.SOC_index_to] - SOC_offset - SOC_margin
 
-    def initialize_model(self, project, net_load, SOC_margin, SOC_offset, calls):
+    def initialize_model(self, project, net_load, SOC_margin, SOC_offset, max_export):
         """Select the vehicles that were plugged at controlled chargers and create
         the optimization variables (see inputs of optimization)
 
@@ -272,6 +282,7 @@ def initialize_model(self, project, net_load, SOC_margin, SOC_offset, calls):
         for vehicle in project.vehicles:
             if vehicle.result is not None:
                 # Get SOC init and SOC end
+                print("MADE IT TO SOC ASSIGNMENT")
                 SOC_init = self.get_initial_SOC(vehicle, SOC_offset)
                 SOC_final = self.get_final_SOC(vehicle, SOC_margin, SOC_offset)
 
@@ -308,15 +319,30 @@ def initialize_model(self, project, net_load, SOC_margin, SOC_offset, calls):
                 temp_vehicle_result['emax'] = (temp_vehicle_result.energy * (project.timestep / (60 * self.optimization_timestep)) + 10000 +
                                                (self.maximum_SOC - SOC_init) * vehicle.car_model.battery_capacity *
                                                (60 / self.optimization_timestep))
+
                 self.emax.update(temp_vehicle_result.to_dict()['emax'])
 
+                
+
+
                 # Push efinal with vehicle key
                 self.efinal.update({vehicle.id: (temp_vehicle_result.tail(1).energy.values[0] * (project.timestep / (60 * self.optimization_timestep)) +
                                                  (SOC_final - SOC_init) * vehicle.car_model.battery_capacity *
                                                  (60 / self.optimization_timestep))})
 
-                # Create max_calls for each vehicle for emergency optimization ######################
-                self.max_calls.update({vehicle.id: calls})
+                print("")
+                print(vehicle.id)
+                print('last energy value:')
+                print(temp_vehicle_result.tail(1).energy.values[0])
+                print('SOC_final: '+str(SOC_final))
+                print('SOC_init: '+str(SOC_init))
+                print('HALF 1: '+str((temp_vehicle_result.tail(1).energy.values[0] * (project.timestep / (60 * self.optimization_timestep)))))
+                print('HALF 2: '+str((SOC_final - SOC_init) * vehicle.car_model.battery_capacity *(60 / self.optimization_timestep)))
+                print('EFINAL: '+str(self.efinal[vehicle.id]))
+                print("")
+
+                # Create export_limit for each vehicle for emergency optimization ######################
+                self.export_limit.update({vehicle.id: max_export})
 
         print('There is ' + str(vehicle_to_optimize) + ' vehicle participating in the optimization (' +
               str(vehicle_to_optimize * 100 / len(project.vehicles)) + '%)')
@@ -324,7 +350,7 @@ def initialize_model(self, project, net_load, SOC_margin, SOC_offset, calls):
         print('')
 
     def process(self, times, vehicles, d, pmax, pmin, emin, emax,
-                efinal, peak_shaving, penalization, peak_scalar, peak_subtractor, price, max_calls, solver="gurobi"):
+                efinal, peak_shaving, penalization, peak_scalar, peak_subtractor, price, export_limit, solver="gurobi"):
 
         """The process function creates the pyomo model and solve it.
         Minimize sum( net_load(t) + sum(power_demand(t, v)))**2
@@ -383,7 +409,7 @@ def process(self, times, vehicles, d, pmax, pmin, emin, emax,
             #HYBRID
             if peak_shaving == 'hybrid':
                 # Lambda for hybrid model--1 is all weight on peak shaving
-                model.lamb = Param(initialize = .5)
+                model.lamb = Param(initialize = 1)
                 # Scaling factor for peak-shaving sub-objective for hybrid model
                 model.peak_scalar = Param(initialize = peak_scalar)
                 # Normalization factor for peak-shaving sub-objective
@@ -397,47 +423,47 @@ def process(self, times, vehicles, d, pmax, pmin, emin, emax,
             #EMERGENCY
             if peak_shaving == 'emergency':
                 # Maximum Wh a vehicle can export over course of optimization (battery capacity * given fraction)
-                model.max_calls = Param(model.v, initialize = max_calls)
+                model.export_limit = Param(model.v, initialize = export_limit)
 
 
             # ###### Variable
-            
-
-            #model.b = Var(model.t, model.v, domain=Binary, doc='Power exported at timestep?') #should it be within=Binary?
-            #model.u_neg = Var(model.t, model.v, domain=Integers, doc='Power used')
-            #model.u_pos = Var(model.t, model.v, domain=Integers, doc='Power used') 
+            if peak_shaving == 'emergency':
+                model.u_neg = Var(model.t, model.v, domain=Integers, doc='Power exported')
+                model.u_pos = Var(model.t, model.v, domain=Integers, doc='Power imported')
             model.u = Var(model.t, model.v, domain=Integers, doc='Power used')
 
 
             # ###### Rules
-            def maximum_power_rule(model, t, v):
-                return model.u[t, v] <= model.p_max[t, v]
-            model.power_max_rule = Constraint(model.t, model.v, rule=maximum_power_rule, doc='P max rule')
 
-            def minimum_power_rule(model, t, v):
-                return model.u[t, v] >= model.p_min[t, v]
-            model.power_min_rule = Constraint(model.t, model.v, rule=minimum_power_rule, doc='P min rule')
-            # you'd have four rules instead: 1 for charging (>0, <pmax), 1 for discharging (<0, >pmin)
-            # u is just the sum of those two things. basically the same thing, but now 3 decision vars
+            if peak_shaving == 'emergency':
+                def max_charging_rule(model, t, v):
+                    return model.u_pos[t, v] <= model.p_max[t, v]
+                model.max_charging_rule = Constraint(model.t, model.v, rule=max_charging_rule, doc='u_pos <= p_max')
 
-            #do we need to worry about big neg + big +? constrain one to be 0 if other is not 0?
-            # if i do whole problem w u_neg and u_post instead of u, takes care of the issue
-            # where u_neg is super neg & vv
-            # bc when calculating the amount of energy in the battery, there's an efficiency term
-            # assigned to the in and the out...bc if both were values, there'd be more losses
-            # 
+                def min_charging_rule(model, t, v):
+                    return model.u_pos[t, v] >= 0
+                model.min_charging_rule = Constraint(model.t, model.v, rule=min_charging_rule, doc='u_pos >= 0')
 
-            #def maximum_power_rule(model, t, v):
-            #    return model.u_pos[t, v] <= model.p_max[t, v]
-            #model.power_max_rule = Constraint(model.t, model.v, rule=maximum_power_rule, doc='P max rule')
+                def max_discharging_rule(model, t, v):
+                    return model.u_neg[t, v] >= model.p_min[t, v]
+                model.max_discharging_rule = Constraint(model.t, model.v, rule=max_discharging_rule, doc='u_neg >= p_min')
 
-            #def positive_power_rule(model, t, v):
-             #   return model.u_pos[t, v] >= 0
-            #rederfine, and do for 2 negative ones -- model.power_max_rule = Constraint(model.t, model.v, rule=maximum_power_rule, doc='P max rule')
+                def min_discharging_rule(model, t, v):
+                    return model.u_neg[t, v] <= 0
+                model.min_discharging_rule = Constraint(model.t, model.v, rule=min_discharging_rule, doc='u_neg <= 0')
 
-            #def u_sum(model, t, v):
-             #   return model.u = model.u_neg + model.u_pos #...
+                def total_power_rule(model, t, v):
+                    return model.u[t, v] == model.u_neg[t, v] + model.u_pos[t, v]
+                model.total_power_rule = Constraint(model.t, model.v, rule=total_power_rule, doc='u = u_neg + u_pos')
+            
+            elif peak_shaving != 'emergency':
+                def maximum_power_rule(model, t, v):
+                    return model.u[t, v] <= model.p_max[t, v]
+                model.power_max_rule = Constraint(model.t, model.v, rule=maximum_power_rule, doc='P max rule')
 
+                def minimum_power_rule(model, t, v):
+                    return model.u[t, v] >= model.p_min[t, v]
+                model.power_min_rule = Constraint(model.t, model.v, rule=minimum_power_rule, doc='P min rule')
 
             def minimum_energy_rule(model, t, v):
                 return sum(model.u[i, v] for i in range(0, t + 1)) >= model.e_min[t, v]
@@ -451,22 +477,6 @@ def final_energy_balance(model, v):
                 return sum(model.u[i, v] for i in model.t) >= model.e_final[v]
             model.final_energy_rule = Constraint(model.v, rule=final_energy_balance, doc='E final rule')
 
-            #def not_two_values:
-             #   return minimum( maximum(u_neg, 0), minimum(u_pos, 0))==0 # too hard...can you even put min/max in constraints? no?
-                # have to somehow penalize them with an efficiency thing to prevent simultaneous charging/discharging
-                # need variable that's gonna be relative to the losses...add a little bit of each charge/discharge
-                # to objective function (an additive term...just like adding a constant...doesn't have to be on same scale, but if scale too small, could be ~epsilon). 3 orders of magnitude
-                # def wanna look at u_neg/u_pos afterwards to see if one is 0 and other not*****
-                # check out google pyomo forum; respnsive
-                # if optimization works, it works
-                # could do fake losses terms where you multiply u_neg*0.9 &
-
-            #def power_export_rule(model, t, v):
-             #   return model.b[t, v] for model.u[t, v] < 0 == 1 ###this is the tricky part ## b is min of 0 or u. if else ok in constr
-                ## if u is +...if t is something, but NOT variables
-                # constraint: model.b <0
-                # constraint: model.b > u
-            
 
             # Set the objective to be either peak shaving or ramp mitigation
             if peak_shaving == 'peak_shaving':
@@ -498,36 +508,15 @@ def objective_rule(model):
             # Definition: Emergency V2G; mostly V1G but each vehicle can get called for V2G a certain number of times/year
             # Can be used either for peak shaving or ramp mitigation; for now just implementing peak shaving
             elif peak_shaving == 'emergency':
+                # Peak shaving objective
                 def objective_rule(model):
                     return sum([(model.d[t] + sum([model.u[t, v] for v in model.v]))**2 for t in model.t])
-                    # minimize obj function that has model.b term in it
-                    # u should be cut into two parts--charging and discharging
                 model.objective = Objective(rule=objective_rule, sense=minimize, doc='Define objective function')
 
-                
-                ### Put constraint on negative part!!
-                # define another variable: model.c (counting) 
-
-                #model.c = Var(...)
-
-              #  rule count(): # try to force count to be 0 unless it can be below 0...model.c > model.u...so when u is positive
-                ### make it follow u, and set clear limits (can't be below 0), and ensure that optimization will
-                # try to make it 0 most of the time, except when it has the ability to not make it 0, only when
-                # you want to count one of the session. only time optimization can minimize this term (ie, be something /= 0)
-                # is when u is negative
-
-#                def negative_power(model, v):
-             #       return sum([model.u_neg for ])
-
-                # Could use mod to constrain sessions of discharging...as soon as you have X timesteps
-                # consecutive that you're discharging...harder
-                # if you want to count how many times u_neg is below zero, easy--if < 0, +1 this var
-                # this would include counting variable in objective function...
-                
-                def max_calls_rule(model, v):
-                    return sum([model.b[t, v] for t in model.t]) <= model.max_calls[v] 
-                    #return len([x for x in [model.u[t, v] for t in model.t] if x < 0]) <= model.max_calls[v] 
-                model.max_calls_rule = Constraint(model.v, rule=max_calls_rule, doc='Max output rule')
+                # V2G (export) limits
+                def max_export_rule(model, v): 
+                    return sum(model.u_neg[i, v] for i in model.t) >= model.export_limit[v]
+                model.max_export_rule = Constraint(model.v, rule=max_export_rule, doc='Max export rule')
 
             # Cost-based
             # Definition: minimizes charging cost based on given cost stack

From 6f83ff5669460e248ce33f4bac30252429a09dc2 Mon Sep 17 00:00:00 2001
From: Margaret McCall <majemc@gmail.com>
Date: Thu, 15 Oct 2020 12:35:44 -0500
Subject: [PATCH 13/14] opti plus

---
 v2gsim/post_simulation/netload_opti_plus.py | 42 +++++++++------------
 1 file changed, 18 insertions(+), 24 deletions(-)

diff --git a/v2gsim/post_simulation/netload_opti_plus.py b/v2gsim/post_simulation/netload_opti_plus.py
index 01d5e53..bb012b4 100644
--- a/v2gsim/post_simulation/netload_opti_plus.py
+++ b/v2gsim/post_simulation/netload_opti_plus.py
@@ -31,7 +31,7 @@ def __init__(self, project, optimization_timestep, date_from,
         self.SOC_index_from = int((date_from - project.date).total_seconds() / project.timestep)
         self.SOC_index_to = int((date_to - project.date).total_seconds() / project.timestep)
 
-    def solve(self, project, net_load, unoptimized_demand, real_number_of_vehicle, SOC_margin=0.02, 
+    def solve(self, project, net_load, unoptimized_demand, real_number_of_vehicle, SOC_init_opti=None, SOC_end_opti=None, SOC_margin=0.02, ##ADDED
               SOC_offset=0.0, peak_shaving='peak_shaving', penalization=5, beta=None, plot=False, peak_scalar=.01, peak_subtractor=30000, 
               price=pandas.DataFrame(), max_export=-3000): #price can't be array, is dict; can't have non-default arg follow default
         """Launch the optimization and the post_processing fucntion. Results
@@ -58,9 +58,10 @@ def solve(self, project, net_load, unoptimized_demand, real_number_of_vehicle, S
         self.export_limit = {}
 
 
+
         # Set the variables for the optimization
         new_net_load, new_unoptimized_demand = self.initialize_net_load(net_load, real_number_of_vehicle, project, unoptimized_demand)
-        self.initialize_model(project, new_net_load, SOC_margin, SOC_offset, max_export) 
+        self.initialize_model(project, new_net_load, SOC_margin, SOC_offset, SOC_init_opti, SOC_end_opti, max_export) ##ADDED
         
         if peak_shaving=='cost':
             price = self.initialize_price(price, net_load)
@@ -80,7 +81,7 @@ def solve(self, project, net_load, unoptimized_demand, real_number_of_vehicle, S
         print('')
 
         # Post process results
-        return self.post_process(project, net_load, opti_model, result, plot), opti_model ##ADDED OPTI_MODEL
+        return self.post_process(project, net_load, opti_model, result, plot), opti_model, result ##ADDED OPTI_MODEL ad RESULT
         #return opti_model, result
 
     def get_peak_scalar(self, net_load_updated, unoptimized_demand):
@@ -191,10 +192,6 @@ def check_energy_constraints_feasible(self, vehicle, SOC_init, SOC_final, SOC_of
         def print_status(SOC_final, SOC_init, diff):
                 print('Set battery gain: ' + str((SOC_final - SOC_init) * 100) + '%')
                 print('Simulation battery gain: ' + str(diff * 100))
-
-        print('CHECKING SOC for vehicle '+str(vehicle.id))
-        print('SOC_init: '+str(SOC_init))
-        print('SOC_final: '+str(SOC_final))
         
 
         # Check if below minimum SOC at any time
@@ -207,7 +204,7 @@ def print_status(SOC_final, SOC_init, diff):
         # Check SOC difference between date_from and date_to ?
         # Diff represent the minimum loss or the maximum gain
         diff = vehicle.SOC[self.SOC_index_to] - vehicle.SOC[self.SOC_index_from]
-        print('Diff between opti SOC diff and unopti SOC diff: '+str(SOC_final - SOC_init - diff))
+        
         # Simulation show a battery gain
         if diff > 0:
             # Gain should be greater than the one we set up
@@ -264,7 +261,7 @@ def get_final_SOC(self, vehicle, SOC_margin, SOC_offset, SOC_end=None):
         else:
             return vehicle.SOC[self.SOC_index_to] - SOC_offset - SOC_margin
 
-    def initialize_model(self, project, net_load, SOC_margin, SOC_offset, max_export):
+    def initialize_model(self, project, net_load, SOC_margin, SOC_offset, SOC_init_opti, SOC_end_opti, max_export): ##ADDED
         """Select the vehicles that were plugged at controlled chargers and create
         the optimization variables (see inputs of optimization)
 
@@ -274,6 +271,9 @@ def initialize_model(self, project, net_load, SOC_margin, SOC_offset, max_export
         Return:
             times, vehicles, d, pmax, pmin, emin, emax, efinal
         """
+
+
+
         # Create a dict with the net load and get time index in a data frame
         self.d, self.times = self.initialize_time_index(net_load)
         vehicle_to_optimize = 0
@@ -282,9 +282,13 @@ def initialize_model(self, project, net_load, SOC_margin, SOC_offset, max_export
         for vehicle in project.vehicles:
             if vehicle.result is not None:
                 # Get SOC init and SOC end
-                print("MADE IT TO SOC ASSIGNMENT")
-                SOC_init = self.get_initial_SOC(vehicle, SOC_offset)
-                SOC_final = self.get_final_SOC(vehicle, SOC_margin, SOC_offset)
+                
+                if SOC_init_opti is not None: ##ADDED
+                    SOC_init = self.get_initial_SOC(vehicle, SOC_offset, SOC_init=SOC_init_opti[vehicle.id])
+                    SOC_final = self.get_final_SOC(vehicle, SOC_margin, SOC_offset, SOC_end=SOC_end_opti[vehicle.id])
+                elif SOC_init_opti is None:
+                    SOC_init = self.get_initial_SOC(vehicle, SOC_offset)
+                    SOC_final = self.get_final_SOC(vehicle, SOC_margin, SOC_offset)
 
                 # Find out if vehicle itinerary is feasible
                 if not self.check_energy_constraints_feasible(vehicle, SOC_init, SOC_final, SOC_offset):
@@ -330,16 +334,6 @@ def initialize_model(self, project, net_load, SOC_margin, SOC_offset, max_export
                                                  (SOC_final - SOC_init) * vehicle.car_model.battery_capacity *
                                                  (60 / self.optimization_timestep))})
 
-                print("")
-                print(vehicle.id)
-                print('last energy value:')
-                print(temp_vehicle_result.tail(1).energy.values[0])
-                print('SOC_final: '+str(SOC_final))
-                print('SOC_init: '+str(SOC_init))
-                print('HALF 1: '+str((temp_vehicle_result.tail(1).energy.values[0] * (project.timestep / (60 * self.optimization_timestep)))))
-                print('HALF 2: '+str((SOC_final - SOC_init) * vehicle.car_model.battery_capacity *(60 / self.optimization_timestep)))
-                print('EFINAL: '+str(self.efinal[vehicle.id]))
-                print("")
 
                 # Create export_limit for each vehicle for emergency optimization ######################
                 self.export_limit.update({vehicle.id: max_export})
@@ -528,6 +522,7 @@ def objective_rule(model):
                     #return sum((model.u[t, v] * model.price[t] for v in model.v) for t in model.t)
                 model.objective = Objective(rule=objective_rule, sense=minimize, doc='Define objective function')
             
+        
             solver_parameters = "ResultFile=model.ilp"
             results = opt.solve(model, options_string=solver_parameters, symbolic_solver_labels=True)
             # results.write()
@@ -698,7 +693,6 @@ def save_vehicle_state_for_optimization(vehicle, timestep, date_from,
                         # Energy is 0.0 because it's already accounted in power_demand
                         vehicle.result['energy'][location_index1:location_index2] -= (
                             [0.0] * (activity_index2 - activity_index1))
-
     elif init:
         vehicle.SOC = [vehicle.SOC[0]]
         vehicle.result = None
@@ -718,4 +712,4 @@ def save_vehicle_state_for_optimization(vehicle, timestep, date_from,
             i = pandas.date_range(start=date_from, end=date_to,
                                   freq=str(timestep) + 's', closed='left')
             vehicle.result = pandas.DataFrame(index=i, data=vehicle.result)
-            vehicle.result['energy'] = vehicle.result['energy'].cumsum()
+            vehicle.result['energy'] = vehicle.result['energy'].cumsum()
\ No newline at end of file

From afa11b852ec538a3e01ad98ad8a11ecf1ba41db4 Mon Sep 17 00:00:00 2001
From: Margaret McCall <majemc@gmail.com>
Date: Fri, 18 Dec 2020 11:58:59 -0600
Subject: [PATCH 14/14] working vers of new optis

---
 v2gsim/post_simulation/netload_opti_emergency.py | 10 ++++++----
 v2gsim/post_simulation/netload_opti_plus.py      |  5 +++--
 2 files changed, 9 insertions(+), 6 deletions(-)

diff --git a/v2gsim/post_simulation/netload_opti_emergency.py b/v2gsim/post_simulation/netload_opti_emergency.py
index 9890d66..d01d263 100644
--- a/v2gsim/post_simulation/netload_opti_emergency.py
+++ b/v2gsim/post_simulation/netload_opti_emergency.py
@@ -33,7 +33,7 @@ def __init__(self, project, optimization_timestep, date_from,
 
     def solve(self, project, net_load, real_number_of_vehicle, SOC_margin=0.02,
               SOC_offset=0.0, peak_shaving='peak_shaving', penalization=5, beta=None, plot=False, peak_scalar=.01, peak_subtractor=30000, 
-              unoptimized_demand = pandas.DataFrame(), price=pandas.DataFrame(), max_export=-3000): #price can't be array, is dict; can't have non-default arg follow default
+              unoptimized_demand = pandas.DataFrame(), price=pandas.DataFrame(), max_export=.1): #price can't be array, is dict; can't have non-default arg follow default
         """Launch the optimization and the post_processing fucntion. Results
         and assumptions are appended to a data frame.
 
@@ -73,13 +73,13 @@ def solve(self, project, net_load, real_number_of_vehicle, SOC_margin=0.02,
         timer = time.time()
         opti_model, result = self.process(self.times, self.vehicles, self.d, self.pmax,
                                           self.pmin, self.emin, self.emax,
-                                          self.efinal, peak_shaving, penalization, peak_scalar, peak_subtractor, price, max_export)
+                                          self.efinal, peak_shaving, penalization, peak_scalar, peak_subtractor, price, self.export_limit)
         timer2 = time.time()
         print('The optimization duration was ' + str((timer2 - timer) / 60) + ' minutes')
         print('')
 
         # Post process results
-        return self.post_process(project, net_load, opti_model, result, plot)
+        return self.post_process(project, net_load, opti_model, result, plot), opti_model, result ##ADDED OPTI_MODEL ad RESULT
         #return opti_model, result
 
     def get_peak_scalar(self, net_load_updated, unoptimized_demand):
@@ -133,11 +133,13 @@ def initialize_net_load(self, net_load, real_number_of_vehicle, project, unoptim
         # Check against the actual lenght it should have
         diff = (len(new_net_load) -
                 int((self.date_to - self.date_from).total_seconds() / (60 * self.optimization_timestep)))
+        
         if diff > 0:
             # We should trim the net load with diff elements (normaly 1 because of slicing inclusion)
             new_net_load.drop(new_net_load.tail(diff).index, inplace=True)
         elif diff < 0:
             print('The net load does not contain enough data points')
+            
 
         if real_number_of_vehicle:
             # Set scaling factor
@@ -317,7 +319,7 @@ def initialize_model(self, project, net_load, SOC_margin, SOC_offset, max_export
                                                  (60 / self.optimization_timestep))})
 
                 # Create export_limit for each vehicle for emergency optimization ######################
-                self.export_limit.update({vehicle.id: max_export})
+                self.export_limit.update({vehicle.id: -1 * vehicle.car_model.battery_capacity * max_export})
 
         print('There is ' + str(vehicle_to_optimize) + ' vehicle participating in the optimization (' +
               str(vehicle_to_optimize * 100 / len(project.vehicles)) + '%)')
diff --git a/v2gsim/post_simulation/netload_opti_plus.py b/v2gsim/post_simulation/netload_opti_plus.py
index bb012b4..4f3603f 100644
--- a/v2gsim/post_simulation/netload_opti_plus.py
+++ b/v2gsim/post_simulation/netload_opti_plus.py
@@ -403,7 +403,7 @@ def process(self, times, vehicles, d, pmax, pmin, emin, emax,
             #HYBRID
             if peak_shaving == 'hybrid':
                 # Lambda for hybrid model--1 is all weight on peak shaving
-                model.lamb = Param(initialize = 1)
+                model.lamb = Param(initialize = 0)
                 # Scaling factor for peak-shaving sub-objective for hybrid model
                 model.peak_scalar = Param(initialize = peak_scalar)
                 # Normalization factor for peak-shaving sub-objective
@@ -523,7 +523,8 @@ def objective_rule(model):
                 model.objective = Objective(rule=objective_rule, sense=minimize, doc='Define objective function')
             
         
-            solver_parameters = "ResultFile=model.ilp"
+            solver_parameters = "ResultFile=model.ilp, LogFile=model_log.log"
+            print('actually got to the point of running solve')
             results = opt.solve(model, options_string=solver_parameters, symbolic_solver_labels=True)
             # results.write()