Commit 32880363 authored by Jaber Ibne Mahboob's avatar Jaber Ibne Mahboob
Browse files

Only availability of wind

Only availability of wind is integrated to HyDesign tool
parent 071ee3bf
......@@ -3,6 +3,7 @@ from __future__ import print_function, division, absolute_import, unicode_litera
# Copyright (C) 2019 DTU Wind
import os
from termcolor import cprint
import logging
import sys
......@@ -256,7 +257,130 @@ class HPP(object):
wind_ts, solar_ts, price_ts)
return hpp_wind_capacity, hpp_solar_capacity, hpp_battery_power_rating, hpp_battery_energy_capacity, P_RES_available_t, P_HPP_t, P_curtailment_t, P_charge_discharge_t, E_SOC_t, hpp_investment_cost, hpp_maintenance_cost, LCOE, NPV, IRR
def sizing_Wind_Solar(self, wind_ts, wt_av_ts, solar_ts, pv_av_ts, price_ts):
def sizing_Wind_Solar(self, wind_ts, solar_ts, price_ts):
"""
Method to calculate sizing of wind and solar
Returns
-------
Capacity of Wind Power
Capacity of Solar Power
HPP power output timeseries
HPP power curtailment timeseries
HPP total CAPEX
HPP total OPEX
Levelised cost of energy
"""
# extract parameters into the variable space
globals().update(self.__dict__)
mdl = Model(name='Sizing')
time = price_ts.index
# Variables definition
P_HPP_t = mdl.continuous_var_dict(
time, lb=0, ub=hpp_grid_connection,
name='HPP power output')
P_curtailment_t = mdl.continuous_var_dict(
time, lb=0, name='Curtailment')
Wind_MW = mdl.integer_var(lb=0, name='Wind capacity')
Solar_MW = mdl.integer_var(lb=0, name='Solar capacity')
for t in time:
# Constraints
#mdl.add_constraint(((turbine_spacing * rotor_diameter * 0.001 * (nWT_per_string-1)) * (turbine_row_spacing * rotor_diameter * 0.001 * ((X*(1/(rating_WT*nWT_per_string)))-1)) ) + (land_SOL*Y) <= land_area_available)
mdl.add_constraint(
P_curtailment_t[t] >= wind_ts[t] * Wind_MW +
solar_ts[t] * Solar_MW -
hpp_grid_connection)
mdl.add_constraint(
P_HPP_t[t] == wind_ts[t] * Wind_MW +
solar_ts[t] * Solar_MW -
P_curtailment_t[t])
# Objective function
mdl.maximize(
# E_u(NPV)
#mdl.sum(NPV[u]*weights[u] for u in unc)
# Q_u(NPV, q) approximate by a N
#mdl.sum(NPV[u]**N *weights[u] for u in unc)**(1/N)
#
#E_u(NPV) - 2*S_u(NPV)
# CAPEX
-(wind_turbine_cost * Wind_MW + \
solar_PV_cost * Solar_MW + \
hpp_BOS_soft_cost * (Wind_MW + Solar_MW) + \
hpp_grid_connection_cost * hpp_grid_connection + \
solar_hardware_installation_cost * Solar_MW + \
wind_civil_works_cost * Wind_MW) + \
# revenues and OPEX
mdl.sum(
(mdl.sum(
price_ts[t] * (wind_ts[t] * Wind_MW + \
solar_ts[t] * Solar_MW - \
P_curtailment_t[t])
for t in time) + \
- (wind_fixed_onm_cost * Wind_MW + \
solar_fixed_onm_cost * Solar_MW)
) / np.power(1 + hpp_discount_factor, i)
for i in range(1, hpp_lifetime + 1)
)
)
# Solving the problem
sol = mdl.solve(log_output=False)
# print(mdl.export_to_string())
# sol.display()
P_curtailment_ts = pd.DataFrame.from_dict(
sol.get_value_dict(P_curtailment_t), orient='index')
P_HPP_ts = pd.DataFrame.from_dict(
sol.get_value_dict(P_HPP_t), orient='index')
AEP = P_HPP_ts[0].sum()
AEP_per_year = np.ones(hpp_lifetime) * AEP
# Investment cost
hpp_investment_cost = \
wind_turbine_cost * sol.get_value(Wind_MW) + \
solar_PV_cost * sol.get_value(Solar_MW) + \
hpp_grid_connection_cost * hpp_grid_connection + \
hpp_BOS_soft_cost * (sol.get_value(Wind_MW) + sol.get_value(Solar_MW)) + \
wind_civil_works_cost * sol.get_value(Wind_MW) + \
solar_hardware_installation_cost * sol.get_value(Solar_MW)
# Maintenance cost
hpp_maintenance_cost = np.ones(hpp_lifetime) * (wind_fixed_onm_cost * sol.get_value(
Wind_MW) + solar_fixed_onm_cost * sol.get_value(Solar_MW))
LCOE = self.calculate_LCOE(
hpp_investment_cost,
hpp_maintenance_cost,
hpp_discount_factor,
hpp_lifetime,
AEP_per_year)
P_HPP_ts_2 = pd.DataFrame.from_dict(
sol.get_value_dict(P_HPP_t), orient='index')
IRR = self.calculate_IRR(
P_HPP_ts_2,
price_ts,
hpp_investment_cost,
hpp_maintenance_cost)
NPV = self.calculate_NPV(
P_HPP_ts_2,
price_ts,
hpp_investment_cost,
hpp_maintenance_cost,
hpp_discount_factor)
return sol.get_value(Wind_MW), sol.get_value(
Solar_MW), P_HPP_ts[0], P_curtailment_ts[0], hpp_investment_cost, hpp_maintenance_cost, LCOE, NPV, IRR
def sizing_Wind_Solar_av(self, wind_ts, wt_av_ts, solar_ts, pv_av_ts, price_ts):
"""
Method to calculate sizing of wind and solar
......@@ -408,6 +532,265 @@ class HPP(object):
hpp_discount_factor)
return sol.get_value(Wind_MW), sol.get_value(
Solar_MW), P_HPP_ts[0], P_curtailment_ts[0], hpp_investment_cost, hpp_maintenance_cost, LCOE, NPV, IRR
def sizing_Wind_Solar_wtav2(self, wind_ts, wt_av_ts, solar_ts, price_ts):
def sizing_Wind_Solar_wtav(self, wind_ts, wt_av_ts, solar_ts, price_ts):
"""
Method to calculate sizing of wind and solar
Returns
-------
Capacity of Wind Power
Capacity of Solar Power
HPP power output timeseries
HPP power curtailment timeseries
HPP total CAPEX
HPP total OPEX
Levelised cost of energy
"""
# extract parameters into the variable space
globals().update(self.__dict__)
iIter = 0
nWT = -1 # give a random value # nWT = random.randint(1, len(wt_av_ts.columns)) # this will give wrong result if by random number is same as the number found from Alessandra's solution
nWT_last = -2 # make sure, it's not same value of nWT
nWT_diff_last = -1
count_mdl_over_damping = 0
toler_mdl_over_damping = 3
err_mdl_over_damping = False
count_mdl_critically_damping = 0
toler_mdl_critically_damping = -1
mdl_critically_damping_nWT = -1
err_mdl_critically_damping = False
NPV_last = -1
LCOE_last = -1
AEPi_last = -1
# print(f'initial number of wind turbines = {nWT}')
while ( nWT != nWT_last ):
iIter = iIter + 1
#if nWT>len(wt_av_ts.columns) or nPV_unit>len(pv_av_ts.columns):
# ********* call the reliability model with more components **********
# for the first iter, consider 100% availability, which will then become equivalent of Alessandra's model
if iIter == 1: wind_av_ts = (wt_av_ts.iloc[:, 0] * 0) + 1
elif nWT > 0: wind_av_ts = wt_av_ts.iloc[:, 0:nWT].sum(axis=1)/nWT
else: wind_av_ts = wt_av_ts.iloc[:, 0] * 0
# else: wind_av_ts = (wt_av_ts.iloc[:, 0] * 0) + 1 # nWT = 0, so availability 1 or, 0 doesn't effect the model
mdl = Model(name='Sizing')
time = price_ts.index
# Variables definition
P_HPP_t = mdl.continuous_var_dict(
time, lb=0, ub=hpp_grid_connection,
name='HPP power output')
P_curtailment_t = mdl.continuous_var_dict(
time, lb=0, name='Curtailment')
Wind_MW = mdl.integer_var(lb=0, name='Wind capacity') * wind_rating_WT
Solar_MW = mdl.integer_var(lb=0, name='Solar capacity') * solar_rating_PV_unit
for t in time:
# Constraints
#mdl.add_constraint(((turbine_spacing * rotor_diameter * 0.001 * (nWT_per_string-1)) * (turbine_row_spacing * rotor_diameter * 0.001 * ((X*(1/(rating_WT*nWT_per_string)))-1)) ) + (land_SOL*Y) <= land_area_available)
mdl.add_constraint(
P_curtailment_t[t] >= wind_ts[t] * wind_av_ts[t] * Wind_MW +
solar_ts[t] * Solar_MW -
hpp_grid_connection)
mdl.add_constraint(
P_HPP_t[t] == wind_ts[t] * wind_av_ts[t] * Wind_MW +
solar_ts[t] * Solar_MW -
P_curtailment_t[t])
# Objective function
mdl.maximize(
# E_u(NPV)
#mdl.sum(NPV[u]*weights[u] for u in unc)
# Q_u(NPV, q) approximate by a N
#mdl.sum(NPV[u]**N *weights[u] for u in unc)**(1/N)
#
#E_u(NPV) - 2*S_u(NPV)
# CAPEX
-(wind_turbine_cost * Wind_MW + \
solar_PV_cost * Solar_MW + \
hpp_BOS_soft_cost * (Wind_MW + Solar_MW) + \
hpp_grid_connection_cost * hpp_grid_connection + \
solar_hardware_installation_cost * Solar_MW + \
wind_civil_works_cost * Wind_MW) + \
# revenues and OPEX
mdl.sum(
(mdl.sum(
price_ts[t] * (wind_ts[t] * wind_av_ts[t] * Wind_MW + \
solar_ts[t] * Solar_MW - \
P_curtailment_t[t])
for t in time) + \
- (wind_fixed_onm_cost * Wind_MW + \
solar_fixed_onm_cost * Solar_MW)
) / np.power(1 + hpp_discount_factor, i)
for i in range(1, hpp_lifetime + 1)
)
)
# Solving the problem
sol = mdl.solve(log_output=False)
Wind_MW_i = sol.get_value(Wind_MW)
Solar_MW_i = sol.get_value(Solar_MW)
nWT_last = nWT
nWT = math.ceil(Wind_MW_i/wind_rating_WT)
nPV_unit = math.ceil(Solar_MW_i/solar_rating_PV_unit)
# print(mdl.export_to_string())
# sol.display()
availability = wind_av_ts.mean(axis=0)
if iIter == 1: print(f"Iteration =====> {iIter} (Alessandra's solution)")
else: print(f'Iteration =====> {iIter} (Used availability calculated for {nWT_last} WTs)')
print(pd.DataFrame([[nWT, nPV_unit, availability]], columns = ['nWT', 'nPV_unit','Availability']).to_string(index=False))
"""
Here is the twist:
This through the while loop, the model tries to reach a point of solution
In a single word - it tries to under damped to the solution point
Now, availability is an important factor for HPP size optimization
If availability become very sensetive, i.e. if the number of WTs reduced to 50 or so
- overall availability fluctuate too much to cause the model is getting Over Damped from solution point
- or, may be get stuck in a loop that keep running without Over damping or, under damping towards solution point
If such situation happens, we need to figure it out and find a solution for that case
- 1st attempt: track the difference of nWT & nWT_last,
if the difference is not same, no need to evaluate any parameter, re-run the loop
elif difference becomes constant, means no Over damping or, conerging
- if the difference > 1: take nWT = avg(nWT & nWT_last)
- else (the difference == 1): track NPV of both solution and find which is lowest
if NPV is same, track LCOE of both solution, check which is higher
if LCOE is same, track AEP of both solution, check which is higher
if AEP is same, take the lower nWT solution as the final solution as it will give lower investment cost but same output
"""
"""
Well besed on the discussion of the above, now we need to evaluate all the parameters in each loop :(
"""
P_curtailment_ts = pd.DataFrame.from_dict(
sol.get_value_dict(P_curtailment_t), orient='index')
P_HPP_ts = pd.DataFrame.from_dict(
sol.get_value_dict(P_HPP_t), orient='index')
AEP = P_HPP_ts[0].sum()
AEP_per_year = np.ones(hpp_lifetime) * AEP
# Investment cost
hpp_investment_cost = \
wind_turbine_cost * sol.get_value(Wind_MW) + \
solar_PV_cost * sol.get_value(Solar_MW) + \
hpp_grid_connection_cost * hpp_grid_connection + \
hpp_BOS_soft_cost * (sol.get_value(Wind_MW) + sol.get_value(Solar_MW)) + \
wind_civil_works_cost * sol.get_value(Wind_MW) + \
solar_hardware_installation_cost * sol.get_value(Solar_MW)
# Maintenance cost
hpp_maintenance_cost = np.ones(hpp_lifetime) * (wind_fixed_onm_cost * sol.get_value(
Wind_MW) + solar_fixed_onm_cost * sol.get_value(Solar_MW))
LCOE = self.calculate_LCOE(
hpp_investment_cost,
hpp_maintenance_cost,
hpp_discount_factor,
hpp_lifetime,
AEP_per_year)
P_HPP_ts_2 = pd.DataFrame.from_dict(
sol.get_value_dict(P_HPP_t), orient='index')
IRR = self.calculate_IRR(
P_HPP_ts_2,
price_ts,
hpp_investment_cost,
hpp_maintenance_cost)
NPV = self.calculate_NPV(
P_HPP_ts_2,
price_ts,
hpp_investment_cost,
hpp_maintenance_cost,
hpp_discount_factor)
"""
Now start checking *****
"""
nWT_diff = abs((nWT - nWT_last))
AEPi = P_HPP_ts[0].mean() * 8760 / 1000
if (toler_mdl_critically_damping != -1) or ((nWT_diff == nWT_diff_last) and (nWT_diff > 1)):
##### MODEL IS NOT RUNNING EITHER over damped OR under damped WAY (INFINITE LOOP)#####
##### MODEL IS RUNNING IN A critically_damping WAY #####
if (toler_mdl_critically_damping == -1): # initial value is -1
toler_mdl_critically_damping = abs((nWT - nWT_last)) - 1
if nWT < nWT_last:
mdl_critically_damping_nWT = nWT
else:
mdl_critically_damping_nWT = nWT_last
if (count_mdl_critically_damping >= toler_mdl_critically_damping) and (toler_mdl_critically_damping != -1):
##### NO SOLUTION WILL BE FOUND ####
err_mdl_critically_damping = True
print(f'nWT_diff: {nWT_diff}')
print(f'nWT_diff_last: {nWT_diff_last}')
cprint("**MODEL critically damping Error**: Due to small number of WTs, the fluctuation of the value of availability of WTs increased to a range that makes the Model running in a critically_damping way and no solution is found. Need a better model to fix this issue.","red","on_grey",attrs=['underline'])
break #break the model
else:
if nWT == nWT_last :
continue # solution found :)
else:
count_mdl_critically_damping = count_mdl_critically_damping + 1
nWT = mdl_critically_damping_nWT + count_mdl_critically_damping # value of count_mdl_critically_damping is increament from 1 to tolar_mdl_critically_damping
print(f'mdl critically damping tolerance -> {toler_mdl_critically_damping} iterations')
print(f'mdl critically damping count -> {count_mdl_critically_damping}')
print(f'nWT -> {nWT}')
continue # go for next loop directly to find a way out of this infinit loop
elif (nWT_diff == nWT_diff_last) and (nWT_diff == 1):
##### MODEL IS NOT RUNNING EITHER over damped OR under damped WAY (INFINITE LOOP)#####
##### SOLUTION SHOULD BE LIED IN BETWEEN nWT & nWT_solution #####
if NPV > NPV_last:
break # we already got the optimal solution in this loop
elif NPV == NPV_last:
if LCOE < LCOE_last:
break # we already got the optimal solution in this loop
elif LCOE == LCOE_last:
if AEPi > AEPi_last:
break # we already got the optimal solution in this loop
elif AEPi == AEPi_last:
if nWT < nWT_last:
break # we already got the optimal solution in this loop
elif (nWT_diff > nWT_diff_last) and (iIter != 1): # don't do this in first loop
##### MODEL TRYING TO FIND A SOLUTION over damped WAY (INFINITE LOOP)#####
if count_mdl_over_damping >= toler_mdl_over_damping:
##### NO SOLUTION WILL BE FOUND ####
err_mdl_over_damping = True
print(f'nWT_diff: {nWT_diff}')
print(f'nWT_diff_last: {nWT_diff_last}')
cprint("**MODEL Over damping Error**: Due to small number of WTs, the fluctuation of the value of availability of WTs increased to a range that makes the Model running in a over damped way and no solution is found. Need a better model to fix this issue.","red","on_grey",attrs=['underline'])
break #break the model
else:
count_mdl_over_damping = count_mdl_over_damping + 1 # lets try again
print(f'mdl Over damping count -> {count_mdl_over_damping}')
nWT_diff_last = nWT_diff
NPV_last = NPV
LCOE_last = LCOE
AEPi_last = AEPi
if err_mdl_over_damping == True or err_mdl_critically_damping == True:
return 0, 0, (P_HPP_ts[0] * 0), (P_curtailment_ts[0] * 0), 0, (hpp_maintenance_cost * 0), 0, 0, 0 # failed to generate a result
else:
return sol.get_value(Wind_MW), sol.get_value(
Solar_MW), P_HPP_ts[0], P_curtailment_ts[0], hpp_investment_cost, hpp_maintenance_cost, LCOE, NPV, IRR
def sizing_Wind_Solar_Battery(self, wind_ts, solar_ts, price_ts):
"""
......
......@@ -8,7 +8,6 @@ Created on Sun Nov 10 12:05:52 2019
# -*- coding: utf-8 -*-
import pandas as pd
import numpy as np
from HyDesign import HPP, read_csv_Data
......@@ -54,9 +53,9 @@ if __name__ == "__main__":
'wind_hub_height': 120, # in m
'wind_turbine_spacing': 5, # in terms of Rotor Diameter
'wind_turbine_row_spacing': 10, # in terms of Rotor Diameter
'wind_turbine_cost': 851000, # [EUR/MW]
'wind_civil_works_cost': 116986, # [Eur/MW]
'wind_fixed_onm_cost': 12805, # Wind fixed O&M cost per year [Eur/MW]
'wind_turbine_cost': round(851000 * 1.00, 0), # [EUR/MW]
'wind_civil_works_cost': round(116986 * 1.00, 0), # [Eur/MW]
'wind_fixed_onm_cost': round(12805 * 1.00, 0), # Wind fixed O&M cost per year [Eur/MW]
# hpp solar parameters
'solar_rating_PV_unit': 1, # [MW]
......@@ -138,11 +137,11 @@ if __name__ == "__main__":
# --------------------------------------------------
print()
print("HPP size optimization for a wind and solar with availability")
print("HPP size optimization by Alessandra's model:")
[hpp_wind_capacity, hpp_solar_capacity, P_HPP_t,
P_curtailment_t, hpp_investment_cost, hpp_maintenance_cost,
LCOE, NPV, IRR] = ExampleHPP.sizing_Wind_Solar(
wind_power_t, wind_av_t, solar_power_t, solar_av_t, spot_price_t)
P_curtailment_t, hpp_investment_cost, hpp_maintenance_cost,
LCOE, NPV, IRR] = ExampleHPP.sizing_Wind_Solar(
wind_power_t, solar_power_t, spot_price_t)
# %%
......@@ -165,12 +164,11 @@ if __name__ == "__main__":
# --------------------------------------------------
print()
print("Evaluation of a non-optimal HPP for wind and solar with availability")
print("HPP size optimization with availability factor of Wind turbines:")
[hpp_wind_capacity, hpp_solar_capacity, P_HPP_t,
P_curtailment_t, hpp_investment_cost, hpp_maintenance_cost,
LCOE, NPV, IRR] = ExampleHPP.eval_Wind_Solar(
wind_power_t, solar_power_t, spot_price_t,
nWT=80, wt_av_ts=wind_av_t, nPV_unit=300, pv_av_ts=solar_av_t)
LCOE, NPV, IRR] = ExampleHPP.sizing_Wind_Solar_wtav(
wind_power_t, wind_av_t, solar_power_t, spot_price_t)
# %%
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment