#= ercot_stoch_reliability.jl
solves the a stochastic two stage unit-commitment model in which slow generators
are committed in the first stage and fast generators can be committed in real time.
Uses formulation based on Papavasiliou and Oren (2013).
Demand response can function as a slow generator with an advance commmitment
and an advane production schedule if desired.
Uncertainty can come from total demand and from the reliabiltiy of DR. To do
the latter, DR must receive a first-stage production schedule (i.e. be type
3.).

cmd line format should be:
julia ercot_stoch_reliability.jl <date> <inputs_file_name.csv> <multi-runTF> <period_name>

multi-runTF is usually true: it indicates that the current simulation is one time
period of a group of periods that make up a single simulation. This just affects
the file structure.

Written under Julia 0.6.4
Patricia Levi
pjlevi@stanford.edu
NOV 2018
=#

### Packages ###
using JuMP
using Gurobi
using DataFrames
using CSV
include("convert3dto2d.jl")
include("make_scenarios_rejection_sample.jl")
include("writecsvmulti.jl")
include("combine_scenarios.jl")

test = Gurobi.Env() # test that gurobi is working

# -------------------------------------------
# USER CONTROLS & CMD LINE ARGS ##########
# -------------------------------------------

# USER PARAMS #
no_vars = false #stops execution before making variables
debug = false # stops execution before solving model

sherlock_fol = "~/dr_stoch_uc/"
sherlock_input_file = "inputs/"
sherlock_output_file = "outputs/"
sherlock_params_fol = "code/"
laptop_fol = "~/Documents/dr_stoch_uc/"
laptop_input_file = "inputs/"
laptop_output_file = "outputs/"
laptop_params_fol = "code/"


### DR is added to slow/fast categories depending on value of
## DRtype (in inputs csv)
dr_gens = ["DR"]
notdr_gens = ["BIOMASS","COAL","GAS","GAS_CC","GAS_CT","GAS_ICE","GAS_ST","HYDRO",
            "LANDFILL_GAS","NUCLEAR","OIL","SOLAR","WIND"    ]

####### END USER CONTROLS ##########

# --------------------------------------------------------------------------------------
# FILEPATH SETUP
# -------------------------------------------
if split(pwd(),"/")[2] == "Users" #returns true if we're on Patricia's laptop
    Sherlock = false
else
    Sherlock = true
end
@show Sherlock

if Sherlock
    base_fol = sherlock_fol
    input_fol = string(sherlock_fol,sherlock_input_file)
    params_fol = string(sherlock_fol, sherlock_params_fol)
else
    base_fol = laptop_fol
    input_fol = string(laptop_fol,laptop_input_file)
    params_fol = string(laptop_fol, laptop_params_fol)
end
default_data_fol = string(input_fol,"ercot_default/")


# PARSE CMD LINE ARGS #
ARGNAMES = ["date" ,"inputs_file_name","input_verion" ,"multi-runTF", "period_name" ]
defaultARGS = [Dates.format(Dates.now(),"Y-m-d"),"inputs_ercot.csv","rand_o25_00pp_keyDays2","true","periods_1_468_588.csv"]
localARGS = length(ARGS) > 0 ? ARGS : defaultARGS #if ARGS supplied, use those. otherwise, use default
nargs = length(localARGS)
@show localARGS

if nargs == 5
    submitdate = localARGS[1]
    input_file_name = localARGS[2]
    input_version = localARGS[3]
    multiTF = parse(Bool,lowercase(localARGS[4]))
    periodID = localARGS[5]
elseif nargs > 5
    error(string("Too many arguments supplied. Need ",join(ARGNAMES[1,:]," ")))
elseif nargs <5
    warn("not enough arguments supplied. Need ", join(ARGNAMES[1,:]," "))
end


read_inputs = CSV.read(string(params_fol,input_file_name))
read_inputs = read_inputs[:,[:input_name, Symbol(input_version)]]
# add cmd line args to read_inputs file
newdf = DataFrame(input_name = ARGNAMES[1:length(localARGS)], args = localARGS)
names!(newdf.colindex,map(parse,["input_name",input_version]))
read_inputs = vcat(read_inputs, newdf)
@show read_inputs
# ~ transpose read_inputs so that values can be referenced by name
read_inputs[:,:rowkey] = 1
inputs = unstack(read_inputs,:rowkey, :input_name, Symbol(input_version))

# parse out non-string inputs
startlim = parse(Float64,inputs[1,:startlim])
hourlim = parse(Float64,inputs[1,:hourlim])
energylim = parse(Float64,inputs[1,:energylim])
ramplims = parse(Float64,inputs[1,:ramplims])
durationlim = parse(Float64, inputs[1,:durationlim])
dr_override = parse(Bool,lowercase(inputs[1,:dr_override]))
dr_varcost = parse(Float64,inputs[1,:dr_varcost])
randScenarioSel = parse(Bool,lowercase(inputs[1,:randScenarioSel]))
trueBinaryStartup = parse(Bool,lowercase(inputs[1,:trueBinaryStartup]))
DRtype = parse(Int64,inputs[1,:DRtype])
DRrand = parse(Bool,lowercase(inputs[1,:DRrand]))
nd_nrand_o = parse(Int64,inputs[1,:nd_nrand_o])
int_length = parse(Int64,inputs[1,:intlength])
debug = !parse(Bool, lowercase(inputs[1,:solve_model]))
MIPFocusParam = parse(Int64,inputs[1,:MIPFocusParam])
MIPGapParam = parse(Float64,inputs[1,:MIPGapParam])
NodefileStartParam = parse(Float64,inputs[1,:NodefileStartParam])
MethodParam = parse(Int64,inputs[1,:MethodParam])
ThreadsParam = parse(Int64,inputs[1,:ThreadsParam])

if Sherlock
    output_fol = string(sherlock_fol, sherlock_output_file, input_version,"_",submitdate,"/")
else
    output_fol = string(laptop_fol, laptop_output_file, input_version,"_",submitdate,"/")
end
subsel_data_fol = string(input_fol,inputs[1,:timeseriesID],"/")

# setup unique ID for period identification in output
if multiTF
    pbits = split(periodID,r"_|\.|-")
    periodnum = pbits[2]
    periodfirst = pbits[3]
    periodlast = pbits[4]
    periodsave = string("p",periodnum,"_",periodfirst,"_",periodlast)
end

# --------------------------------------------------------------------------------------
### TIME PERIOD DATA
# -------------------------------------------
if multiTF
    hours = CSV.read(string(subsel_data_fol,periodID),datarow=1,types=[Int])[1]
else
    hours = CSV.read(string(subsel_data_fol,"periods_",timeseriesID,".csv"),datarow=1,types=[Int])[1]
end
n_days = convert(Int32,ceil(length(hours)/24))
n_periods = length(hours)

# convert first and notfirst periods to indices
t_firsts = 1
t_notfirst = collect(2:length(hours))

dem2 = CSV.read(string(default_data_fol, inputs[1,:demandFile]),datarow=2,missingstring="NA")
# subselect for just the rows corresponding to 'hours'
dem = dem2[hours,2]

# --------------------------------------------------------------------------------------
# STOCHASTIC VARIABLE DATA
# -------------------------------------------
# TESTING
# vdr = [0.9,1,1.1]
# pro = [0.25,0.5,0.25]
# probs = CSV.read(string(default_data_fol , "dist_input_",stochID,".csv"))
# vdr_in = convert(Array,probs[1,:]) # converts the first row of probs to an Array
# pro_in = rationalize.(convert(Array,probs[2,:])) #to avoid rounding issues later
    # if this becomes a problem, look into https://github.com/JuliaMath/DecFP.jl

## Net Demand uncertainty ##
# if _vdem files do not exist already, assume that the distribution data
# is given by a file describing the probability distribution of deviations
if !isfile(string(subsel_data_fol,"demandScenarios_vdem","_",inputs[1,:stochID],"_",periodsave,".csv"))
    ndprobs = CSV.read(string(default_data_fol , "dist_input_",inputs[1,:stochID],"_nd.csv"))
    ndv_in = convert(Array,ndprobs[1,:]) # converts the first row of probs to an Array
    ndpro_in = rationalize.(convert(Array,ndprobs[2,:])) #to avoid rounding issues later

    demandScenarios = make_scenarios(n_periods, ndv_in, ndpro_in, int_length; randsel = randScenarioSel, nrand = nd_nrand_o)
    vdem_in = demandScenarios[1]
    vdem_p = demandScenarios[2]
    writecsvmulti(DataFrame(vdem_in),subsel_data_fol,string("demandScenarios_vdem","_",inputs[1,:stochID]),multiTF,periodsave)
    writecsvmulti(DataFrame(vdem_p),subsel_data_fol,string("demandScenarios_prob","_",inputs[1,:stochID]),multiTF,periodsave)
else
    vdem_in = convert(Array,CSV.read(string(subsel_data_fol,"demandScenarios_vdem","_",inputs[1,:stochID],"_",periodsave,".csv")))
    vdem_p = convert(Array,CSV.read(string(subsel_data_fol,"demandScenarios_prob","_",inputs[1,:stochID],"_",periodsave,".csv")))
end

## DR Response uncertainty ##
if !isfile(string(subsel_data_fol,"drResponseScenarios_vdr","_",inputs[1,:DRrand_ID],"_",periodsave,".csv"))
    dr_int_length = parse(Int64,inputs[1,:dr_int_length])
    dr_nrand_o = parse(Int64,inputs[1,:dr_nrand_o])

    drprobs = CSV.read(string(default_data_fol , "dist_input_",inputs[1,:DRrand_ID],"_dr.csv"))
    drv_in = convert(Array,drprobs[2,:]) # converts the first row of probs to an Array
    drpro_in = rationalize.(convert(Array,drprobs[1,:])) #to avoid rounding issues later

    demandScenarios = make_scenarios(n_periods, drv_in, drpro_in, dr_int_length; randsel = randScenarioSel, nrand = dr_nrand_o)
    vdr_in = demandScenarios[1]
    vdr_p = demandScenarios[2]
    writecsvmulti(DataFrame(vdr_in),subsel_data_fol,string("drResponseScenarios_vdr","_",inputs[1,:DRrand_ID]),multiTF,periodsave)
    writecsvmulti(DataFrame(vdr_p),subsel_data_fol,string("drResponseScenarios_vdr_p","_",inputs[1,:DRrand_ID]),multiTF,periodsave)
else
    vdr_in = convert(Array,CSV.read(string(subsel_data_fol,"drResponseScenarios_vdr_",inputs[1,:DRrand_ID],"_",periodsave,".csv")))
    vdr_p = convert(Array,CSV.read(string(subsel_data_fol,"drResponseScenarios_vdr_p_",inputs[1,:DRrand_ID],"_",periodsave,".csv")))
end


## Create new vdr, vdem, pro representing combined scenarios from demand
## and DR reliability
println(vdr_p)
println(vdem_p)
scenario_info = combine_scenarios(vdr_in,vdr_p,vdem_in,vdem_p)
vdr = scenario_info[1];
vnd = scenario_info[2];
pro = scenario_info[3];
writecsvmulti(DataFrame(vdr),subsel_data_fol,string("vdr","_",inputs[1,:DRrand_ID]),multiTF,periodsave)
writecsvmulti(DataFrame(vnd),subsel_data_fol,string("vnd","_",inputs[1,:DRrand_ID]),multiTF,periodsave)
writecsvmulti(DataFrame(pro),subsel_data_fol,string("pro","_",inputs[1,:DRrand_ID]),multiTF,periodsave) ## WARNING: passing columns argument with non-AbstractVector entries is deprecated



println("sum of probabilities is ", sum(pro))
n_omega = length(pro) #redefine for new number of scenarios

# Set up demand realizations matrix dem_real[t,o] #
dem_real = dem .* vnd


# --------------------------------------------------------------------------------------
# GENERATOR DATA
# -------------------------------------------

genset = CSV.read(string(default_data_fol,inputs[1,:genFile]), missingstring ="NA")
# genset[Symbol("Plant Name")] # this is how to access by column name if there are spaces
# names(genset) # this is how to get the column names
# anscombe[:,[:X3, :Y1]]  #how to grab several columns by colname

dr_ind = findin(genset[:Fuel],dr_gens)
slow_ind = findin(genset[:Speed],["SLOW"])
fast_ind = findin(genset[:Speed],["FAST"])

# adjust slow/fast defn based on DRtype
if DRtype == 1
    fast_ind = vcat(fast_ind_, dr_ind)
elseif (DRtype == 2) | (DRtype == 3)
    slow_ind = vcat(slow_ind, dr_ind)
else
    error("DRtype must be 1, 2 or 3")
end

#generator min and max
pmin = genset[:PMin]
pmax = genset[:Capacity]
startup = genset[:StartCost]
varcost = genset[:VCost]
rampmax = genset[:ramprate]

# for manual override of DR variable cost
if dr_override
    varcost[dr_ind] = dr_varcost
end

# Generator PARAMS ###
n_gsl= length(slow_ind)# number of slow generators
n_g =nrow(genset)# number of generators
n_gf = length(fast_ind)
n_gdr = length(dr_ind) #number of DR generators
n_t = n_periods # number of timesteps

## GENERATOR AVAILABILITY ##
pf = repeat([1.0], inner = [n_g, n_t])

### Wind and solar
# load wind, solar info
solar_avail = CSV.read(string(default_data_fol,"solar_availability_factors_2016.csv"))
wind_avail = CSV.read(string(default_data_fol,"wind_availability_factors_2016.csv"))
# remove last col of each
solar_avail = solar_avail[:,1:(ncol(solar_avail)-1)]
wind_avail = wind_avail[:,1:(ncol(wind_avail)-1)]
# loop through all colnames, use findin(genset[:plantUnique],XX) to get row
# sub in new info
for i in 1:length(names(solar_avail))
    col = names(solar_avail)[i]
    ind = findin(genset[:plantUnique],[convert(String, col)])
    pf[ind,:] = solar_avail[hours,i]
end

for i in 1:length(names(wind_avail))
    col = names(wind_avail)[i]
    ind = findin(genset[:plantUnique],[convert(String, col)])
    pf[ind,:] = wind_avail[hours,i]
end

### Demand Response
if parse(inputs[1,:availID])!=0
    dr_avail = CSV.read(string(default_data_fol,inputs[1,:availID],".csv"))
    # remove first col
    dr_avail = dr_avail[:,2:ncol(dr_avail)]
    # sub in new info
    for i in 1:length(names(dr_avail))
        col = names(dr_avail)[i]
        ind = findin(genset[:plantUnique],[convert(String, col)])
        if length(ind) > 0
            pf[ind,:] = dr_avail[hours,i]
        else
            error("No corresponding DR unit found for ", names(dr_avail)[i])
        end
    end
end


# --------------------------------------------------------------------------------------
# MODEL SETUP
# -------------------------------------------
### SETS ###
TIME = 1:n_t
SCENARIOS = 1:n_omega
GENERATORS = 1:n_g #all generators
GEN_NODR = findin(genset[:Fuel],notdr_gens)
GF = fast_ind
GSL = slow_ind #slow generators
GDR = dr_ind #DR generators
# need an index for where the DR is in the slow generators
GDR_SL_ind = findin(GSL,GDR)

# -------------------------------------------
### MODEL ###
# m = Model(solver = ClpSolver())
m = Model(solver=GurobiSolver(Threads = ThreadsParam,ImproveStartGap = 0.0002, Method=MethodParam,#1,
                              MIPFocus=MIPFocusParam,#3,
                              MIPGap=MIPGapParam,#0.0003,
                              NodefileStart = NodefileStartParam,#0.05,
                              NodefileDir = "/scratch/users/pjlevi/gurobi_solving_outputs/",
                              Seed = convert(Int64,abs(floor(rand(Float64)*2000000000)))))
                # can also try reducing threadcount, reducing nodefilestart
if no_vars
    error("just testing model so we are stopping here")
else
    println("begin building model")
end

# -----------------------------------------------------------------------------------------------------
### VARIABLES ###
# @variable(m, 0 <= x <= 2 )
@variable(m, z[1:n_gsl,1:n_t]) # slow generator startup
if trueBinaryStartup
    @variable(m, w[1:n_gsl,1:n_t], Bin) # slow generator commitment TRUE BINARY
else
    @variable(m, 0 <= w[1:n_gsl,1:n_t] <= 1) # slow generator commitment RELAXED BINARY
end

#real-time commitment, startup
@variable(m, v[1:n_g,1:n_t,1:n_omega]) # generator startup
if trueBinaryStartup
    @variable(m, u[1:n_g,1:n_t,1:n_omega], Bin) # generator commitment TRUE BINARY
else
    @variable(m, 0 <= u[1:n_g,1:n_t,1:n_omega] <= 1) # generator commitment RELAXED BINARY
end

#production variables
@variable(m, p[1:n_g,1:n_t,1:n_omega] >= 0) #generator production

# if DR requires advance production schedule
if DRtype == 3
    @variable(m, p_dr[1:n_gdr, 1:n_t] >= 0) # DR day-ahead production commitment
end

# track number of startups to determine costs
@variable(m, start_num[1:n_g, 1:n_t, 1:n_omega] >= 0)


# -----------------------------------------------------------------------------------------------------
### CONSTRAINTS ###
# @constraint(m, 1x + 5y <= 3.0 )

#SUPPLY-DEMAND
@constraint(m,supplydemand[t=1:n_t, o=1:n_omega],
    sum( p[g,t,o] for g=1:n_g) >= dem_real[t,o])

#GENMIN
@constraint(m, mingen[g= GEN_NODR, t= 1:n_t, o=1:n_omega ],
    p[g,t,o] >= pmin[g] * u[g,t,o])
#GENMAX - only for GEN_NODR if dr is random
@constraint(m,[g = GEN_NODR, t = 1:n_t, o=1:n_omega],
    p[g,t,o] <= pmax[g] * u[g,t,o] * pf[g,t] )

if DRtype == 3
#GENMAXDR
    @constraint(m,[g = 1:n_gdr, t = 1:n_t],
        p_dr[g,t] <= pmax[GDR[g]] * w[GDR_SL_ind[g],t] * pf[GDR_SL_ind[g],t]) #needed for p_dr
#GENMINDR
    @constraint(m,[g = 1:n_gdr, t = 1:n_t],
        p_dr[g,t] >= pmin[GDR[g]] * w[GDR_SL_ind[g],t])
end

#START_S
@constraint(m,[g = 1:n_gsl, t=t_notfirst],
    z[g,t] == w[g,t] - w[g,t-1])
#START_F
@constraint(m,[g=GENERATORS ,t=t_notfirst, o = SCENARIOS],
    v[g,t,o] == u[g,t,o] - u[g,t-1,o])
#INIT_S
@constraint(m,[g=1:n_gsl,t=t_firsts],
    z[g,t] == w[g,t])
#INIT_F
@constraint(m,[g=GENERATORS,t=t_firsts, o = SCENARIOS],
    v[g,t,o] == u[g,t,o])
#NAN_ST
@constraint(m,[g = 1:n_gsl, t = TIME, o = SCENARIOS],
    v[GSL[g],t,o] == z[g,t])
#NAN_CM
@constraint(m,[g = 1:n_gsl, t = TIME, o = SCENARIOS],
    u[GSL[g],t,o] == w[g,t])

#DR_RAND
if DRtype == 3 && DRrand
    @constraint(m, dr_rand[g=1:n_gdr,t = TIME, o = SCENARIOS],
         p[GDR[g],t,o] == p_dr[g,t] * vdr[t,o])
elseif DRtype != 3 && DRrand
    error("DRrand can only be used with a DRtype of 3")
end
#TODO: make an elseif for DRype ==3 and !DRrand

#STARTUP COUNT
@constraint(m, [g=GENERATORS, t = TIME, o = SCENARIOS],
    start_num[g,t,o] >= v[g,t,o])

#RAMP RATE -
if ramplims !=0
    @constraint(m,ramplim[g=GENERATORS,t=t_notfirst, o = SCENARIOS],
        p[g,t,o] - p[g,t-1,o] <= (rampmax[g] * pmax[g]))
end

# ------ DR USAGE LIMITS
# Require that each period starts at the beginning of a new day
# to simplify computation

#STARTLIM
# number of startups per period
if startlim !=0
    # # per day
    # @constraint(m,[g = 1:n_gdr,d = 1:n_days, o = SCENARIOS],
    #     sum(v[GDR[g],t,o] for t = (24*(d-1)+1):(24*d)) <= startlim)
    # per period
    @constraint(m,limstart[g = 1:n_gdr, o = SCENARIOS],
        sum(start_num[GDR[g],t,o] for t = TIME) <= startlim)
end

#HOURLIM
#number of hours used per period
if hourlim != 0
    # # per day
    # @constraint(m,[g = 1:n_gdr,d = 1:n_days, o = SCENARIOS],
    #     sum(u[GDR[g],t,o] for t = (24*(d-1)+1):(24*d)) <= hourlim)
    # per period
    @constraint(m,limhrs[g = 1:n_gdr, o = SCENARIOS],
        sum(u[GDR[g],t,o] for t = TIME) <= hourlim)
end

#ENERGYLIM
# amount of energy used per period
if energylim != 0
    # per day
    # @constraint(m,[g = 1:n_gdr,d = 1:n_days ,o=SCENARIOS],
    #     sum(p[GDR[g],t,o] for t = (24*(d-1)+1):(24*d)) <= energylim)
    # per period
    @constraint(m,limenergy[g = 1:n_gdr, o=SCENARIOS],
        sum(p[GDR[g],t,o] for t = TIME) <= energylim)
end

#DURATIONLIM
# number of hours in a row that DR can be committed
if durationlim != 0
    l = Int(durationlim + 1)
    s = 1
    WINDOWS = [TIME[i:i+l-1] for i in range(1,s,Int(length(TIME)-l+1))]
    for i in range(1,s,Int(length(TIME)-l+1))
        print(TIME[i:i+l-1])
    end
    for W in WINDOWS
        @constraint(m,[g = 1:n_gdr, o = SCENARIOS],
            sum(u[GDR[g],t,o] for t = W) <= durationlim)
    end
end

### OBJECTIVE ###
@objective(m, Min, sum(pro[o] *
    sum(start_num[g,t,o]*startup[g] + p[g,t,o]*varcost[g] for g = GENERATORS, t = TIME)
    for o = SCENARIOS))

if debug
    error("just testing model so we are stopping here")
end

# -------------------------------------------
# SOLVE MODEL
# -------------------------------------------
@printf("\nSolving:\n")
status = solve(m)
@printf("Status: %s\n", status)

if !isdir(output_fol)
    mkdir(output_fol)
end

## save a copy of inputs & demand scenarios to the output file
writecsvmulti(read_inputs,output_fol,"inputfile",multiTF,periodsave)

# --------------------------------------------------------------------------------------
# SAVE OUTPUT
# -------------------------------------------

# production of DR
y = getvalue(p[GDR,:,:])

y_out = convert3dto2d(y,1, 3,  2,
    vcat([String("o$i") for i in 1:n_omega],"DR_IND","t"),
     genset[dr_ind,:plantUnique])
writecsvmulti(y_out,output_fol,"DR_production",multiTF,periodsave)

# production of slow generators
zs = getvalue(p[GSL,:,:])

y_out = convert3dto2d(zs,1, 3, 2,
    vcat([String("o$i") for i in 1:n_omega],"GEN_IND","t"),
     genset[slow_ind,:plantUnique])
writecsvmulti(y_out,output_fol,"slow_production",multiTF,periodsave)

# production of fast generators
zf = getvalue(p[GF,:,:])
y_out = convert3dto2d(zf,1, 3, 2,
    vcat([String("o$i") for i in 1:n_omega],"GEN_IND","t"),
     genset[fast_ind,:plantUnique])
writecsvmulti(y_out,output_fol,"fast_production",multiTF,periodsave)

# commitment of slow generators
w_out = getvalue(w)
wdf = DataFrame(transpose(w_out))
names!(wdf,[Symbol("$input") for input in genset[slow_ind,:plantUnique]])
writecsvmulti(wdf,output_fol,"slow_commitment",multiTF,periodsave)

# commitment of all generators
u_out = getvalue(u)
y_out = convert3dto2d(u_out,1, 3, 2,
    vcat([String("o$i") for i in 1:n_omega],"GEN_IND","t"),
     genset[:plantUnique])
writecsvmulti(y_out,output_fol,"u_commitment",multiTF,periodsave)

# startup of all generators
v_out = getvalue(v)
y_out = convert3dto2d(v_out,1, 3, 2,
    vcat([String("o$i") for i in 1:n_omega],"GEN_IND","t"),
     genset[:plantUnique])
writecsvmulti(y_out,output_fol,"v_startup",multiTF,periodsave)

# calculate costs
totcost = getvalue(sum(pro[o] *
    sum(start_num[g,t,o]*startup[g] + p[g,t,o]*varcost[g] for g = GENERATORS, t = TIME)
    for o = SCENARIOS))
totstartupcost = getvalue(sum(pro[o] *
    sum(start_num[g,t,o]*startup[g] for g = GENERATORS, t = TIME)
    for o = SCENARIOS))
totvarcost = getvalue(sum(pro[o] *
    sum(p[g,t,o]*varcost[g] for g = GENERATORS, t = TIME)
    for o = SCENARIOS))

output_summary = DataFrame(TotalCost = totcost, TotStartupCst = totstartupcost,
                            TotVarCst = totvarcost,
                            FracStartupCost = totstartupcost/totcost,
                            FracVarCost = totvarcost/totcost,
                            timeseries = inputs[1,:timeseriesID],
                            stoch_scenario = inputs[1,:stochID])

writecsvmulti(output_summary,output_fol,"summary_stats",multiTF,periodsave)

# Get dual variables and save
# only works if not using true binary startup
if !trueBinaryStartup
    #supplydemand #2D
    sd_shadow = DataFrame(getdual(supplydemand))
    writecsvmulti(sd_shadow,output_fol,"supplydemand_shadow",multiTF,periodsave)

    ### mingen ###
    mingen_shadow = getdual(mingen)
    mingen_sf = convert3dto2d(mingen_shadow,1, 3,  2,
        vcat([String("o$i") for i in 1:n_omega],"GEN_IND","t"),
         genset[:,:plantUnique])
    # only save rows where there is a nonzero shadow price
    shadowsum = sum(convert(Array{Float64},mingen_sf[:,1:n_omega]),2)
    saveind = find(shadowsum .!= 0)
    writecsvmulti(mingen_sf[saveind,:],output_fol,"mingen_shadow",multiTF,periodsave)

    ### type3dr ###
    ## TODO: This does not work.
    if DRtype == 3
        type3dr_shadow = getdual(type3dr)
        type3dr_sf = convert3dto2d(type3dr_shadow,1, 3,  2,
            vcat([String("o$i") for i in 1:n_omega],"DR_IND","t"),
             genset[dr_ind,:plantUnique])
        # only save rows where there is a nonzero shadow price
        shadowsum = sum(convert(Array{Float64},type3dr_sf[:,1:n_omega]),2)
        saveind = find(shadowsum .!= 0)
        writecsvmulti(type3dr_sf[saveind,:],output_fol,"type3dr_shadow",multiTF,periodsave)
    end
end
