diff --git a/src/calculate_emissions.jl b/src/calculate_emissions.jl index 49c4e3f..5ebbb19 100644 --- a/src/calculate_emissions.jl +++ b/src/calculate_emissions.jl @@ -10,12 +10,18 @@ end function calculate_ACE(data::Dict{String,<:Any}, optimizer) # data checks @assert _contains_generator_emissions(data) + @assert _all_gen_costs_single_type(data) - #PowerModels.standardize_cost_terms!(data, order=2) PowerModels.calc_thermal_limits!(data) + cost_type = collect(values(data["gen"]))[1]["model"] + cost_type = ["PWL","POLY"][cost_type] + + #build ref, solve dcopf ref = PowerModels.build_ref(data)[:it][:pm][:nw][0] - model = _solve_dcopf(ref, optimizer) + model = _solve_dcopf(ref, optimizer,cost_type=cost_type) + update_ref!(ref,model) + total_system_emissions = sum(JuMP.value(model[:pg][id])*gen["emissions"] for (id,gen) in ref[:gen]) #get sum of all system emissions avg_system_emissions = total_system_emissions/sum(JuMP.value(model[:pg][id]) for (id,gen) in ref[:gen]) #get average system emissions @@ -94,12 +100,16 @@ end function calculate_LACE(data::Dict{String,<:Any}, optimizer) @assert _contains_generator_emissions(data) @assert _no_duplicate_costs_at_bus(data) + @assert _all_gen_costs_single_type(data) - #PowerModels.standardize_cost_terms!(data, order=2) PowerModels.calc_thermal_limits!(data) + #get generator cost type + cost_type = collect(values(data["gen"]))[1]["model"] + cost_type = ["PWL","POLY"][cost_type] + #build ref, solve dcopf ref = PowerModels.build_ref(data)[:it][:pm][:nw][0] - model = _solve_dcopf(ref, optimizer) + model = _solve_dcopf(ref, optimizer,cost_type=cost_type) update_ref!(ref,model) @@ -221,105 +231,206 @@ function calculate_LMCE(file::String, optimizer) calculate_LMCE(data, optimizer) end - -# function calc_LMCE(data;return_delta_pg=false,adjusted_LMCE = false,account=false) function calculate_LMCE(data::Dict{String,<:Any}, optimizer) # data checks @assert _contains_generator_emissions(data) @assert _no_duplicate_costs_at_bus(data) + @assert _all_gen_costs_single_type(data) - #PowerModels.standardize_cost_terms!(data, order=2) - PowerModels.calc_thermal_limits!(data) + #grab first gen cost type, asserted to be same as all others + gencost_type = collect(values(data["gen"]))[1]["model"] - ref = PowerModels.build_ref(data)[:it][:pm][:nw][0] - model = _solve_dcopf(ref, optimizer) + if gencost_type == 1 #PWL cost + PowerModels.calc_thermal_limits!(data) - ## Calculate ΔPg ## - #initialize A matrix to be NxN where n = N_b + N_g (number of buses plus number of generators) - A = assemble_A(model) - A_inv = inv(A) + ref = PowerModels.build_ref(data)[:it][:pm][:nw][0] + model = _solve_dcopf(ref, optimizer,cost_type="PWL") - #grab relevant portion of A_inv matrix (Last num_gens rows, first num_bus columns of A_inv) - Nb = length(model[:va]) - Ng = length(model[:pg]) - B = A_inv[Nb+1:Ng+Nb,1:Nb] + ## Calculate ΔPg ## + #initialize A matrix to be NxN where n = N_b + N_g (number of buses plus number of generators) + A = assemble_A_pwl(model) + A_inv = inv(A) - ## Calculate lmce by multiplying B (generator response) by generator emissions intensity (emissions/MWh) - sorted_gen_emission_factors = [ref[:gen][id]["emissions"] for id in sort(collect(keys(ref[:gen])))] - for r in eachcol(B) - r .= r .* sorted_gen_emission_factors + #grab relevant portion of A_inv matrix (Last num_gens rows, first num_bus columns of A_inv) + Nb = length(model[:va]) + Ng = length(model[:pg]) + B = A_inv[Nb+1:Ng+Nb,1:Nb] + + ## Calculate lmce by multiplying B (generator response) by generator emissions intensity (emissions/MWh) + sorted_gen_emission_factors = [ref[:gen][id]["emissions"] for id in sort(collect(keys(ref[:gen])))] + for r in eachcol(B) + r .= r .* sorted_gen_emission_factors + end + + bus_ids = sort(collect(keys(ref[:bus]))) + bus_lmce = vec(sum(B,dims=1)) + return LMCE(bus_ids, bus_lmce) + + elseif gencost_type==2 #polynomial cost + PowerModels.calc_thermal_limits!(data) + ref = PowerModels.build_ref(data)[:it][:pm][:nw][0] + model = _solve_dcopf(ref, optimizer,cost_type="POLY") + + ## Calculate ΔPg ## + #initialize A matrix to be NxN where n = N_b + N_g (number of buses plus number of generators) + A = assemble_A_poly(model) + A_inv = inv(A) + + #grab relevant portion of A_inv matrix (Last num_gens rows, first num_bus columns of A_inv) + Nb = length(model[:va]) + Ng = length(model[:pg]) + B = A_inv[end-Ng+1:end,1:Nb] + + ## Calculate lmce by multiplying B (generator response) by generator emissions intensity (emissions/MWh) + sorted_gen_emission_factors = [ref[:gen][id]["emissions"] for id in sort(collect(keys(ref[:gen])))] + for r in eachcol(B) + r .= r .* sorted_gen_emission_factors + end + + bus_ids = sort(collect(keys(ref[:bus]))) + bus_lmce = vec(sum(B,dims=1)) + return LMCE(bus_ids, bus_lmce) + else + error("Unsupported generator cost type $(gencost_type) detected") + return end - bus_ids = sort(collect(keys(ref[:bus]))) - bus_lmce = vec(sum(B,dims=1)) - return LMCE(bus_ids, bus_lmce) + end + function calculate_ALMCE(file::String, optimizer) data = PowerModels.parse_file(file) calculate_ALMCE(data, optimizer) end -# function calc_LMCE(data;return_delta_pg=false,adjusted_LMCE = false,account=false) + function calculate_ALMCE(data::Dict{String,<:Any}, optimizer) # data checks @assert _contains_generator_emissions(data) @assert _no_duplicate_costs_at_bus(data) + @assert _all_gen_costs_single_type(data) - #PowerModels.standardize_cost_terms!(data, order=2) - PowerModels.calc_thermal_limits!(data) + #grab first gen cost type, asserted to be same as all others + gencost_type = collect(values(data["gen"]))[1]["model"] - ref = PowerModels.build_ref(data)[:it][:pm][:nw][0] - model = _solve_dcopf(ref, optimizer) + if gencost_type==1 #PWL + PowerModels.calc_thermal_limits!(data) - ## Calculate ΔPg ## - #initialize A matrix to be NxN where n = N_b + N_g (number of buses plus number of generators) - A = assemble_A(model) - A_inv = inv(A) + ref = PowerModels.build_ref(data)[:it][:pm][:nw][0] + model = _solve_dcopf(ref, optimizer) - #grab relevant portion of A_inv matrix (Last num_gens rows, first num_bus columns of A_inv) - Nb = length(model[:va]) - Ng = length(model[:pg]) - B = B = A_inv[Nb+1:Ng+Nb,1:Nb] + ## Calculate ΔPg ## + #initialize A matrix to be NxN where n = N_b + N_g (number of buses plus number of generators) + A = assemble_A_pwl(model) + A_inv = inv(A) - ## Calculate emissions/MWh - sorted_gen_emission_factors = [ref[:gen][id]["emissions"] for id in sort(collect(keys(ref[:gen])))] - for r in eachcol(B) - r .= r .* sorted_gen_emission_factors - end + #grab relevant portion of A_inv matrix (Last num_gens rows, first num_bus columns of A_inv) + Nb = length(model[:va]) + Ng = length(model[:pg]) + B = A_inv[Nb+1:Ng+Nb,1:Nb] - bus_ids = sort(collect(keys(ref[:bus]))) - bus_lmce = vec(sum(B,dims=1)) - lmce = LMCE(bus_ids, bus_lmce) - - # adusted LMCE - sorted_bus_list = sort(collect(keys(ref[:bus_loads]))) #grab list of buses - sorted_bus_loads = zeros(length(sorted_bus_list)) #initialize zeros for total load at each bus - - #loop through buses and their constituent loads to get total system load - for (bus_idx,bus) in enumerate(sorted_bus_list) - bus_loads = ref[:bus_loads][bus] #grab list of loads at bus - for load in bus_loads - sorted_bus_loads[bus_idx] += ref[:load][load]["pd"] #add real power demand for load to corresponding bus + ## Calculate emissions/MWh + sorted_gen_emission_factors = [ref[:gen][id]["emissions"] for id in sort(collect(keys(ref[:gen])))] + for r in eachcol(B) + r .= r .* sorted_gen_emission_factors + end + + bus_ids = sort(collect(keys(ref[:bus]))) + bus_lmce = vec(sum(B,dims=1)) + lmce = LMCE(bus_ids, bus_lmce) + + # adusted LMCE + sorted_bus_list = sort(collect(keys(ref[:bus_loads]))) #grab list of buses + sorted_bus_loads = zeros(length(sorted_bus_list)) #initialize zeros for total load at each bus + + #loop through buses and their constituent loads to get total system load + for (bus_idx,bus) in enumerate(sorted_bus_list) + bus_loads = ref[:bus_loads][bus] #grab list of loads at bus + for load in bus_loads + sorted_bus_loads[bus_idx] += ref[:load][load]["pd"] #add real power demand for load to corresponding bus + end end - end - #grab lmce-assigned emissions and total system load: - lmce_bus_emissions = sorted_bus_loads.*lmce.emissions_intensity - total_system_load = sum(sorted_bus_loads) + #grab lmce-assigned emissions and total system load: + lmce_bus_emissions = sorted_bus_loads.*lmce.emissions_intensity + total_system_load = sum(sorted_bus_loads) - #now need to get current generator set points and calculate total system emissions - sorted_pg_idx = sort(axes(model[:pg])[1]) - sorted_pg_vars = Array(model[:pg][sorted_pg_idx]) - gen_set_points = JuMP.value.(sorted_pg_vars) - total_system_emissions = sum(gen_set_points.*sorted_gen_emission_factors) #get sum of all system emissions + #now need to get current generator set points and calculate total system emissions + sorted_pg_idx = sort(axes(model[:pg])[1]) + sorted_pg_vars = Array(model[:pg][sorted_pg_idx]) + gen_set_points = JuMP.value.(sorted_pg_vars) + total_system_emissions = sum(gen_set_points.*sorted_gen_emission_factors) #get sum of all system emissions - #calculate emissions discrepancy with LMCE-assigned emissions - emissions_discrep = total_system_emissions - sum(lmce_bus_emissions) + #calculate emissions discrepancy with LMCE-assigned emissions + emissions_discrep = total_system_emissions - sum(lmce_bus_emissions) + + LMCE_adjustment_factor = emissions_discrep/total_system_load + + adjusted_LMCE = lmce.emissions_intensity .+LMCE_adjustment_factor + return ALMCE(sorted_bus_list, adjusted_LMCE) + + elseif gencost_type==2 #Polynomial cost + PowerModels.calc_thermal_limits!(data) + ref = PowerModels.build_ref(data)[:it][:pm][:nw][0] + model = _solve_dcopf(ref, optimizer,cost_type="POLY") + + ## Calculate ΔPg ## + #initialize A matrix to be NxN where n = N_b + N_g (number of buses plus number of generators) + A = assemble_A_poly(model) + A_inv = inv(A) + + #grab relevant portion of A_inv matrix (Last num_gens rows, first num_bus columns of A_inv) + Nb = length(model[:va]) + Ng = length(model[:pg]) + B = A_inv[end-Ng+1:end,1:Nb] + + ## Calculate lmce by multiplying B (generator response) by generator emissions intensity (emissions/MWh) + sorted_gen_emission_factors = [ref[:gen][id]["emissions"] for id in sort(collect(keys(ref[:gen])))] + for r in eachcol(B) + r .= r .* sorted_gen_emission_factors + end - LMCE_adjustment_factor = emissions_discrep/total_system_load + bus_ids = sort(collect(keys(ref[:bus]))) + bus_lmce = vec(sum(B,dims=1)) + lmce = LMCE(bus_ids, bus_lmce) + + # adusted LMCE + sorted_bus_list = sort(collect(keys(ref[:bus_loads]))) #grab list of buses + sorted_bus_loads = zeros(length(sorted_bus_list)) #initialize zeros for total load at each bus + + #loop through buses and their constituent loads to get total system load + for (bus_idx,bus) in enumerate(sorted_bus_list) + bus_loads = ref[:bus_loads][bus] #grab list of loads at bus + for load in bus_loads + sorted_bus_loads[bus_idx] += ref[:load][load]["pd"] #add real power demand for load to corresponding bus + end + end + + #grab lmce-assigned emissions and total system load: + lmce_bus_emissions = sorted_bus_loads.*lmce.emissions_intensity + total_system_load = sum(sorted_bus_loads) + + #now need to get current generator set points and calculate total system emissions + sorted_pg_idx = sort(axes(model[:pg])[1]) + sorted_pg_vars = Array(model[:pg][sorted_pg_idx]) + gen_set_points = JuMP.value.(sorted_pg_vars) + total_system_emissions = sum(gen_set_points.*sorted_gen_emission_factors) #get sum of all system emissions + + #calculate emissions discrepancy with LMCE-assigned emissions + emissions_discrep = total_system_emissions - sum(lmce_bus_emissions) + + LMCE_adjustment_factor = emissions_discrep/total_system_load + + adjusted_LMCE = lmce.emissions_intensity .+LMCE_adjustment_factor + return ALMCE(sorted_bus_list, adjusted_LMCE) + + + else + error("Unsupported generator cost type $(gencost_type) detected") + return + + end - adjusted_LMCE = lmce.emissions_intensity .+LMCE_adjustment_factor - return ALMCE(sorted_bus_list, adjusted_LMCE) end \ No newline at end of file diff --git a/src/utils.jl b/src/utils.jl index 3bead94..10ba0d3 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -29,8 +29,22 @@ function _no_duplicate_costs_at_bus(data) return no_duplicate end +function _all_gen_costs_single_type(data) + same_type = true + first_type = collect(values(data["gen"]))[1]["model"] + for (id,gen) in data["gen"] + if !(gen["model"]==first_type) + print("Generators have multiple cost function types, not supported.") + same_type=false + end + end + + return same_type + +end + -function _build_dcopf(ref::Dict{Symbol,<:Any}, optimizer,pwl=true) +function _build_dcopf(ref::Dict{Symbol,<:Any}, optimizer;cost_type="PWL") model = JuMP.Model() JuMP.set_optimizer(model, optimizer) @@ -45,7 +59,7 @@ function _build_dcopf(ref::Dict{Symbol,<:Any}, optimizer,pwl=true) # sum(gen["cost"][2]*pg[i] for (i,gen) in ref[:gen]) # ) - if pwl + if cost_type=="PWL" #create dummy cost variable for each generator. Will be constrained such that cg is higher than all lines formed by extending PWL segments (assumes convexity) JuMP.@variable(model,cg[i in keys(ref[:gen])]) #constrain cg variables @@ -78,14 +92,16 @@ function _build_dcopf(ref::Dict{Symbol,<:Any}, optimizer,pwl=true) #and finally, set objective to be minimization of sum of cg variables JuMP.@objective(model, Min, sum(cg)) - else - #approximation of PWL, drawing line from first point to last point, right now tuned for RTS-GMLC system - #..assumes there are at least two points that define the PWL cost + + elseif cost_type=="POLY" + #no support for quadratic terms, use only linear and constant + #this assumes that these are the final two terms in a polynomial cost fn JuMP.@objective(model, Min, - sum( - ((gen["cost"][end] - gen["cost"][2])/(gen["cost"][end-1]-gen["cost"][1]))*pg[i] for (i,gen) in ref[:gen] - ) + sum(gen["cost"][end-1]*pg[i] + gen["cost"][end] for (i,gen) in ref[:gen]) ) + + else + error("Unsupported cost function type $(cost_type) provided") end # Add Constraints @@ -131,8 +147,8 @@ function _build_dcopf(ref::Dict{Symbol,<:Any}, optimizer,pwl=true) end -function _solve_dcopf(data, optimizer) - model = _build_dcopf(data, optimizer) +function _solve_dcopf(data, optimizer;cost_type="PWL") + model = _build_dcopf(data, optimizer,cost_type=cost_type) JuMP.optimize!(model) return model end @@ -219,7 +235,54 @@ function collect_binding_inequality_constraints(model::JuMP.AbstractModel) end -function assemble_A(model::JuMP.AbstractModel) +function assemble_A_poly(model::JuMP.AbstractModel) + Nb = length(model[:va]) + Ng = length(model[:pg]) + n = Nb + Ng #voltage angle and pg variables + A = zeros(n,n) + + all_binding_constraints = collect_binding_constraints(model, Nb) #Nb passed in in order to identify the nodal power balance constraints + + #println("There are $(length(all_binding_constraints)) total binding constraints") + + i = 1 + for con in all_binding_constraints + #initialize row + temp_row = zeros(n) + + #grab the voltage angle and pg vars in order + sorted_va_idx = sort(axes(model[:va])[1]) + sorted_pg_idx = sort(axes(model[:pg])[1]) + + va_vars = Array(model[:va][sorted_va_idx]) + pg_vars = Array(model[:pg][sorted_pg_idx]) + + #and combine them into one list + va_pg_vars = [va_vars;pg_vars] + + if typeof(con)==JuMP.VariableRef #if constraint is just variable ref, then row in A should just have a 1 where that variable is + var_idx = findfirst(==(con),va_pg_vars) + temp_row[var_idx]=1 + else #otherwise if constraint is proper funtion (equality constraint or non-variable-definition inequality) + #now grab the constraint function that we're encoding + con_func_terms = con.terms + #and loop through the variables in va_pg_vars, checking whether they're in the con_func (and if so, grabbing their coefficient) + for (j,var) in enumerate(va_pg_vars) + if var in keys(con_func_terms) + temp_row[j] = con_func_terms[var] + end + end + end + + #make temp row permanent + A[i,:] = temp_row + i = i+1 + end + return A +end + + +function assemble_A_pwl(model::JuMP.AbstractModel) Nb = length(model[:va]) Ng = length(model[:pg]) n = Nb + Ng + Ng #2x Ng because one cg variable for each gen as well as one pg variable