Commit 9a0920ba authored by brugeron's avatar brugeron

Create a linearization module

Implementation of the module in the utils.py
Test on the storage design file
parent dd9f7717
......@@ -10,6 +10,8 @@
##########################
*.bak
*.lp
*.csv
*.mps
#Eclipse Settings #
#####################
......
This source diff could not be displayed because it is too large. You can view the blob instead.
#! usr/bin/env python3
# coding=utf-8 #
"""
** This Module is an example of storage capacity optimisation.**
This example describes a simple microgrid with :
- A load profile known in advance
- An adjustable production unit, with a maximum power limit
- A storage system with power in charge and power in discharge
The objective consists on minimizing the storage capacity.
"""
from matplotlib import pyplot as plt
from pulp import LpStatus, GUROBI_CMD, COIN_CMD
import pickle
import os
import cProfile
from omegalpes.energy.energy_nodes import EnergyNode
from omegalpes.energy.units.consumption_units import FixedConsumptionUnit
from omegalpes.energy.units.production_units import FixedProductionUnit, VariableProductionUnit, IdentificationVariableProductionUnit
from omegalpes.energy.units.storage_units import StorageUnit, WaterTank, IdentificationWaterTank
from omegalpes.general.optimisation.model import OptimisationModel
from omegalpes.general.time import TimeUnit
from omegalpes.general.plots import plt, plot_quantity, plot_node_energetic_flows
from omegalpes.general.utils import save_energy_flows
from omegalpes.general.optimisation.elements import Objective
def ballon_solaire_backup(load_profile, loss_profile, production_p, backup_p,
emax_storage, emin_storage, einit_storage,
horizon, time_step, alpha, u_load):
"""
:param production_pmax: Maximal power delivered by the production unit [kW]
:type production_pmax: float or int
:param storage_pcharge_max: Maximal charging power for the storage [kW]
:type storage_pcharge_max: float or int
:param storage_pdischarge_max: Maximal discharging power for the storage [kW]
:type storage_pdischarge_max: float or int
:param load_profile: hourly load profile during a day
:type load_profile: list of 24 float or int
"""
# global model, time, load, loss, production, storage, back_up, node
# Create an empty model
model = OptimisationModel(name='example') # Optimisation model
time = TimeUnit(start='01/01/2018', periods=horizon, dt=time_step/60) # Study on a day, hour
# Create the load - The load profile is known
load = FixedConsumptionUnit(time, 'load', p=load_profile)
load.u.value = u_load
# Create the losses - modelled as a consumption
loss = FixedConsumptionUnit(time, 'loss', p=loss_profile)
# Create the production unit - The production profile is unknown
production = FixedProductionUnit(time, 'production', p=production_p)
# Create back-up
back_up = IdentificationVariableProductionUnit(time, 'back_up', measurement = backup_p, pmin = 0)
back_up.min_power_abs()
# Create the storage
storage = IdentificationWaterTank(time, name='water_tank', soc_min = 0, measurement=backup_p,
soc_max= emax_storage/emax_storage, e_0=einit_storage, ef_is_e0 = False)#,
# Create the energy node and connect units
node = EnergyNode(time, 'energy_node')
# Add the energy node to the model
node.connect_units(load, loss, production, storage, back_up) # Connect all units on the same energy node
model.add_nodes(node)
# Optimisation process
model.writeLP('C:/Users/brugeron/Dimosim/OMEGAlpes/examples/optim_models/optimMPC.lp')
model.solve_and_update()
return model, time, load, loss, production, storage, back_up, node
def print_results():
"""
*** This function print the optimisation result:
- The optimal capacity
And plot the power curves :
Figure 1 :
- Power consumed by the storage, labelled 'Storage'
- Power consumed by the load, labelled 'Load'
- Power delivered by the production unit, labelled 'Production'
Figure 2 :
- The state of charge of the storage unit
"""
if LpStatus[MODEL.status] == 'Optimal':
print("\n - - - - - OPTIMISATION RESULTS - - - - - ")
# Show the graph
# Power curves
plot_node_energetic_flows(NODE)
# SOC curve
plot_quantity(TIME, STORAGE.e)
plt.xlabel('Time (h)')
plt.ylabel('Energy (kWh)')
plt.title('Storage temperature')
plt.show()
plot_quantity(TIME, BACK_UP.p)
plt.xlabel('Time (h)')
plt.ylabel('Power (kW)')
plt.title('Backup usage')
plt.show()
elif LpStatus[MODEL.status] == 'Infeasible':
print("Sorry, the optimisation problem has no feasible solution !")
elif LpStatus[MODEL.status] == 'Unbounded':
print("The cost function of the optimisation problem is unbounded !")
elif LpStatus[MODEL.status] == 'Undefined':
print("Sorry, a feasible solution has not been found (but may exist). "
"PuLP does not manage to interpret the solver's output, "
"the infeasibility of the MILP problem may have been "
"detected during presolve.")
else:
print("Sorry, the optimisation problem has not been solved.")
if __name__ == '__main__':
# OPTIMIZATION PARAMETERS #
WORK_PATH = os.getcwd()
# SIMULATION PARAMETERS #
with open("C:/Users/brugeron/Dimosim/dimosim_on_branch/dimosim/tests/test_external_control/Database.txt",
'rb') as file_:
stored = pickle.load(file_)
simu_time_step = stored['SIMU_TIME_STEP']
optim_time_step = stored['OPTIM_STEP_TIME']
horizon = stored['HORIZON']
tank_loss_factor = stored['THERMAL_LOSS_FACTOR']
temperature_loss = stored['TEMPERATURE_LOSS']
inter_load = stored['INTER_LOAD']
inter_prod = stored['INTER_PROD']
inter_losses = stored['INTER_LOSS']
real_load = stored['REAL_LOAD']
real_production = stored['REAL_PROD']
real_thermal_losses = stored['REAL_LOSS']
tank_temperature = stored['STATE_OF_CHARGE']
BACK_UP_POWER = stored['BACKUP_POWER']
backup_measurement = stored['BACKUP_MEASUREMENT']
volume = stored['VOLUME']
u_load = {}
for i in range(len(real_load)):
if real_load[i] > 0.:
u_load[i] = 1
else:
u_load[i] = 0
# horizon = round(24*3 * 60/optim_time_step)
horizon = len(real_production)
rho_water = 1000 # kg*/m3
cp_water = 4186 / 3600 # [Wh/kg.K]
storage_min_temperature = 57.5 + 273.15 # K
storage_max_temperature = 90 + 273.15 # K
(storage_Emin, storage_Emax) = round(storage_min_temperature * rho_water * cp_water * volume / 1000, 3), \
round(storage_max_temperature * rho_water * cp_water * volume / 1000, 3)
backup_normalized = 97649.0049
availability_normalized = 402
# à modifier à chaque reprise
new_load_signal = 3795
new_production_signal = 0
thermal_loss = 175
# optim_load = step_load.copy()
# optim_load[0] = new_load_signal
# optim_production = step_production.copy()
# optim_production[0] = new_production_signal
# optim_losses = step_thermal_losses.copy()
# optim_losses[0] = thermal_loss
# tank_temperature_float = float(np.array(received[obj]['tank_temperature']['value']))
# tank_temperature = round(tank_temperature_float, 0)
storage_init_temperature = tank_temperature + 273.15
storage_Einit = rho_water * volume * cp_water * storage_init_temperature / 1000
storage_Einit = round(storage_Einit, 3)
alpha = 0.7
MODEL, TIME, LOAD, LOSS, PRODUCTION, STORAGE, BACK_UP, NODE \
= ballon_solaire_backup(load_profile=real_load,
loss_profile=real_thermal_losses,
production_p=real_production,
backup_p=backup_measurement,
emin_storage=storage_Emin,
emax_storage=storage_Emax,
einit_storage=storage_Einit,
horizon=horizon,
time_step=optim_time_step,
alpha = alpha,
u_load = u_load)
new_signal_float = (STORAGE.e.value[1]) / rho_water / volume / cp_water - 273.15
is_backup = BACK_UP.p.value[0]
# Save energy flows into a CSV file
save_energy_flows(NODE, file_name=WORK_PATH + r'\results\storage_design')
# # print_results()
#
# test = BACK_UP.p.value
# below_area = STORAGE.below_area.value
# is_load = STORAGE.is_load.value
# back_up_usage = 0.
# below_value = 0.
# is_load_value = 0.
# thr_usage = 0.
#
# for i in range(len(test)):
# back_up_usage += test[i]/6
#
# for i in range(len(test)):
# below_value += below_area[i]/6
#
# for i in range(len(test)):
# is_load_value += is_load[i]/6
#
# print('The back-up consumption value is ' + str(back_up_usage) + ' kWh')
# print('The energy in the water tank is going ' + str(below_value) + ' times below the minimal value of discomfort')
# print('is_load ' + str(is_load_value))
\ No newline at end of file
#! usr/bin/env python3
# coding=utf-8 #
"""
** This Module is an example of storage capacity optimisation.**
This example describes a simple microgrid with :
- A load profile known in advance
- An adjustable production unit, with a maximum power limit
- A storage system with power in charge and power in discharge
The objective consists on minimizing the storage capacity.
"""
from matplotlib import pyplot as plt
from pulp import LpStatus, GUROBI_CMD, COIN_CMD
import pickle
import os
import cProfile
from omegalpes.energy.energy_nodes import EnergyNode
from omegalpes.energy.units.consumption_units import FixedConsumptionUnit
from omegalpes.energy.units.production_units import FixedProductionUnit, VariableProductionUnit, IdentificationVariableProductionUnit
from omegalpes.energy.units.storage_units import StorageUnit, WaterTank, IdentificationWaterTank, InitializationWaterTank
from omegalpes.general.optimisation.model import OptimisationModel
from omegalpes.general.time import TimeUnit
from omegalpes.general.plots import plt, plot_quantity, plot_node_energetic_flows
from omegalpes.general.utils import save_energy_flows
from omegalpes.general.optimisation.elements import Objective
def ballon_solaire_backup(load_profile, loss_profile, production_p, energy_profile,
emax_storage, emin_storage,
horizon, time_step, alpha):
"""
:param production_pmax: Maximal power delivered by the production unit [kW]
:type production_pmax: float or int
:param storage_pcharge_max: Maximal charging power for the storage [kW]
:type storage_pcharge_max: float or int
:param storage_pdischarge_max: Maximal discharging power for the storage [kW]
:type storage_pdischarge_max: float or int
:param load_profile: hourly load profile during a day
:type load_profile: list of 24 float or int
"""
# global model, time, load, loss, production, storage, back_up, node
# Create an empty model
model = OptimisationModel(name='example') # Optimisation model
time = TimeUnit(start='01/01/2018', periods=horizon, dt=time_step/60) # Study on a day, hour
# Create the load - The load profile is known
load = FixedConsumptionUnit(time, 'load', p=load_profile)
# load.u.value = u_load
# Create the losses - modelled as a consumption
loss = FixedConsumptionUnit(time, 'loss', p=loss_profile)
# Create the production unit - The production profile is unknown
production = FixedProductionUnit(time, 'production', p=production_p)
# Create back-up
back_up = VariableProductionUnit(time, 'back_up', pmin = 0, pmax=1.)
# Create the storage
storage = InitializationWaterTank(time, name='water_tank', soc_min = 0, measurement=energy_profile,
soc_max= emax_storage/emax_storage, capacity= 90.512, ef_is_e0 = False)#,
storage.min_energy_abs()
# Create the energy node and connect units
node = EnergyNode(time, 'energy_node')
# Add the energy node to the model
node.connect_units(load, loss, production, storage, back_up) # Connect all units on the same energy node
model.add_nodes(node)
# Optimisation process
model.writeLP('C:/Users/brugeron/Dimosim/OMEGAlpes/examples/optim_models/Initialization_Test.lp')
model.solve_and_update()
return model, time, load, loss, production, storage, back_up, node
def print_results():
"""
*** This function print the optimisation result:
- The optimal capacity
And plot the power curves :
Figure 1 :
- Power consumed by the storage, labelled 'Storage'
- Power consumed by the load, labelled 'Load'
- Power delivered by the production unit, labelled 'Production'
Figure 2 :
- The state of charge of the storage unit
"""
if LpStatus[MODEL.status] == 'Optimal':
print("\n - - - - - OPTIMISATION RESULTS - - - - - ")
# Show the graph
# Power curves
plot_node_energetic_flows(NODE)
# SOC curve
plot_quantity(TIME, STORAGE.e)
plt.xlabel('Time (h)')
plt.ylabel('Energy (kWh)')
plt.title('Storage temperature')
plt.show()
plot_quantity(TIME, BACK_UP.p)
plt.xlabel('Time (h)')
plt.ylabel('Power (kW)')
plt.title('Backup usage')
plt.show()
elif LpStatus[MODEL.status] == 'Infeasible':
print("Sorry, the optimisation problem has no feasible solution !")
elif LpStatus[MODEL.status] == 'Unbounded':
print("The cost function of the optimisation problem is unbounded !")
elif LpStatus[MODEL.status] == 'Undefined':
print("Sorry, a feasible solution has not been found (but may exist). "
"PuLP does not manage to interpret the solver's output, "
"the infeasibility of the MILP problem may have been "
"detected during presolve.")
else:
print("Sorry, the optimisation problem has not been solved.")
if __name__ == '__main__':
# OPTIMIZATION PARAMETERS #
WORK_PATH = os.getcwd()
import pandas as pd
path = 'C:/Users/brugeron/Dimosim/dimosim_on_branch/dimosim/tests/test_external_control'
optim_time_step = 60
database = pd.DataFrame.from_csv(path + '/database.csv', sep=';')
rho_water = 1000 # kg*/m3
cp_water = 4186 / 3600 # [Wh/kg.K]
storage_min_temperature = 57.5 + 273.15 # K
storage_max_temperature = 90 + 273.15 # K
(storage_Emin, storage_Emax) = round(storage_min_temperature * rho_water * cp_water * 0.2 / 1000, 3), \
round(storage_max_temperature * rho_water * cp_water * 0.2 / 1000, 3)
backup_normalized = 97649.0049
availability_normalized = 402
# à modifier à chaque reprise
new_load_signal = 3795
new_production_signal = 0
thermal_loss = 175
residu = 0
path = 'C:/Users/brugeron/Dimosim/OMEGAlpes/examples/results'
d = {}
df = pd.DataFrame(data = d)
df.to_csv(path + '/results.csv', sep=';')
for i in range(int(len(database.values)/24)):
print(i)
j1= i*24
j2 = (i+1)*24
real_load = list(database.values[j1:j2, 0])
real_production = list(database.values[j1:j2, 1])
energy_profile = list(database.values[j1:j2, 2])
real_thermal_losses = list(database.values[j1:j2, 3])
# horizon = round(24*3 * 60/optim_time_step)
horizon = len(real_production)
alpha = 0.7
MODEL, TIME, LOAD, LOSS, PRODUCTION, STORAGE, BACK_UP, NODE \
= ballon_solaire_backup(load_profile=real_load,
loss_profile=real_thermal_losses,
production_p=real_production,
energy_profile=energy_profile,
emin_storage=storage_Emin,
emax_storage=storage_Emax,
horizon=horizon,
time_step=optim_time_step,
alpha = alpha)
optim_result = list(STORAGE.e.value.values())[-1]
measurement = energy_profile[-1]
residu += (optim_result - measurement) ** (2)
residu = (residu / len(database.values)/24) ** (1/2)
print('test')
......@@ -46,29 +46,24 @@ def prediction_unit(measurements, prediction_mode, step, horizon):
if prediction_mode == 'perfect':
len_meas = len(measurements)
for i in range(len_meas):
temp_tup = (measurements[i][step:step + horizon],)
prediction = prediction + temp_tup
predictions = dict()
elif prediction_mode == 'naive':
predicted_solar_production = measurements[1].copy()[step-1:step + horizon-1]
predicted_DHW_load = measurements[0].copy()[step-1:step+horizon-1]
predictions["Solar_Production"] = predicted_solar_production
predictions["DHW_load"] = predicted_DHW_load
len_meas = len(measurements)
for i in range(len_meas):
if i % 2 == 0:
elif prediction_mode == 'naive':
step_data = measurements[i].copy()
step_data = step_data[0:len(measurements[i])-1][::-1]
new_signal_float = measurements[i + 1]
if new_signal_float < 0:
new_signal = 0
else:
new_signal = round(new_signal_float, 3)
predictions = dict()
step_data.append(new_signal)
step_data = step_data[::-1]
temp_tup = (step_data,)
prediction = prediction + temp_tup
predicted_solar_production = measurements[2]["thermal_output_solar"][step - horizon: step].copy()
predicted_solar_production[0] = predicted_solar_production[-1]
predicted_DHW_load = measurements[0]["draw_off_load"][step - horizon: step].copy()
predicted_DHW_load[0] = predicted_DHW_load[-1]
predictions["Solar_Production"] = predicted_solar_production
predictions["DHW_load"] = predicted_DHW_load
return prediction
\ No newline at end of file
return predictions
\ No newline at end of file
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
......@@ -169,10 +169,10 @@ class StorageUnit(VariableEnergyUnit):
name='set_soc_max', parent=self)
# --- Charging and discharging powers ---
self.pc = Quantity(name='pc', opt=True, unit='kW', vlen=time.LEN, lb=0,
self.pc = Quantity(name='pc', opt=True, unit='kW', vlen=time.LEN, lb=pc_min,
ub=pc_max,
parent=self)
self.pd = Quantity(name='pd', opt=True, unit='kW', vlen=time.LEN, lb=0,
self.pd = Quantity(name='pd', opt=True, unit='kW', vlen=time.LEN, lb=pd_min,
ub=pd_max,
parent=self)
......@@ -468,9 +468,9 @@ class WaterTank(StorageUnit):
once at its maximal value during the period Tcycl.
"""
def __init__(self, time, name='WaterTank', pc_min=1e-5, pc_max=1e+5,
pd_min=1e-5, pd_max=1e+5, consumption_unit = None, ref_min = None,
capacity=None, e_0=None, e_f=None, soc_min = None,
def __init__(self, time, name='WaterTank', pc_min=0., pc_max=1e+5,
pd_min=0., pd_max=1e+5, consumption_unit = None, ref_min = None,
capacity=None, e_0=None, e_f=None, soc_min = 0.,
soc_max=None, eff_c=1, eff_d=1, self_disch=0, Tcycl=120,
ef_is_e0=False, owner=None):
"""
......@@ -500,7 +500,7 @@ class WaterTank(StorageUnit):
capacity=capacity, e_0=e_0, e_f=e_f,
soc_min=soc_min, soc_max=soc_max, eff_c=eff_c,
eff_d=eff_d, self_disch=self_disch,
ef_is_e0=ef_is_e0, energy_type='Heat', owner=owner)
ef_is_e0=ef_is_e0, energy_type='Heat')
# DECISION VARIABLES AND PARAMETERS
# --- Checking threshold behavior ---
......@@ -519,11 +519,6 @@ class WaterTank(StorageUnit):
unit='kWh', vlen=time.LEN,
lb=0, ub=capacity, parent=self)
self.is_load = Quantity(name='is_load', opt=True,
description='energy at t in the storage below the minimal value of discomfort',
unit='kWh', vlen=time.LEN,
lb=0, ub=capacity, parent=self)
# CONSTRAINTS
# Defines the threshold evaluation variable
......@@ -557,10 +552,6 @@ class WaterTank(StorageUnit):
'- (1 - {0}_thr[t]) * {0}_capacity'.format(self.name, ref_min),
t_range='for t in time.I', name='def_below_area_one_sup', parent=self)
self.def_is_load = DynamicConstraint(
exp_t='{0}_is_load[t] == {0}_below_area[t] * {1}[t]'.format(self.name, consumption_unit),
t_range='for t in time.I', name='def_is_load', parent=self)
class InitializationWaterTank(StorageUnit):
"""
**Description**
......
......@@ -488,8 +488,8 @@ def compute_gurobi_IIS(gurobi_exe_path=r'C:\gurobi800\win64\bin',
raise ValueError('You should provide either the OptimisationModel '
'or the MPS model')
opt_model.solve()
opt_model.writeMPS('{}.mps'.format(opt_model.name))
MPS_file = glob.glob("*.mps")
# opt_model.writeMPS('{}.mps'.format(opt_model.name))
# MPS_file = glob.glob("*.mps")
else:
if MPS_model[-4:] == '.mps':
......
......@@ -256,35 +256,26 @@ def linearize_square(quantity, nd):
i_nd = list(i_nd_array)
parent.i_nd = Quantity(name=q_name + '_i_nd', value=i_nd, opt=False, vlen=len(i_nd), parent=parent)
'eval({0}[0])[t] + eval({0}[1])[t] + eval({0}[2])[t] + eval({0}[3])[t] + eval({0}[4])[t] + eval({0}[5])[t] == 1'
exp1 = 'lpSum(eval({0}[i])[t] for i in range({1})) == 1'.format(list_u, nd-1)
exp2 = 'lpSum(eval({0})[i][t] for i in range({1})) == 1'.format(list_w[nd], nd+1)
exp3 = 'lpSum(eval({2})[i][t]*{0}_Q[i] for i in range({3})) == {0}_{1}_square_quantity[t]'.format(parent.name, q_name,
list_w, nd)
exp4 = 'lpSum(eval({2})[i][t]*{5}[i] for i in range({3}))== {0}_{4}[t]'.format(parent.name, q_name, list_w, nd,
square_name, Q)
''
def_square_1 = DynamicConstraint(name='def_square_1_' + q_name,
exp_t=exp1,
t_range='for t in time.I',
parent=parent)
def_square_2 = DynamicConstraint(name='def_square_2_' + q_name,
exp_t=exp2,
t_range='for t in time.I',
parent=parent)
def_square_3 = DynamicConstraint(name='def_square_3_' + q_name,
exp_t=exp3,
t_range='for t in time.I',
parent=parent)
def_square_4 = DynamicConstraint(name='def_square_4_' + q_name,
exp_t=exp4,
t_range='for t in time.I',
parent=parent)
setattr(parent, 'def_square_1_' + q_name, def_square_1)
# setattr(parent, 'def_square_2_' + q_name, def_square_2)
# setattr(parent, 'def_square_3_' + q_name, def_square_3)
# setattr(parent, 'def_square_4_' + q_name, def_square_4)
for t in time.I:
exp1 = 'lpSum(eval({0}[i])[{2}] for i in range({1})) == 1'.format(list_u, nd-1, t)
setattr(parent, 'def_square_1_' + q_name + '_' + str(t), Constraint(name='def_square_1_' + q_name + '_' + str(t),
exp=exp1,
parent=parent))
exp2 = 'lpSum(eval({0}[i])[{2}] for i in range({1})) == 1'.format(list_w, nd, t)
setattr(parent, 'def_square_2_'+ q_name + '_' + str(t), Constraint(name='def_square_2_'+ q_name + '_' + str(t),
exp=exp2,
parent=parent))
exp3 = 'lpSum(eval({2}[i])[{4}]*{0}_{1}_Q[i] for i in range({3})) == {0}_{1}[{4}]'.format(parent.name, q_name,
list_w, nd, t)
setattr(parent, 'def_square_3_'+ q_name + '_' + str(t), Constraint(name='def_square_3_'+ q_name + '_' + str(t),
exp=exp3,
parent=parent))
exp4 = 'lpSum(eval({2}[i])[{6}]*{5}[i]*{5}[i] for i in range({3}))== {0}_{4}[{6}]'.format(parent.name, q_name, list_w, nd,
square_name, Q, t)
setattr(parent, 'def_square_4_'+ q_name + '_' + str(t), Constraint(name='def_square_4_'+ q_name + '_' + str(t),
exp=exp4,
parent=parent))
return parent.square_quantity
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!