This optimsation problem is the same as in the deterministic one (see separate notebook), but now we consider uncertainty
We should cater for the likely future wind farm, whilst at the same time avoid investing in standed assets
%matplotlib inline
%load_ext autoreload
%autoreload 2
import ipywidgets
import IPython.display
import pyomo.pysp.util.rapper as rapper
import pyomo.pysp.plugins.csvsolutionwriter as csvw
import pyomo.environ as pyo
import pandas as pd
import copy
import powergama
import powergama.powergim as pgim
import powergama.plots
grid_data = powergama.GridData()
grid_data.readSipData(nodes = "data/dog_nodes.csv",
branches = "data/dog_branches_min.csv",
generators = "data/dog_generators.csv",
consumers = "data/dog_consumers.csv")
# Profiles:
samplesize = 100
grid_data.readProfileData(
filename= "data/timeseries_sample_100_rnd2016.csv",
timerange=range(samplesize), timedelta=1.0)
sip = pgim.SipModel()
dict_data = sip.createModelData(grid_data,
datafile='data/dog_data_irpwind.xml',
maxNewBranchNum=5,maxNewBranchCap=5000)
with open('data/dog_data_irpwind.xml',"r") as file:
parameters_xml = file.read()
Computing B and DA matrices... Creating B and DA coefficients...
@ipywidgets.interact(
data=['','node','branch','generator','consumer','PARAMETERS'])
def showdata(data):
if data=='': print("Select data to display")
elif data=='PARAMETERS': print(parameters_xml)
else: display(getattr(grid_data,data))
interactive(children=(Dropdown(description='data', options=('', 'node', 'branch', 'generator', 'consumer', 'PA…
powergama.plots.plotMap(pg_data=grid_data,pg_res=None,
nodetype='powergim_type',branchtype='type',spread_nodes_r=0.04,zoom_start=7)
This problem has three scenarios. The uncertainty is related to what will happen with wind farm Teeside D that is a planned phase 2 development.
scenario_numberof=3
scenario_probability = [0.50,0.25,0.25]
stm = sip.createScenarioTreeModel(num_scenarios=scenario_numberof,
probabilities=scenario_probability)
def _scenario_data(scenario_name,return_var='dict_data'):
"""Modify data according to scenario"""
global dict_data
global grid_data
# Uncertain wind farm capacity (coming in stage 2)
# deterministic: 1200
# Deterministicv solution connects NO-GB via Teeside C (indxC)
indxD = grid_data.generator[grid_data.generator['desc']=='Teeside D'].index[0]
capD = {'Scenario1': 1200,
'Scenario2': 600,
'Scenario3': 0}
if return_var=='dict_data':
dict_data['powergim']['genCapacity2'][indxD] = capD[scenario_name]
return dict_data
elif return_var=='grid_data':
grid_data2 = copy.deepcopy(grid_data)
grid_data2.generator.loc[indxD,'pmax'] = capD[scenario_name]
return grid_data2
else:
raise Exception("Wrong return_var")
def pysp_instance_creation_callback(xxx,scenario_name, node_names):
'''call-back function to create model instance'''
#print("Creating model instance for scenario={}".format(scenario_name))
#print("Arguments:",xxx,scenario_name,node_names)
dict_data = _scenario_data(scenario_name)
instance = sip.createConcreteModel(dict_data=dict_data)
#instance.write('sto_model_{}.lp'.format(scenario_name))
#instance.pprint('sto_prob_{}.txt'.format(scenario_name))
return instance
# test call-back function:
#inst=pysp_instance_creation_callback('hei','Scenario1', ['root','Scenario1'])
solvername='cbc'
sopts = {}
sopts['--drop-proximal-terms'] = None
sopts['--linearize-nonbinary-penalty-terms'] = 5
sopts = None
phopts = {}
phopts['--output-solver-log'] = None
phopts['--max-iterations'] = '3'
stsolver = rapper.StochSolver(
fsfile=None,fsfct=pysp_instance_creation_callback,
tree_model = stm,
phopts=phopts)
ef_sol = stsolver.solve_ef(subsolver=solvername,sopts=sopts)
print(stsolver.root_E_obj()/1e9)
-5.6002680225
csvw.write_csv_soln(stsolver.scenario_tree, "sto_doggerbank_smaller")
Scenario tree solution written to file=sto_doggerbank_smaller.csv Scenario stage costs written to file=sto_doggerbank_smaller_StageCostDetail.csv
Ref: https://pyomo.readthedocs.io/en/stable/advanced_topics/pysp_rapper/demorapper.html#ph
#ph_sol = stsolver.solve_ph(subsolver=solvername, default_rho=1,phopts=phopts,sopts=sopts)
# With PH, it is important to be careful to distinguish x-bar from x-hat.
#obj, xhat = rapper.xhat_from_ph(ph_sol)
m_opt
dictionary)allres
dictionary)# spread coordinates for better plotting
grid_data.spreadNodeCoordinates(radius=0.04,inplace=True)
m_opt={}
grid_res={}
allres={}
costs=pd.DataFrame()
m_opt['input'] = powergama.plots.plotMap(pg_data=grid_data,pg_res=None,
nodetype='powergim_type',branchtype='powergim_type')
for scen in stsolver.scenario_tree.subproblems:
scen_name=scen.name
inst = scen.instance
# Update data according to scenario (affects stage 2)
grid_data = _scenario_data(scen_name,return_var='grid_data')
for stage in stm.Stages:
res = sip.extractResultingGridData(
grid_data=grid_data,
model=inst,stage=stage,scenario=scen_name)
m = powergama.plots.plotMap(
pg_data=res,res=None,
nodetype='powergim_type',branchtype='powergim_type')
k = '{}_{}'.format(scen_name,stage)
m_opt[k] = m
grid_res[k] = res
# res.writeGridDataToFiles("sto_{}_".format(k))
print('Scenario: {}'.format(k))
# Storing all results in dataframes
for v in inst.component_objects(pyo.Var,active=True):
#print(v,end=": ")
df = pd.DataFrame.from_dict(v.get_values(),orient="index")
df = df.reset_index()
vi = v.index_set()
if hasattr(vi,'_sets'):
cols=[str(k) for k in vi._sets]
else:
cols=[str(v.index_set())]
#print(cols)
df1 = pd.DataFrame(df['index'].tolist(),index=df.index,columns=cols)
df1['value'] = df[0]
allres[v.name] = df1
# Compute costs (investment+operation)
df=allres['{}.generation'.format(scen_name)]
df.columns=['gen','time','stage','value']
generation = [ df[df['stage']==1], df[df['stage']==2]]
cost = pgim.computeSTOcosts(
grid_res[k],dict_data,generation=generation)
costs = pd.concat([costs,pd.DataFrame.from_dict(
cost,orient="index",columns=[scen_name])],axis=1)
Scenario: Scenario1_1 Scenario: Scenario1_2 Warning! powergim.computeSTOcosts is not fully implemented yet. Scenario: Scenario2_1 Scenario: Scenario2_2 Warning! powergim.computeSTOcosts is not fully implemented yet. Scenario: Scenario3_1 Scenario: Scenario3_2 Warning! powergim.computeSTOcosts is not fully implemented yet.
costs
Scenario1 | Scenario2 | Scenario3 | |
---|---|---|---|
(1, invest) | 9.038486e+09 | 9.038486e+09 | 9.038486e+09 |
(1, op) | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 |
(2, invest) | 1.071958e+09 | 9.835806e+08 | 1.497586e+08 |
(2, op) | -1.555478e+10 | -1.559269e+10 | -1.512568e+10 |
costs.sum()
Scenario1 -5.444332e+09 Scenario2 -5.570623e+09 Scenario3 -5.937439e+09 dtype: float64
dfcost=pd.concat([pd.DataFrame.from_dict(inst.opCost.get_values(),orient="index"),
pd.DataFrame.from_dict(inst.investmentCost.get_values(),orient="index")],axis=1)
dfcost.columns=['opex','capex']
dfcost
opex | capex | |
---|---|---|
1 | -9.370928e+08 | 9.037656e+09 |
2 | -1.418859e+10 | 1.497665e+08 |
allres.keys()
dict_keys(['Scenario1.branchNewCapacity', 'Scenario1.branchNewCables', 'Scenario1.newNodes', 'Scenario1.genNewCapacity', 'Scenario1.branchFlow12', 'Scenario1.branchFlow21', 'Scenario1.voltageAngle', 'Scenario1.generation', 'Scenario1.loadShed', 'Scenario1.investmentCost', 'Scenario1.opCost', 'Scenario2.branchNewCapacity', 'Scenario2.branchNewCables', 'Scenario2.newNodes', 'Scenario2.genNewCapacity', 'Scenario2.branchFlow12', 'Scenario2.branchFlow21', 'Scenario2.voltageAngle', 'Scenario2.generation', 'Scenario2.loadShed', 'Scenario2.investmentCost', 'Scenario2.opCost', 'Scenario3.branchNewCapacity', 'Scenario3.branchNewCables', 'Scenario3.newNodes', 'Scenario3.genNewCapacity', 'Scenario3.branchFlow12', 'Scenario3.branchFlow21', 'Scenario3.voltageAngle', 'Scenario3.generation', 'Scenario3.loadShed', 'Scenario3.investmentCost', 'Scenario3.opCost'])
@ipywidgets.interact(x= m_opt.keys())
def plotmap(x):
print("Offshore grid - {}".format(x))
display(m_opt[x])
interactive(children=(Dropdown(description='x', options=('input', 'stage1', 0, 1, 2), value='input'), Output()…
This is useful when wanting to re-processing previous results without running the optimisation again
res_file="sto_doggerbank_smaller.csv"
grid_res = sip.extractResultingGridData(grid_data=grid_data,
file_ph=res_file,stage=1)
grid_res.writeGridDataToFiles("sto_smaller_")
m_opt={}
m_opt['input'] = powergama.plots.plotMap(pg_data=grid_data,pg_res=None,
nodetype='powergim_type',branchtype='type',spread_nodes_r=0.04,zoom_start=7)
m_opt['stage1'] = powergama.plots.plotMap(pg_data=grid_res,pg_res=None,
nodetype='powergim_type',branchtype='type',spread_nodes_r=0.04,zoom_start=7)
for s in range(scenario_numberof):
scen_name = 'Scenario{}'.format(s+1)
dict_data = _scenario_data(scen_name,return_var='grid_data')
grid_res = sip.extractResultingGridData(grid_data=grid_data,
file_ph=res_file,stage=2,
scenario=s+1)
grid_res.writeGridDataToFiles("sto_{}_".format(s+1))
m_opt[s]=powergama.plots.plotMap(pg_data=grid_res,res=None,
nodetype='powergim_type',branchtype='type',spread_nodes_r=0.04,zoom_start=7)
print('Scenario {}'.format(s+1))
Scenario 1 Scenario 2 Scenario 3