diff --git a/src/core/optimization_container.jl b/src/core/optimization_container.jl index a0a549d9d..6ae04ecd8 100644 --- a/src/core/optimization_container.jl +++ b/src/core/optimization_container.jl @@ -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 + # TODO: why not add_to...??? container.objective_function.invariant_terms += cost_expr else JuMP.add_to_expression!(container.objective_function.invariant_terms, cost_expr) diff --git a/src/devices_models/devices/AC_branches.jl b/src/devices_models/devices/AC_branches.jl index 48767cd79..d830e6535 100644 --- a/src/devices_models/devices/AC_branches.jl +++ b/src/devices_models/devices/AC_branches.jl @@ -387,6 +387,7 @@ function add_constraints!( slack_lb = get_variable(container, FlowActivePowerSlackLowerBound(), T) end + # TODO: loop order for device in devices ci_name = PSY.get_name(device) if ci_name ∈ PNM.get_radial_branches(radial_network_reduction) @@ -394,6 +395,7 @@ function add_constraints!( end limits = get_min_max_limits(device, RateLimitConstraint, U) # depends on constraint type and formulation type for t in time_steps + # TODO: ternary con_ub[ci_name, t] = JuMP.@constraint(get_jump_model(container), array[ci_name, t] - (use_slacks ? slack_ub[ci_name, t] : 0.0) <= @@ -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? ) ) end @@ -648,7 +651,7 @@ function _make_flow_expressions!( #= Leaving serial code commented out for debugging purposes in the future for name in branches - ptdf_col = ptdf[name, :] + ptdf_col = @view ptdf[name, :] branch_flow_expr[name, :] .= _make_flow_expressions!( jump_model, name, @@ -733,7 +736,7 @@ function add_constraints!( jump_model = get_jump_model(container) for br in devices name = PSY.get_name(br) - ptdf_col = ptdf[name, :] + ptdf_col = @view ptdf[name, :] inv_x = 1 / PSY.get_x(br) for t in time_steps branch_flow[name, t] = JuMP.@constraint( @@ -914,12 +917,12 @@ function add_constraints!( for br in devices name = PSY.get_name(br) inv_x = 1.0 / PSY.get_x(br) - flow_variables_ = flow_variables[name, :] + 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, :] @assert inv_x > 0.0 for t in time_steps branch_flow[name, t] = JuMP.@constraint( diff --git a/src/devices_models/devices/TwoTerminalDC_branches.jl b/src/devices_models/devices/TwoTerminalDC_branches.jl index 29685dbb8..1e8694729 100644 --- a/src/devices_models/devices/TwoTerminalDC_branches.jl +++ b/src/devices_models/devices/TwoTerminalDC_branches.jl @@ -467,14 +467,15 @@ function _add_hvdc_flow_constraints!( check_hvdc_line_limits_consistency(d) max_rate = get_variable_upper_bound(var, d, HVDCTwoTerminalDispatch()) min_rate = get_variable_lower_bound(var, d, HVDCTwoTerminalDispatch()) + name = PSY.get_name(d) for t in time_steps - constraint_ub[PSY.get_name(d), t] = JuMP.@constraint( + constraint_ub[name, t] = JuMP.@constraint( get_jump_model(container), - variable[PSY.get_name(d), t] <= max_rate + variable[name, t] <= max_rate ) - constraint_lb[PSY.get_name(d), t] = JuMP.@constraint( + constraint_lb[name, t] = JuMP.@constraint( get_jump_model(container), - min_rate <= variable[PSY.get_name(d), t] + min_rate <= variable[name, t] ) end end @@ -684,32 +685,32 @@ function add_constraints!( R_min_from, R_max_from = PSY.get_active_power_limits_from(d) R_min_to, R_max_to = PSY.get_active_power_limits_to(d) for t in get_time_steps(container) - constraint_tf_ub[PSY.get_name(d), t] = JuMP.@constraint( + constraint_tf_ub[name, t] = JuMP.@constraint( get_jump_model(container), tf_var[name, t] <= R_max_to * direction_var[name, t] ) - constraint_tf_lb[PSY.get_name(d), t] = JuMP.@constraint( + constraint_tf_lb[name, t] = JuMP.@constraint( get_jump_model(container), tf_var[name, t] >= R_min_to * (1 - direction_var[name, t]) ) - constraint_ft_ub[PSY.get_name(d), t] = JuMP.@constraint( + constraint_ft_ub[name, t] = JuMP.@constraint( get_jump_model(container), ft_var[name, t] <= R_max_from * (1 - direction_var[name, t]) ) - constraint_ft_lb[PSY.get_name(d), t] = JuMP.@constraint( + constraint_ft_lb[name, t] = JuMP.@constraint( get_jump_model(container), ft_var[name, t] >= R_min_from * direction_var[name, t] ) - constraint_loss[PSY.get_name(d), t] = JuMP.@constraint( + constraint_loss[name, t] = JuMP.@constraint( get_jump_model(container), tf_var[name, t] + ft_var[name, t] == losses[name, t] ) - constraint_loss_aux1[PSY.get_name(d), t] = JuMP.@constraint( + constraint_loss_aux1[name, t] = JuMP.@constraint( get_jump_model(container), losses[name, t] >= l1 * ft_var[name, t] + l0 - M_VALUE * direction_var[name, t] ) - constraint_loss_aux2[PSY.get_name(d), t] = JuMP.@constraint( + constraint_loss_aux2[name, t] = JuMP.@constraint( get_jump_model(container), losses[name, t] >= l1 * tf_var[name, t] + l0 - M_VALUE * (1 - direction_var[name, t]) diff --git a/src/devices_models/devices/common/add_to_expression.jl b/src/devices_models/devices/common/add_to_expression.jl index 386d27ca2..3d16b4a14 100644 --- a/src/devices_models/devices/common/add_to_expression.jl +++ b/src/devices_models/devices/common/add_to_expression.jl @@ -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], - fuel_expr, + prop_term_per_unit * dt, + variable[name, t], ) end elseif value_curve isa PSY.QuadraticCurve @@ -1709,22 +1709,33 @@ function add_to_expression!( for t in time_steps sos_status = _get_sos_value(container, W, d) if sos_status == SOSStatusVariable.NO_VARIABLE - bin = 1.0 + JuMP.add_to_expression!( + expression[name, t], + P_min * prop_term_per_unit * dt, + ) elseif sos_status == SOSStatusVariable.PARAMETER param = get_default_on_parameter(d) bin = get_parameter(container, param, V).parameter_array[name, t] + JuMP.add_to_expression!( + expression[name, t], + P_min * prop_term_per_unit * dt, + bin + ) elseif sos_status == SOSStatusVariable.VARIABLE var = get_default_on_variable(d) bin = get_variable(container, var, V)[name, t] + JuMP.add_to_expression!( + expression[name, t], + P_min * prop_term_per_unit * dt, + bin + ) else @assert false end - fuel_expr = - variable[name, t] * prop_term_per_unit * dt + - P_min * bin * prop_term_per_unit * dt JuMP.add_to_expression!( expression[name, t], - fuel_expr, + prop_term_per_unit * dt, + variable[name, t], ) end elseif value_curve isa PSY.QuadraticCurve diff --git a/src/devices_models/devices/common/duration_constraints.jl b/src/devices_models/devices/common/duration_constraints.jl index 12a1f8d7b..105df0df1 100644 --- a/src/devices_models/devices/common/duration_constraints.jl +++ b/src/devices_models/devices/common/duration_constraints.jl @@ -181,7 +181,7 @@ function device_duration_look_ahead!( end end if t <= duration_data[ix].up - lhs_on += get_value(ic) + JuMP.add_to_expression!(lhs_on, get_value(ic)) end con_up[name, t] = JuMP.@constraint( get_jump_model(container), @@ -196,11 +196,12 @@ function device_duration_look_ahead!( lhs_off = JuMP.GenericAffExpr{Float64, JuMP.VariableRef}(0) for i in UnitRange{Int}(Int(t - duration_data[ix].down + 1), t) if i in time_steps - 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]) end end if t <= duration_data[ix].down - lhs_off += get_value(ic) + JuMP.add_to_expression!(lhs_off, get_value(ic)) end con_down[name, t] = JuMP.@constraint( get_jump_model(container), @@ -300,7 +301,7 @@ function device_duration_parameters!( end end if t <= duration_data[ix].up - lhs_on += get_value(ic) + JuMP.add_to_expression!(lhs_on, get_value(ic)) con_up[name, t] = JuMP.@constraint( get_jump_model(container), varstop[name, t] * duration_data[ix].up - lhs_on <= 0.0 @@ -322,14 +323,15 @@ function device_duration_parameters!( for i in UnitRange{Int}(Int(t - duration_data[ix].down + 1), t) if t <= duration_data[ix].down if in(i, time_steps) - 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]) end else JuMP.add_to_expression!(lhs_off, varstop[name, i]) end end if t <= duration_data[ix].down - lhs_off += get_value(ic) + JuMP.add_to_expression!(lhs_off, get_value(ic)) con_down[name, t] = JuMP.@constraint( get_jump_model(container), varstart[name, t] * duration_data[ix].down - lhs_off <= 0.0 diff --git a/src/devices_models/devices/common/objective_function/market_bid.jl b/src/devices_models/devices/common/objective_function/market_bid.jl index 635f59e07..f075c47a1 100644 --- a/src/devices_models/devices/common/objective_function/market_bid.jl +++ b/src/devices_models/devices/common/objective_function/market_bid.jl @@ -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) + for (i, cost) in enumerate(_cost_data) JuMP.add_to_expression!( gen_cost, - cost * multiplier * pwl_var_container[(name, i, time_period)], + cost * multiplier, + pwl_var_container[(name, i, time_period)], ) end return gen_cost @@ -196,7 +197,8 @@ function _get_pwl_cost_expression( for i in 1:length(slopes) JuMP.add_to_expression!( ordc_cost, - slopes[i] * multiplier * pwl_var_container[(name, i, time_period)], + slopes[i] * multiplier, + pwl_var_container[(name, i, time_period)], ) end return ordc_cost @@ -241,7 +243,8 @@ function _get_pwl_cost_expression( for (i, cost) in enumerate(cost_data) JuMP.add_to_expression!( gen_cost, - cost * multiplier * pwl_var_container[(name, i, time_period)], + cost * multiplier, + pwl_var_container[(name, i, time_period)], ) end return gen_cost diff --git a/src/devices_models/devices/common/objective_function/piecewise_linear.jl b/src/devices_models/devices/common/objective_function/piecewise_linear.jl index aff4260c7..268f2fd42 100644 --- a/src/devices_models/devices/common/objective_function/piecewise_linear.jl +++ b/src/devices_models/devices/common/objective_function/piecewise_linear.jl @@ -260,11 +260,12 @@ function _get_pwl_cost_expression( name = PSY.get_name(component) pwl_var_container = get_variable(container, PieceWiseLinearCostVariable(), 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) + for (i, cost) in enumerate(_cost_data) JuMP.add_to_expression!( gen_cost, - cost * multiplier * pwl_var_container[(name, i, time_period)], + cost * multiplier, + pwl_var_container[(name, i, time_period)], ) end return gen_cost diff --git a/src/devices_models/devices/electric_loads.jl b/src/devices_models/devices/electric_loads.jl index 6725b2bd3..dc17fa1f7 100644 --- a/src/devices_models/devices/electric_loads.jl +++ b/src/devices_models/devices/electric_loads.jl @@ -87,8 +87,8 @@ function add_constraints!( name = PSY.get_name(d) pf = sin(atan((PSY.get_max_reactive_power(d) / PSY.get_max_active_power(d)))) reactive = get_variable(container, U(), V)[name, t] - 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) end end diff --git a/src/devices_models/devices/regulation_device.jl b/src/devices_models/devices/regulation_device.jl index aba9df5a6..1fdcc9e97 100644 --- a/src/devices_models/devices/regulation_device.jl +++ b/src/devices_models/devices/regulation_device.jl @@ -227,8 +227,8 @@ function add_constraints!( R_dn[name, t] == (p_factor.dn * area_reserve_dn[agc_name, t]) + R_dn_emergency[name, t] ) - JuMP.add_to_expression!(expr_up[area_name, t], -1 * R_up_emergency[name, t]) - JuMP.add_to_expression!(expr_dn[area_name, t], -1 * R_dn_emergency[name, t]) + JuMP.add_to_expression!(expr_up[area_name, t], -1, R_up_emergency[name, t]) + JuMP.add_to_expression!(expr_dn[area_name, t], -1, R_dn_emergency[name, t]) end end diff --git a/src/parameters/update_cost_parameters.jl b/src/parameters/update_cost_parameters.jl index 259141da6..10f97f05d 100644 --- a/src/parameters/update_cost_parameters.jl +++ b/src/parameters/update_cost_parameters.jl @@ -107,7 +107,8 @@ function _update_pwl_cost_expression( for i in 1:length(cost_data) JuMP.add_to_expression!( gen_cost, - slopes[i] * upb[i] * dt * pwl_var_container[(component_name, i, time_period)], + slopes[i] * upb[i] * dt, + pwl_var_container[(component_name, i, time_period)], ) end return gen_cost @@ -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) add_to_objective_variant_expression!(container, gen_cost) set_expression!(container, ProductionCostExpression, gen_cost, component, time_period) return @@ -160,6 +161,7 @@ function update_variable_cost!( PSY.PiecewiseLinearData(cost_data), ) add_to_objective_variant_expression!(container, mult_ * gen_cost) + # TODO: missing _mult below? set_expression!(container, ProductionCostExpression, gen_cost, component, time_period) return end @@ -180,8 +182,9 @@ function update_variable_cost!( end mult_ = parameter_multiplier[component_name, time_period] expression = get_expression(container, FuelConsumptionExpression(), T) - cost_expr = expression[component_name, time_period] * fuel_cost * mult_ + cost_expr = expression[component_name, time_period] * (fuel_cost * mult_) add_to_objective_variant_expression!(container, cost_expr) + # TODO: missing _mult below? set_expression!(container, ProductionCostExpression, cost_expr, component, time_period) return end diff --git a/src/services_models/agc.jl b/src/services_models/agc.jl index f8e37598d..1cfe0b053 100644 --- a/src/services_models/agc.jl +++ b/src/services_models/agc.jl @@ -152,14 +152,14 @@ function add_constraints!( const_container = add_constraints_container!(container, T(), PSY.System, time_steps) for t in time_steps - system_balance = sum(area_balance.data[:, t]) + system_balance = JuMP.@expression(container.JuMPmodel, sum(@view area_balance.data[:, t])) for agc in agcs a = PSY.get_name(agc) area_name = PSY.get_name(PSY.get_area(agc)) JuMP.add_to_expression!(system_balance, R_up[a, t]) - JuMP.add_to_expression!(system_balance, -1 * R_dn[a, t]) + JuMP.add_to_expression!(system_balance, -1, R_dn[a, t]) JuMP.add_to_expression!(system_balance, R_up_emergency[area_name, t]) - JuMP.add_to_expression!(system_balance, -1 * R_dn_emergency[area_name, t]) + JuMP.add_to_expression!(system_balance, -1, R_dn_emergency[area_name, t]) end const_container[t] = JuMP.@constraint( container.JuMPmodel, @@ -303,6 +303,7 @@ function add_proportional_cost!( ) where {T <: PSY.AGC, U <: LiftVariable} lift_variable = get_variable(container, U(), T) for index in Iterators.product(axes(lift_variable)...) + # TODO add_to_objective_invariant_expression!( container, SERVICES_SLACK_COST * lift_variable[index...], diff --git a/src/services_models/reserve_group.jl b/src/services_models/reserve_group.jl index f02d99af4..dd63e7f18 100644 --- a/src/services_models/reserve_group.jl +++ b/src/services_models/reserve_group.jl @@ -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, + JuMP.@expression(container.JuMPmodel, sum(@view reserve_variable[:, t]))) + # consider a for loop to add the reserve variables end if use_slacks - resource_expression += slack_vars[t] + JuMP.add_to_expression!(resource_expression, slack_vars[t]) end constraint[service_name, t] = JuMP.@constraint(container.JuMPmodel, resource_expression >= requirement) diff --git a/src/services_models/reserves.jl b/src/services_models/reserves.jl index af42f5231..e00ff1315 100644 --- a/src/services_models/reserves.jl +++ b/src/services_models/reserves.jl @@ -161,9 +161,11 @@ 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( + jump_model, sum(@view reserve_variable[:, t]) + slack_vars[t]) else - resource_expression = sum(reserve_variable[:, t]) + resource_expression = JuMP.@expression( + jump_model, sum(@view reserve_variable[:, t])) end constraint[service_name, t] = JuMP.@constraint(jump_model, resource_expression >= param[t] * requirement) @@ -171,9 +173,11 @@ function add_constraints!( else for t in time_steps if use_slacks - resource_expression = sum(reserve_variable[:, t]) + slack_vars[t] + resource_expression = JuMP.@expression( + jump_model, sum(@view reserve_variable[:, t]) + slack_vars[t]) else - resource_expression = sum(reserve_variable[:, t]) + resource_expression = JuMP.@expression( + jump_model, sum(@view reserve_variable[:, t])) end constraint[service_name, t] = JuMP.@constraint( jump_model, @@ -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] ) else cons[name, t] = JuMP.@constraint( jump_model, - var_r[name, t] <= ts_vector[t] * requirement * max_participation_factor + var_r[name, t] <= (requirement * max_participation_factor) * ts_vector[t] ) end end @@ -268,9 +272,11 @@ function add_constraints!( jump_model = get_jump_model(container) for t in time_steps resource_expression = JuMP.GenericAffExpr{Float64, JuMP.VariableRef}() - JuMP.add_to_expression!(resource_expression, sum(reserve_variable[:, t])) + JuMP.add_to_expression!(resource_expression, + JuMP.@expression(jump_model, sum(@view reserve_variable[:, t]))) + # consider a for loop if use_slacks - resource_expression += slack_vars[t] + JuMP.add_to_expression!(resource_expression, slack_vars[t]) end constraint[service_name, t] = JuMP.@constraint(jump_model, resource_expression >= requirement) @@ -316,7 +322,7 @@ function add_constraints!( for t in time_steps constraint[service_name, t] = JuMP.@constraint( jump_model, - sum(reserve_variable[:, t]) >= requirement_variable[service_name, t] + sum(@view reserve_variable[:, t]) >= requirement_variable[service_name, t] ) end @@ -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 end for t in time_steps cons[name, t] = JuMP.@constraint( @@ -553,6 +559,7 @@ function add_proportional_cost!( for index in Iterators.product(axes(reserve_variable)...) add_to_objective_invariant_expression!( container, + # possibly decouple DEFAULT_RESERVE_COST / base_p * reserve_variable[index...], ) end