1from .centrals import PowerCentral as pc 2from .centrals import HydroCentral as hc 3from .Demand import Demand 4import nevergrad as ng 5from .nevergradBased.Optimizer import Optimizer 6import numpy as np # type: ignore 7import pandas as pd # type: ignore 8import pkgutil 9import csv 10import os 11import warnings 12from math import ceil 13#import time 14from typing import List, Any, Type, Dict 15from datetime import datetime 16import matplotlib.pyplot as plt # type: ignore 17 18class MixSimulator: 19 """ 20 The simulator Base 21 """ 22 def __init__(self, carbon_cost: float = 0, penalisation_cost: float = 1000000000000) -> None: 23 self.__centrals : List[Any] = [] 24 self.__reset_centrals() 25 self.__demand = Demand(20, 0.2, 0.2) 26 self.__lost = 0. 27 self.__penalisation_cost = penalisation_cost 28 self.__optimizer = Optimizer() 29 self.__carbon_cost = carbon_cost 30 self.__carbon_quota = 800139. # (g/MWh) 31 32 # for reuse and get_params() 33 self.time_index = 24*7 34 self.step = 1 35 self.time_interval = 1 36 self.plot = "default" 37 self.average_wide = 0 38 39 def __reset_centrals(self): 40 self.__centrals = [] 41 42 ## DATA ## 43 def set_variation_csv(self, bind = None, delimiter: str=";") -> None: 44 if self.__centrals == []: 45 warnings.warn("Please load a original dataset by using MixSimulator.set_data_csv(...)") 46 raise 47 else : 48 try : 49 data = pd.DataFrame(pd.read_csv(bind,delimiter=delimiter)) 50 except FileNotFoundError as e : 51 print("Error occured on pandas.read_csv : ",e) 52 print("Please check your file") 53 raise 54 except Exception as e: 55 print("Error occured on pandas.read_csv : ",e) 56 raise 57 58 for central_index in range(len(self.__centrals)): 59 if self.__centrals[central_index].is_tuneable(): 60 for i in range (0,data.shape[0]): 61 if self.__centrals[central_index].get_id() == data["centrals"][i]: 62 tmp_list=[] 63 upper = str(data["upper"][i]).split(":") 64 upper = [float(numeric_string) for numeric_string in upper] 65 lower = str(data["lower"][i]).split(":") 66 lower = [float(numeric_string) for numeric_string in lower] 67 discrete = str(data["discrete"][i]).split(":") 68 discrete = [float(numeric_string) for numeric_string in discrete] 69 self.__centrals[central_index].set_variation_params(lower = lower, upper = upper, choices = discrete) 70 71 def set_data_csv(self, bind = None, raw_data = None, delimiter: str=";"): 72 if raw_data is not None : 73 data = pd.DataFrame(raw_data) 74 #set columns & index 75 header = data.iloc[0] 76 data = data[1:] 77 data.columns = header 78 data = data.reset_index(drop=True) 79 for column in data.columns.tolist(): 80 try: 81 # convert numeric values 82 data[column] = pd.to_numeric(data[column]) 83 except: 84 pass 85 else : 86 try : 87 data = pd.DataFrame(pd.read_csv(bind,delimiter=delimiter)) 88 except FileNotFoundError as e : 89 print("Error occured on pandas.read_csv : ",e) 90 print("Please check your file") 91 raise 92 except Exception as e: 93 print("Error occured on pandas.read_csv : ",e) 94 raise 95 96 self.__reset_centrals() 97 centrale = pc.PowerCentral() 98 try : 99 for i in range (0,data.shape[0]): 100 isHydro = data["hydro"][i] 101 if isHydro == True : 102 centrale = hc.HydroCentral(data["height"][i],data["flow"][i],data["capacity"][i],data["stock_available"][i],0.1,0.8) 103 else : 104 centrale = pc.PowerCentral() 105 centrale.set_tuneable(data["tuneable"][i]) 106 centrale.set_id(str(data["centrals"][i])) 107 centrale.set_fuel_consumption(data["fuel_consumption"][i]) 108 centrale.set_availability(data["availability"][i]) 109 centrale.set_fuel_cost(data["fuel_cost"][i]) 110 centrale.set_initial_value(data["init_value"][i]) 111 centrale.set_lifetime(data["lifetime"][i]) 112 centrale.set_carbon_prod(data["carbon_production"][i]) 113 centrale.set_raw_power(data["raw_power"][i]) 114 centrale.set_nb_employees(data["nb_employees"][i]) 115 centrale.set_mean_employees_salary(data["mean_salary"][i]) 116 centrale.set_max_var(data["max_var"][i]) 117 self.__centrals.append(centrale) 118 self.__demand.set_mean_demand(data["Demand"][0]) 119 self.__lost=data["lost"][0] 120 except KeyError: 121 print("Columns must be in: tuneable, centrals, fuel_consumption, availability, fuel_cost, init_value, lifetime, carbon_cost, raw_power, nb_employees, mean_salary, demand, lost, height, flow, capacity, stock_available") 122 raise 123 124 def set_data_to(self, dataset, delimiter: str=";"): 125 if dataset == "Toamasina": 126 #by defaut we keep it "Toamasina" 127 data = pkgutil.get_data('mixsimulator', '/data/RIToamasina/dataset_RI_Toamasina_v2.csv') 128 data = csv.reader(data.decode('utf-8').splitlines(), delimiter = delimiter) 129 self.set_data_csv(raw_data=data) 130 else : 131 #by defaut we keep it "Toamasina" 132 data = pkgutil.get_data('mixsimulator', '/data/RIToamasina/dataset_RI_Toamasina_v2.csv') 133 data = csv.reader(data.decode('utf-8').splitlines(), delimiter = delimiter) 134 self.set_data_csv(raw_data=data) 135 136 def set_demand(self, demand: Demand): 137 self.__demand = demand 138 139 def get_demand(self): 140 return self.__demand 141 142 def set_lost(self, lost: float): 143 self.__lost = lost 144 145 def set_optimizer(self, optimizer: Optimizer): 146 self.__optimizer = optimizer 147 148 def get_optimizer(self) -> Optimizer: 149 return self.__optimizer 150 151 def set_carbon_quota(self, cb_quota: float ): 152 self.__carbon_quota = cb_quota 153 154 def set_penalisation_cost(self, k): 155 self.__penalisation_cost = k 156 157 def get_penalisation_cost(self): 158 return self.__penalisation_cost 159 160 def set_carbon_cost(self, carbon_cost): 161 self.__carbon_cost = carbon_cost 162 163 def get_carbon_cost(self): 164 return self.__carbon_cost 165 166 ## EVALUATION FONCTION ## 167 def get_production_cost_at_t(self, usage_coef, time_index, time_interval): 168 production_cost = 0 169 for centrale_index in range (0, len(self.__centrals)): 170 central = self.__centrals[centrale_index] 171 production_cost += ( (central.get_fuel_cost() * central.get_fuel_consumption() 172 * central.get_raw_power() * usage_coef[centrale_index]) + ( central.get_employees_salary() 173 + central.get_amortized_cost(time_index) ) ) * time_interval 174 return production_cost 175 176 def get_production_at_t(self, usage_coef, time_interval): 177 production = 0 178 for centrale_index in range (0, len(self.__centrals)): 179 central = self.__centrals[centrale_index] 180 production += (central.get_raw_power() * usage_coef[centrale_index]) * time_interval 181 return production 182 183 def get_unsatisfied_demand_at_t(self, usage_coef, time_index, time_interval): 184 #return ( self.__demand.get_demand_approxima(time_index, time_interval) - self.get_production_at_t(usage_coef, time_interval)) 185 return ( self.__demand.get_demand_monthly(time_index, time_interval) - self.get_production_at_t(usage_coef, time_interval)) 186 187 def get_carbon_production_at_t(self, usage_coef, time_interval): 188 carbon_prod = 0 189 for centrale_index in range (0, len(self.__centrals)): 190 central = self.__centrals[centrale_index] 191 carbon_prod += (central.get_carbon_production() * usage_coef[centrale_index] * central.get_raw_power()) * time_interval 192 return carbon_prod 193 194 def get_carbon_over_production(self, usage_coef, time_interval): 195 emited_carbon = 0 # (g) 196 total_production = 1 197 for t in range(0, len(usage_coef)): 198 emited_carbon += self.get_carbon_production_at_t(usage_coef[t], time_interval) 199 total_production += self.get_production_at_t(usage_coef[t], time_interval) 200 carbon_production = emited_carbon/total_production 201 return max(0, carbon_production - self.__carbon_quota) # (g/MWh) 202 203 def get_weighted_coef(self, usage_coef, time_interval): 204 for central_index in range(0, len(usage_coef[0])): 205 self.__centrals[central_index].reset_central() 206 weighted_coef = usage_coef.copy() 207 for t in range(0, len(weighted_coef)): 208 for central_index in range(0, len(weighted_coef[t])): 209 min_av = self.__centrals[central_index].get_min_availability(t) 210 max_av = self.__centrals[central_index].get_max_availability(t) 211 if max_av < min_av: 212 min_av = 0 # a verifier 213 weighted_coef[t][central_index] = min_av + weighted_coef[t][central_index]*(max_av-min_av) 214 self.__centrals[central_index].back_propagate(weighted_coef[t][central_index], t, time_interval) 215 return weighted_coef 216 217 def loss_function(self, usage_coef, time_interval : int = 1) -> float : 218 usage_coef = self.arrange_coef_as_array_of_array(usage_coef) 219 weighted_coef = self.get_weighted_coef(usage_coef, time_interval=time_interval) 220 loss = 0 221 for t in range(0, len(weighted_coef)): 222 loss += self.get_production_cost_at_t(weighted_coef[t], t, time_interval) + ( self.get_penalisation_cost() * np.abs( self.get_unsatisfied_demand_at_t(weighted_coef[t], t, time_interval)) ) 223 loss += self.get_carbon_cost() * (self.get_carbon_over_production(weighted_coef, time_interval) ) 224 return loss 225 226 def arrange_coef_as_array_of_array(self, raw_usage_coef): 227 ordered_coef = [] 228 cur_time_coef = [] 229 for coef_index in range(len(raw_usage_coef)): 230 cur_time_coef.append(raw_usage_coef[coef_index]) 231 ## indice de la premiere centrale a t+1 232 if (coef_index+1)%len(self.__centrals) == 0: 233 ordered_coef.append(cur_time_coef) 234 cur_time_coef = [] 235 return ordered_coef 236 237 238 ## OPTiMiZATION ## 239 240 def __opt_params(self, time_index): 241 self.__optimizer.set_parametrization(self.get_opt_params(time_index)) 242 243 def get_opt_params(self, time_index): 244 variable_parametrization = [] 245 for _ in range(time_index): 246 for central_index in range(len(self.__centrals)): 247 if not self.__centrals[central_index].is_tuneable(): 248 variable_parametrization += [ng.p.Choice([0.,1.])] 249 else: 250 #check the params by 251 #print(self.__centrals[central_index].get_variation_params()) 252 variable_parametrization += [self.__centrals[central_index].get_variation_params()] 253 return ng.p.Tuple(*variable_parametrization) 254 255 def optimizeMix(self, carbon_quota: float = None, demand: Demand = None, lost: float = None, 256 optimizer: Optimizer = None, step : int = 1, 257 time_index: int = 24*7, time_interval: float = 1, 258 penalisation : float = None, carbon_cost : float = None, plot : str = "default", average_wide : int = 0): 259 260 # init params 261 self.time_index = time_index 262 self.step = step 263 self.time_interval = time_interval 264 self.plot = plot 265 self.average_wide = average_wide 266 if demand is not None : self.set_demand(demand) 267 if lost is not None : self.set_lost(lost) 268 if penalisation is not None : self.set_penalisation_cost(penalisation) 269 if carbon_cost is not None : self.set_carbon_cost(carbon_cost) 270 if carbon_quota is not None : self.set_carbon_quota(carbon_quota) 271 if optimizer is not None : self.set_optimizer(optimizer) 272 273 # tune optimizer parametrization 274 self.__opt_params(self.time_index) 275 276 #init constraints 277 constraints = {} 278 constraints.update({"time_interval":self.time_interval}) 279 280 #let's optimize 281 results = self.__optimizer.optimize(mix = self , func_to_optimize = self.loss_function, constraints=constraints, step = self.step) 282 283 results = self.__reshape_results(results, self.time_interval) 284 285 self.plotResults(results, mode = self.plot , time_interval = self.time_interval, average_wide = self.average_wide) 286 287 return results 288 289 def __reshape_results(self, results, time_interval): 290 for tmp in results: 291 usage_coef = self.arrange_coef_as_array_of_array(tmp['coef']) 292 tmp.update({"coef":self.get_weighted_coef(usage_coef, time_interval)}) 293 294 # for central_index in range(0, len(self.__centrals)): 295 # try: 296 # print(self.__centrals[central_index].get_stock_evolution()) 297 # # Not a hydro power plant, so the methode does not exist 298 # except: 299 # pass 300 return results 301 302 def get_params(self) -> dict: 303 return {"centrals" : self.__centrals, "optimizer" : self.get_optimizer(), 304 "penalisation_cost" : self.get_penalisation_cost(), "carbon_cost" : self.get_carbon_cost(), 305 "demand" : self.__demand, "lost" : self.__lost, "carbon_quota" : self.__carbon_quota, 306 "step" : self.step, "time_interval" : self.time_interval, "time_index" : self.time_index, 307 "plot" : self.plot, "moving average_wide" : self.average_wide} 308 309 ## PLOT ## 310 def moving_average(self, x, w): 311 return np.convolve(x, np.ones(w), 'valid') / w 312 313 def plotResults(self, optimum : dict = {} , mode : str = "default", time_interval : float = 1, average_wide : int = 0): 314 #set the moving average wide 315 if average_wide == 0 : 316 average_wide = ceil(len(optimum[-1]["coef"])/4) 317 318 if mode == "default" or mode == "save": 319 #set Y 320 Y: Dict[str,List[float]] ={} 321 label_y: List[str]=[] 322 for c in self.__centrals : 323 label_y.append(c.get_id()) 324 Y.update({c.get_id():[]}) 325 for array in optimum[-1]["coef"] : 326 for enum, value in enumerate(array) : 327 Y[label_y[enum]].append(value) 328 329 fig, axs = plt.subplots(1, 1, figsize=(6, 6)) 330 331 # data integration 332 X = [i for i in range(0,len(optimum[-1]["coef"]))] 333 for n_axs in range(0,1) : 334 for central, value in Y.items(): 335 smooth_value = self.moving_average(value,average_wide) 336 axs.plot(X[(average_wide - 1):], smooth_value, '.-' ,alpha=0.5, lw=2, label=central) 337 338 339 # Add execution_time and loss information 340 info = "production_cost: "+ "{:.3f}".format(optimum[-1]["loss"])+" ; execution_time: "+"{:.3f}".format(optimum[-1]["elapsed_time"]) 341 plt.annotate(info, 342 xy=(0.5, 0), xytext=(0, 10), 343 xycoords=('axes fraction', 'figure fraction'), 344 textcoords='offset points', 345 size=10, ha='center', va='bottom') 346 347 # plots parametrizations 348 axs.grid() 349 axs.yaxis.set_tick_params(length=0) 350 axs.xaxis.set_tick_params(length=0) 351 axs.set_xlabel('hours') 352 #axs[n_axs].yaxis.set_major_formatter(StrMethodFormatter("{x}"+units[0])) 353 axs.set_ylabel('coef (%)') 354 axs.legend() 355 356 fig.tight_layout() 357 try : 358 path = "results_"+datetime.now().strftime("%d_%m_%Y") 359 os.makedirs(path) 360 name = path+"/"+"opt_"+str(self.get_optimizer().get_optimizers())+"_"+datetime.now().strftime("%H%M%S")+".png" 361 fig.savefig(name) 362 if mode == "default" : 363 plt.show() 364 365 except OSError: 366 warnings.warn("Can't create the directory "+path+" or already exists") 367 try : 368 name = path+"/"+"opt_"+str(self.get_optimizer().get_optimizers())+"_"+datetime.now().strftime("%H%M%S")+".png" 369 fig.savefig(name) 370 if mode == "default" : 371 plt.show() 372 except FileNotFoundError: 373 name = "opt_"+str(self.get_optimizer().get_optimizers())+"_"+datetime.now().strftime("%H%M%S")+".png" 374 fig.savefig(name) 375 if mode == "default" : 376 plt.show() 377 378 elif mode == "None" : 379 pass 380 else : 381 warnings.warn("Choose an available option : default, save or None") 382 #plt.show() 383