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
|
# 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 matplotlib
import matplotlib.pyplot as plt
import microGrid.parameters as bv
import microGrid.simulateState as simX
import microGrid.valueNext as fv
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
X0 = np.zeros(outSampleSimul)
X0[:] = 0 #initial state for residual demand
Xt, _ = simX.residualDemand(param, param.nstep, outSampleSimul, param.maturity, X0, param.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(39,"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)
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()
|
