1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
|
# Copyright (C) 2018 The Regents of the University of California, Michael Ludkovski and Aditya Maheshwari
# All Rights Reserved
# This code is published under the GNU Lesser General Public License (GNU LGPL)
from __future__ import division
import numpy as np
import math
import sys
import matplotlib
import matplotlib.pyplot as plt
import microgridDEA.parameters as bv
import microgridDEA.simulateState_new as simX
import microgridDEA.valueNext_new as fv
if (sys.version_info > (3, 0)):
import pickle
else:
import cPickle as pickle
import StOptGrids
import StOptReg
import StOptGeners
def currentCost(param,St, dt):
return (param.c1*(dt**param.a))*param.dt, param.c2*np.where(St>0,St,0)*param.dt + param.c3*np.where(St<0,-St,0)*param.dt
def forwardSimulations():
outSampleSimul = 10000
maturityOutsample = 10
param = bv.basicVariables(maturityOutsample)
condExpFilePath = param.filetoDump
with open(param.regParamDump, 'rb') as input:
regParameters = pickle.load(input)
X0 = np.zeros(outSampleSimul)
X0[:] = 0 #initial state for residual demand
Xt = simX.residualDemand(param, param.nstep, outSampleSimul, param.maturity, X0, regParameters[0].rmctype, "outsample")
r,c = Xt.shape
It = np.zeros((r,c)) # inventory evolution
Dt = np.zeros((r,c), dtype=np.int) # regime
Dt[:,0] = 0 #initial regime
St = np.zeros((r,c)) # stores Xt - dt - Bt
Bt = np.zeros((r,c)) # battery output
dt = np.zeros((r,c)) # diesel output
costTrue = np.zeros((r,c)) # cost
It[:,0] = param.I0 #initial state for battery inventory
archiveToRead = StOptGeners.BinaryFileArchive(condExpFilePath,"r")
for iSteps in range(param.nstep):
contValues = archiveToRead.readGridAndRegressedValue(param.nstep-iSteps,"toStore")
for sim in range(outSampleSimul):
_, dt[sim, iSteps], St[sim, iSteps], It[sim, iSteps+1] ,Bt[sim, iSteps] = fv.findOptimalControl(Xt[sim,iSteps], It[sim,iSteps], Dt[sim,iSteps], contValues[0], param, regParameters[iSteps])
Dt[:, iSteps+1] = np.where(dt[:, iSteps] > 0, 1, 0)
costTrue[:,iSteps] = np.where((Dt[:,iSteps]==0) & (Dt[:,iSteps+1]==1),param.K,0)
DieselCost, ArtificalCost = currentCost(param,St, dt)
costTrue = costTrue + DieselCost
costTrue[:,-1] = fv.penalty(param,It[:,-1])
cumCostTrue = np.cumsum(costTrue,1)
print("Cost: Total", np.mean(np.sum(costTrue+ArtificalCost,1)), "true", np.mean(np.sum(costTrue,1)), "artifical", np.mean(np.sum(ArtificalCost,1)))
print("std: Total", np.std(np.sum(costTrue+ArtificalCost,1))/np.sqrt(outSampleSimul), "true", np.std(np.sum(costTrue,1))/np.sqrt(outSampleSimul), "artifical", np.std(np.sum(ArtificalCost,1))/np.sqrt(outSampleSimul))
# plots to show curtailed energy and proportion of blackouts
# St1 = np.where(St<-10**-9,1,0)
# plt.plot(np.sum(St1,0)[:-2]/St1.shape[0])
# plt.show()
# St1 = np.where(St>10**-9,1,0)
# plt.plot(np.sum(St1,0)[:-2]/St1.shape[0])
# plt.show()
if __name__ == "__main__":
forwardSimulations()
|
