Skip to content

Commit

Permalink
small change to _update_system!
Browse files Browse the repository at this point in the history
Roman Bolgaryn committed Jan 16, 2025
1 parent 09189c3 commit 691a73a
Showing 6 changed files with 126 additions and 45 deletions.
6 changes: 4 additions & 2 deletions src/ac_power_flow.jl
Original file line number Diff line number Diff line change
@@ -61,9 +61,11 @@ function PolarPowerFlow(data::ACPowerFlowData; time_step::Int64 = 1)
Q_net = zeros(n_buses)
for ix in 1:n_buses
P_net[ix] =
data.bus_activepower_injection[ix, time_step] - data.bus_activepower_withdrawals[ix, time_step]
data.bus_activepower_injection[ix, time_step] -
data.bus_activepower_withdrawals[ix, time_step]
Q_net[ix] =
data.bus_reactivepower_injection[ix, time_step] - data.bus_reactivepower_withdrawals[ix, time_step]
data.bus_reactivepower_injection[ix, time_step] -
data.bus_reactivepower_withdrawals[ix, time_step]
end
x0 = _calculate_x0(1,
data.bus_type[:, time_step],
74 changes: 49 additions & 25 deletions src/newton_ac_powerflow.jl
Original file line number Diff line number Diff line change
@@ -33,7 +33,7 @@ solve_ac_powerflow!(sys, method=:newton)
function solve_powerflow!(
pf::ACPowerFlow{<:ACPowerFlowSolverType},
system::PSY.System;
time_step::Int64=1,
time_step::Int64 = 1,
kwargs...,
)
#Save per-unit flag
@@ -47,11 +47,16 @@ function solve_powerflow!(
check_connectivity = get(kwargs, :check_connectivity, true),
)

converged, V, Sbus_result = _ac_powereflow(data, pf; time_step=time_step, kwargs...)
converged, V, Sbus_result = _ac_powereflow(data, pf; time_step = time_step, kwargs...)
x = _calc_x(data, V, Sbus_result)

if converged
write_powerflow_solution!(system, x, data, get(kwargs, :maxIter, DEFAULT_NR_MAX_ITER))
write_powerflow_solution!(
system,
x,
data,
get(kwargs, :maxIter, DEFAULT_NR_MAX_ITER),
)
@info("PowerFlow solve converged, the results have been stored in the system")
else
@error("The powerflow solver returned convergence = $(converged)")
@@ -163,15 +168,13 @@ function solve_powerflow(
return results
end



# Multiperiod power flow - work in progress
function solve_powerflow!(
data::ACPowerFlowData;
kwargs...,
)
pf = ACPowerFlow() # todo: somehow store in data which PF to use (see issue #50)

sorted_time_steps = sort(collect(keys(data.timestep_map)))
# preallocate results
ts_converged = zeros(Bool, 1, length(sorted_time_steps))
@@ -191,38 +194,48 @@ function solve_powerflow!(
ts_V[:, t] .= V
ts_S[:, t] .= Sbus_result

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])
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
# write results for REF
data.bus_activepower_injection[ref, t] .= real.(Sbus_result[ref]) .+ data.bus_activepower_withdrawals[ref, t]
data.bus_reactivepower_injection[ref, t] .= imag.(Sbus_result[ref]) .+ data.bus_reactivepower_withdrawals[ref, t]
data.bus_activepower_injection[ref, t] .=
real.(Sbus_result[ref]) .+ data.bus_activepower_withdrawals[ref, t]
data.bus_reactivepower_injection[ref, t] .=
imag.(Sbus_result[ref]) .+ data.bus_reactivepower_withdrawals[ref, t]
# write Q results for PV
data.bus_reactivepower_injection[pv, t] .= imag.(Sbus_result[pv]) .+ data.bus_reactivepower_withdrawals[pv, t]
data.bus_reactivepower_injection[pv, t] .=
imag.(Sbus_result[pv]) .+ data.bus_reactivepower_withdrawals[pv, t]
# results for PQ buses do not need to be updated -> already consistent with inputs

# write bus bus_types
# todo

# write voltage results
data.bus_magnitude[pq, t] .= abs.(V[pq])
data.bus_angles[pq, t] .= angle.(V[pq])
data.bus_angles[pv, t] .= angle.(V[pv])



else
# todo
1+2
1 + 2
end
end

# write branch flows
Sft = ts_V[fb,:] .* conj.(Yft * ts_V)
Stf = ts_V[tb,:] .* conj.(Ytf * 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)
@@ -235,7 +248,6 @@ function solve_powerflow!(
return
end


function _ac_powereflow(
data::ACPowerFlowData,
pf::ACPowerFlow{<:ACPowerFlowSolverType};
@@ -275,12 +287,14 @@ function _check_q_limit_bounds!(
@info "Bus $(bus_names[ix]) changed to PSY.ACBusTypes.PQ"
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]
data.bus_reactivepower_injection[ix, time_step] =
data.bus_reactivepower_bounds[ix, time_step][1]
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[ix, time_step] = PSY.ACBusTypes.PQ
data.bus_reactivepower_injection[ix, time_step] = data.bus_reactivepower_bounds[ix, time_step][2]
data.bus_reactivepower_injection[ix, time_step] =
data.bus_reactivepower_bounds[ix, time_step][2]
else
@debug "Within Limits"
end
@@ -296,7 +310,8 @@ function _solve_powerflow!(
kwargs...,
)
for _ in 1:MAX_REACTIVE_POWER_ITERATIONS
converged, V, Sbus_result = _newton_powerflow(pf, data; time_step=time_step, kwargs...)
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, time_step)
return converged, V, Sbus_result
@@ -525,9 +540,18 @@ 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[:, 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])
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)
6 changes: 4 additions & 2 deletions src/nlsolve_powerflow.jl
Original file line number Diff line number Diff line change
@@ -1,4 +1,5 @@
const _NLSOLVE_AC_POWERFLOW_KWARGS = Set([:check_reactive_power_limits, :check_connectivity])
const _NLSOLVE_AC_POWERFLOW_KWARGS =
Set([:check_reactive_power_limits, :check_connectivity])
function _newton_powerflow(
pf::ACPowerFlow{NLSolveACPowerFlow},
data::ACPowerFlowData;
@@ -13,7 +14,8 @@ function _newton_powerflow(
)
end

nlsolve_solver_kwargs = filter(p -> !(p.first in _NLSOLVE_AC_POWERFLOW_KWARGS), nlsolve_kwargs)
nlsolve_solver_kwargs =
filter(p -> !(p.first in _NLSOLVE_AC_POWERFLOW_KWARGS), nlsolve_kwargs)

pf = PolarPowerFlow(data)
J = PowerFlows.PolarPowerFlowJacobian(data, pf.x0)
58 changes: 48 additions & 10 deletions src/post_processing.jl
Original file line number Diff line number Diff line change
@@ -709,25 +709,63 @@ if a new `PowerFlowData` is constructed from the resulting system it is the same
See also `write_powerflow_solution!`. NOTE that this assumes that `data` was initialized
from `sys` and then solved with no further modifications.
"""
function update_system!(sys::PSY.System, data::PowerFlowData; time_step=1)
function update_system!(sys::PSY.System, data::PowerFlowData; time_step = 1)
for bus in PSY.get_components(PSY.Bus, sys)
if bus.bustype == PSY.ACBusTypes.REF
bus_number = data.bus_lookup[PSY.get_number(bus)]
bus_type = data.bus_type[bus_number, time_step] # use this instead of bus.bustype to account for PV -> PQ
if bus_type == PSY.ACBusTypes.REF
# For REF bus, voltage and angle are fixed; update active and reactive
P_gen = data.bus_activepower_injection[data.bus_lookup[PSY.get_number(bus)], time_step]
Q_gen = data.bus_reactivepower_injection[data.bus_lookup[PSY.get_number(bus)], time_step]
P_gen = data.bus_activepower_injection[bus_number, time_step]
Q_gen = data.bus_reactivepower_injection[bus_number, time_step]
_power_redistribution_ref(sys, P_gen, Q_gen, bus,
DEFAULT_MAX_REDISTRIBUTION_ITERATIONS)
elseif bus.bustype == PSY.ACBusTypes.PV
elseif bus_type == PSY.ACBusTypes.PV
# For PV bus, active and voltage are fixed; update reactive and angle
Q_gen = data.bus_reactivepower_injection[data.bus_lookup[PSY.get_number(bus)], time_step]
Q_gen = data.bus_reactivepower_injection[bus_number, time_step]
_reactive_power_redistribution_pv(sys, Q_gen, bus,
DEFAULT_MAX_REDISTRIBUTION_ITERATIONS)
PSY.set_angle!(bus, data.bus_angles[data.bus_lookup[PSY.get_number(bus)], time_step])
elseif bus.bustype == PSY.ACBusTypes.PQ
PSY.set_angle!(bus, data.bus_angles[bus_number, time_step])
elseif bus_type == PSY.ACBusTypes.PQ
# For PQ bus, active and reactive are fixed; update voltage and angle
Vm = data.bus_magnitude[data.bus_lookup[PSY.get_number(bus)], time_step]
Vm = data.bus_magnitude[bus_number, time_step]
PSY.set_magnitude!(bus, Vm)
PSY.set_angle!(bus, data.bus_angles[data.bus_lookup[PSY.get_number(bus)], time_step])
PSY.set_angle!(bus, data.bus_angles[bus_number, time_step])
# if it used to be a PV bus, also set the Q value:
if bus.bustype == PSY.ACBusTypes.PV
Q_gen = data.bus_reactivepower_injection[bus_number, time_step]
_reactive_power_redistribution_pv(sys, Q_gen, bus,
DEFAULT_MAX_REDISTRIBUTION_ITERATIONS)
# now both the Q and the Vm, Va are correct for this kind of buses
end
end
end
end

# This cannot work because we do not know the redistribution keys from the power values of devices
# to the power values of buses. Leaving this here anyways to document this issue.
function set_system_time_step(sys::PSY.System, data::PowerFlowData; time_step = 1)
@error("This function cannot work")
for bus in PSY.get_components(PSY.Bus, sys)
bus_number = data.bus_lookup[PSY.get_number(bus)]
bus_type = data.bus_type[bus_number, time_step]
sys_bus_type = sys.bustype # here we should be using the system bus type to know which inputs to modify

# for REF bus, we don't need to set anything
if sys_bus_type == PSY.ACBusTypes.PV
# we need to modify the P value
# TODO if we also modify time-series data for voltage set-points, also set that value.
# However, we cannot do this easily because it can happen that the PV -> PQ transformation
# caused a different Vm value as the result, while the set-point Vm is still the same.
# We need a solution for this.
P_gen = data.bus_activepower_injection[bus_number, time_step]
Q_gen = 0.0
# TODO write time-series value of the bus back to devices (not possible)
elseif sys_bus_type == PSY.ACBusTypes.PQ
P_gen = data.bus_activepower_injection[bus_number, time_step]
Q_gen = data.bus_reactivepower_injection[bus_number, time_step]
P_load = data.bus_activepower_withdrawals[bus_number, time_step]
Q_load = data.bus_reactivepower_withdrawals[bus_number, time_step]
# TODO write time-series value of the bus back to devices (not possible)
end
end
end
17 changes: 13 additions & 4 deletions test/test_newton_ac_powerflow.jl
Original file line number Diff line number Diff line change
@@ -505,12 +505,21 @@ end
pf_default,
sys;
check_connectivity = true)

time_step = 1

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])
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)
10 changes: 8 additions & 2 deletions test/test_utils/legacy_pf.jl
Original file line number Diff line number Diff line change
@@ -47,8 +47,14 @@ 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[:, time_step])
pq = findall(x -> x == PowerSystems.ACBusTypesModule.ACBusTypes.PQ, 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)

0 comments on commit 691a73a

Please sign in to comment.