# vim: set wrap
#: by David Thierry 2021
#: Set all the eqns to hslice - 1
using JuMP
using SCIP
using DataFrames
using CSV

#: Data frames section
#: Load parameters

df_gas = DataFrame(CSV.File("../reg/gas_coeffs.csv"))
df_steam_full_power = DataFrame(CSV.File("../reg/steam_coeffs.csv"))
df_steam_full_steam = DataFrame(CSV.File("../reg/steam_coeffs_v3.csv"))

#: Load Prices
df_pow_c = DataFrame(CSV.File("../resources/FLECCSPriceSeriesData.csv"))
df_ng_c = DataFrame(CSV.File("../resources/natural_gas_price.csv"))


#: Assign parameters

#: Gas parameters
bPowGasTeLoad = df_gas[1, 2]
aPowGasTeLoad = df_gas[1, 3]

bFuelEload = df_gas[2, 2]
aFuelEload = df_gas[2, 3]
lbcoToTonneco = 0.4535924 / 1000
bEmissFactEload = df_gas[3, 2] * lbcoToTonneco
aEmissFactEload = df_gas[3, 3] * lbcoToTonneco

bPowHpEload = df_gas[4, 2] / 1000  #: To scale the kW to MW
aPowHpEload = df_gas[4, 3] / 1000

bPowIpEload = df_gas[5, 2] / 1000
aPowIpEload = df_gas[5, 3] / 1000

bAuxRateGas = df_gas[6, 2] / 1000
aAuxRateGas = df_gas[6, 3] / 1000

#: Steam params
bCcRebDutyEload = df_steam_full_steam[2, 2]
aCcRebDutyEload = df_steam_full_steam[2, 3]

#: Full power gives you the min steam
bDacSteaBaseEload = df_steam_full_power[3, 2]  
aDacSteaBaseEload = df_steam_full_power[3, 3]

bSideSteaEload = df_steam_full_steam[3, 2] - df_steam_full_power[3, 2]
aSideSteaEload = df_steam_full_steam[3, 3] - df_steam_full_power[3, 3]

bAuxRateStea = df_steam_full_power[4, 2] / 1000
aAuxRateStea = df_steam_full_power[4, 3] / 1000

aLpSteaToPow = 78.60233832  # MMBtu to kwh

kwhToMmbtu = 3412.1416416 / 1e+06
#aSteaUseRateDacAir = 1944 * kwhToMmbtu
#aSteaUseRateDacFlue = 1944 * kwhToMmbtu
# 7 GJ/tonneCO
aSteaUseRateDacAir = 5 * (7 * 1e+06 / 3600) * kwhToMmbtu
aSteaUseRateDacFlue = 5 * (7 * 1e+06 / 3600) * kwhToMmbtu

aPowUseRateDacAir = 500 / 1000
aPowUseRateDacFlue = 250 / 1000
# 1 mmol/gSorb #
# per gCo/gSorb
gCogSorbRatio = 1e-03 * 44.0095

aSorbCo2CapFlue = 1 * gCogSorbRatio
aSorbCo2CapAir = gCogSorbRatio
#aSorbAmountFreshFlue = 176. * 10  # Tonne sorb
#aSorbAmountFreshAir = 176. * 10  # Tonne sorb

aSorbAmountFreshFlue = 176. * 10   # Tonne sorb (Max. heat basis)
aSorbAmountFreshAir = 3162.18 - 176. * 10  # Tonne sorb (Max. heat basis)


aCapRatePcc = 0.97
# 2.4 MJ/kg (1,050 Btu/lb) CO2 page 379/
#aSteaUseRatePcc = aSteaUseRateDacFlue * 0.2
#aSteaUseRatePcc = 2.4 * 1000 * 1000 / 3600 * kwhToMmbtu 
#println(aSteaUseRatePcc)
# aPowUseRatePcc = 0.173514487  # MWh/tonneCoi2 (old)
aSteaUseRatePcc = 2.69 + 0.0218 + 0.00127 # MMBTU/tonne CO2 (trimeric)
println(aSteaUseRatePcc)
aPowUseRatePcc = 0.047 # MWh/tonne CO2 (trimeric)

#: Horizon Lenght
tHorz = 24 * 10 - 1


#: Slices per hour
hSlice = 4  # the number of slices of a given hour
# There are tHorz - 1 slices
# Each slice has hSlice points, but only states have the 0th

# If there's several slices in an hour we kind of need to divide the
# hourly-based quantities :(



# Correction for low load
slopeFactor =  (0.01/ 0.1) / ((.95 - .75) / (0.9 - 0.6))

# Dictionaries with the disjunction parameters

## Slope
aPowGt = Dict(0 => 1, 1 => aPowGasTeLoad * slopeFactor, 2 => aPowGasTeLoad)
aFuel = Dict(0 => 1, 1 => aFuelEload, 2 => aFuelEload)
aEmis = Dict(0 => 1, 1 => aEmissFactEload, 2 => aEmissFactEload)
aAuxGt = Dict(0 => 1, 
    1 => aAuxRateGas * slopeFactor, 
    2 => aAuxRateGas)
aPowHp = Dict(0 => 1, 1 => aPowHpEload * slopeFactor, 
    2 => aPowHpEload)
aPowIp = Dict(0 => 1, 1 => aPowIpEload * slopeFactor, 
    2 => aPowIpEload)
aCcReb = Dict(
    0 => 1, 
    1 => aCcRebDutyEload * slopeFactor, 
    2 => aCcRebDutyEload)
aDacSb = Dict(0 => 1, 1 => aDacSteaBaseEload * slopeFactor, 
    2 => aDacSteaBaseEload)
aSideS = Dict(0 => 1, 1 => aSideSteaEload * slopeFactor, 2 => aSideSteaEload)
aAuxSt = Dict(0 => 1, 1 => aAuxRateStea * slopeFactor, 2 => aAuxRateStea)

ldics = ["aPowGt", "aFuel", "aEmis", "aAuxGt", "aPowHp", "aPowIp", "aCcReb", "aDacSb", "aSideS", "aAuxSt"]
println("Slopes")
for v in ldics
    s = Symbol(v)
    println(v, eval(s))
end

## b
bPowGt = Dict(0 => 0, 
    1 => 60 * (-aPowGt[1] + aPowGt[2]) + bPowGasTeLoad, 
    2 => bPowGasTeLoad)
bFuel = Dict(0 => 0, 1 => bFuelEload, 2 => bFuelEload)
bEmis = Dict(0 => 0, 1 => bEmissFactEload, 2 => bEmissFactEload)
bAuxGt = Dict(0 => 0, 
    1 => 60 * (-aAuxGt[1] + aAuxGt[2]) + bAuxRateGas, 
    2 => bAuxRateGas)
bPowHp = Dict(0 => 0, 
    1 => 60 * (-aPowHp[1] + aPowHp[2]) + bPowHpEload,
    2 => bPowHpEload)
bPowIp = Dict(0 => 0,
    1 => 60 * (-aPowIp[1] + aPowIp[2]) + bPowIpEload,
    2 => bPowIpEload
    )
bCcReb = Dict(0 => 0,
    1 => 60 * (-aCcReb[1] + aCcReb[2]) + bCcRebDutyEload,
    2 => bCcRebDutyEload)
bDacSb = Dict(0 => 0,
    1 => 60 * (-aDacSb[1] + aDacSb[2]) + bDacSteaBaseEload,
    2 => bDacSteaBaseEload)
bSideS = Dict(0 => 0,
    1 => 60 * (-aSideS[1] + aSideS[2]) + bSideSteaEload,
    2 => bSideSteaEload)
bAuxSt = Dict(0 => 0,
    1 => 60 * (-aAuxSt[1] + aAuxSt[2]) + bAuxRateStea,
    2 => bAuxRateStea)


ldics = ["bPowGt", "bFuel", "bEmis", "bAuxGt", "bPowHp", "bPowIp", "bCcReb", "bDacSb", "bSideS", "bAuxSt"]
println("InterC")
for v in ldics
    s = Symbol(v)
    println(v, eval(s))
end


# USD/MWh
pow_price =(df_pow_c[!, "MiNg_150_ERCOT"])  # USD/MWh

# pow_price =(df_pow_c[!, "MiNg_150_PJM-W"])  # USD/MWh
#: Natural gas price
# 0.056 lb/cuft STP
#std_w_ng1000cuft = 0.056 * 1000
#cNgPerLbUsd = (3.5 / 1000) / 0.056

# Cost of natural gas
cNgPerMmbtu = 3.5

m = Model()

# aPowUseRateComp = 0.279751187  # MWh/tonneCo2
aPowUseRateComp = 0.076 # MWh/tonneCo2 (Trimeric)

# Other costs
cCostInvCombTurb = 1e+02
cCostInvSteaTurb = 1e+02
cCostInvTransInter = 1e+02
cCostInvPcc = 1e+02
cCostInvDac = 1e+03
cCostInvComp = 1e+01

# Cost parameters.
cEmissionPrice = 1.5e+02 # USD/tonne CO2
cCo2TranspPrice = 1e+00
pCo2Credit = 1e+00


#vCapCombTurb = 3.
vCapSteaTurb = 2.
vCapTransInter = 5.
vCapPcc = 20.
vCapComp = 1000.
# Capital Cost DAC USD/tCo2/yr 
cCostInvDacUsdtCo2yr = 750
cCostFixedDacUsdtCo2yr = 25
cCostVariableDacUsdtCo2yr = 12

nMod = 2

# 0 -> off
# 1 -> warm-up
# 2 -> on

# 0 for off 2 for other ones.
extrPoint = Dict(0 => 0, 1 => 0:1, 2 => 0:1)
@variable(m, 0 <= lambda[0:tHorz, i=0:nMod, extrPoint[i]] <= 1)
# On/Off
@variable(m, y[0:tHorz, 0:nMod], Bin, start = 0)  # On and off
# Transition
@variable(m, z[0:tHorz, j1 = 0:nMod, j2 = 0:nMod], Bin, start = 0)

for i in 0:tHorz
    set_start_value(y[i, 1], 1)
end

@variable(m, xLoad[0:tHorz, 0:nMod])


@variable(m, 0 <= xeLoad[0:tHorz] <= 100)


@variable(m, 0 <= xPowGasTur[0:tHorz])
@variable(m, 0 <= xPowGross[0:tHorz])
@variable(m, 0 <= xPowOut[0:tHorz])

@variable(m, 0 <= xAuxPowGasT[0:tHorz])

# Steam Turbine
@variable(m, 0 <= xPowHp[0:tHorz])
@variable(m, 0 <= xPowIp[0:tHorz])
@variable(m, 0 <= xPowLp[0:tHorz])

@variable(m, 0 <= xFuel[0:tHorz])
@variable(m, 0 <= xCo2Fuel[0:tHorz])
@variable(m, 0 <= xDacSteaDuty[0:tHorz])


@variable(m, 0 <= xCcRebDuty[0:tHorz])
@variable(m, 0 <= xDacSteaBaseDuty[0:tHorz])

@variable(m, 0 <= xAllocSteam[0:tHorz])
@variable(m, 0 <= xSteaPowLp[0:tHorz])
@variable(m, 0 <= xSideSteaDac[0:tHorz])

#
@variable(m, 0 <= xPowSteaTur[0:tHorz])
@variable(m, 0 <= xAuxPowSteaT[0:tHorz])

# Disagretated variables

@variable(m, 0 <= xPowGt[0:tHorz, 0:nMod])
@variable(m, 0 <= xFuelNu[0:tHorz, 0:nMod])
@variable(m, 0 <= xEmis[0:tHorz, 0:nMod])
@variable(m, 0 <= xAuxGt[0:tHorz, 0:nMod])
@variable(m, 0 <= xPowHpNu[0:tHorz, 0:nMod])
@variable(m, 0 <= xPowIpNu[0:tHorz, 0:nMod])
@variable(m, 0 <= xPccReb[0:tHorz, 0:nMod])
@variable(m, 0 <= xDacSb[0:tHorz, 0:nMod])
@variable(m, 0 <= xAllocS[0:tHorz, 0:nMod])
@variable(m, 0 <= xAuxSt[0:tHorz, 0:nMod])



#@constraint(m, powgtDisEq[i = 0:tHorz, k = 0:nMod],
#    xPowGt[i, k] = aPowGt[k] * xLoad[i, k] + bPowGt[k])
#@constraint(m, powGasTur[i = 0:tHorz], 
#            xPowGasTur[i] == sum(xPowGt[i, k] for k in 0:nMod)
#           )

# Pcc
#@variable(m, 0 <= xCo2CapPcc[0:tHorz - 1] <= vCapPcc)
@variable(m, 0 <= xCo2CapPcc[0:tHorz])
@variable(m, 0 <= xSteaUsePcc[0:tHorz])
@variable(m, 0 <= xPowUsePcc[0:tHorz])
@variable(m, 0 <= xCo2PccOut[0:tHorz])
#@variable(m, 0 <= vCo2PccVent[0:tHorz - 1] <= 0.1)
@variable(m, 0 <= vCo2PccVent[0:tHorz])
#vCo2PccVent = 0.0
@variable(m, 0 <= xCo2DacFlueIn[0:tHorz])
@variable(m, 0 <= xPccSteaSlack[0:tHorz])

# Dac-Flue
@variable(m, 0 <= xA0Flue[0:tHorz, 0:hSlice-1]) # Kind of input 

@variable(m, 0 <= xA1Flue[0:tHorz, 0:hSlice]) # State


@variable(m, 0 <= xR0Flue[0:tHorz, 0:hSlice-1])  # Kind of input

@variable(m, 0 <= xR1Flue[0:tHorz, 0:hSlice])  # State
@variable(m, 0 <= xFflue[0:tHorz, 0:hSlice])  # State
@variable(m, 0 <= xSflue[0:tHorz, 0:hSlice])  # State


@variable(m, 0 <= xCo2CapDacFlue[0:tHorz, 0:hSlice-1])
@variable(m, 0 <= xSteaUseDacFlue[0:tHorz, 0:hSlice-1])
@variable(m, 0 <= xPowUseDacFlue[0:tHorz, 0:hSlice-1])
@variable(m, 0 <= xCo2DacVentFlue[0:tHorz])

# Dac-Air
@variable(m, 0 <= xA0Air[0:tHorz, 0:hSlice-1]) # Kind of input

@variable(m, 0 <= xA1Air[0:tHorz, 0:hSlice]) # State
@variable(m, 0 <= xA2Air[0:tHorz, 0:hSlice]) # State

@variable(m, 0 <= xR0Air[0:tHorz, 0:hSlice-1])  # Kind of input

@variable(m, 0 <= xR1Air[0:tHorz, 0:hSlice])  # State
@variable(m, 0 <= xFair[0:tHorz, 0:hSlice])  # State
@variable(m, 0 <= xSair[0:tHorz, 0:hSlice])  # State


@variable(m, 0 <= xCo2CapDacAir[0:tHorz, 0:hSlice-1])
@variable(m, 0 <= xSteaUseDacAir[0:tHorz, 0:hSlice-1])
@variable(m, 0 <= xPowUseDacAir[0:tHorz, 0:hSlice-1])

@variable(m, 0 <= xDacSteaSlack[0:tHorz])
# DAC hourly capacity
#



# CO2 compression
@variable(m, 0 <= xCo2Comp[0:tHorz])
@variable(m, 0 <= xPowUseComp[0:tHorz])
#@variable(m, 0 <= vCapComp)
@variable(m, xCo2Vent[0:tHorz])  # This used to be only positive.

@variable(m, 0 <= xAuxPow[0:tHorz])

# Constraints
# Op mode
# Down times
# Up times

# Disjunction 0 (off)

extreme_d_0 = [0]
extreme_d_1 = [20.0, 50.0]
extreme_d_2 = [60.0, 100.0]

ep_list = [extreme_d_0, extreme_d_1, extreme_d_2]

# Convex conbination
@constraint(m, cConvxEq[i = 0:tHorz, mod = 0:nMod],
    sum(lambda[i, mod, k] * ep_list[mod + 1][k + 1] for k in extrPoint[mod]) ==
    xLoad[i, mod]  # There is only a single extreme
    )

# jth mode
# Lambda constraint
@constraint(m, xLambdaEq[i = 0:tHorz, mod = 0:nMod],
    sum(lambda[i, mod, k] for k in extrPoint[mod]) == y[i, mod]
    )

# Big M
@constraint(m, bmconEq[i = 0:tHorz, mod = 0:nMod],
    xLoad[i, mod] <= maximum(ep_list[mod + 1])* y[i, mod]
    )


# Overall Disjunction

@constraint(m, gasLeq[i = 0:tHorz],
    xeLoad[i] == sum(xLoad[i, mod] for mod in 0:nMod)
    )


@constraint(m, sumyEq[i = 0:tHorz],
    sum(y[i, mod] for mod in 0:nMod) == 1
    )

# Switch

@constraint(m, switchCon[i = 1:tHorz, mod = 0:nMod],
    sum(z[i, k, mod] for k in 0:nMod) - 
    sum(z[i, mod, k] for k in 0:nMod) 
    == y[i, mod] - y[i - 1, mod])

# Forbidden
# from off to full
@constraint(m, forbConOffFullEq[i = 0:tHorz],
    z[i, 0, 2] == 0)
# from full to warm-up
@constraint(m, forbConEq2[i = 0:tHorz],
    z[i, 2, 1] == 0)
# from warm-up to off
#@constraint(m, forbConwarmpuoff[i = 0:tHorz],
#    z[i, 1, 0] == 0)

# Ramping Constraints
rup = [0, 15, 30]
rdown = [0, 15, 30]

@constraints(m, begin 
    rampup1eq[i=1:tHorz], xLoad[i, 1] - xLoad[i-1, 1] <= rup[2] * y[i, 1]
    rampup2eq[i=1:tHorz], xLoad[i, 2] - xLoad[i-1, 2] <= rup[3] * y[i, 2]
end)

#@constraints(m, begin 
    #rampdo1eq[i=1:tHorz], xLoad[i - 1, 1] - xLoad[i, 1] <= rdown[2] * y[i, 1]
    #rampdo2eq[i=1:tHorz], xLoad[i - 1, 2] - xLoad[i, 2] <= rdown[3] * y[i, 2]
#end)

# Minimum stay 
KminOff = [[0, 36, 0], [0, 0, 0], [48, 0, 0]]

# off to on(1)
@constraint(m, minstay01[i = 1:tHorz],
    y[i, 1] >= sum(z[i - k, 0, 1] for k in 0:(24-1) if (i-k) >= 0)
    )
# on(2) to off
@constraint(m, minstay20[i = 1:tHorz],
    y[i, 0] >= sum(z[i - k, 2, 0] for k in 0:(48-1) if (i-k) >= 0)
    )


# Disagregated variables constraints
@constraint(m, powgtdeq[i = 0:tHorz, mod = 0:nMod],
    xPowGt[i, mod] == aPowGt[mod] * xLoad[i, mod] + bPowGt[mod] * y[i, mod]
    )
@constraint(m, fueldeq[i = 0:tHorz, mod = 0:nMod],
    xFuelNu[i, mod] == aFuel[mod] * xLoad[i, mod] + bFuel[mod] * y[i, mod]
    )
@constraint(m, emissdeq[i = 0:tHorz, mod = 0:nMod],
    xEmis[i, mod] == aEmis[mod] * xLoad[i, mod] + bEmis[mod] * y[i, mod]
    )
@constraint(m, auxgtdeq[i = 0:tHorz, mod = 0:nMod],
    xAuxGt[i, mod] == aAuxGt[mod] * xLoad[i, mod] + bAuxGt[mod] * y[i, mod]
    )
@constraint(m, powhpdeq[i = 0:tHorz, mod = 0:nMod],
    xPowHpNu[i, mod] == aPowHp[mod] * xLoad[i, mod] + bPowHp[mod] * y[i, mod]
    )
@constraint(m, powipdeq[i = 0:tHorz, mod = 0:nMod],
    xPowIpNu[i, mod] == aPowIp[mod] * xLoad[i, mod] + bPowIp[mod] * y[i, mod]
    )
@constraint(m, pccrebdeq[i = 0:tHorz, mod = 0:nMod],
    xPccReb[i, mod] == aCcReb[mod] * xLoad[i, mod] + bCcReb[mod] * y[i, mod]
    )
@constraint(m, dacsbtdeq[i = 0:tHorz, mod = 0:nMod],
    xDacSb[i, mod] == aDacSb[mod] * xLoad[i, mod] + bDacSb[mod] * y[i, mod]
    )
@constraint(m, allocstdeq[i = 0:tHorz, mod = 0:nMod],
    xAllocS[i, mod] == aSideS[mod] * xLoad[i, mod] + bSideS[mod] * y[i, mod]
    )
@constraint(m, auxstdeq[i = 0:tHorz, mod = 0:nMod],
    xAuxSt[i, mod] == aAuxSt[mod] * xLoad[i, mod] + bAuxSt[mod] * y[i, mod]
    )


# Gas Turbine
@constraint(m, powGasTur[i = 0:tHorz], 
            xPowGasTur[i] == sum(xPowGt[i, mod] for mod in 0:nMod)
           )
# 
@constraint(m, fuelEq[i = 0:tHorz], 
            xFuel[i] == sum(xFuelNu[i, mod] for mod in 0:nMod)
           )
# 
@constraint(m, co2FuelEq[i = 0:tHorz], 
            xCo2Fuel[i] == sum(xEmis[i, mod] for mod in 0:nMod)
           )

@constraint(m, auxPowGasT[i = 0:tHorz],
            xAuxPowGasT[i] == sum(xAuxGt[i, mod] for mod in 0:nMod)
           )
# 
# Steam
# 
@constraint(m, powHpEq[i = 0:tHorz], 
            xPowHp[i] == sum(xPowHpNu[i, mod] for mod in 0:nMod)
           )
# 
@constraint(m, powIpEq[i = 0:tHorz], 
            xPowIp[i] == sum(xPowIpNu[i, mod] for mod in 0:nMod)
           )

# 
@constraint(m, powLpEq[i = 0:tHorz], 
            xPowLp[i] == xSteaPowLp[i] * aLpSteaToPow / 1000
           )
# 
@constraint(m, powerSteaEq[i = 0:tHorz], 
            xPowSteaTur[i] == 
            xPowHp[i] + xPowIp[i] + xPowLp[i]
           )

@constraint(m, ccRebDutyEq[i = 0:tHorz],
            xCcRebDuty[i] == sum(xPccReb[i, mod] for mod in 0:nMod)
           )

@constraint(m, dacSteaDutyEq[i = 0:tHorz],
            xDacSteaBaseDuty[i] == sum(xDacSb[i, mod] for mod in 0:nMod)
           )


@constraint(m, sideSteaEloadEq[i = 0:tHorz],
            xAllocSteam[i] == sum(xAllocS[i, mod] for mod in 0:nMod)
           )

@constraint(m, sideSteaRatioEq[i = 0:tHorz],
            xAllocSteam[i] == xSideSteaDac[i] + xSteaPowLp[i]
           )

@constraint(m, availSteaDacEq[i = 0:tHorz],
            xDacSteaDuty[i] == xDacSteaBaseDuty[i] + xSideSteaDac[i]
           )

@constraint(m, auxPowSteaTEq[i = 0:tHorz],
            xAuxPowSteaT[i] == sum(xAuxSt[i, mod] for mod in 0:nMod)
           )

# PCC
# 
#@constraint(m, co2CapPccEq[i = 0:tHorz - 1], 
#xCo2CapPcc[i] == aCo2PccCapRate * xCo2Fuel[i])
@constraint(m, co2CapPccEq[i = 0:tHorz], 
            xCo2CapPcc[i] == aCapRatePcc * xCo2Fuel[i])
# 
@constraint(m, co2PccOutEq[i = 0:tHorz], 
            xCo2PccOut[i] == xCo2Fuel[i] - xCo2CapPcc[i])
# 
@constraint(m, co2DacFlueInEq[i = 0:tHorz], 
            xCo2DacFlueIn[i] == xCo2PccOut[i] - vCo2PccVent[i])
# 
# @constraint(m, co2CapPccIn[i = 0:tHorz - 1], xCo2CapPcc[i] <= vCapPcc)
# Dav: Sometimes there is not enough steam, so we have to relax this constraint 
@constraint(m, steamUsePccEq[i = 0:tHorz], 
            xSteaUsePcc[i] <= aSteaUseRatePcc * xCo2CapPcc[i])
# 
@constraint(m, powerUsePccEq[i = 0:tHorz], 
            xPowUsePcc[i] == aPowUseRatePcc * xCo2CapPcc[i]
           )

@constraint(m, pccSteaSlack[i = 0:tHorz], 
            xPccSteaSlack[i] == xCcRebDuty[i] - xSteaUsePcc[i])

# DAC-Flue
# Flue gas takes 15 minutes to saturation?
#: "State equation"
@constraint(m, a1dFlueEq[i = 0:tHorz, j=1:hSlice], 
            xA1Flue[i, j] == xA0Flue[i, j-1]
           )

#: "State equation"
@constraint(m, aRdFlueEq[i = 0:tHorz, j=1:hSlice], 
            xR1Flue[i, j] == xR0Flue[i, j-1]
           )
#: "State equation"
@constraint(m, storeFflueeq[i = 0:tHorz, j = 1:hSlice], 
            xFflue[i, j] == xFflue[i, j-1] - xA0Flue[i, j-1] + xR1Flue[i, j-1]
           )
#: "State equation"
@constraint(m, storeSflueeq[i = 0:tHorz, j = 1:hSlice], 
            xSflue[i, j] == xSflue[i, j-1] - xR0Flue[i, j-1] + xA1Flue[i, j-1]
           )
# Initial conditions
@constraint(m, icXa1FlueEq, xA1Flue[0, 0] == 0.)
@constraint(m, icAR1FlueEq, xR1Flue[0, 0] == 0.)
@constraint(m, capDacFlueEq, xFflue[0, 0] == aSorbAmountFreshFlue)
@constraint(m, icSsFlueEq, xSflue[0, 0] == 0.)
# End-point constraints we need to get rid of them and then put them back
@constraint(m, endXa1FlueEq, xA1Flue[tHorz, hSlice] == 0.)
@constraint(m, endAR1FlueEq, xR1Flue[tHorz, hSlice] == 0.)
@constraint(m, endDacFlueEq, xFflue[tHorz, hSlice] == aSorbAmountFreshFlue)
@constraint(m, endSsFlueEq, xSflue[tHorz, hSlice] == 0.)
#
#These dac related variables must start at 0 and end at 1-hslice
@constraint(m, co2CapDacFlueEq[i = 0:tHorz, j = 0:hSlice-1], 
            xCo2CapDacFlue[i, j] == 
            #aSorbCo2CapFlue * xR1Flue[i, j]
            aSorbCo2CapFlue * xA1Flue[i, j]
           )
#
@constraint(m, steamUseDacFlueEq[i = 0:tHorz, j = 0:hSlice-1], 
            xSteaUseDacFlue[i, j] == 
            #aSteaUseRateDacFlue * xCo2CapDacFlue[i, j]
            aSteaUseRateDacFlue * aSorbCo2CapFlue * xR1Flue[i, j]
           )
#
@constraint(m, powUseDacFlueEq[i = 0:tHorz, j = 0:hSlice-1], 
            xPowUseDacFlue[i, j] == aPowUseRateDacFlue * xCo2CapDacFlue[i, j]
           )

# DAC-Air
# Bluntly assume we can just take CO2 from air in pure form.
# "State equation"
@constraint(m, a1dAirEq[i = 0:tHorz, j = 1:hSlice], 
            xA1Air[i, j] == xA0Air[i, j - 1]
           )
# "State equation"
@constraint(m, a2dAirEq[i = 0:tHorz, j = 1:hSlice], 
            xA2Air[i, j] == xA1Air[i, j - 1]
           )
# "State equation"
@constraint(m, aRdAirEq[i = 0:tHorz, j = 1:hSlice], 
            xR1Air[i, j] == xR0Air[i, j - 1]
           )
# "State equation"
@constraint(m, storeFairEq[i = 0:tHorz, j = 1:hSlice], 
            xFair[i, j] == xFair[i, j-1] - xA0Air[i, j-1] + xR1Air[i, j-1]
           )
# "State equation"
@constraint(m, storeSaireq[i = 0:tHorz, j = 1:hSlice], 
            xSair[i, j] == xSair[i, j-1] - xR0Air[i, j-1] + xA2Air[i, j-1]
           )

# Initial conditions - Air
#@constraint(m, capDacAirEq, xFair[0] == xSorbFreshAir)
@constraint(m, capDacAirEq, xFair[0, 0] == aSorbAmountFreshAir)
@constraint(m, icA1AirEq, xA1Air[0, 0] == 0.)
@constraint(m, icA2AirEq, xA2Air[0, 0] == 0.)
@constraint(m, icAR1AirEq, xR1Air[0, 0] == 0.)
@constraint(m, icSsAirEq, xSair[0, 0] == 0.)
# End-point conditions - Air

@constraint(m, endDacAirEq, xFair[tHorz, hSlice] == aSorbAmountFreshAir)
@constraint(m, endA1AirEq, xA1Air[tHorz, hSlice] == 0.)
@constraint(m, endA2AirEq, xA2Air[tHorz, hSlice] == 0.)
@constraint(m, endAR1AirEq, xR1Air[tHorz, hSlice] == 0.)
@constraint(m, endSsAirEq, xSair[tHorz, hSlice] == 0.)

#@constraint(m, endDacAirEq, xFair[0] == xSorbFreshAir)
#@constraint(m, endDacAirEq, xFair[0, 0] == aSorbAmountFreshAir)
#@constraint(m, endA1AirEq, xA1Air[0, 0] == 0.)
#@constraint(m, endA2AirEq, xA2Air[0, 0] == 0.)
#@constraint(m, endAR1AirEq, xR1Air[0, 0] == 0.)
#@constraint(m, endSsAirEq, xSair[0, 0] == 0.)

#
@constraint(m, co2CapDacAirEq[i = 0:tHorz, j=0:hSlice-1], 
            xCo2CapDacAir[i, j] == 
            #aSorbCo2CapAir * xR1Air[i, j]
            (aSorbCo2CapAir * xA1Air[i, j]/2 + aSorbCo2CapAir * xA2Air[i, j]/2)
           )
# 
@constraint(m, steamUseDacAirEq[i = 0:tHorz, j=0:hSlice-1], 
            xSteaUseDacAir[i, j] == 
            #aSteaUseRateDacAir * xCo2CapDacAir[i, j]
            aSteaUseRateDacAir * aSorbCo2CapAir * xR1Air[i, j]
           )
# 
@constraint(m, powUseDacAirEq[i = 0:tHorz, j=0:hSlice-1], 
            xPowUseDacAir[i, j] == aPowUseRateDacAir * xCo2CapDacAir[i, j]
           )
# We need the integral over the whole time for DAC
@constraint(m, dacSteaSlackEq[i = 0:tHorz], 
            xDacSteaSlack[i] == xDacSteaDuty[i] 
            - sum(xSteaUseDacFlue[i, j] for j in 0:hSlice - 1)
            - sum(xSteaUseDacAir[i, j] for j in 0:hSlice - 1)
           )

# Equal to the amount vented, at least in flue mode.
# Integral for DAC
@constraint(m, co2DacFlueVentEq[i = 0:tHorz], 
            xCo2DacVentFlue[i] == xCo2DacFlueIn[i] 
            - sum(xCo2CapDacFlue[i, j] for j in  0:hSlice-1)
           )


# Co2 Compression
# 
# We need the integral over the whole time for DAC
@constraint(m, co2CompEq[i = 0:tHorz], 
            xCo2Comp[i] == xCo2CapPcc[i] 
            + sum(xCo2CapDacFlue[i, j] for j in 0:hSlice-1)
            + sum(xCo2CapDacAir[i, j] for j in 0:hSlice-1)
           )
# 
@constraint(m, powUseCompEq[i = 0:tHorz], 
            xPowUseComp[i] == aPowUseRateComp * xCo2Comp[i]
           )
# 
@constraint(m, co2VentEq[i = 0:tHorz], 
            xCo2Vent[i] == vCo2PccVent[i] 
            + xCo2DacVentFlue[i] - sum(xCo2CapDacAir[i, j] 
                for j in 0:hSlice-1)
           )

## Overall

#
@constraint(m, powGrossEq[i = 0:tHorz], 
            xPowGross[i] == xPowGasTur[i] + xPowSteaTur[i]
           )
@constraint(m, auxPowEq[i = 0:tHorz],
            xAuxPow[i] == xAuxPowGasT[i] + xAuxPowSteaT[i])

@constraint(m, powOutEq[i = 0:tHorz], 
            xPowOut[i] == xPowGross[i] 
            - xPowUsePcc[i]
            - sum(xPowUseDacFlue[i, j] for j in 0:hSlice-1)
            - sum(xPowUseDacAir[i, j] for j in 0:hSlice-1)
            - xPowUseComp[i] 
            - xAuxPow[i]
           )


# Continuity of states
@constraint(m, 
            contxfflue[i = 1:tHorz], xFflue[i, 0] == xFflue[i - 1, hSlice])
@constraint(m, 
            conta1flue[i = 1:tHorz], xA1Flue[i, 0] == xA1Flue[i - 1, hSlice])

@constraint(m, 
            contcxsflue[i = 1:tHorz], xSflue[i, 0] == xSflue[i - 1, hSlice])
@constraint(m, 
            contr1flue[i = 1:tHorz], xR1Flue[i, 0] == xR1Flue[i - 1, hSlice])

@constraint(m, 
            contxfair[i = 1:tHorz], xFair[i, 0] == xFair[i - 1, hSlice])
@constraint(m, 
            conta1air[i = 1:tHorz], xA1Air[i, 0] == xA1Air[i - 1, hSlice])
@constraint(m, 
            conta2air[i = 1:tHorz], xA2Air[i, 0] == xA2Air[i - 1, hSlice])
@constraint(m, 
            contcxsair[i = 1:tHorz], xSair[i, 0] == xSair[i - 1, hSlice])
@constraint(m, 
            contr1air[i = 1:tHorz], xR1Air[i, 0] == xR1Air[i - 1, hSlice])


shutdowncost = cNgPerMmbtu * (bFuelEload + 60 * aFuelEload)
shutdownfromwarm = cNgPerMmbtu * (bFuelEload + 60 * aFuelEload) * 24
startupcost = cNgPerMmbtu * (bFuelEload + 60 * aFuelEload) / 2
# Objective function expression
@expression(m, eObjfExpr, 
    sum(cNgPerMmbtu * xFuel[i]
        + cEmissionPrice * xCo2Vent[i] 
        + cCo2TranspPrice * xCo2Comp[i]
        - pow_price[i + 1] * xPowOut[i]
        for i in 0:tHorz)
        +
    sum(shutdowncost * z[i, 2, 0] + startupcost * z[i, 0, 1]
        + shutdownfromwarm * z[i, 1, 0]
        for i in 0:tHorz)
    )

@objective(m, Min, eObjfExpr)

print("The number of variables\t")
println(num_variables(m))
print("The number of constraints\n")

n = 0
for i in list_of_constraint_types(m)
    global n
    println(num_constraints(m, i[1], i[2]))
    n += num_constraints(m, i[1], i[2])
end


println()

# Set optimizer options
set_optimizer(m, SCIP.Optimizer)
#set_optimizer_attribute(m, "LogLevel", 3)
#set_optimizer_attribute(m, "PresolveType", 1)

optimize!(m)
println(termination_status(m))

#f = open("model.lp", "w")
#print(f, m)
#close(f)

#write_to_file(m, "lp_mk0.mps")
#write_to_file(m, "lp_mk10.lp", format=MOI.FileFormats.FORMAT_LP)

#format::MOI.FileFormats.FileFormat = MOI.FileFormats.FORMAT_AUTOMATIC

# Co2 Data Frame
df_co = DataFrame(Symbol("Co2Fuel") => Float64[], # Pairs.
                  Symbol("Co2CapPcc") => Float64[],
                  Symbol("Co2PccOut") => Float64[],
                  Symbol("vCo2PccVent") => Float64[],
                  Symbol("Co2DacFlueIn") => Float64[],
                  Symbol("Co2CapDacFlue") => Float64[],
                  Symbol("Co2CapDacAir") => Float64[],
                  Symbol("Co2DacVentFlue") => Float64[],
                  Symbol("Co2Vent") => Float64[],
                 )

# Co2 / hSlice
for i in 0:tHorz
    co2fuel = value(xCo2Fuel[i])
    co2pcc = value(xCo2CapPcc[i])
    co2pccout = value(xCo2PccOut[i])
    co2pccvent = value(vCo2PccVent[i])
    co2dacfluein = value(xCo2DacFlueIn[i])
    co2dacflue = sum(value(xCo2CapDacFlue[i, j]) for j in 0:hSlice-1)
    co2dacair = sum(value(xCo2CapDacAir[i, j]) for j in 0:hSlice-1)
    co2dacventflue = value(xCo2DacVentFlue[i])
    co2vent = value(xCo2Vent[i])
    push!(df_co,(
        value(co2fuel),
        value(co2pcc),
        value(co2pccout),
        value(co2pccvent),
        value(co2dacfluein),
        value(co2dacflue),
        value(co2dacair),
        value(co2dacventflue),
        value(co2vent),
        ))
end

# Power Data Frame.
df_pow = DataFrame(
                  Symbol("PowGasTur") => Float64[], # Pairs.
                  Symbol("PowSteaTurb") => Float64[],
                  Symbol("PowHp") => Float64[],
                  Symbol("PowIp") => Float64[],
                  Symbol("PowLp") => Float64[],
                  Symbol("PowUsePcc") => Float64[],
                  Symbol("PowUseDacFlue") => Float64[],
                  Symbol("PowUseDacAir") => Float64[],
                  Symbol("PowUseComp") => Float64[],
                  Symbol("AuxPowGasT") => Float64[],
                  Symbol("AuxPowSteaT") => Float64[],
                  Symbol("PowGross") => Float64[],
                  Symbol("PowOut") => Float64[],
                  Symbol("xeLoad") => Float64[],
                 )


# Pow / hSlice
for i in 0:tHorz
    ygastelecload = value(xeLoad[i])
    powgastur = value(xPowGasTur[i])
    powsteatur = value(xPowSteaTur[i])
    powhp = value(xPowHp[i])
    powip = value(xPowIp[i])
    powlp = value(xPowLp[i])
    powusepcc = value(xPowUsePcc[i])
    powusedacflue = sum(value(xPowUseDacFlue[i, j]) for j in 0:hSlice-1)
    powusedacair = sum(value(xPowUseDacAir[i, j]) for j in 0:hSlice-1)
    powusecomp = value(xPowUseComp[i])
    auxpowgast = value(xAuxPowGasT[i])
    auxpowsteat = value(xAuxPowSteaT[i])
    powgross = value(xPowGross[i])
    powout = value(xPowOut[i])
    push!(df_pow, (
        powgastur,
        powsteatur,
        powhp,
        powip,
        powlp,
        powusepcc,
        powusedacflue,
        powusedacair,
        powusecomp,
        auxpowgast,
        auxpowsteat,
        powgross,
        powout,
        ygastelecload
        ))
end


# Steam DataFrame
df_steam = DataFrame(
                     Symbol("CcRebDuty") => Float64[],
                     Symbol("SteaUsePcc") => Float64[],
                     Symbol("PccSteaSlack") => Float64[],
                     Symbol("DacSteaDuty") => Float64[],
                     Symbol("SteaUseDacFlue") => Float64[],
                     Symbol("SteaUseDacAir") => Float64[],
                     Symbol("DacSteaSlack") => Float64[],
                     Symbol("SideStea") => Float64[],
                     Symbol("DacSteaBaseDuty") => Float64[],
                     Symbol("SideSteaDac") => Float64[],
                     Symbol("Fuel") => Float64[]
                    )

# Steam / hSlice
for i in 0:tHorz
    ccrebduty = value(xCcRebDuty[i])
    steausepcc = value(xSteaUsePcc[i])
    pccsteaslack = value(xPccSteaSlack[i])
    dacsteaduty = value(xDacSteaDuty[i])
    steausedacflue = sum(value(xSteaUseDacFlue[i, j]) for j in 0:hSlice-1)
    steausedacair = sum(value(xSteaUseDacAir[i, j]) for j in 0:hSlice-1)
    dacsteaslack = value(xDacSteaSlack[i])
    sidesteam = value(xAllocSteam[i])
    dacsteabaseduty = value(xDacSteaBaseDuty[i])
    sidesteadac = value(xSideSteaDac[i])
    xfuel = value(xFuel[i])
    push!(df_steam, 
        (
            ccrebduty,
            steausepcc,
            pccsteaslack,
            dacsteaduty,
            steausedacflue,
            steausedacair,
            dacsteaslack,
            sidesteam,
            dacsteabaseduty,
            sidesteadac,
            xfuel
            ))
end

# DAC-flue DataFrame
df_dac_flue = DataFrame(
    :time => Float64[],
    :xFflue => Float64[],
    :xSflue => Float64[],
    :xA0Flue => Float64[],
    :xA1Flue => Float64[],
    :xR0Flue => Float64[],
    :xR1Flue => Float64[]
    )

sliceFact = 1/hSlice
for i in 0:tHorz
    for j in 0:hSlice-1
        currtime = i + j * sliceFact
        push!(df_dac_flue,(
            currtime,
            value(xFflue[i, j]), 
            value(xSflue[i, j]),
            value(xA0Flue[i, j]), 
            value(xA1Flue[i, j]), 
            value(xR0Flue[i, j]), 
            value(xR1Flue[i, j]))
        )
    end
end

# DAC-air DataFrame
df_dac_air = DataFrame(
    :time => Float64[],
    :xFair => Float64[],
    :xSair => Float64[],
    :xA0Air => Float64[],
    :xA1Air => Float64[],
    :xA2Air => Float64[],
    :xR0Air => Float64[],
    :xR1Air => Float64[])

for i in 0:tHorz
    for j in 0:hSlice-1
        currtime = i + j * sliceFact
        push!(df_dac_air,(
            currtime,
            value(xFair[i, j]), 
            value(xSair[i, j]),
            value(xA0Air[i, j]), 
            value(xA1Air[i, j]), 
            value(xA2Air[i, j]),
            value(xR0Air[i, j]), 
            value(xR1Air[i, j])))
    end
end


df_pow_price = DataFrame(
                         price = Float64[]
                        )
for i in 0:tHorz
    push!(df_pow_price, (pow_price[i + 1],))
end



# Cost DataFrame
df_cost = DataFrame(
                    cNG = Float64[],
                    cCo2Em = Float64[],
                    cTransp = Float64[],
                    PowSales = Float64[]
                   )

for i in 0:tHorz
    cng = cNgPerMmbtu * value(xFuel[i])
    cco = cEmissionPrice * value(xCo2Vent[i])
    ctr = cCo2TranspPrice * value(xCo2Comp[i])
    cpow = pow_price[i + 1] * value(xPowOut[i])
    push!(df_cost, (cng, cco, ctr, cpow))
end

df_time_slice = DataFrame(:time => Float64[])

for i in 0:tHorz
    for j in 0:hSlice-1
        currtime = i + j * sliceFact
        push!(df_time_slice, (currtime, ))
    end
end


df_binary = DataFrame(
    :y0 => Float64[],
    :y1 => Float64[],
    :y2 => Float64[],
    :z00 => Float64[],
    :z01 => Float64[],
    :z02 => Float64[],
    :z10 => Float64[],
    :z11 => Float64[],
    :z12 => Float64[],
    :z20 => Float64[],
    :z21 => Float64[],
    :z22 => Float64[],
    )

for i in 0:tHorz
    push!(df_binary, 
        (
            value(y[i, 0]),
            value(y[i, 1]),
            value(y[i, 2]),
            value(z[i, 0, 0]),
            value(z[i, 0, 1]),
            value(z[i, 0, 2]),
            value(z[i, 1, 0]),
            value(z[i, 1, 1]),
            value(z[i, 1, 2]),
            value(z[i, 2, 0]),
            value(z[i, 2, 1]),
            value(z[i, 2, 2]),
        ))
end

df_pr = DataFrame(
    :load0 => Float64[], 
    :load1 => Float64[], 
    :load2 => Float64[], 
    :lambda0 => Float64[], 
    :lambda10 => Float64[],
    :lambda11 => Float64[],
    :lambda20 => Float64[],
    :lambda21 => Float64[]
    )

for i in 0:tHorz
    push!(df_pr, (
        value(xLoad[i, 0]), 
        value(xLoad[i, 1]), 
        value(xLoad[i, 2]),
        value(lambda[i, 0, 0]), 
        value(lambda[i, 1, 0]),
        value(lambda[i, 1, 1]),
        value(lambda[i, 2, 0]),
        value(lambda[i, 2, 1]),
        ))
end

# Write CSV
CSV.write("df_co.csv", df_co)
CSV.write("df_pow.csv", df_pow)
CSV.write("df_steam.csv", df_steam)
CSV.write("df_dac_flue.csv", df_dac_flue)
CSV.write("df_dac_air.csv", df_dac_air)
CSV.write("df_pow_price.csv", df_pow_price)
CSV.write("df_cost.csv", df_cost)
CSV.write("df_time_slice.csv", df_time_slice)
CSV.write("df_binary.csv", df_binary)
CSV.write("df_pr.csv", df_pr)

df_diss = DataFrame(
    "load0" => Float64[],
    "load1" => Float64[],
    "load2" => Float64[],
    "powgt0" => Float64[],
    "powgt1" => Float64[],
    "powgt2" => Float64[],
    "fuel0" => Float64[],
    "fuel1" => Float64[],
    "fuel2" => Float64[],
    "emiss0" => Float64[],
    "emiss1" => Float64[],
    "emiss2" => Float64[],
    "auxgt0" => Float64[],
    "auxgt1" => Float64[],
    "auxgt2" => Float64[],
    "powhp0" => Float64[],
    "powhp1" => Float64[],
    "powhp2" => Float64[],
    "powip0" => Float64[],
    "powip1" => Float64[],
    "powip2" => Float64[],
    "pccreb0" => Float64[],
    "pccreb1" => Float64[],
    "pccreb2" => Float64[],
    "dacsb0" => Float64[],
    "dacsb1" => Float64[],
    "dacsb2" => Float64[],
    "allocst0" => Float64[],
    "allocst1" => Float64[],
    "allocst2" => Float64[],
    "auxst0" => Float64[],
    "auxst1" => Float64[],
    "auxst2" => Float64[])

dvars = ["xLoad", "xPowGt", "xFuelNu", "xEmis", "xAuxGt", "xPowHpNu", "xPowIpNu", "xPccReb", "xDacSb", "xAllocS", "xAuxSt"]
for i in 0:tHorz
    l0 = []
    for v in dvars
        s = Symbol(v)
        l = [value(m[s][i, mod]) for mod in 0:nMod]
        l0 = vcat(l0, l)
    end
    push!(df_diss, l0)
end

CSV.write("df_diss.csv", df_diss)
