Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Feature/nr_pf_solver #54

Open
wants to merge 17 commits into
base: hrgks/psse_exporter_psy4
Choose a base branch
from

Conversation

rbolgaryn
Copy link

Integrate the implementation of the Newton-Raphson Power Flow calculation leveraging KLU factorization

Roman Bolgaryn added 4 commits December 5, 2024 08:55
…CPowerFlow, introduce KLUACPowerFlow for the use in the new NR power flow function
…nction supports KLU PF; implement tests for KLUACPowerFlow
@rbolgaryn rbolgaryn added the enhancement New feature or request label Dec 6, 2024
@rbolgaryn rbolgaryn requested a review from GabrielKS December 6, 2024 21:54
Copy link
Contributor

@github-actions github-actions bot left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Remaining comments which cannot be posted as a review comment to avoid GitHub Rate Limit

JuliaFormatter

[JuliaFormatter] reported by reviewdog 🐶

if isapprox(Psources, P_gen; atol=0.001) &&
isapprox(Qsources, Q_gen; atol=0.001)


[JuliaFormatter] reported by reviewdog 🐶

sort(collect(devices_); by=x -> _get_limits_for_power_distribution(x).max)


[JuliaFormatter] reported by reviewdog 🐶

if !isapprox(p_residual, 0.0; atol=ISAPPROX_ZERO_TOLERANCE)


[JuliaFormatter] reported by reviewdog 🐶

while !isapprox(p_residual, 0.0; atol=ISAPPROX_ZERO_TOLERANCE)


[JuliaFormatter] reported by reviewdog 🐶

if isapprox(p_residual, 0; atol=ISAPPROX_ZERO_TOLERANCE)


[JuliaFormatter] reported by reviewdog 🐶

if !isapprox(p_residual, 0.0; atol=ISAPPROX_ZERO_TOLERANCE)


[JuliaFormatter] reported by reviewdog 🐶

if isapprox(Qsources, Q_gen; atol=0.001)


[JuliaFormatter] reported by reviewdog 🐶

devices = sort(collect(devices_); by=x -> PSY.get_max_reactive_power(x))


[JuliaFormatter] reported by reviewdog 🐶

if isapprox(total_active_power, 0.0; atol=ISAPPROX_ZERO_TOLERANCE)


[JuliaFormatter] reported by reviewdog 🐶

if isapprox(q_limits.max, 0.0; atol=BOUNDS_TOLERANCE) &&
isapprox(q_limits.min, 0.0; atol=BOUNDS_TOLERANCE)


[JuliaFormatter] reported by reviewdog 🐶

if isapprox(q_residual, 0.0; atol=ISAPPROX_ZERO_TOLERANCE)


[JuliaFormatter] reported by reviewdog 🐶

if !isapprox(q_residual, 0.0; atol=ISAPPROX_ZERO_TOLERANCE)


[JuliaFormatter] reported by reviewdog 🐶

while !isapprox(q_residual, 0.0; atol=ISAPPROX_ZERO_TOLERANCE)


[JuliaFormatter] reported by reviewdog 🐶

if isapprox(q_residual, 0; atol=ISAPPROX_ZERO_TOLERANCE)


[JuliaFormatter] reported by reviewdog 🐶

if !isapprox(q_residual, 0.0; atol=ISAPPROX_ZERO_TOLERANCE)


[JuliaFormatter] reported by reviewdog 🐶

atol=ISAPPROX_ZERO_TOLERANCE,


[JuliaFormatter] reported by reviewdog 🐶

sort!(collect(PSY.get_components(PSY.Bus, sys)); by=x -> PSY.get_number(x)),


[JuliaFormatter] reported by reviewdog 🐶

P_gen = result[2*ix-1]
Q_gen = result[2*ix]


[JuliaFormatter] reported by reviewdog 🐶

Q_gen = result[2*ix-1]
bus.angle = result[2*ix]


[JuliaFormatter] reported by reviewdog 🐶

Vm = result[2*ix-1]
θ = result[2*ix]


[JuliaFormatter] reported by reviewdog 🐶

function _get_branches_buses(data::Union{PTDFPowerFlowData,vPTDFPowerFlowData})


[JuliaFormatter] reported by reviewdog 🐶

bus_number=buses,
Vm=bus_magnitude,
θ=bus_angles,
P_gen=P_gen_vect,
P_load=P_load_vect,
P_net=P_gen_vect - P_load_vect,
Q_gen=Q_gen_vect,
Q_load=Q_load_vect,
Q_net=Q_gen_vect - Q_load_vect,


[JuliaFormatter] reported by reviewdog 🐶

line_name=branches,
bus_from=from_bus,
bus_to=to_bus,
P_from_to=P_from_to_vect,
Q_from_to=Q_from_to_vect,
P_to_from=P_to_from_vect,
Q_to_from=Q_to_from_vect,
P_losses=zeros(length(branches)),
Q_losses=zeros(length(branches)),


[JuliaFormatter] reported by reviewdog 🐶

data::Union{PTDFPowerFlowData,vPTDFPowerFlowData,ABAPowerFlowData},


[JuliaFormatter] reported by reviewdog 🐶

result_dict = Dict{Union{String,Char},Dict{String,DataFrames.DataFrame}}()


[JuliaFormatter] reported by reviewdog 🐶

buses = sort!(collect(PSY.get_components(PSY.Bus, sys)); by=x -> PSY.get_number(x))


[JuliaFormatter] reported by reviewdog 🐶

P_gen_vect[ix] = result[2*ix-1] * sys_basepower
Q_gen_vect[ix] = result[2*ix] * sys_basepower


[JuliaFormatter] reported by reviewdog 🐶

θ_vect[ix] = result[2*ix]


[JuliaFormatter] reported by reviewdog 🐶

Q_gen_vect[ix] = result[2*ix-1] * sys_basepower


[JuliaFormatter] reported by reviewdog 🐶

Vm_vect[ix] = result[2*ix-1]
θ_vect[ix] = result[2*ix]


[JuliaFormatter] reported by reviewdog 🐶

bus_number=PSY.get_number.(buses),
Vm=Vm_vect,
θ=θ_vect,
P_gen=P_gen_vect,
P_load=P_load_vect,
P_net=P_gen_vect - P_load_vect,
Q_gen=Q_gen_vect,
Q_load=Q_load_vect,
Q_net=Q_gen_vect - Q_load_vect,


[JuliaFormatter] reported by reviewdog 🐶

line_name=PSY.get_name.(branches),
bus_from=PSY.get_number.(PSY.get_from.(PSY.get_arc.(branches))),
bus_to=PSY.get_number.(PSY.get_to.(PSY.get_arc.(branches))),
P_from_to=P_from_to_vect,
Q_from_to=Q_from_to_vect,
P_to_from=P_to_from_vect,
Q_to_from=Q_to_from_vect,
P_losses=P_from_to_vect + P_to_from_vect,
Q_losses=Q_from_to_vect + Q_to_from_vect,


[JuliaFormatter] reported by reviewdog 🐶

sys = PSB.build_system(PSB.PSITestSystems, "c_sys14"; add_forecasts=false)
data = PowerFlows.PowerFlowData(ACPowerFlow(), sys; check_connectivity=true)


[JuliaFormatter] reported by reviewdog 🐶

@test solve_ac_powerflow!(ACPowerFlow(), sys; method=:newton)


[JuliaFormatter] reported by reviewdog 🐶

data = PowerFlows.PowerFlowData(ACPowerFlow(), sys; check_connectivity=true)


[JuliaFormatter] reported by reviewdog 🐶

@testset "AC Power Flow 14-Bus Line Configurations" for ACPowerFlow in (NLSolveACPowerFlow, KLUACPowerFlow)
sys = PSB.build_system(PSB.PSITestSystems, "c_sys14"; add_forecasts=false)


[JuliaFormatter] reported by reviewdog 🐶

sys = PSB.build_system(PSB.PSITestSystems, "c_sys14"; add_forecasts=false)


[JuliaFormatter] reported by reviewdog 🐶

@test isapprox(PSY.get_magnitude(test_bus), 1.002; atol=1e-3, rtol=0)


[JuliaFormatter] reported by reviewdog 🐶

sys = PSB.build_system(PSB.PSITestSystems, "c_sys14"; add_forecasts=false)


[JuliaFormatter] reported by reviewdog 🐶

@test isapprox(df["bus_results"].P_gen, p_gen_matpower_3bus, atol=1e-4)
@test isapprox(df["bus_results"].Q_gen, q_gen_matpower_3bus, atol=1e-4)


[JuliaFormatter] reported by reviewdog 🐶

@testset "AC Power Flow convergence fail testing" for ACPowerFlow in (NLSolveACPowerFlow, KLUACPowerFlow)
pf_sys5_re = PSB.build_system(PSB.PSITestSystems, "c_sys5_re"; add_forecasts=false)


[JuliaFormatter] reported by reviewdog 🐶

bus_name_formatter=x -> strip(string(x["name"])) * "-" * string(x["index"]),
runchecks=false,


[JuliaFormatter] reported by reviewdog 🐶

number=1,
name="01",
bustype=ACBusTypes.REF,
angle=0.0,
magnitude=1.1,
voltage_limits=(0.0, 2.0),
base_voltage=230,


[JuliaFormatter] reported by reviewdog 🐶

name="source_1",
available=true,
bus=b,
active_power=0.5,
reactive_power=0.1,
R_th=1e-5,
X_th=1e-5,


[JuliaFormatter] reported by reviewdog 🐶

name="source_2",
available=true,
bus=b,
active_power=-0.5,
reactive_power=-0.1,
R_th=1e-5,
X_th=1e-5,


[JuliaFormatter] reported by reviewdog 🐶

number=1,
name="01",
bustype=ACBusTypes.REF,
angle=0.0,
magnitude=1.1,
voltage_limits=(0.0, 2.0),
base_voltage=230,


[JuliaFormatter] reported by reviewdog 🐶

number=2,
name="02",
bustype=ACBusTypes.PV,
angle=0.0,
magnitude=1.1,
voltage_limits=(0.0, 2.0),
base_voltage=230,


[JuliaFormatter] reported by reviewdog 🐶

a = Arc(; from=b1, to=b2)


[JuliaFormatter] reported by reviewdog 🐶

name="l1",
available=true,
active_power_flow=0.0,
reactive_power_flow=0.0,
arc=a,
r=1e-3,
x=1e-3,
b=(from=0.0, to=0.0),
rating=0.0,
angle_limits=(min=-pi / 2, max=pi / 2),


[JuliaFormatter] reported by reviewdog 🐶

name="source_1",
available=true,
bus=b1,
active_power=0.5,
reactive_power=0.1,
R_th=1e-5,
X_th=1e-5,


[JuliaFormatter] reported by reviewdog 🐶

name="source_2",
available=true,
bus=b2,
active_power=0.5,
reactive_power=1.1,
R_th=1e-5,
X_th=1e-5,


[JuliaFormatter] reported by reviewdog 🐶

name="source_3",
available=true,
bus=b2,
active_power=-0.5,
reactive_power=-1.1,
R_th=1e-5,
X_th=1e-5,


[JuliaFormatter] reported by reviewdog 🐶

number=1,
name="01",
bustype=ACBusTypes.REF,
angle=0.0,
magnitude=1.1,
voltage_limits=(0.0, 2.0),
base_voltage=230,


[JuliaFormatter] reported by reviewdog 🐶

name="source_1",
available=true,
bus=b,
active_power=0.5,
reactive_power=0.1,
R_th=1e-5,
X_th=1e-5,


[JuliaFormatter] reported by reviewdog 🐶

name="init",
available=true,
status=false,
bus=b,
active_power=0.1,
reactive_power=0.1,
rating=0.0,
active_power_limits=(min=0.0, max=0.0),
reactive_power_limits=nothing,
ramp_limits=nothing,
operation_cost=ThermalGenerationCost(nothing),
base_power=100.0,
time_limits=nothing,
prime_mover_type=PrimeMovers.OT,
fuel=ThermalFuels.OTHER,
services=Device[],
dynamic_injector=nothing,
ext=Dict{String,Any}(),


[JuliaFormatter] reported by reviewdog 🐶


[JuliaFormatter] reported by reviewdog 🐶


[JuliaFormatter] reported by reviewdog 🐶

number=1,
name="01",
bustype=ACBusTypes.REF,
angle=0.0,
magnitude=1.1,
voltage_limits=(0.0, 2.0),
base_voltage=230,


[JuliaFormatter] reported by reviewdog 🐶

number=2,
name="02",
bustype=ACBusTypes.PV,
angle=0.0,
magnitude=1.1,
voltage_limits=(0.0, 2.0),
base_voltage=230,


[JuliaFormatter] reported by reviewdog 🐶

a = Arc(; from=b1, to=b2)


[JuliaFormatter] reported by reviewdog 🐶

name="l1",
available=true,
active_power_flow=0.0,
reactive_power_flow=0.0,
arc=a,
r=1e-3,
x=1e-3,
b=(from=0.0, to=0.0),
rating=0.0,
angle_limits=(min=-pi / 2, max=pi / 2),


[JuliaFormatter] reported by reviewdog 🐶

name="source_1",
available=true,
bus=b1,
active_power=0.5,
reactive_power=0.1,
R_th=1e-5,
X_th=1e-5,


[JuliaFormatter] reported by reviewdog 🐶

name="source_2",
available=true,
bus=b2,
active_power=0.5,
reactive_power=1.1,
R_th=1e-5,
X_th=1e-5,


[JuliaFormatter] reported by reviewdog 🐶

name="init",
available=true,
status=false,
bus=b2,
active_power=0.1,
reactive_power=0.1,
rating=0.0,
active_power_limits=(min=0.0, max=0.0),
reactive_power_limits=nothing,
ramp_limits=nothing,
operation_cost=ThermalGenerationCost(nothing),
base_power=100.0,
time_limits=nothing,
prime_mover_type=PrimeMovers.OT,
fuel=ThermalFuels.OTHER,
services=Device[],
dynamic_injector=nothing,
ext=Dict{String,Any}(),


[JuliaFormatter] reported by reviewdog 🐶


[JuliaFormatter] reported by reviewdog 🐶


[JuliaFormatter] reported by reviewdog 🐶

sys = PSB.build_system(PSB.PSITestSystems, "c_sys14"; add_forecasts=false)


[JuliaFormatter] reported by reviewdog 🐶

sys = PSB.build_system(PSB.PSITestSystems, "c_sys14"; add_forecasts=false)


[JuliaFormatter] reported by reviewdog 🐶

@test PowerFlowData(DCPowerFlow(), sys; time_steps=time_steps) isa PF.ABAPowerFlowData
@test PowerFlowData(PTDFDCPowerFlow(), sys; time_steps=time_steps) isa


[JuliaFormatter] reported by reviewdog 🐶

@test PowerFlowData(vPTDFDCPowerFlow(), sys; time_steps=time_steps) isa


[JuliaFormatter] reported by reviewdog 🐶

isdir(test_psse_export_dir) && rm(test_psse_export_dir; recursive=true)


[JuliaFormatter] reported by reviewdog 🐶

macro log_assert(ex, comparison_name=nothing)


[JuliaFormatter] reported by reviewdog 🐶

default_tol=SYSTEM_REIMPORT_COMPARISON_TOLERANCE;


[JuliaFormatter] reported by reviewdog 🐶

all(isapprox.(col1, col2; atol=my_tol))


[JuliaFormatter] reported by reviewdog 🐶

default_tol=SYSTEM_REIMPORT_COMPARISON_TOLERANCE;


[JuliaFormatter] reported by reviewdog 🐶

isapprox(a, b; atol=SYSTEM_REIMPORT_COMPARISON_TOLERANCE) || IS.isequivalent(a, b)


[JuliaFormatter] reported by reviewdog 🐶

bus_name_mapping=Dict{String,String}(),


[JuliaFormatter] reported by reviewdog 🐶

include_types=[


[JuliaFormatter] reported by reviewdog 🐶

exclude_fields=Set([


[JuliaFormatter] reported by reviewdog 🐶

exclude_fields_for_type=Dict(


[JuliaFormatter] reported by reviewdog 🐶

generator_comparison_fns=[ # TODO rating


[JuliaFormatter] reported by reviewdog 🐶

ignore_name_order=true,
ignore_extra_of_type=Union{PSY.ThermalStandard,PSY.StaticLoad},
exclude_reactive_power=false)


[JuliaFormatter] reported by reviewdog 🐶

result &= IS.compare_values(sys1, sys2; exclude=[:data])


[JuliaFormatter] reported by reviewdog 🐶

exclude=my_excludes,


[JuliaFormatter] reported by reviewdog 🐶

GenSource = Union{Generator,Source}


[JuliaFormatter] reported by reviewdog 🐶

any(Union{typeof(gen1),typeof(gen2)} .<: include_types) && continue


[JuliaFormatter] reported by reviewdog 🐶

exclude=exclude_fields,


[JuliaFormatter] reported by reviewdog 🐶

function test_power_flow(sys1::System, sys2::System; exclude_reactive_flow=false)


[JuliaFormatter] reported by reviewdog 🐶

POWERFLOW_COMPARISON_TOLERANCE; line_name=nothing, Q_to_from=reactive_power_tol,
Q_from_to=reactive_power_tol, Q_losses=reactive_power_tol)


[JuliaFormatter] reported by reviewdog 🐶

do_power_flow_test=true,
exclude_reactive_flow=false,


[JuliaFormatter] reported by reviewdog 🐶

test_power_flow(sys, sys2; exclude_reactive_flow=exclude_reactive_flow)


[JuliaFormatter] reported by reviewdog 🐶

exclude_metadata_keys=["case_name"],


[JuliaFormatter] reported by reviewdog 🐶

exclude_reactive_flow=true) # TODO why is reactive flow not matching?


[JuliaFormatter] reported by reviewdog 🐶

test_power_flow(sys2, reread_sys2; exclude_reactive_flow=true) # TODO why is reactive flow not matching?


[JuliaFormatter] reported by reviewdog 🐶

exclude_reactive_flow=true) # TODO why is reactive flow not matching?


[JuliaFormatter] reported by reviewdog 🐶

test_power_flow(sys2, reread_sys2; exclude_reactive_flow=true) # TODO why is reactive flow not matching?


[JuliaFormatter] reported by reviewdog 🐶

exclude_reactive_power=true) # TODO why is reactive power not matching?


[JuliaFormatter] reported by reviewdog 🐶

test_power_flow(sys2, reread_sys3; exclude_reactive_flow=true) # TODO why is reactive flow not matching?


[JuliaFormatter] reported by reviewdog 🐶

exporter = PSSEExporter(sys, :v33, export_location; write_comments=true)


[JuliaFormatter] reported by reviewdog 🐶

exclude_reactive_flow=true) # TODO why is reactive flow not matching?

Copy link
Contributor

@GabrielKS GabrielKS left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

A few initial comments. I'll dig into the math more later. Thanks!

Project.toml Outdated Show resolved Hide resolved
src/PowerFlows.jl Show resolved Hide resolved
src/nlsolve_ac_powerflow.jl Show resolved Hide resolved
src/nlsolve_ac_powerflow.jl Outdated Show resolved Hide resolved
write_powerflow_solution!(system, res.zero, max_iterations)
converged, x = _solve_powerflow!(pf, data, check_reactive_power_limits; kwargs...)
if converged
write_powerflow_solution!(system, x, max_iterations)
@info("PowerFlow solve converged, the results have been stored in the system")
#Restore original per unit base
PSY.set_units_base_system!(system, settings_unit_cache)
Copy link
Contributor

@GabrielKS GabrielKS Dec 10, 2024

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Not necessarily your responsibility but just flagging that this is a good use case for

function with_units(f::Function, sys::System, units::Union{PSY.UnitSystem, String})
(and a good argument for moving it someplace more central)

src/nlsolve_ac_powerflow.jl Outdated Show resolved Hide resolved
src/nlsolve_ac_powerflow.jl Outdated Show resolved Hide resolved
# 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)
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Probably no significant difference, but note that these three lines could be combined to only iterate over the bus types once

src/nlsolve_ac_powerflow.jl Outdated Show resolved Hide resolved
src/nlsolve_ac_powerflow.jl Outdated Show resolved Hide resolved
J = [j11 j12; j21 j22]

factor_J = KLU.klu(J)
dx = -(factor_J \ F)
Copy link
Contributor

@GabrielKS GabrielKS Dec 10, 2024

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Comment from José: do this and previous two lines in place using functor approach

src/powerflow_types.jl Outdated Show resolved Hide resolved
src/newton_ac_powerflow.jl Outdated Show resolved Hide resolved
src/newton_ac_powerflow.jl Outdated Show resolved Hide resolved
src/newton_ac_powerflow.jl Outdated Show resolved Hide resolved
src/newton_ac_powerflow.jl Outdated Show resolved Hide resolved
F_imag = view(F, npvpq + 1:npvpq + npq)

mis .= V .* conj.(Ybus * V) .- Sbus
@. F_real = real(mis_pvpq) # In-place assignment to the real part, using views
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@rbolgaryn I think you are overwriting F_real here. You might want to try F_real[:] =

@rbolgaryn rbolgaryn linked an issue Jan 9, 2025 that may be closed by this pull request
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
enhancement New feature or request
Projects
None yet
Development

Successfully merging this pull request may close these issues.

Add solve_powerflow! for AC power flow
3 participants