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

[DO NOT MERGE] Notes os JuMP issues #1218

Open
wants to merge 2 commits into
base: main
Choose a base branch
from

Conversation

joaquimg
Copy link

@joaquimg joaquimg commented Jan 17, 2025

main take aways:

1 - use add_to_expression! instead of += (saves allocations of extra expressions)
2 - use add_to_expression!(exp, number, var) instead of add_to_expression!(exp, number * var) because number * var will allocate.
3 - prefer (number * number) * var or var * (number * number) to force computing product of floats first (very very fast) than number * number * var or var * number * number which will depend of the JuMP order of operations (one will be fast, the other slow). use style to keep all coefs in the left and variables in the right. Alsoo use parenthesis when there are multiple coefs.
4 - avoid operations (+, -, /, *, sum) between variables (or varaibles and constants) outside macros, because macros will try to allocate less.
5 - array[:, i] allocates new vectors, use @view to avoid this
6 - attention with type stability when: a) assining values in if-else blocks, b) re-using names inside a function
7 - attention with loop order: julia is column major so inner loops should change the leftmost index.

@@ -1627,6 +1627,7 @@ function add_to_objective_invariant_expression!(
) where {T <: JuMP.AbstractJuMPScalar}
T_cf = typeof(container.objective_function.invariant_terms)
if T_cf <: JuMP.GenericAffExpr && T <: JuMP.GenericQuadExpr
# why not add_to...???
Copy link
Member

Choose a reason for hiding this comment

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

before add didn't work for mixed ExpressionTypes. We can re-test

@@ -1633,10 +1633,10 @@ function add_to_expression!(
device_base_power,
)
for t in time_steps
fuel_expr = variable[name, t] * prop_term_per_unit * dt
JuMP.add_to_expression!(
expression[name, t],
Copy link
Member

Choose a reason for hiding this comment

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

Copy link
Author

Choose a reason for hiding this comment

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

for more info:
JuMP.add_to_expression!(exp1, float_coef * jump_var) is worse because it creates an extra temporary expression float_coef * jump_var

JuMP.add_to_expression!(
gen_cost,
cost * multiplier * pwl_var_container[(name, i, time_period)],
cost * multiplier,
Copy link
Member

Choose a reason for hiding this comment

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

Copy link
Member

@jd-lara jd-lara Jan 17, 2025

Choose a reason for hiding this comment

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

The AGC functionality is offline for now. We need LHS parameters to revive it

Copy link
Author

Choose a reason for hiding this comment

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

POI is much better nowadays.
Moreover, there is room for improvement, so let me know whenever you want to test stuff and bring it back

@@ -60,10 +60,12 @@ function add_constraints!(
for t in time_steps
resource_expression = JuMP.GenericAffExpr{Float64, JuMP.VariableRef}()
for reserve_variable in reserve_variables
JuMP.add_to_expression!(resource_expression, sum(reserve_variable[:, t]))
JuMP.add_to_expression!(resource_expression,
Copy link
Member

Choose a reason for hiding this comment

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

Is this more effective than pre-allocating the expression in 61?

Copy link
Author

Choose a reason for hiding this comment

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

not sure I got what you mean

end
if use_slacks
resource_expression += slack_vars[t]
JuMP.add_to_expression!(resource_expression, slack_vars[t])
Copy link
Member

Choose a reason for hiding this comment

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

good catch

@@ -161,19 +161,23 @@ function add_constraints!(
param = get_parameter_column_refs(param_container, service_name)
for t in time_steps
if use_slacks
resource_expression = sum(reserve_variable[:, t]) + slack_vars[t]
resource_expression = JuMP.@expression(
Copy link
Member

Choose a reason for hiding this comment

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

this could be related to #1000

Copy link
Member

Choose a reason for hiding this comment

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

These improvements should make it to main

Copy link
Member

Choose a reason for hiding this comment

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

These ones are quite important too.

Copy link
Member

@jd-lara jd-lara left a comment

Choose a reason for hiding this comment

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

These are great catches @joaquimg I would like to put some of these into main sooner

@@ -1627,7 +1627,7 @@ function add_to_objective_invariant_expression!(
) where {T <: JuMP.AbstractJuMPScalar}
T_cf = typeof(container.objective_function.invariant_terms)
if T_cf <: JuMP.GenericAffExpr && T <: JuMP.GenericQuadExpr
# why not add_to...???
# TODO: why not add_to...???
Copy link
Member

Choose a reason for hiding this comment

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

The add_to here usually fails is not allowed by JuMP

Copy link
Author

Choose a reason for hiding this comment

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

This is weird, do you have an example? It does not happen in RTS. Maybe it was an old JuMP bug?

Copy link
Member

Choose a reason for hiding this comment

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

This only happens when we have data that mixes PWL and Quadratic cost functions.

Copy link
Author

Choose a reason for hiding this comment

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

If there is any case I can test with those features I'd be happy to.

add_to_objective_variant_expression!(container, cost_expr)
# TODO: missing _mult below?
Copy link
Member

Choose a reason for hiding this comment

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

no, one expression tracks the cost and the other the fuel consumption

Copy link
Author

Choose a reason for hiding this comment

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

yup not much to do here, as the two need to be kept different

@@ -61,7 +61,7 @@ function add_constraints!(
resource_expression = JuMP.GenericAffExpr{Float64, JuMP.VariableRef}()
for reserve_variable in reserve_variables
JuMP.add_to_expression!(resource_expression,
JuMP.@expression(container.JuMPmodel, sum(reserve_variable[:, t])))
JuMP.@expression(container.JuMPmodel, sum(@view reserve_variable[:, t])))
Copy link
Member

Choose a reason for hiding this comment

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

nice trick. @GabrielKS take note

@joaquimg
Copy link
Author

joaquimg commented Jan 28, 2025

Feel free to add them in main before I do.
I was trying to catch as many possible issues as I could.

Many of these are general suggestions and should also make it to the JuMP docs.

The gain is not big in the RTS case, about 5% less allocations, but these should make a difference in larger cases with more details.

My plan was to make cleaner PR's after this pass.

for device in devices
ci_name = PSY.get_name(device)
if ci_name ∈ PNM.get_radial_branches(radial_network_reduction)
continue
end
limits = get_min_max_limits(device, RateLimitConstraint, U) # depends on constraint type and formulation type
for t in time_steps
# TODO: ternary
Copy link
Author

Choose a reason for hiding this comment

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

ignore this: issue solved recently in MutableArithmetics

@@ -606,6 +608,7 @@ function _make_flow_expressions!(
sum(
ptdf_col[i] * nodal_balance_expressions[i, t] for
i in 1:length(ptdf_col)
# maybe only sum if nonempy?
Copy link
Author

Choose a reason for hiding this comment

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

ignore

Comment on lines -199 to +200
JuMP.add_to_expression!(lhs_off, (1 - varon[name, i]))
JuMP.add_to_expression!(lhs_off, 1)
JuMP.add_to_expression!(lhs_off, -1, varon[name, i])
Copy link
Author

Choose a reason for hiding this comment

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

the previous allocates more because (1 - varon[name, i]) becomes a new expression

@@ -137,11 +137,12 @@ function _get_pwl_cost_expression(
name = PSY.get_name(component)
pwl_var_container = get_variable(container, PieceWiseLinearBlockOffer(), T)
gen_cost = JuMP.AffExpr(0.0)
cost_data = PSY.get_y_coords(cost_data)
for (i, cost) in enumerate(cost_data)
_cost_data = PSY.get_y_coords(cost_data)
Copy link
Author

Choose a reason for hiding this comment

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

this is for type stability

Comment on lines -90 to +91
real = get_variable(container, ActivePowerVariable(), V)[name, t] * pf
constraint[name, t] = JuMP.@constraint(jump_model, reactive == real)
real = get_variable(container, ActivePowerVariable(), V)[name, t]
constraint[name, t] = JuMP.@constraint(jump_model, reactive == real * pf)
Copy link
Author

Choose a reason for hiding this comment

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

the product outside the macro lead to additional allocations

@@ -131,7 +132,7 @@ function update_variable_cost!(
end
mult_ = parameter_multiplier[component_name, time_period]
variable = get_variable(container, get_variable_type(attributes)(), T)
gen_cost = variable[component_name, time_period] * mult_ * cost_data * base_power * dt
gen_cost = variable[component_name, time_period] * (mult_ * cost_data * base_power * dt)
Copy link
Author

Choose a reason for hiding this comment

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

this is interesting:
the parenthesis forces the floats the be multiplied before and simplified into a single float. the previous one requires a bunch of operations in an expression. not an allocation issue, but slower than multiplying floats.

@@ -224,12 +228,12 @@ function add_constraints!(
cons[name, t] =
JuMP.@constraint(
jump_model,
var_r[name, t] <= param[t] * requirement * max_participation_factor
var_r[name, t] <= (requirement * max_participation_factor) * param[t]
Copy link
Author

Choose a reason for hiding this comment

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

rule of thumb: first multiply floats and multiply them with jump variables

@@ -473,7 +479,7 @@ function add_constraints!(
PSY.get_active_power_limits(d).min +
(reserve_response_time - startup_time) * minutes_per_period * ramp_limits.up
else
reserve_limit = 0
reserve_limit = 0.0
Copy link
Author

Choose a reason for hiding this comment

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

type stability as the if block would return a float and 0 is ana Int

@@ -387,13 +387,15 @@ function add_constraints!(
slack_lb = get_variable(container, FlowActivePowerSlackLowerBound(), T)
end

# TODO: loop order
Copy link
Author

Choose a reason for hiding this comment

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

this is interesting, and PSI uses mixed patterns.

Julia is column based, so if you vector is array[name, t]
you should iterate as:

for t in T
    for name in names
        # do something with array[name, t]
    end
end

which is equivalent to:

for t in T, name in names
    # do something with array[name, t]
end

In the case of sums:

sum(array[name, t] for t in T, name in names)

equivalent to:

sum(array[name, t] for t in T for name in names)

almost equivalent to:

sum(sum(array[name, t] for name in names) for t in T )

Copy link
Author

Choose a reason for hiding this comment

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

so names should typically be the inner loop and t the outer loop.
but if the operation inside the loop is complex both are almost the same.
I would keep the pattern so style would lead you to the correct choices in thigh loops.

Copy link
Member

Choose a reason for hiding this comment

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

Yeah, I originally wrote all the loops in the right order and then had to modify some depending on the context.

Comment on lines +920 to +925
flow_variables_ = @view flow_variables[name, :]
from_bus = PSY.get_name(PSY.get_from(PSY.get_arc(br)))
to_bus = PSY.get_name(PSY.get_to(PSY.get_arc(br)))
angle_variables_ = ps_angle_variables[name, :]
bus_angle_from = bus_angle_variables[from_bus, :]
bus_angle_to = bus_angle_variables[to_bus, :]
angle_variables_ = @view ps_angle_variables[name, :]
bus_angle_from = @view bus_angle_variables[from_bus, :]
bus_angle_to = @view bus_angle_variables[to_bus, :]
Copy link
Author

Choose a reason for hiding this comment

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

not sure about these. the @view saves many allocations but the the iteration is worst in the loop below as there are many cache misses.

@jd-lara
Copy link
Member

jd-lara commented Jan 28, 2025

Feel free to add them in main before I do. I was trying to catch as many possible issues as I could.

Many of these are general suggestions and should also make it to the JuMP docs.

The gain is not big in the RTS case, about 5% less allocations, but these should make a difference in larger cases with more details.

My plan was to make cleaner PR's after this pass.

I can test this branch on the super large systems and report back w.r.t. to main.

@joaquimg
Copy link
Author

That would be nice, just build time benchmarking should be enough. no need to run solvers.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

2 participants