
#Pkg.add("ExcelFiles")
# Info: Using the Python distribution in the Conda package by default.
# To use a different Python version, set ENV["PYTHON"]="pythoncommand" and re-run Pkg.build("PyCall").

using ExcelFiles, DataFrames

years = [Symbol("$i") for i in 2006:2015];

demand = DataFrame(load("Data.xlsx", "Load!A1:K8761"));

head(demand)

demand = demand[2:end];

demand.colindex

size(demand)

windSouthEast = DataFrame(load("Data.xlsx", "SouthEast!A1:K8761"));


head(windSouthEast)

windNorth = DataFrame(load("Data.xlsx", "North!A1:K8761"));

windNorth =windNorth[2:end];

windSouthWest = DataFrame(load("Data.xlsx", "SouthWest!A1:K8761"));

windSouthWest = windSouthWest[2:end];

#SouthsolarInput rng=PINCHER-CREEK
solarSouth = DataFrame(load("Data.xlsx", "PINCHER-CREEK"));

solarSouth = solarSouth[2:end-1];

# par= MiddleSolarInput rng=EDMONTON-STONY-PLAIN-CS!A1:K8761 Rdim=1 Cdim=1
# par= NorthSolarInput rng=GRANDE-PRAIRIE-A!A1:K8761 Rdim=1 Cdim=1
solarMiddle = DataFrame(load("Data.xlsx", "EDMONTON-STONY-PLAIN-CS!A1:K8761"));

#solarMiddle

solarMiddle = solarMiddle[2:end]; 

solarNorth = DataFrame(load("Data.xlsx", "GRANDE-PRAIRIE-A!A1:K8761"));

solarNorth = solarNorth[2:end];

#solarNorth[Symbol("2006")]

size(solarSouth)

#solarSouth

#solarSouth

typeof(solarSouth)

#solarSouth1 = similar(solarSouth)

#* the unit of output in data table is w. transfor to MW. A solar farm with 50000 PV panel? may cause floating point problem, number is too small
 
for (i,j) in eachcol(solarSouth) # a tuple with symbol and array
  solarSouth[i] = j*50000/1000000;
end

for (i,j) in eachcol(solarMiddle) # a tuple with symbol and array
  solarMiddle[i] = j*50000/1000000;
end

for (i,j) in eachcol(solarNorth) # a tuple with symbol and array
  solarNorth[i] = j*50000/1000000;
end

#Para.csv
using CSV
para =  CSV.read("Para.csv")

para

para[1].== "wind"

turbine = parse(para[ para[1].== "wind",:UnitCap][1]);

turbine

panel =  parse(para[ para[1].== "solar",:UnitCap][1]);
panel

generation_assets = [Symbol(i) for i in para[:,1][2:end]]  

fossil_fuel = [Symbol(i) for i in para[:,1][end-1:end]]  

renewables = [Symbol(i) for i in para[:,1][2:3]]  

windLoc = [:se, :sw, :n]
solarLoc = [:s, :m, :n]

size(para)

#para.columns[1]

para.colindex

FixOMCosts = Dict( zip( generation_assets    , ([parse(i) for i in para[:FixOMCost][2:end]]) )   )

VarOMCosts = Dict( zip( generation_assets    , ([parse(i) for i in para[:VarOMCost][2:end]]) )   )

ONights= Dict( zip( generation_assets    , ([parse(i) for i in para[:ONight][2:end]]) )   )

UnitCaps  = Dict( zip( generation_assets    , ([parse(i) for i in para[:UnitCap][2:end]]) )   )

Maxis  = Dict( zip( generation_assets    , ([parse(i) for i in para[:Max][2:end]]) )   )

Emits  = Dict( zip( generation_assets    , ([parse(i) for i in para[:Emit][2:end]]) )   )

Ramps  = Dict( zip( generation_assets    , ([parse(i) for i in para[:Ramp][2:end]]) )   )

LTimes= Dict( zip( generation_assets    , ([parse(i) for i in para[:LTime][2:end]]) )   )

# # parameters
# ''''''
# round = Para("battery", "Round") ;
# duration = Para("battery", "duration") ;
# taxEm = Para("CT","CTax");
# turbine = Para("wind","UnitCap");
# maxturbine = Para("wind","Max");
# panel = Para("solar","UnitCap");
# maxPanel = Para("solar","Max");
# rate = Para("CT","Discount");
# ''''''






maximums = Dict( zip( generation_assets    , ([parse(i) for i in para[:Max][2:end]]) )   )

roundEfficiency = Dict( zip( generation_assets    , ([parse(i) for i in para[:Round][2:end]]) )   )

duration= Dict( zip( generation_assets    , ([i for i in para[:duration][2:end]]) )   )

taxEm = Dict( zip( generation_assets    , ([parse(i) for i in para[:Ctax][2:end]]) )   )

rate = Dict( zip( generation_assets    , ([parse(i) for i in para[:Discount][2:end]]) )   )

Vcost = Dict( zip( generation_assets    , ( VarOMCosts[i] + taxEm[i]* Emits[i]  for i in generation_assets) )   ) ;



Vcost

Fcost =   Dict( zip(generation_assets, (ONights[i]*rate[i]*((1+rate[i])^LTimes[i])/(((1+rate[i])^LTimes[i])-1) + 
                                            FixOMCosts[i] for i in generation_assets) 
                    )  
                );

Fcost

using JuMP
using Gurobi

m = Model(solver=GurobiSolver())

T = 8760

@variable(m, batteryCap>=0)

# energy in battery
@variable(m, v[1:T]>=0)

@variable(m, charge[1:T]>=0)
@variable(m, discharge[1:T]>=0)

@variable(m, fuelCap[fossil_fuel] >=0)


@variable(m, gen[ 1:T,fossil_fuel] >=0)

@variable(m, nSolarPlant[solarLoc] >=0)
@variable(m, nWindTurbine[windLoc] >=0)



wind = DataFrame(s = 1.0:T,m = 1.0:T, n = 1.0:T);

solar = DataFrame(s = 1.0:T,m = 1.0:T, n = 1.0:T);

year = 2006

Symbol("$year")

wind[:se] = windSouthEast[Symbol("$year")];
wind[:sw] = windSouthWest[Symbol("$year")];
wind[:n] = windNorth[ Symbol("$year")];

#* solar unit output per hour  (MW), number of panel is in 50000, so total solar ouput is MW
solar[:s] = solarSouth[Symbol("$year")];
solar[:m] = solarMiddle[Symbol("$year")];
solar[:n] = solarNorth[Symbol("$year")];

demand0 = demand[Symbol("$year")];

@objective(m, Min, sum(sum(wind[t,L]*nWindTurbine[L]*Vcost[:wind] for L in windLoc) +
                            (sum(solar[t,Ls]*nSolarPlant[Ls]*Vcost[:solar] for Ls in solarLoc))  +  
                            (sum(gen[t,f] *Vcost[f] for f in fossil_fuel )) +   
                            discharge[t]*Vcost[:battery] 
                            for t in 1:T
                            ) + 
                    sum(Fcost[:wind]*nWindTurbine[L]*turbine for L in windLoc) +
                    sum(Fcost[:solar]*nSolarPlant[Ls]*turbine for Ls in solarLoc) +
                    sum(Fcost[f]*fuelCap[f] for f in fossil_fuel) +
                    Fcost[:battery]*batteryCap
            );

@constraint(m,  sum(wind[1,L]*nWindTurbine[L] for L in windLoc) +
                                sum(solar[1,Ls]*nSolarPlant[Ls] for Ls in solarLoc) +
                                  sum(gen[1,f] for f in fossil_fuel) +
                                   discharge[1] -charge[1] >= demand0[1]);

   



# for t in 1:T    
#     @constraint(m, Balance[t], sum(wind[t,L].*nWindTurbine[L] for L in windLoc) +
#                                 sum(solar[t,Ls].*nSolarPlant[Ls] for Ls in solarLoc) +
#                                   sum(gen[t,f] for f in fossil_fuel) +
#                                    discharge[t] -charge[t] >= demand0[t]);
# end


@constraint(m, Balance0[t=1:T], sum(wind[t,L].*nWindTurbine[L] for L in windLoc) +
                            sum(solar[t,Ls].*nSolarPlant[Ls] for Ls in solarLoc) +
                              sum(gen[t,f] for f in fossil_fuel) +
                               discharge[t] -charge[t] >= demand0[t]);

 @constraint(m, capacityFuel[t=1:T, f = fossil_fuel], gen[t,f] <= fuelCap[f] );

 @constraint(m, capacityWind[L=windLoc], nWindTurbine[L] <= maximums[:wind] );

 @constraint(m, capacitySolar[Ls=solarLoc], nSolarPlant[Ls] <= maximums[:solar] );

@constraint(m, RampUp[t=2:T, f = fossil_fuel], gen[t, f] <= gen[t-1,f]+Ramps[f]*fuelCap[f] );

@constraint(m, RampDn[t=2:T, f = fossil_fuel], gen[t, f] >= gen[t-1,f]-Ramps[f]*fuelCap[f] );

@constraint(m, Storage[t=2:T], v[t] ==   v[t-1]+(roundEfficiency[:battery]*charge[t-1]-discharge[t-1]));


@constraint(m, MaxStorage[t=1:T], v[t] <= duration[:battery] * batteryCap);

@constraint(m, MaxCharge[t=1:T], charge[t] <= duration[:battery] * batteryCap);

@constraint(m, MaxDischarge[t=1:T], discharge[t] <= batteryCap);

solve(m)

# # include("examples_julia/2_uc.jl");

# n_machines = 3;
# time_steps = 12;

# fixed_costs = [5, 10, 1000];
# variable_costs = [800, 700, 20];
# max_power = [500, 700, 1500];
# startup_costs = [10, 5, 850];
# demands = [1500, 1200, 1800, 2000, 2200, 1750, 1600, 1500, 1200, 1500, 1600, 1800];


# using JuMP
# using Gurobi

# m = Model(solver=GurobiSolver())

# @variable(m, on[1:time_steps, 1:n_machines], Bin)
# @variable(m, power[1:time_steps, 1:n_machines] >= 0)
# @variable(m, start[1:time_steps, 1:n_machines], Bin)

# @objective(m, Min, sum(dot(fixed_costs, vec(on[t, :])) + 
#         dot(variable_costs, vec(power[t, :])) + dot(startup_costs, vec(start[t, :])) for t in 1:time_steps))


# for t in 1:time_steps
#     for i in 1:n_machines
#         @constraint(m, power[t, i] <= max_power[i] * on[t, i])
#     end
#     @constraint(m, sum(power[t, :]) == demands[t])
# end


# for i in 1:n_machines
#     @constraint(m, start[1, i] == on[1, i])
# end

# for t in 2:time_steps
#     for i in 1:n_machines
#         @constraint(m, start[t, i] >= on[t, i] - on[t - 1, i])
#     end
# end


# for t in 1:time_steps
#     for i in 1:n_machines
#         @constraint(m, on[t, i] >= start[t, i])
#         if t <= time_steps - 1
#             @constraint(m, on[t + 1, i] >= start[t, i])
#         end
#         if t <= time_steps - 2
#             @constraint(m, on[t + 2, i] >= start[t, i])
#         end
#     end
# end

# print(m)
# solve(m)
# getvalue(on)


Gurobi.Env()
