-
Notifications
You must be signed in to change notification settings - Fork 2
/
optimODS.py
366 lines (335 loc) · 10.9 KB
/
optimODS.py
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
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
#!/usr/bin/env python
# coding: utf-8
# In[1]:
# optimODS.py
# see README.txt file for details.
#----------------------------------
import numpy as np
from scipy.optimize import minimize
# In[2]:
#--------------------------------------------
# Necesary functions
#-----------------------------
def readOD(dir):
# Read all the files required for the problem.
# All files must be located in directory: dir
# files to read:
# namelistOD.txt, file-cf1.txt, file-ref.txt,file-cf2.txt,file-min.txt
# See README.txt for details.
#-------------------------------------------------
# Read file namelistOD.txt
# lines starting with # are comments
# In the first line: Ns,Nw,rs2w,Ssc,Swc,S_Sscmin,Sscmax,Swcmin,Swcmax
# In the second line: SscminR (vector of length Ns)
# In the third line: SscmaxR (vector of length Ns)
# In the fourth line: SwcminR (vector of length Nw)
# In the fifth line: SwcmaxR (vector of length Nw)
global Ns,Nw,rs2w,Ssc,Swc,Sscmin,Sscmax,Swcmin,Swcmax
global SscminR,SscmaxR,SwcminR,SwcmaxR
file_param=dir+'/namelistOD.txt'
with open(file_param,'r') as file:
param={}
icount=-1
for line in file.readlines():
t=line.split()
icount=icount+1
#print(icount,len(t))
if(len(t)==0 or t[0].strip()[0]=='#'):
icount=icount-1
else:
if(icount==0):
Ns=int(t[0].strip())
Nw=int(t[1].strip())
rs2w=float((t[2].strip()))
Ssc=float((t[3].strip()))
Swc=float((t[4].strip()))
Sscmin=float((t[5].strip()))
Sscmax=float((t[6].strip()))
Swcmin=float((t[7].strip()))
Swcmax=float((t[8].strip()))
if(icount==1):
SscminR=[]
for i in range(len(t)):
SscminR=np.append(SscminR,float(t[i].strip()))
if(icount==2):
SscmaxR=[]
for i in range(len(t)):
SscmaxR=np.append(SscmaxR,float(t[i].strip()))
if(icount==3):
SwcminR=[]
for i in range(len(t)):
SwcminR=np.append(SwcminR,float(t[i].strip()))
if(icount==4):
SwcmaxR=[]
for i in range(len(t)):
SwcmaxR=np.append(SwcmaxR,float(t[i].strip()))
if(Ssc <0 or Swc <0 or Ssc>1 or Swc>1 or abs(Ssc+Swc-1)> 0.01):
print('Check Ssc or Swc values')
if(Sscmin>1 or Sscmax<=0 or Swcmin>1 or Swcmax<=0):
print('Check Sscmin, Sscmax, Swcmin, Swcmax')
# put Sscmin/max and Swcmin/max into a single vector
global SminR,SmaxR
SminR=np.append(SscminR,SwcminR)
SmaxR=np.append(SscmaxR,SwcmaxR)
# guarantee that Ssc and Swc add upto 1.
tot=Ssc+Swc
Ssc=Ssc/tot
Swc=Swc/tot
# read the file: file-cf1.txt
# AA : matrix of shape NTTx(Ns+Nw)
# first NS columns correspond to the matrix As and last Nw columns corrspond to Aw.
# See README.txt for details.
global NTT,AA
AA=np.loadtxt(dir+'/file-cf1.txt')
NTT=AA.shape[0] # number of row
l1=AA.shape[1] # number of columns
if(l1 != Ns+Nw):
print('wrong dimensions in file-cf1.txt')
# read the file: file-ref.txt. One column and NTT rows.
# BB : vector of length NTT.
# See README.txt for details.
global BB
BB=np.loadtxt(dir+'/file-ref.txt')
if(BB.shape[0] != NTT):
print('wrong dimensions in file-ref.txt')
# read the file: file-cf2.txt
# CC : matrix of shape 12x(Ns+Nw)
# first NS columns correspond to the matrix As and last Nw columns corrspond to Aw.
# each row correspods to one month
# See README.txt for details.
global NMP,CC
CC=np.loadtxt(dir+'/file-cf2.txt')
NMP=CC.shape[0] # number of row
l1=CC.shape[1] #number of columns
if(l1 != Ns+Nw):
print('wrong dimensions in file-cf2.txt')
# read the file: file-min.txt. One column and NMP rows
# MM : vector of length NMP.
# See README.txt for details.
global MM
MM=np.loadtxt(dir+'/file-min.txt')
if(MM.shape[0] != NMP):
print('wrong dimensions in file-min.txt')
for j in range(NMP):
ccmax=max(CC[j,:])
if(ccmax < MM[j]):
print('wrong construction of file-min.txt')
return
def P_anom(S): # Total production anomalies.
p=np.dot(AA,S)-BB
# print('p:',p.shape,AA.shape,S.shape)
return p
def fun1(S): # Defines the minimizing function. Eq. (1) in README.txt.
p=P_anom(S)
# f=np.var(p) not the same if the mean is not zero.
f=sum(p*p/(len(p)))
return f
def fconstr1(S):
# defines the constrain of minimun capacity of production. Eq. (2) in README.txt.
# See section 1.4 of READM.txt for explanation.
Prod=np.dot(CC,S) # produccion
fcons=0
for i in range(NMP):
if(Prod[i] < MM[i]):
fcons=fcons+(Prod[i]-MM[i])**2
return fcons
def fconstr2(S):
# defines the constrain of minimun/maximun threshold per sub-region. Eq. (3) in README.txt.
fcons=0
for i in range(len(SminR)):
if(S[i] < SminR[i]):
fcons=fcons+(S[i]-SminR[i])**2
if(S[i] > SmaxR[i]):
fcons=fcons+(S[i]-SmaxR[i])**2
return fcons
def fconstr3(S):
# defines the constrain total shares of each technology. Eq. (4) in README.txt.
fcons=(sum(S[0:Ns])-Ssc)**2+(sum(S[Ns:Ns+Nw])-Swc)**2
return fcons
def fconstr4(S):
#defines the constrain of the ratio solar to wind technology: rs2w, eq. (5).
Ss=sum(S[:Ns])
Sw=sum(S[Ns:])
fcons=0
if(rs2w>1 and Ss<Sw):
fcons=fcons+(Ss-Sw)**2
elif(0<rs2w and rs2w<1 and Ss>Sw):
fcons=fcons+(Ss-Sw)**2
return fcons
def fconstr5(S):
#defines the constrain of maximum and minimum threshold for toal share, eq. (6).
Ss=sum(S[:Ns])
Sw=sum(S[Ns:])
fcons=0
if(Ss<Sscmin):
fcons=fcons+(Ss-Sscmin)**2
elif(Ss>Sscmax):
fcons=fcons+(Ss-Sscmax)**2
if(Sw<Swcmin):
fcons=fcons+(Sw-Swcmin)**2
elif(Sw>Swcmax):
fcons=fcons+(Sw-Swcmax)**2
return fcons
def funOD(X,C1,C2,C3):
# defines function to minimize for the OD problem.
Ctot=np.dot(X,X)
S=X*X/Ctot
f0=fun1(S)
fcons1=fconstr1(S)
fcons2=fconstr2(S)
fcons3=fconstr3(S)
f=f0+C1*fcons1+C2*fcons2+C3*fcons3
# print('f:',f0,fcons1,fcons2,fcons3)
return f
def funODS(X,C1,C2,C4,C5):
# defines function to minimize for the OD problem.
Ctot=np.dot(X,X)
S=X*X/Ctot
f0=fun1(S)
fcons1=fconstr1(S)
fcons2=fconstr2(S)
fcons4=fconstr4(S)
fcons5=fconstr5(S)
f=f0+C1*fcons1+C2*fcons2+C4*fcons4+C5*fcons5
# print('f:',f0,fcons1,fcons2,fcons4)
return f
def funtryOD(X):
# used to try to verify possible solutions within the constrains
Ctot=np.dot(X,X)
S=X*X/Ctot
fcons1=fconstr1(S)
fcons2=fconstr2(S)
fcons3=fconstr3(S)
f=fcons1+fcons2+fcons3
# print(f0,fcons1,fcons2,fcons3)
return f
def funtryODS(X):
# used to try to verify possible solutions within the constrains
Ctot=np.dot(X,X)
S=X*X/Ctot
fcons1=fconstr1(S)
fcons2=fconstr2(S)
fcons4=fconstr4(S)
fcons5=fconstr5(S)
f=fcons1+fcons2+fcons4+fcons5
# print(f0,fcons1,fcons2,fcons4,fcons5)
return f
# In[3]:
# Read all data of the OD or ODS problem.
dir='.'
readOD(dir)
print('Ns=',Ns)
print('Nw=',Nw)
print('rs2w=',rs2w)
print('Ssc=',Ssc)
print('Swc=',Swc)
print('Sscmin=',Sscmin)
print('Sscmax=',Sscmax)
print('Swcmin=',Swcmin)
print('Swcmax=',Swcmax)
print('SscminR:',SscminR)
print('SscmaxR:',SscmaxR)
print('SwcminR:',SwcminR)
print('SwcmaxR:',SwcmaxR)
print('NTT=',NTT)
print('NMP=',NMP)
#------------------------------------------------------------------------------------
# Check constrains
# inicializa X
X=np.random.rand(Ns+Nw)
# search for a possible solution of the constrains with tolerance=tol
# If a solution is found it is used as initial guess for the full minimization problem
tol=1.e-8
for i in range(Ns+Nw):
X[i]=np.random.uniform(np.sqrt(max(SminR[i],0)),np.sqrt(min(SmaxR[i],1)))
res=minimize(funtryODS,X,tol=tol)
X=res.x
ff=funtryODS(X)
if(ff > 1.e-11):
print('It seems that constrains are not compatible')
print('____________________________________________')
print('Fisrt Constrain:')
Ctot=np.dot(X,X)
S=X*X/Ctot
Prod=np.dot(CC,S)
for i in range(NMP):
if(Prod[i]<MM[i]):
print(i,Prod[i],'<',MM[i])
print()
print('Second Constrain:')
for i in range(len(SminR)):
if(i<Ns):
Ty='Solar'
ii=i+1
else:
Ty='Wind'
ii=i+1-Ns
if(S[i] < SminR[i]):
print(Ty+' region:',ii,' S[i]=%6.3f < Smin[i]=%6.3f' %(S[i],SminR[i]))
if(S[i] > SmaxR[i]):
print(Ty+' region:',ii,' S[i]=%6.3f > Smax[i]=%6.3f' %(S[i],SmaxR[i]))
print()
print('Third Constrain:')
Ss=sum(S[:Ns])
Sw=sum(S[Ns:])
if(rs2w>1 and Ss<Sw):
print('Ss=',Ss,'should be > Sw=',Sw)
elif(0<rs2w and rs2w<1 and Ss>Sw):
print('Ss=',Ss,'should be < Sw=',Sw)
print()
print('Fourth Constrain:')
if(Ss<Sscmin):
print('Ss=',Ss,'should be > Sscmin=',Sscmin)
elif(Ss>Sscmax):
print('Ss=',Ss,'should be < Sscmax=',Sscmax)
if(Sw<Swcmin):
print('Sw=',Sw,'should be > Swcmin=',Swcmin)
elif(Sw>Swcmax):
print('Sw=',Sw,'should be < Swcmax=',Swcmax)
# In[5]:
#Optimize ODS problem
tol=1.e-6 # tolerance related to the constrain functions.
C1=10.0
C2=100.0
C4=100.0
C5=100
fc1=100
fc2=100
fc4=100
fc5=100
while(fc1>tol or fc2>tol or fc4>tol or fc5>tol):
# print('C:',C1,C2,C4,C5)
# print('f:',fc1,fc2,fc4,fc5)
if(fc1>tol):
C1=2.0*C1
if(fc2>tol):
C2=2.0*C2
if(fc4>tol):
C4=2.0*C4
if(fc5>tol):
C5=2.0*C5
res=minimize(funODS,X,args=(C1,C2,C4,C5))
X=res.x
Ctot=np.dot(X,X)
S=X*X/Ctot
fc1=fconstr1(S)
fc2=fconstr2(S)
fc4=fconstr4(S)
fc5=fconstr5(S)
f0=fun1(S)
print('fmin=%5.3f fcons1=%9.3e fcons1=%9.3e fcons4=%9.3e fcons5=%9.3e' % (f0,fc1,fc2,fc4,fc5))
Ss=sum(S[:Ns])
Sw=sum(S[Ns:])
print(' Sscmin=%5.3f Sum(Ss)=%5.3f Sscmax=%5.3f' %(Sscmin,Ss,Sscmax))
print(' Swcmin=%5.3f Sum(Sw)=%5.3f Swcmax=%5.3f' %(Swcmin,Sw,Swcmax))
print('_____________')
for i in range(Ns):
print(' Ss(%2d)= %5.3f' %(i+1,S[i]))
print('_____________')
for i in range(Nw):
print(' Sw(%2d)= %5.3f' %(i+1,S[i+Ns]))
print('_____________')
Prod=np.dot(CC,S) # produccion
for i in range(NMP):
print(' Prod[%2d]=%8.3f, MM[%2d]=%8.3f' %(i+1,Prod[i],i+1,MM[i]))
# In[51]: