From bc9b3719a7b6b2e1b09f5650de068d1aecc9dbfa Mon Sep 17 00:00:00 2001 From: Jose Daniel Lara Date: Mon, 30 Oct 2023 17:40:23 -0600 Subject: [PATCH] enforce reactive power limits --- src/PowerFlowData.jl | 38 +++++++++++++----------- src/nlsolve_ac_powerflow.jl | 59 ++++++++++++++++++++++++++++++++++--- 2 files changed, 75 insertions(+), 22 deletions(-) diff --git a/src/PowerFlowData.jl b/src/PowerFlowData.jl index 344fcd97..79737bb6 100644 --- a/src/PowerFlowData.jl +++ b/src/PowerFlowData.jl @@ -48,19 +48,19 @@ flows and angles, as well as these ones. depending on the method considered. - `neighbors::Vector{Set{Int}}`: Vector with the sets of adjacent buses. """ -struct PowerFlowData{M <: PNM.PowerNetworkMatrix, N, S <: Union{String, String}} +struct PowerFlowData{M <: PNM.PowerNetworkMatrix, N <: Union{PNM.PowerNetworkMatrix, Nothing}} bus_lookup::Dict{Int, Int} branch_lookup::Dict{String, Int} bus_activepower_injection::Matrix{Float64} bus_reactivepower_injection::Matrix{Float64} bus_activepower_withdrawals::Matrix{Float64} bus_reactivepower_withdrawals::Matrix{Float64} - bus_reactivepower_bounds::Vector{Float64} + bus_reactivepower_bounds::Vector{Vector{Float64}} bus_type::Vector{PSY.ACBusTypes} bus_magnitude::Matrix{Float64} bus_angles::Matrix{Float64} branch_flow_values::Matrix{Float64} - timestep_map::Dict{Int, S} + timestep_map::Dict{Int, String} valid_ix::Vector{Int} power_network_matrix::M aux_network_matrix::N @@ -178,18 +178,20 @@ function PowerFlowData( sys, ) - # initialize data - init_1 = zeros(n_buses, timesteps) - init_2 = zeros(n_branches, timesteps) - # define fields as matrices whose number of columns is eqault to the number of timesteps - bus_activepower_injection_1 = deepcopy(init_1) - bus_reactivepower_injection_1 = deepcopy(init_1) - bus_activepower_withdrawals_1 = deepcopy(init_1) - bus_reactivepower_withdrawals_1 = deepcopy(init_1) - bus_magnitude_1 = deepcopy(init_1) - bus_angles_1 = deepcopy(init_1) - branch_flow_values_1 = deepcopy(init_2) + bus_activepower_injection_1 = zeros(n_buses, timesteps) + bus_reactivepower_injection_1 = zeros(n_buses, timesteps) + bus_activepower_withdrawals_1 = zeros(n_buses, timesteps) + bus_reactivepower_withdrawals_1 = zeros(n_buses, timesteps) + bus_magnitude_1 = zeros(n_buses, timesteps) + bus_angles_1 = zeros(n_buses, timesteps) + branch_flow_values_1 = zeros(n_branches, timesteps) + + bus_reactivepower_bounds = Vector{Vector{Float64}}(undef, n_buses) + for i in 1:n_buses + bus_reactivepower_bounds[i] = [0.0, 0.0] + end + _get_reactive_power_bound!(bus_reactivepower_bounds, bus_lookup, sys) # initial values related to first timestep allocated in the first column bus_activepower_injection_1[:, 1] .= bus_activepower_injection @@ -207,7 +209,7 @@ function PowerFlowData( bus_reactivepower_injection_1, bus_activepower_withdrawals_1, bus_reactivepower_withdrawals_1, - Vector{Float64}(), + bus_reactivepower_bounds, bus_type, bus_magnitude_1, bus_angles_1, @@ -336,7 +338,7 @@ function PowerFlowData( bus_reactivepower_injection_1, bus_activepower_withdrawals_1, bus_reactivepower_withdrawals_1, - Vector{Float64}(), + Vector{Vector{Float64}}(), bus_type, bus_magnitude_1, bus_angles_1, @@ -465,7 +467,7 @@ function PowerFlowData( bus_reactivepower_injection_1, bus_activepower_withdrawals_1, bus_reactivepower_withdrawals_1, - Vector{Float64}(), + Vector{Vector{Float64}}(), bus_type, bus_magnitude_1, bus_angles_1, @@ -594,7 +596,7 @@ function PowerFlowData( bus_reactivepower_injection_1, bus_activepower_withdrawals_1, bus_reactivepower_withdrawals_1, - Vector{Float64}(), + Vector{Vector{Float64}}(), bus_type, bus_magnitude_1, bus_angles_1, diff --git a/src/nlsolve_ac_powerflow.jl b/src/nlsolve_ac_powerflow.jl index 5c028628..cfc3bb7b 100644 --- a/src/nlsolve_ac_powerflow.jl +++ b/src/nlsolve_ac_powerflow.jl @@ -66,7 +66,7 @@ res = solve_powerflow(sys, method=:newton) ``` """ function solve_powerflow( - ::ACPowerFlow, + pf::ACPowerFlow, system::PSY.System; kwargs..., ) @@ -79,7 +79,9 @@ function solve_powerflow( system; check_connectivity = get(kwargs, :check_connectivity, true), ) - res = _solve_powerflow(data; kwargs...) + + res = _solve_powerflow!(data, pf.check_reactive_power_limits; kwargs...) + if res.f_converged @info("PowerFlow solve converged, the results are exported in DataFrames") df_results = write_results(ACPowerFlow(), system, res.zero) @@ -92,12 +94,61 @@ function solve_powerflow( return res.f_converged end -function _solve_powerflow(data::ACPowerFlowData; kwargs...) - nlsolve_kwargs = (k for k in kwargs if first(k) ∉ AC_PF_KW) +function _check_q_limit_bounds!(data::ACPowerFlowData, zero::Vector{Float64}) + within_limits = true + for (ix, b) in enumerate(data.bus_type) + if b == PSY.ACBusTypes.PV + Q_gen = zero[2 * ix - 1] + else + continue + end + + if Q_gen <= data.bus_reactivepower_bounds[ix][1] + @show "here under" + within_limits = false + data.bus_type[ix] = PSY.ACBusTypes.PQ + data.bus_reactivepower_injection[ix] = data.bus_reactivepower_bounds[ix][1] + elseif Q_gen >= data.bus_reactivepower_bounds[ix][2] + @show "here over" + within_limits = false + data.bus_type[ix] = PSY.ACBusTypes.PQ + data.bus_reactivepower_injection[ix] = data.bus_reactivepower_bounds[ix][2] + else + @debug "Within Limits" + end + end + return within_limits +end + +function _solve_powerflow!( + data::ACPowerFlowData, + check_reactive_power_limits; + nlsolve_kwargs..., +) + if check_reactive_power_limits + for _ in 1:MAX_REACTIVE_POWER_ITERATIONS + res = _nlsolve_powerflow(data; nlsolve_kwargs...) + if res.f_converged + if _check_q_limit_bounds!(data, res.zero) + return res + end + else + return res + end + end + else + return _nlsolve_powerflow(data; nlsolve_kwargs...) + end +end + +function _nlsolve_powerflow(data::ACPowerFlowData; 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...) + if !res.f_converged + @error("The powerflow solver returned convergence = $(res.f_converged)") + end return res end