diff --git a/src/ac_power_flow.jl b/src/ac_power_flow.jl index 617683cf..6daf57c7 100644 --- a/src/ac_power_flow.jl +++ b/src/ac_power_flow.jl @@ -55,16 +55,15 @@ function _calculate_x0(n::Int, return x0 end -function PolarPowerFlow(data::ACPowerFlowData) - time_step = 1 # TODO placeholder time_step +function PolarPowerFlow(data::ACPowerFlowData; time_step::Int64 = 1) n_buses = first(size(data.bus_type)) P_net = zeros(n_buses) Q_net = zeros(n_buses) for ix in 1:n_buses P_net[ix] = - data.bus_activepower_injection[ix] - data.bus_activepower_withdrawals[ix] + data.bus_activepower_injection[ix, time_step] - data.bus_activepower_withdrawals[ix, time_step] Q_net[ix] = - data.bus_reactivepower_injection[ix] - data.bus_reactivepower_withdrawals[ix] + data.bus_reactivepower_injection[ix, time_step] - data.bus_reactivepower_withdrawals[ix, time_step] end x0 = _calculate_x0(1, data.bus_type[:, time_step], diff --git a/src/common.jl b/src/common.jl index 28a095ec..cd2d1fb3 100644 --- a/src/common.jl +++ b/src/common.jl @@ -199,7 +199,7 @@ function make_powerflowdata( bus_reactivepower_injection_1 = zeros_bus_time() bus_activepower_withdrawals_1 = zeros_bus_time() bus_reactivepower_withdrawals_1 = zeros_bus_time() - bus_reactivepower_bounds_1 = Matrix{Vector{Float64}}(undef, n_buses, time_stepsm) + bus_reactivepower_bounds_1 = Matrix{Vector{Float64}}(undef, n_buses, time_steps) bus_magnitude_1 = ones_bus_time() bus_angles_1 = zeros_bus_time() diff --git a/src/newton_ac_powerflow.jl b/src/newton_ac_powerflow.jl index 1c2a1156..f2eaefc5 100644 --- a/src/newton_ac_powerflow.jl +++ b/src/newton_ac_powerflow.jl @@ -1,4 +1,3 @@ -const _SOLVE_AC_POWERFLOW_KWARGS = Set([:check_reactive_power_limits, :check_connectivity]) """ Solves a the power flow into the system and writes the solution into the relevant structs. Updates generators active and reactive power setpoints and branches active and reactive @@ -34,6 +33,7 @@ solve_ac_powerflow!(sys, method=:newton) function solve_powerflow!( pf::ACPowerFlow{<:ACPowerFlowSolverType}, system::PSY.System; + time_step::Int64=1, kwargs..., ) #Save per-unit flag @@ -47,8 +47,7 @@ function solve_powerflow!( check_connectivity = get(kwargs, :check_connectivity, true), ) - solver_kwargs = filter(p -> !(p.first in _SOLVE_AC_POWERFLOW_KWARGS), kwargs) - converged, V, Sbus_result = _ac_powereflow(data, pf; solver_kwargs...) + converged, V, Sbus_result = _ac_powereflow(data, pf; time_step=time_step, kwargs...) x = _calc_x(data, V, Sbus_result) if converged @@ -124,8 +123,8 @@ function solve_powerflow( sorted_time_steps = sort(collect(keys(data.timestep_map))) # preallocate results ts_converged = zeros(Bool, 1, length(sorted_time_steps)) - ts_V = zeros(Complex{Float64}, length(data.bus_type), length(sorted_time_steps)) - ts_S = zeros(Complex{Float64}, length(data.bus_type), length(sorted_time_steps)) + ts_V = zeros(Complex{Float64}, length(data.bus_type[:, 1]), length(sorted_time_steps)) + ts_S = zeros(Complex{Float64}, length(data.bus_type[:, 1]), length(sorted_time_steps)) results = Dict() @@ -176,8 +175,8 @@ function solve_powerflow!( sorted_time_steps = sort(collect(keys(data.timestep_map))) # preallocate results ts_converged = zeros(Bool, 1, length(sorted_time_steps)) - ts_V = zeros(Complex{Float64}, length(data.bus_type), length(sorted_time_steps)) - ts_S = zeros(Complex{Float64}, length(data.bus_type), length(sorted_time_steps)) + ts_V = zeros(Complex{Float64}, length(data.bus_type[:, 1]), length(sorted_time_steps)) + ts_S = zeros(Complex{Float64}, length(data.bus_type[:, 1]), length(sorted_time_steps)) Yft = data.power_network_matrix.data_ft Ytf = data.power_network_matrix.data_tf @@ -192,9 +191,9 @@ function solve_powerflow!( ts_V[:, t] .= V ts_S[:, t] .= Sbus_result - ref = findall(x -> x == PowerSystems.ACBusTypesModule.ACBusTypes.REF, data.bus_type) - pv = findall(x -> x == PowerSystems.ACBusTypesModule.ACBusTypes.PV, data.bus_type) - pq = findall(x -> x == PowerSystems.ACBusTypesModule.ACBusTypes.PQ, data.bus_type) + ref = findall(x -> x == PowerSystems.ACBusTypesModule.ACBusTypes.REF, data.bus_type[:, t]) + pv = findall(x -> x == PowerSystems.ACBusTypesModule.ACBusTypes.PV, data.bus_type[:, t]) + pq = findall(x -> x == PowerSystems.ACBusTypesModule.ACBusTypes.PQ, data.bus_type[:, t]) # temporary implementation that will need to be improved: if converged @@ -222,7 +221,13 @@ function solve_powerflow!( end # write branch flows - data.branch_flow_values .= real.(ts_V[fb,:] .* conj.(Yft * ts_V)) + Sft = ts_V[fb,:] .* conj.(Yft * ts_V) + Stf = ts_V[tb,:] .* conj.(Ytf * ts_V) + + data.branch_activepower_flow_from_to .= real.(Sft) + data.branch_reactivepower_flow_from_to .= imag.(Sft) + data.branch_activepower_flow_to_from .= real.(Stf) + data.branch_reactivepower_flow_to_from .= imag.(Stf) # todo: # return df_results @@ -243,7 +248,7 @@ function _ac_powereflow( converged, V, Sbus_result = _newton_powerflow(pf, data; time_step = time_step, kwargs...) if !converged || !check_reactive_power_limits || - _check_q_limit_bounds!(data, Sbus_result) + _check_q_limit_bounds!(data, Sbus_result, time_step) return (converged, V, Sbus_result) end end @@ -259,8 +264,8 @@ function _check_q_limit_bounds!( ) bus_names = data.power_network_matrix.axes[1] within_limits = true - for (ix, b) in enumerate(data.bus_type) - if b == PSY.ACBusTypes.PV + for (ix, bt) in enumerate(data.bus_type[:, time_step]) + if bt == PSY.ACBusTypes.PV Q_gen = imag(Sbus_result[ix]) else continue @@ -271,10 +276,10 @@ function _check_q_limit_bounds!( within_limits = false data.bus_type[ix, time_step] = PSY.ACBusTypes.PQ data.bus_reactivepower_injection[ix, time_step] = data.bus_reactivepower_bounds[ix, time_step][1] - elseif Q_gen >= data.bus_reactivepower_bounds[ix][2] + elseif Q_gen >= data.bus_reactivepower_bounds[ix, time_step][2] @info "Bus $(bus_names[ix]) changed to PSY.ACBusTypes.PQ" within_limits = false - data.bus_type[ixm, time_step] = PSY.ACBusTypes.PQ + data.bus_type[ix, time_step] = PSY.ACBusTypes.PQ data.bus_reactivepower_injection[ix, time_step] = data.bus_reactivepower_bounds[ix, time_step][2] else @debug "Within Limits" @@ -469,23 +474,24 @@ end function _calc_x( data::ACPowerFlowData, V::Vector{Complex{Float64}}, - Sbus_result::Vector{Complex{Float64}}, + Sbus_result::Vector{Complex{Float64}}; + time_step::Int64 = 1, ) Vm = abs.(V) Va = angle.(V) n_buses = length(V) # mock the expected x format, where the values depend on the type of the bus: x = zeros(Float64, 2 * n_buses) - for (ix, b) in enumerate(data.bus_type) - if b == PSY.ACBusTypes.REF + for (ix, bt) in enumerate(data.bus_type[:, time_step]) + if bt == PSY.ACBusTypes.REF # When bustype == REFERENCE PSY.Bus, state variables are Active and Reactive Power Generated x[2 * ix - 1] = real(Sbus_result[ix]) + data.bus_activepower_withdrawals[ix] x[2 * ix] = imag(Sbus_result[ix]) + data.bus_reactivepower_withdrawals[ix] - elseif b == PSY.ACBusTypes.PV + elseif bt == PSY.ACBusTypes.PV # When bustype == PV PSY.Bus, state variables are Reactive Power Generated and Voltage Angle x[2 * ix - 1] = imag(Sbus_result[ix]) + data.bus_reactivepower_withdrawals[ix] x[2 * ix] = Va[ix] - elseif b == PSY.ACBusTypes.PQ + elseif bt == PSY.ACBusTypes.PQ # When bustype == PQ PSY.Bus, state variables are Voltage Magnitude and Voltage Angle x[2 * ix - 1] = Vm[ix] x[2 * ix] = Va[ix] @@ -519,16 +525,15 @@ function _newton_powerflow( Ybus = data.power_network_matrix.data # Find indices for each bus type - ref = findall(x -> x == PowerSystems.ACBusTypesModule.ACBusTypes.REF, data.bus_type) - pv = findall(x -> x == PowerSystems.ACBusTypesModule.ACBusTypes.PV, data.bus_type) - pq = findall(x -> x == PowerSystems.ACBusTypesModule.ACBusTypes.PQ, data.bus_type) + ref = findall(x -> x == PowerSystems.ACBusTypesModule.ACBusTypes.REF, data.bus_type[:, time_step]) + pv = findall(x -> x == PowerSystems.ACBusTypesModule.ACBusTypes.PV, data.bus_type[:, time_step]) + pq = findall(x -> x == PowerSystems.ACBusTypesModule.ACBusTypes.PQ, data.bus_type[:, time_step]) pvpq = [pv; pq] # nref = length(ref) npv = length(pv) npq = length(pq) npvpq = npv + npq - n_buses = length(data.bus_type) Vm = data.bus_magnitude[:, time_step] # prevent unfeasible starting values for Vm; for pv and ref buses we cannot do this: diff --git a/src/nlsolve_powerflow.jl b/src/nlsolve_powerflow.jl index b4537780..77274ade 100644 --- a/src/nlsolve_powerflow.jl +++ b/src/nlsolve_powerflow.jl @@ -1,4 +1,4 @@ - +const _NLSOLVE_AC_POWERFLOW_KWARGS = Set([:check_reactive_power_limits, :check_connectivity]) function _newton_powerflow( pf::ACPowerFlow{NLSolveACPowerFlow}, data::ACPowerFlowData; @@ -13,11 +13,13 @@ function _newton_powerflow( ) end + nlsolve_solver_kwargs = filter(p -> !(p.first in _NLSOLVE_AC_POWERFLOW_KWARGS), nlsolve_kwargs) + pf = PolarPowerFlow(data) J = PowerFlows.PolarPowerFlowJacobian(data, pf.x0) df = NLsolve.OnceDifferentiable(pf, J, pf.x0, pf.residual, J.Jv) - res = NLsolve.nlsolve(df, pf.x0; nlsolve_kwargs...) + res = NLsolve.nlsolve(df, pf.x0; nlsolve_solver_kwargs...) if !res.f_converged @error( "The powerflow solver NLSolve did not converge (returned convergence = $(res.f_converged))" diff --git a/src/post_processing.jl b/src/post_processing.jl index aa23c7d6..4bb88cf0 100644 --- a/src/post_processing.jl +++ b/src/post_processing.jl @@ -461,6 +461,7 @@ function write_powerflow_solution!( system_bustype = PSY.get_bustype(bus) data_bustype = data.bus_type[ix] (system_bustype == data_bustype) && continue + println("sys bustype: $system_bustype, data bustype: $data_bustype") @assert system_bustype == PSY.ACBusTypes.PV @assert data_bustype == PSY.ACBusTypes.PQ @debug "Updating bus $(PSY.get_name(bus)) reactive power and type to PQ due to check_reactive_power_limits" @@ -620,7 +621,8 @@ function write_results( ::ACPowerFlow{<:ACPowerFlowSolverType}, sys::PSY.System, data::ACPowerFlowData, - result::Union{Vector{Float64}, Matrix{Float64}}, #todo + result::Union{Vector{Float64}, Matrix{Float64}}; + time_step::Int64 = 1, ) @info("Voltages are exported in pu. Powers are exported in MW/MVAr.") buses = sort!(collect(PSY.get_components(PSY.Bus, sys)); by = x -> PSY.get_number(x)) @@ -634,27 +636,27 @@ function write_results( P_load_vect = fill(0.0, N_BUS) Q_load_vect = fill(0.0, N_BUS) - for (ix, bustype) in enumerate(data.bus_type) - P_load_vect[ix] = data.bus_activepower_withdrawals[ix] * sys_basepower - Q_load_vect[ix] = data.bus_reactivepower_withdrawals[ix] * sys_basepower + for (ix, bustype) in enumerate(data.bus_type[:, time_step]) + P_load_vect[ix] = data.bus_activepower_withdrawals[ix, time_step] * sys_basepower + Q_load_vect[ix] = data.bus_reactivepower_withdrawals[ix, time_step] * sys_basepower P_admittance, Q_admittance = _get_fixed_admittance_power(sys, buses[ix], result, ix) P_load_vect[ix] += P_admittance * sys_basepower Q_load_vect[ix] += Q_admittance * sys_basepower if bustype == PSY.ACBusTypes.REF - Vm_vect[ix] = data.bus_magnitude[ix] - θ_vect[ix] = data.bus_angles[ix] + Vm_vect[ix] = data.bus_magnitude[ix, time_step] + θ_vect[ix] = data.bus_angles[ix, time_step] P_gen_vect[ix] = result[2 * ix - 1] * sys_basepower Q_gen_vect[ix] = result[2 * ix] * sys_basepower elseif bustype == PSY.ACBusTypes.PV - Vm_vect[ix] = data.bus_magnitude[ix] + Vm_vect[ix] = data.bus_magnitude[ix, time_step] θ_vect[ix] = result[2 * ix] - P_gen_vect[ix] = data.bus_activepower_injection[ix] * sys_basepower + P_gen_vect[ix] = data.bus_activepower_injection[ix, time_step] * sys_basepower Q_gen_vect[ix] = result[2 * ix - 1] * sys_basepower elseif bustype == PSY.ACBusTypes.PQ Vm_vect[ix] = result[2 * ix - 1] θ_vect[ix] = result[2 * ix] - P_gen_vect[ix] = data.bus_activepower_injection[ix] * sys_basepower - Q_gen_vect[ix] = data.bus_reactivepower_injection[ix] * sys_basepower + P_gen_vect[ix] = data.bus_activepower_injection[ix, time_step] * sys_basepower + Q_gen_vect[ix] = data.bus_reactivepower_injection[ix, time_step] * sys_basepower end end diff --git a/test/test_newton_ac_powerflow.jl b/test/test_newton_ac_powerflow.jl index 9c4b6717..d029a49b 100644 --- a/test/test_newton_ac_powerflow.jl +++ b/test/test_newton_ac_powerflow.jl @@ -45,25 +45,25 @@ x1 = PowerFlows._calc_x(data, V1, S1) @test LinearAlgebra.norm(result_14 - x1, Inf) <= 1e-6 - # Test that solve_ac_powerflow! succeeds + # Test that solve_powerflow! succeeds solved1 = deepcopy(sys) @test solve_powerflow!(pf, solved1) # Test that passing check_reactive_power_limits=false is the default and violates limits solved2 = deepcopy(sys) - @test solve_ac_powerflow!(solved2; check_reactive_power_limits = false) + @test solve_powerflow!(pf, solved2; check_reactive_power_limits = false) @test IS.compare_values(solved1, solved2) @test get_reactive_power(get_component(ThermalStandard, solved2, "Bus8")) > get_reactive_power_limits(get_component(ThermalStandard, solved2, "Bus8")).max # Test that passing check_reactive_power_limits=true fixes that solved3 = deepcopy(sys) - @test solve_ac_powerflow!(solved3; check_reactive_power_limits = true) + @test solve_powerflow!(pf, solved3; check_reactive_power_limits = true) @test get_reactive_power(get_component(ThermalStandard, solved3, "Bus8")) <= get_reactive_power_limits(get_component(ThermalStandard, solved3, "Bus8")).max # Test Newton method - @test solve_ac_powerflow!(deepcopy(sys); method = :newton) + @test solve_powerflow!(pf, deepcopy(sys); method = :newton) # Test enforcing the reactive power limits in closer detail set_reactive_power!(get_component(PowerLoad, sys, "Bus4"), 0.0) @@ -505,10 +505,12 @@ end pf_default, sys; check_connectivity = true) + + time_step = 1 - ref = findall(x -> x == PowerSystems.ACBusTypesModule.ACBusTypes.REF, data.bus_type) - pv = findall(x -> x == PowerSystems.ACBusTypesModule.ACBusTypes.PV, data.bus_type) - pq = findall(x -> x == PowerSystems.ACBusTypesModule.ACBusTypes.PQ, data.bus_type) + ref = findall(x -> x == PowerSystems.ACBusTypesModule.ACBusTypes.REF, data.bus_type[:, time_step]) + pv = findall(x -> x == PowerSystems.ACBusTypesModule.ACBusTypes.PV, data.bus_type[:, time_step]) + pq = findall(x -> x == PowerSystems.ACBusTypesModule.ACBusTypes.PQ, data.bus_type[:, time_step]) pvpq = [pv; pq] npvpq = length(pvpq) diff --git a/test/test_utils/legacy_pf.jl b/test/test_utils/legacy_pf.jl index a6ba33df..eafb8077 100644 --- a/test/test_utils/legacy_pf.jl +++ b/test/test_utils/legacy_pf.jl @@ -47,8 +47,8 @@ function _newton_powerflow( # Find indices for each bus type #ref = findall(x -> x == PowerSystems.ACBusTypesModule.ACBusTypes.REF, data.bus_type) - pv = findall(x -> x == PowerSystems.ACBusTypesModule.ACBusTypes.PV, data.bus_type) - pq = findall(x -> x == PowerSystems.ACBusTypesModule.ACBusTypes.PQ, data.bus_type) + pv = findall(x -> x == PowerSystems.ACBusTypesModule.ACBusTypes.PV, data.bus_type[:, time_step]) + pq = findall(x -> x == PowerSystems.ACBusTypesModule.ACBusTypes.PQ, data.bus_type[:, time_step]) pvpq = [pv; pq] #nref = length(ref)