-
Notifications
You must be signed in to change notification settings - Fork 1
/
Cascata_sDemanda.jl
494 lines (405 loc) · 21.9 KB
/
Cascata_sDemanda.jl
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
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
#o if a seguir verifica se o Julia está na pasta correta para buscar os arquivos necessários
if ~occursin("Finardi", pwd())
try
cd("Dropbox/LAC/Alocação de unidades/Curtissimo Cascata Finardi") #direciona o Julia para a pasta correta
catch
nothing
end
end
#os if's a seguir verificam a existência dos arquivos necessários para rodar o modelo
if ~isfile("Dados/aflu_h.csv")
include("aflu_h.jl")
end
if ~isfile("Dados/entrada.csv") || ~isfile("Dados/potencia_instalada_hidr.csv")
include("hidr.jl")
end
if ~isfile("Dados/Dessem/deflu_h.csv")
include("deflu_h.jl")
end
if ~isfile("Dados/Dessem/demanda_sub_h.csv") || ~isfile("Dados/Dessem/tempo_viagem.csv")
include("Entrada_dados.jl")
end
if ~isfile("Dados/demanda_demanda.csv") || ~isfile("Dados/demanda_hidr.csv") #precisa ser o último dos if's
include("entrada.jl")
end
#---------------------
using JuMP, Juniper, Ipopt, DataFrames
using CSV
using Dates
import XLSX
using Interpolations
using Plots
using LinearAlgebra
pyplot() #utilizando o pyplot
dia, mes, ano = 18,10,2019 #inserindo o dia para obtenção dos dados
#t1 = now()
optimizer = Juniper.Optimizer
nl_solver = optimizer_with_attributes(Ipopt.Optimizer, "print_level" => 0)
model = Model(optimizer_with_attributes(optimizer, "nl_solver" => nl_solver))#,"mip_solver"=>mip_solver))
grau = "2" #grau para coeficientes da curva colina
coefs_cc = CSV.File(string("Dados/Curva Colina/coefs_",grau,"grau.csv")) |> Tables.matrix
Hidr = CSV.File(string("Dados/entrada.csv")) |> DataFrame!
DefAnt = CSV.File(string("Dados/Dessem/deflu_h.csv")) |> DataFrame! #defluência anterior
#Dem_sub = CSV.File(string("Dados/Dessem/demanda_sub_h.csv")) |> DataFrame! #demanda
dem_dem = CSV.File(string("Dados/demanda_demanda.csv")) |> DataFrame! #demanda tipo demanda
dem_hidr = CSV.File(string("Dados/demanda_hidr.csv")) |> DataFrame! #demanda tipo hidr
temp_viag = CSV.File(string("Dados/Dessem/tempo_viagem.csv")) |> DataFrame! #tempo de viagem
Aflu = CSV.File(string("Dados/aflu_h.csv")) |> DataFrame!
cod_usinas = collect(Iterators.flatten((71:74, 76:77))) #códigos da ONS das usinas presentes neste código
usinas = length(cod_usinas)
DefAnt_usinas = DefAnt[findall(x -> x in cod_usinas, DefAnt[!,:Jus]), 1:end]
temp_viag_usinas = temp_viag[findall(x -> x in cod_usinas, temp_viag[!,:Jus]), 1:end]
# 71 - Pot, liquida # Sta Clara PR
# 72 - cc Constante/Fio d'agua # Fundão
# 73 - N temos cc/Reservatorio # Jordão
# 74 - Pot, liquida # Foz do Areia/G.B. Munhoz
# 76 - Vaz, liquida # Segredo
# 77 - Pot, bruta # S. Santiago
usinas_geradoras = collect(Iterators.flatten((1:4, 6:usinas))) #pega todas usinas que geram algo, retira as que são somente reservatório
#usi_mon_natural = [[-1] [71] [72] [-1] [74] [73 76]] #usinas a montante no rio
#usi_mon_natural = Array{Any}(ones(6))#cod. usinas a montante no rio
#usi_mon_natural[6] = [3 5]
#usi_mon_natural[1:5] = [[0] [1] [2] [0] [4]]
##-------
#VARIÁVEIS MANUAIS
temp = 3 #períodos que o modelo irá rodar
tipo_aflu = "Media Movel" # = "Media Movel Scaled" #variável que seleciona entre média móvel ou média móvel escalonada
tipo_demanda = "hidr" # = dem #define o timpo de demanda utilizada durante o código
##-------
Vaz_aflu = Aflu[Aflu[!,:x3] .== tipo_aflu, 1:end]
v_vert = Hidr[!,"v_vert"] #vol min necessário para verter
v_min = Hidr[!,"v_min"] #
v_max = Hidr[!,"v_max"] #
H = 60.0 * 60.0 #contante de transformação de vazão [m³/s] para volume [m³]
n_unidades = Hidr[!,"n_unidades"]
#L é um dataframe ou seja pra pegar a demanda horária de uma usina com código X (Inteiro) se usa L[!,string(X)] ou L[!,Symbol(X)]
if tipo_demanda == "hidr"
L = dem_hidr
else
L = dem_dem
end
tipo_perda = Hidr[!,"tipo_perda"] #tipo de perda para altura líquida
perda_valor = Hidr[!,"perda_valor"] #valor utilizado na conta da altura líquida
w_min = Hidr[!,"w_min"]#limites da vazão turbinada
w_max = Hidr[!,"w_max"]#limites da vazão turbinada
w_min[3] = 10 #jordão turbina 10 por todo momento para não secar o resto do rio
w_max[3] = 10 #jordão turbina 10 por todo momento para não secar o resto do rio
g_min = Hidr[!,"g_min"]#limites da geração
g_max = Hidr[!,"g_max"]#limites da geração
I_max = ones(6)*3 #CHUTE
T_up = ones(6)*6 #tempo mínimo da máquina ligada CHUTE
T_down = ones(6)*2 #tempo mínimo da máquina desligada CHUTE
jus_0 = Hidr[!,"jus_0"]
jus_1 = Hidr[!,"jus_1"]
jus_2 = Hidr[!,"jus_2"]
jus_3 = Hidr[!,"jus_3"]
jus_4 = Hidr[!,"jus_4"]
mon_0 = Hidr[!,"mon_0"]
mon_1 = Hidr[!,"mon_1"]
mon_2 = Hidr[!,"mon_2"]
mon_3 = Hidr[!,"mon_3"]
mon_4 = Hidr[!,"mon_4"]
#---------------------
#max_h = Hidr[!,"max_h"]
#min_h = Hidr[!,"min_h"]
rendimento = Hidr[!,"rendi_nominal"] #"RENDIMENTO NOMINAL'
periodos_aflu = size(Aflu)[2] - 3 #tiro as 2 colunas de texto
Y = zeros(usinas, periodos_aflu)
O = zeros(usinas)
indice_usinas_montante = Dict() #dicionário que entra com o código da usina e retorna o índice i/k dela
δ = zeros(usinas)
for i in 1:usinas
cod_aux = cod_usinas[i]
if cod_aux in temp_viag_usinas[!,:Jus]
δ[i] = temp_viag_usinas[temp_viag_usinas[!,:Jus] .== cod_aux,:hora][1]
end
posic_usina = findall(x -> x == cod_aux, Vaz_aflu[!,:x2])[1]
Y[i, 1:end] = Array(Vaz_aflu[posic_usina, 4:end])
O[i] = sum(DefAnt_usinas[DefAnt_usinas[!,:Jus] .== cod_aux, :defluencia]) #vazões defluêntes anteriores
get!(indice_usinas_montante, cod_aux, findall(x -> x in DefAnt_usinas[findall(x -> x == cod_aux, DefAnt_usinas[!,:Jus]),:Mont], cod_usinas))
end
#@variable(model, 1 ≥ x[ativo = 1:6] ≥ 0)
#OBS: Todas variáveis que contem n_unidades[5]*, não existem na entrada 5, logo o código não rodaria
@variable(model, q[i =1:usinas, t = 1:temp] ≥ 0)#, start = w_max[1]*n_unidades[1])#vazão turbinada total usina
@variable(model, s[i =1:usinas, t = 1:temp] ≥ 0)#, start = 0)#vazão vertida
@variable(model, v_max[i] ≥ v[i = 1:usinas, t = 0:temp] ≥ v_min[i])#, start = v_max[1])#volume
@variable(model, g[i =1:usinas, t = 1:temp, j = 1:n_unidades[i]]≥ 0)#, start = g_max[i]) #geração por turbina
@variable(model, hb[i=1:usinas, t = 1:temp] ≥ 0) #altura bruta, no código dá pra usar hl já que ele faz a transformação para hl ou deidxa em hb dependendo da cc
@variable(model, 1 ≥ η[i=1:usinas,t = 1:temp, j = 1:n_unidades[i]]≥ 0) #rendimento
@variable(model, w_max[i] ≥ w[i =1:usinas, t = 1:temp, j = 1:n_unidades[i]]≥ 0) #,start = w_max[i]) #vazão da turbina
@variable(model, u[i =1:usinas, t = 0:temp, j = 1:n_unidades[i]], Bin) #estado on off da turbina
@variable(model, ç[i =1:usinas, t = 0:temp, j = 1:n_unidades[i]], Bin) #ação de ligar
@variable(model, o[i =1:usinas, t = 0:temp, j = 1:n_unidades[i]], Bin) #ação de desligar
@variable(model, def[i=1:usinas, t = 1:temp] ≥ 0) #deficit
@variable(model, vaz_t[i = 1:usinas, t = 1:temp], start = 0.) #vazão transferida, positiva se a usina recebeu água, negativa se a usina passou água
#sta clara e foz do areia previsão de afluência por todo período
#fundão previsão de 2h, I_aflu = [10, 20, q+s[sta clara], ...]
#jordão previsão de 1h, I_aflu = [15, q+s[fundão], ...]
@NLexpression(model, mon[i = 1:usinas,
t = 1:temp],
mon_4[i]*(v[i,t])^4 + mon_3[i]*(v[i,t])^3 +mon_2[i]*(v[i,t])^2 + mon_1[i]*(v[i,t]) +mon_0[i])#se der estranho pode ser aqui JORDAO
@NLexpression(model, jus[i = 1:usinas,
t = 1:temp],
jus_4[i]*(q[i,t])^4 + jus_3[i]*(q[i,t])^3 +jus_2[i]*(q[i,t])^2 + jus_1[i]*(q[i,t]) +jus_0[i]) #se der estranho pode ser aqui JORDAO
#Restrição diz exatamente
#se i ∈ [1,4,5] então hl = (hb - perda_valor) ou (1-perda_valor/100)*hb
#caso contrário (i ∈ [2,3,6]) então hl = hb
@expression(model, hl[i = 1:usinas,
t = 1:temp],
i in collect(Iterators.flatten((1:1, 4:5))) ? (2 - tipo_perda[i])*(1 - (perda_valor[i]/100.0))*hb[i,t] + (tipo_perda[i] - 1)*(hb[i,t] - perda_valor[i]) : hb[i,t])
@expression(model, pot[i=collect(Iterators.flatten((1, 4, 6))),
t = 1:temp,
j = 1:n_unidades[i]], #as cc da usina 1,4 e 6 são em pot
9.80665*1e-3*rendimento[i]*hl[i,t]*w[i,t,j]) #usa a perda 1 se tipo perda == 1, análogo para 2
#Expressão que retorna a vazão transferida de Jordão para Segredo
#Se positiva -> Jordão recebe água de Segredo
#Se negativa -> Jordão envia água para Segredo
#Explicação: O túnel que liga as duas usinas funciona através da física dos vasos comunicantes
#@NLexpression(model, jord_segr[t = 1:temp], #vazão transferida de Jorão para Segredo
# if mon[3,t] - mon[5,t] > 0
# -72*sqrt(mon[3,t] - mon[5,t])
# else
# 72*sqrt(-mon[3,t] + mon[5,t])
# end
# )
#Hidr[[3,6],[:max_h, :min_h, :v_min, :v_max]]
#collect(Iterators.flatten((1:4, 6:usinas)))
@expression(model,cc[i = 1:usinas,
t = 1:temp,
j = 1:n_unidades[i]],#[1,P,P*P,P*A,A,A*A]
coefs_cc[i,1]+
i == 5 ? coefs_cc[i,2]*w[i,t,j] : coefs_cc[i,2]*g[i,t,j] + #se i = 5 então w[] senão g[]
i == 5 ? coefs_cc[i,3]*w[i,t,j]^2 : coefs_cc[i,3]*g[i,t,j]^2 +
i == 5 ? coefs_cc[i,4]*w[i,t,j]*hl[t] : coefs_cc[i,4]*g[i,t,j]*hl[t] +
coefs_cc[i,5]*hl[i,t]+
coefs_cc[i,6]*hl[i,t]^2)
# Lembre que Segredo é com vazao_real
#@NLexpression(model,cc_segredo[i=5,t=1:temp, j=1:n_unidades[1]],#[1,P,P*P,P*A,A,A*A]
# coefs_cc[i,1]+
# coefs_cc[i,2]*w[i,t,j] +
# coefs_cc[i,3]*w[i,t,j]^2 +
# coefs_cc[i,4]*w[i,t,j]*hl[t] +
# coefs_cc[i,5]*hl[i,t]+
# coefs_cc[i,6]*hl[i,t]^2)
for k = 1:usinas
for turbi = 1:n_unidades[k]
fix.(ç[k,0,turbi],false , force = true)
fix.(o[k,0,turbi],false , force = true)
fix.(u[k,0,turbi],true , force = true)
end
end
for t = 1:temp
fix.(q[3,t],10,force = true)
end
for k = 1:usinas
if ~(k in [3,6]) # para todas usinas menos Jordão (3) e Segredo (6)
for t = 1:temp #em todos os tempos
fix.(vaz_t[k,t], 0., force = true) #fixo vazão transferida = 0
end
end
fix.(v[k,0], v_max[k] , force = true)## encontrar os dados iniciais de volume
end
@objective(model, Min, sum(q[i,t] + s[i,t] + 1e6*def[i,t] for i = 1:usinas, t = 1:temp))#REVER tem que multiplicar o deficit por alguma coisa!!
#se i = 5 => q[5,t] = 0.
#senão q[i,t] = sum(w[i,t,j] for j = 1:n_unidades[i])
@NLconstraint(model, turbinamento_usina[i = 1:usinas, t = 1:temp],
q[i,t] == sum(w[i,t,j] for j = 1:n_unidades[i]))
#i == 5 ? q[i,t] == 0. : q[i,t] == sum(w[i,t,j] for j = 1:n_unidades[i]))
#restrição do volume de vertimento(soleira), if v <= k, s = 0 elseif v > k -> s ∈ [0, s_max]
#@constraint(model, vol_soleira[i = 1:usinas, t = 1:temp], ########################################################
# v[i,t] <= constante_por_usina ? s[i,t] = 0. : nothing)
#Se i = 3 então vaz_t[3,t] = vazão transferida de Jordão para Segredo
#Se i = 6 então vaz_t[6,t] = - vazão transferida de Jordão para Segredo
#Exemplo: Se Jordão manda 10 m3/s para Segredo então vaz_t[3,t] = -10 e vaz_t[6,t] = 10
vet_aux = zeros(usinas); vet_aux[3] = 1.; vet_aux[5] = -1. # usina 5 negativo pois se 3 recebe água então 5 perde e vice versa
@constraint(model,transferencia_vazao[i = 1:usinas, t = 1:temp],
vaz_t[i,t] == vet_aux[i]*(mon[3,t] - mon[5,t])/abs(mon[3,t]- mon[5,t])*(-72*sqrt(abs(mon[3,t]- mon[5,t])))) #0 se i not in [3,5]
#essa conta resulta na mesma operação que está abaixo:
#if i == 3
# if mon[3,t] - mon[5,t] > 0
# vaz_t[i,t] == -72*sqrt(mon[3,t] - mon[5,t])
# else
# vaz_t[i,t] == 72*sqrt(-mon[3,t] + mon[5,t])
# end
#elseif i == 5
# if mon[3,t] - mon[5,t] > 0
# vaz_t[i,t] == 72*sqrt(mon[3,t] - mon[5,t])
# else
# vaz_t[i,t] == -72*sqrt(-mon[3,t] + mon[5,t])
# end
#else
# vaz_t[i,t] == 0.
#end)
#usinas_geradoras = collect(Iterators.flatten((1:4, 6:usinas)) #pega todas usinas que geram algo, retira as que são somente reservatório
#Separando com if/elseif/else entre os rios
#Rio Jordão -> i = [1,2,3] observação: i = 6 também está no rio mas vamos considerar diferente
#Rio Iguaçu -> i = [4,5,6]
#Se i = 1 -> Estamos no rio Jordão
#Se i = 4 -> Estamos no rio Iguaçu
#Se i = 6 -> 1ᵃ usina que recebe os dois rios
@constraint(model,balanço_hidrico[i = 1:usinas, t = 1:temp],
#v[i,t] == v[i,t-1] - H.*(q[i,t]+s[i,t]-Y[i,t]+vaz_t[i,t])) #o turbinado de jordão entra no incremental de Segredo e o vertido de jordão entra no incremental de S.Santiago
if length(indice_usinas_montante[cod_usinas[i]]) == 0. #se é a primeira usina do rio
v[i,t] == v[i,t-1] - H.*(q[i,t]+s[i,t]-Y[i,t]+vaz_t[i,t]))
else
if t > δ[i]
v[i,t] == v[i,t-1] - H.*(q[i,t]+s[i,t]- sum(q[i,t-Int(t-δ[i])] + s[i,t-Int(t-δ[i])] for i = indice_usinas_montante[cod_usinas[i]]) + vaz_t[i,t]))
else
v[i,t] == v[i,t-1] - H.*(q[i,t]+s[i,t]- O[i] + vaz_t[i,t]))
end
end
#if i == 1 (v[i,t] == v[i,t-1] - H.*(q[i,t]+s[i,t]-Y[i,t]+vaz_t[i,t]))
#elseif i == 4 (v[i,t] == v[i,t-1] - H.*(q[i,t]+s[i,t]-Y[Int(i/2),t]+vaz_t[i,t]))
#elseif i == 6 (v[i,t] == v[i,t-1] - H.*(q[i,t]+s[i,t]-q[i,t]-s[i,t]+vaz_t[i,t]))
#else (v[i,t] == v[i,t-1] - H.*(q[i,t]+s[i,t]-q[i,t]-s[i,t]+vaz_t[i,t]))
# sum(q[k,t] + s[k,t] for k = usi_mon_natural[i])
@constraint(model,balanço_hidrico_qx[ t = 1:temp],
0 == (q[2,t]+s[2,t]-q[1,t]-s[1,t]))
@NLconstraint(model, atend_dem[i = collect(Iterators.flatten((1:2, 4:usinas))),t = 1:temp],
sum(g[i,t,j] for j = 1:n_unidades[i]) == L[t,cod_usinas[i]])#rever o def[i,t] com o Tiago
@NLconstraint(model, geracao[i = collect(Iterators.flatten((1:4, 6:usinas))),t = 1:temp, j = 1:n_unidades[i]],
9.80665*1e-3*η[i,t,j]*(1 - perda_valor[i])*hb[t]*w[i,t,j] == g[i,t,j])
@NLconstraint(model, geracao_segredo[t = 1:temp, j = 1:n_unidades[5]],
9.80665*1e-3*η[5,t,j]*hl*w[5,t,j]==g[5,t,j])
@NLconstraint(model, rendimento[i=collect(Iterators.flatten((1:4, 6:usinas))),t = 1:temp,j = 1:n_unidades[1]],
cc[i,t,j] == η[i,t,j])
@NLconstraint(model, rendimento_segredo[i=5,t = 1:temp,j = 1:n_unidades[1]],
cc_segredo[i,t,j] == η[i,t,j])
@NLconstraint(model, altura_bruta[i=1:usinas,t = 1:temp],
hb[i,t] - mon[i,t] + jus[i,t] == 0)
@constraint(model,lim_turbinamento_max[i = 1:usinas, t = 1:temp, j = 1:n_unidades[i]],
w_max[i]*u[i,t,j] ≥ w[i,t,j])
@constraint(model,lim_turbinamento_min[i = 1:usinas,t = 1:temp, j = 1:n_unidades[i]],
w[i,t,j] ≥ w_min[i]*u[i,t,j])
@constraint(model, lim_geracao_min[i = 1:usinas,t = 1:temp, j = 1:n_unidades[i]],
g[i,t,j] ≥ g_min[i]*u[i,t,j])
@constraint(model, lim_geracao_max[i = 1:usinas, t = 1:temp, j = 1:n_unidades[i]],
g_max[i]*u[i,t,j] ≥ g[i,t,j])
@constraint(model, start_up[i = 1:usinas,t = 1:temp, j = 1:n_unidades[i]],
u[i,t,j]-u[i,t-1,j] == ç[i,t,j]-o[i,t,j])
@constraint(model, n_start_up_max[i = 1:usinas,j = 1:n_unidades[i]],
sum(ç[i,t,j] for t = 1:temp) ≤ I_max[i] )
@constraint(model, n_star_up[i = 1:usinas,t = 1:temp,j = 1:n_unidades[i]],
sum(ç[i,n,j] for n = max(1,t - T_up[i]):t ) ≤ u[i,t,j])
@constraint(model, n_shut_down[i = 1:usinas, j = 1:n_unidades[i], t = 1:temp],
sum(o[i,n,j] for n = max(1,t - T_down[i]):t ) ≤ 1 - u[i,t,j])
optimize!(model)
error("AINDA NÃO")
if objective_value(model) > 0
q1, s1, g1, w1 , v1, hb1 = value.(q), value.(s),value.(g),value.(w), value.(v), value.(hb)#, g, hb, w, n = value.(q), value.(s),value.(v), value.(g), value.(hb),value.(w), value.(η)
#u1, I1, o1 = value.(ç),value.(I),value.(o)
n = value.(η)
println("----")
println("Vazão turbinada usina $(q1)")
println("----")
println("Vazão vertida usina $(s1)")
println("----")
println("Volume usina $(v1)")
println("----")
println("Altura bruta $(hb1)")
println("----")
println("Geração:", g1)
println("----")
println("Vazão turbinada unidade:", w1)
println("----")
horas = 0:temp
aux = []
for k in 0:24
push!(aux,v1[k])
end
volume = plot(horas, aux, xlabel = "Tempo (h)", ylabel ="Volume Útil (m³)" ,
title = string("Volume Útil do Reservatório com Grau ",grau), lw = 3, label = "Volume Útil")
g_esp = ones(24,3)
g_qx = ones(24,4)
unidade = [3 4]
for usi =1:2
for t =1:24
for uni =1:unidade[usi]
if usi ==1
g_esp[t,uni] = g1[usi,t,uni]
else
g_qx[t,uni] = g1[usi,t,uni]
end
end
end
end
horas = 1:24
geracao_esp = plot(horas, g_esp, xlabel = "Tempo (h)", ylabel ="Geração (MW/Médio)" ,
title = string("Geração por Unidade na Usina 1 com Grau ",grau),label = ["Turbina 1" "Turbina 2" "Turbina 3"], lw = 3)
geracao_qx = plot(horas, g_qx, xlabel = "Tempo (h)", ylabel ="Geração (MW/Médio)",
title = string("Geração por Unidade na Usina 2 com Grau ",grau),label = ["Turbina 1" "Turbina 2" "Turbina 3" "Turbina 4"], lw = 3)
w_esp = ones(24,3)
w_qx = ones(24,4)
n_esp = ones(24,3)
n_qx = ones(24,4)
for usi =1:2
for t =1:24
for uni = 1:unidade[usi]
if usi ==1
w_esp[t,uni] = w1[usi,t,uni]
n_esp[t,uni] = n[usi,t,uni]
else
w_qx[t,uni] = w1[usi,t,uni]
n_qx[t,uni] = n[usi,t,uni]
end
end
end
end
vazao_turbi_esp = plot(horas, w_esp, xlabel = "Tempo (h)", ylabel ="Vazão Turbinada (m³/s)" ,
title = string("Vazão turbinada por Unidade na Usina 1 com Grau ",grau),label = ["Turbina 1" "Turbina 2" "Turbina 3" ], lw = 3)
vazao_turbi_qx = plot(horas, w_qx, xlabel = "Tempo (h)", ylabel ="Vazão Turbinada (m³/s)" ,
title = string("Vazão turbinada por Unidade na Usina 2 com Grau ",grau),label = ["Turbina 1" "Turbina 2" "Turbina 3" "Turbina 4"], lw = 3)
rendimento_esp = plot(horas, n_esp*100 , xlabel = "Tempo (h)", ylabel = "Rendimento (%)" ,
title = string("Rendimento por Unidade na Usina 1 com Grau ",grau),label = ["Turbina 1" "Turbina 2" "Turbina 3" ], lw = 3)
#rend_qx = -0.18910719 .+ 0.11302677.*w_qx .- 0.0028301.*(w_qx.^2)
rendimento_qx = plot(horas, n_qx*100 , xlabel = "Tempo (h)", ylabel = "Rendimento (%)" ,
title = string("Rendimento por Unidade na Usina 2 com Grau ",grau),label = ["Turbina 1" "Turbina 2" "Turbina 3" "Turbina 4"], lw = 3)
d_esp = ones(24,2)
d_qx = ones(24,2)
for usi = 1:2
for t = 1:24
if usi ==1
d_esp[t,1] = q1[usi,t]
d_esp[t,2] = s1[usi,t]
else
d_qx[t,2] = s1[usi,t]
d_qx[t,1] = q1[usi,t]
end
end
end
deflu_esp = plot(horas, [d_esp Y], xlabel = "Tempo (h)", ylabel ="Vazão (m³/s)" ,
title = string("Defluência na Usina 1 com Grau ",grau),label = ["Turbinada" "Vertida" "Afluente"], lw = 3)
afluente = d_esp[:,1].+d_esp[:,2] #queixada
deflu_qx = plot(horas, [d_qx afluente], xlabel = "Tempo (h)", ylabel ="Vazão (m³/s)" ,
title = string("Defluência na Usina 2 com Grau ",grau),label = ["Turbinada" "Vertida" "Afluente" ], lw = 3)
dem_esp = plot(horas, L1, xlabel = "Tempo (h)", ylabel ="Geração (MW/Médio)",
title = string("Demanda na Usina 1 com Grau ",grau),label = "Demanda", lw = 3)
texto = string("cascata_espora")
texto2 = string("cascata_queixada")
pasta = string("Finardi/Todos")#/Cascata_sDemanda")
#----
savefig(volume,string(pasta,"/vol_",texto,"_",grau,".svg"))
savefig(geracao_esp,string(pasta,"/ger",texto,"_",grau,".svg"))
savefig(vazao_turbi_esp,string(pasta,"/vaz_turbi_",texto,"_",grau,".svg"))
savefig(deflu_esp,string(pasta,"/deflu_",texto,"_",grau,".svg"))
savefig(rendimento_esp,string(pasta,"/rend_",texto,"_",grau,".svg"))
savefig(dem_esp,string(pasta,"/dem_",texto,"_",grau,".svg"))
savefig(geracao_qx,string(pasta,"/ger",texto2,"_",grau,".svg"))
savefig(vazao_turbi_qx,string(pasta,"/vaz_turbi_",texto2,"_",grau,".svg"))
savefig(deflu_qx,string(pasta,"/deflu_",texto2,"_",grau,".svg"))
savefig(rendimento_qx,string(pasta,"/rend_",texto2,"_",grau,".svg"))
println("Terminou")
end
df = DataFrame(horas = horas, Turb_esp = q1[1,1:end], Turb_qx = q1[2,1:end],
Vert_esp = s1[1,1:end], Vert_qx = s1[2,1:end],
Vol_oti = vol_oti[1:24,1], Vol_real = vol ,Hb_oti = hb1, Hb_real = h_dia,
g1_esp = g_esp[1:end,1], g2_esp = g_esp[1:end,2], g3_esp = g_esp[1:end,3],
g1_qx = g_qx[1:end,1], g2_qx = g_qx[1:end,2], g3_qx = g_qx[1:end,3], g4_qx = g_qx[1:end,4],
w1_esp = w_esp[1:end,1], w2_esp = w_esp[1:end,2], w3_esp = w_esp[1:end,3],
w1_qx = w_qx[1:end,1], w2_qx = w_qx[1:end,2], w3_qx = w_qx[1:end,3], w4_qx = w_qx[1:end,4],
n1_esp = n_esp[1:end,1], n2_esp = n_esp[1:end,2], n3_esp = n_esp[1:end,3],
n1_qx = n_qx[1:end,1], n2_qx = n_qx[1:end,2], n3_qx = n_qx[1:end,3], n4_qx = n_qx[1:end,4],
Dem_esp = L1,
CSV.write(string(pasta,"Dados_rodagem_",dia,"_",mes,"_",ano,".csv"), df)
t2 = now() - t1
#println("Tempo de execução do código: ",canonicalize(Dates.CompoundPeriod(t2)))