diff --git a/src/WaterModels.jl b/src/WaterModels.jl index 2142db73..07787850 100644 --- a/src/WaterModels.jl +++ b/src/WaterModels.jl @@ -99,6 +99,7 @@ include("util/variable_index.jl") include("util/bound_problem.jl") include("util/pairwise_cuts.jl") include("util/obbt.jl") +include("util/fixing_variables.jl") # Deprecated functions. include("deprecated.jl") diff --git a/src/core/objective.jl b/src/core/objective.jl index 70486fec..44a52761 100644 --- a/src/core/objective.jl +++ b/src/core/objective.jl @@ -193,7 +193,7 @@ function objective_max_scaled_demand(wm::AbstractWaterModel)::JuMP.AffExpr # Initialize the objective expression to zero. objective = JuMP.AffExpr(0.0) - scale = 1.0 + scale = 0.0 for n in network_ids_flow # Get the set of dispatchable demands at time index `n`. dispatchable_demands = ref(wm, n, :dispatchable_demand) @@ -205,6 +205,9 @@ function objective_max_scaled_demand(wm::AbstractWaterModel)::JuMP.AffExpr JuMP.add_to_expression!(objective, var(wm, n, :q_demand, i)) end end + if(scale==0.0) + scale = 1.0 + end println("Using this scale for water objective $(scale)*****") # Maximize the total amount of prioritized water volume delivered. return JuMP.@objective(wm.model, JuMP.MAX_SENSE, objective/scale) diff --git a/src/form/crd.jl b/src/form/crd.jl index 3fc67af3..09b56ae9 100644 --- a/src/form/crd.jl +++ b/src/form/crd.jl @@ -163,57 +163,57 @@ function constraint_on_off_pump_head_gain_ne( end -function constraint_on_off_pump_power( - wm::AbstractCRDModel, - n::Int, - a::Int, - q_min_forward::Float64, -) - # Gather pump flow, power, and status variables. - q, P, z = var(wm, n, :qp_pump, a), var(wm, n, :P_pump, a), var(wm, n, :z_pump, a) - - # Compute pump flow and power partitioning. - q_lb, q_ub = q_min_forward, JuMP.upper_bound(q) - - if q_lb == 0.0 && q_ub == 0.0 - c = JuMP.@constraint(wm.model, P == 0.0) - append!(con(wm, n, :on_off_pump_power)[a], [c]) - else - f_ua = _calc_pump_power_ua(wm, n, a, [q_lb, q_ub]) - if f_ua[1] != f_ua[2] - # Build a linear under-approximation of the power. - slope = (f_ua[2] - f_ua[1]) / (q_ub - q_lb) - power_expr = slope * (q - q_lb * z) + f_ua[1] * z - c = JuMP.@constraint(wm.model, power_expr <= P) - append!(con(wm, n, :on_off_pump_power)[a], [c]) - end - end -end - - -function constraint_on_off_pump_power_ne( - wm::AbstractCRDModel, - n::Int, - a::Int, - q_min_forward::Float64, -) - # Gather pump flow, power, and status variables. - q, P, z = var(wm, n, :qp_ne_pump, a), var(wm, n, :P_ne_pump, a), var(wm, n, :z_ne_pump, a) - - # Compute pump flow and power partitioning. - q_lb, q_ub = q_min_forward, JuMP.upper_bound(q) - - if q_lb == 0.0 && q_ub == 0.0 - c = JuMP.@constraint(wm.model, P == 0.0) - append!(con(wm, n, :on_off_pump_power_ne)[a], [c]) - else - f_ua = _calc_pump_power_ua_ne(wm, n, a, [q_lb, q_ub]) - if f_ua[1] != f_ua[2] - # Build a linear under-approximation of the power. - slope = (f_ua[2] - f_ua[1]) / (q_ub - q_lb) - power_expr = slope * (q - q_lb * z) + f_ua[1] * z - c = JuMP.@constraint(wm.model, power_expr <= P) - append!(con(wm, n, :on_off_pump_power_ne)[a], [c]) - end - end -end +# function constraint_on_off_pump_power( +# wm::AbstractCRDModel, +# n::Int, +# a::Int, +# q_min_forward::Float64, +# ) +# # Gather pump flow, power, and status variables. +# q, P, z = var(wm, n, :qp_pump, a), var(wm, n, :P_pump, a), var(wm, n, :z_pump, a) +# +# # Compute pump flow and power partitioning. +# q_lb, q_ub = q_min_forward, JuMP.upper_bound(q) +# +# if q_lb == 0.0 && q_ub == 0.0 +# c = JuMP.@constraint(wm.model, P == 0.0) +# append!(con(wm, n, :on_off_pump_power)[a], [c]) +# else +# f_ua = _calc_pump_power_ua(wm, n, a, [q_lb, q_ub]) +# if f_ua[1] != f_ua[2] +# # Build a linear under-approximation of the power. +# slope = (f_ua[2] - f_ua[1]) / (q_ub - q_lb) +# power_expr = slope * (q - q_lb * z) + f_ua[1] * z +# c = JuMP.@constraint(wm.model, power_expr <= P) +# append!(con(wm, n, :on_off_pump_power)[a], [c]) +# end +# end +# end +# +# +# function constraint_on_off_pump_power_ne( +# wm::AbstractCRDModel, +# n::Int, +# a::Int, +# q_min_forward::Float64, +# ) +# # Gather pump flow, power, and status variables. +# q, P, z = var(wm, n, :qp_ne_pump, a), var(wm, n, :P_ne_pump, a), var(wm, n, :z_ne_pump, a) +# +# # Compute pump flow and power partitioning. +# q_lb, q_ub = q_min_forward, JuMP.upper_bound(q) +# +# if q_lb == 0.0 && q_ub == 0.0 +# c = JuMP.@constraint(wm.model, P == 0.0) +# append!(con(wm, n, :on_off_pump_power_ne)[a], [c]) +# else +# f_ua = _calc_pump_power_ua_ne(wm, n, a, [q_lb, q_ub]) +# if f_ua[1] != f_ua[2] +# # Build a linear under-approximation of the power. +# slope = (f_ua[2] - f_ua[1]) / (q_ub - q_lb) +# power_expr = slope * (q - q_lb * z) + f_ua[1] * z +# c = JuMP.@constraint(wm.model, power_expr <= P) +# append!(con(wm, n, :on_off_pump_power_ne)[a], [c]) +# end +# end +# end diff --git a/src/form/nc.jl b/src/form/nc.jl index 5a3b687f..93f11c18 100644 --- a/src/form/nc.jl +++ b/src/form/nc.jl @@ -190,8 +190,12 @@ function constraint_pipe_head_loss( q, h_i, h_j = var(wm, n, :q_pipe, a), var(wm, n, :h, node_fr), var(wm, n, :h, node_to) # Add nonconvex constraint for the head loss relationship. - c_1 = JuMP.@NLconstraint(wm.model, r * head_loss(q) <= (h_i - h_j) / L) - c_2 = JuMP.@NLconstraint(wm.model, r * head_loss(q) >= (h_i - h_j) / L) + head_loss_form = wm.ref[:it][wm_it_sym][:head_loss] + p = uppercase(head_loss_form) == "H-W" ? 1.852 : 2.0 + alpha = (p-1) + h_l = JuMP.@expression(wm.model, q*(abs(q)^alpha)) + c_1 = JuMP.@NLconstraint(wm.model, r * h_l <= (h_i - h_j) / L) + c_2 = JuMP.@NLconstraint(wm.model, r * h_l >= (h_i - h_j) / L) # Append the :pipe_head_loss constraint array. append!(con(wm, n, :pipe_head_loss)[a], [c_1, c_2]) @@ -604,42 +608,42 @@ subnetwork (or time) index that is considered, `a` is the index of the pump, and `q_min_forward` is the _minimum_ (positive) amount of flow when the pump is active. Note that only a quadratic approximation is used for `AbstractNCModel`. """ -function constraint_on_off_pump_power( - wm::AbstractNCModel, - n::Int, - a::Int, - q_min_forward::Float64, -) - # Gather pump flow, scaled power, and status variables. - q, P, z = var(wm, n, :q_pump, a), var(wm, n, :P_pump, a), var(wm, n, :z_pump, a) - - # Add constraint equating power with respect to the power curve. - power_qa = _calc_pump_power_quadratic_approximation(wm, n, a, z) - c_1 = JuMP.@constraint(wm.model, power_qa(q) <= P) - c_2 = JuMP.@constraint(wm.model, power_qa(q) >= P) - - # Append the :on_off_pump_power constraint array. - append!(con(wm, n, :on_off_pump_power)[a], [c_1, c_2]) -end - - -function constraint_on_off_pump_power_ne( - wm::AbstractNCModel, - n::Int, - a::Int, - q_min_forward::Float64, -) - # Gather pump flow, scaled power, and status variables. - q, P, z = var(wm, n, :q_ne_pump, a), var(wm, n, :P_ne_pump, a), var(wm, n, :z_ne_pump, a) - - # Add constraint equating power with respect to the power curve. - power_qa = _calc_pump_power_quadratic_approximation_ne(wm, n, a, z) - c_1 = JuMP.@constraint(wm.model, power_qa(q) <= P) - c_2 = JuMP.@constraint(wm.model, power_qa(q) >= P) - - # Append the :on_off_pump_power constraint array. - append!(con(wm, n, :on_off_pump_power_ne)[a], [c_1, c_2]) -end +# function constraint_on_off_pump_power( +# wm::AbstractNCModel, +# n::Int, +# a::Int, +# q_min_forward::Float64, +# ) +# # Gather pump flow, scaled power, and status variables. +# q, P, z = var(wm, n, :q_pump, a), var(wm, n, :P_pump, a), var(wm, n, :z_pump, a) +# +# # Add constraint equating power with respect to the power curve. +# power_qa = _calc_pump_power_quadratic_approximation(wm, n, a, z) +# c_1 = JuMP.@constraint(wm.model, power_qa(q) <= P) +# c_2 = JuMP.@constraint(wm.model, power_qa(q) >= P) +# +# # Append the :on_off_pump_power constraint array. +# append!(con(wm, n, :on_off_pump_power)[a], [c_1, c_2]) +# end +# +# +# function constraint_on_off_pump_power_ne( +# wm::AbstractNCModel, +# n::Int, +# a::Int, +# q_min_forward::Float64, +# ) +# # Gather pump flow, scaled power, and status variables. +# q, P, z = var(wm, n, :q_ne_pump, a), var(wm, n, :P_ne_pump, a), var(wm, n, :z_ne_pump, a) +# +# # Add constraint equating power with respect to the power curve. +# power_qa = _calc_pump_power_quadratic_approximation_ne(wm, n, a, z) +# c_1 = JuMP.@constraint(wm.model, power_qa(q) <= P) +# c_2 = JuMP.@constraint(wm.model, power_qa(q) >= P) +# +# # Append the :on_off_pump_power constraint array. +# append!(con(wm, n, :on_off_pump_power_ne)[a], [c_1, c_2]) +# end """ @@ -848,6 +852,7 @@ function constraint_short_pipe_flow_ne( q_max_reverse::Float64, q_min_forward::Float64, ) +# println("Using NC ******* short pipe *********") # Get flow and status variables for the short pipe. q = var(wm, n, :q_ne_short_pipe, a) z = var(wm, n, :z_ne_short_pipe, a) diff --git a/src/form/ncd.jl b/src/form/ncd.jl index 5a39887a..f6881a29 100644 --- a/src/form/ncd.jl +++ b/src/form/ncd.jl @@ -738,7 +738,6 @@ function constraint_on_off_pump_power( # Add constraint equating power with respect to the power curve. power_la = _calc_pump_power_linear_coeff(wm, n, z) - println("power_la = $power_la") # power_la = _calc_pump_power_quadratic_approximation(wm, n, a, z) c_1 = JuMP.@constraint(wm.model, power_la(q) <= P) c_2 = JuMP.@constraint(wm.model, power_la(q) >= P) @@ -904,6 +903,7 @@ function constraint_short_pipe_flow_ne( q_max_reverse::Float64, q_min_forward::Float64, ) +# println("Using NCD ******* short pipe *********") # Get expansion short pipe flow, direction, and status variables. qp, qn = var(wm, n, :qp_ne_short_pipe, a), var(wm, n, :qn_ne_short_pipe, a) y, z = var(wm, n, :y_ne_short_pipe, a), var(wm, n, :z_ne_short_pipe, a) diff --git a/src/form/pwlrd.jl b/src/form/pwlrd.jl index a00ca6d2..a3cdeadc 100644 --- a/src/form/pwlrd.jl +++ b/src/form/pwlrd.jl @@ -43,7 +43,7 @@ function variable_flow(wm::AbstractPWLRDModel; nw::Int=nw_id_default, bounded::B # Create directed head difference (`dhp` and `dhn`) variables for each component. _variable_component_head_difference(wm, name; nw = nw, bounded = bounded, report = report) end - + # Create variables involved in convex combination constraints for pumps. pump_partition = Dict{Int, Vector{Float64}}(a => pump["flow_partition"] for (a, pump) in ref(wm, nw, :pump)) @@ -58,6 +58,20 @@ function variable_flow(wm::AbstractPWLRDModel; nw::Int=nw_id_default, bounded::B base_name = "$(nw)_x", binary = true, start = comp_start_value(ref(wm, nw, :pump, a), "x_start")) + # Create variables involved in convex combination constraints for expansion pumps. + pump_partition = Dict{Int, Vector{Float64}}(a => + pump["flow_partition"] for (a, pump) in ref(wm, nw, :ne_pump)) + + var(wm, nw)[:lambda_ne_pump] = JuMP.@variable(wm.model, + [a in ids(wm, nw, :ne_pump), k in 1:length(pump_partition[a])], + base_name = "$(nw)_ne_lambda", lower_bound = 0.0, upper_bound = 1.0, + start = comp_start_value(ref(wm, nw, :ne_pump, a), "lambda_start", k)) + + var(wm, nw)[:x_con_ne_pump] = JuMP.@variable(wm.model, + [a in ids(wm, nw, :ne_pump), k in 1:length(pump_partition[a]) - 1], + base_name = "$(nw)_ne_x_con", binary = true, + start = comp_start_value(ref(wm, nw, :ne_pump, a), "x_con_ne_start")) + # Create variables involved in convex combination constraints for pipes. pipe_partition_p = Dict{Int, Vector{Float64}}(a => get_pipe_flow_partition_positive(pipe) @@ -110,7 +124,7 @@ function variable_flow(wm::AbstractPWLRDModel; nw::Int=nw_id_default, bounded::B [a in ids(wm, nw, :des_pipe), k in 1:length(des_pipe_partition_n[a])], base_name = "$(nw)_lambda_n", lower_bound = 0.0, upper_bound = 1.0, start = comp_start_value(ref(wm, nw, :des_pipe, a), "lambda_n_start", k)) - + var(wm, nw)[:x_n_des_pipe] = JuMP.@variable(wm.model, [a in ids(wm, nw, :des_pipe), k in 1:length(des_pipe_partition_n[a])-1], base_name = "$(nw)_x_n", binary = true, @@ -205,7 +219,7 @@ function constraint_pipe_head_loss( for flow_value in filter(x -> x > 0.0, partition_p) # Add head loss outer (i.e., lower) approximations. lhs_p = r * _calc_head_loss_oa(qp, y, flow_value, exponent) - + if minimum(abs.(lhs_p.terms.vals)) >= 1.0e-4 scalar = _get_scaling_factor(vcat(lhs_p.terms.vals, [1.0 / L])) c_7 = JuMP.@constraint(wm.model, scalar * lhs_p <= scalar * dhp / L) @@ -243,7 +257,7 @@ function constraint_pipe_head_loss( c_13 = JuMP.@constraint(wm.model, lambda_n[a, bn_range[end]] <= x_n[a, bn_range_m1[end]]) append!(con(wm, n, :pipe_head_loss, a), [c_11, c_12, c_13]) end - + # Add a constraint that upper-bounds the head loss variable. if maximum(partition_n) != 0.0 f_n = r .* partition_n.^exponent @@ -412,7 +426,7 @@ function constraint_on_off_des_pipe_head_loss( lhs_n_2 = r * _calc_head_loss_oa(qn, z, flow_value, exponent) scalar = _get_scaling_factor(vcat(lhs_n_2.terms.vals, [1.0 / L])) c_17 = JuMP.@constraint(wm.model, scalar * lhs_n_2 <= scalar * dhn / L) - + append!(con(wm, n, :on_off_des_pipe_head_loss, a), [c_16, c_17]) end @@ -428,7 +442,7 @@ function constraint_on_off_des_pipe_head_loss( c_19 = JuMP.@constraint(wm.model, x_p_sum + x_n_sum == z) lambda_p_sum = sum(lambda_p[a, k] for k in bp_range) - lambda_n_sum = sum(lambda_n[a, k] for k in bn_range) + lambda_n_sum = sum(lambda_n[a, k] for k in bn_range) c_20 = JuMP.@constraint(wm.model, lambda_p_sum + lambda_n_sum == z) append!(con(wm, n, :on_off_des_pipe_head_loss, a), [c_19, c_20]) end @@ -523,50 +537,135 @@ end """ - constraint_on_off_pump_power( + constraint_on_off_pump_head_gain( wm::AbstractPWLRDModel, n::Int, a::Int, + node_fr::Int, + node_to::Int, q_min_forward::Float64 ) -Adds constraints that model the pump's power consumption, if operating, as -being bounded between piecewise-linear inner and outer approximations of a -(potentially) nonlinear function of a more accurate model. If the pump is -inactive, the power is restricted to a value of zero. Here, `wm` is the -WaterModels object, `n` is the subnetwork (or time) index that is considered, -`a` is the index of the pump, and `q_min_forward` is the _minimum_ (positive) -amount of flow when the pump is active. +Adds constraints that model the pump's head gain, if operating, via outer +and piecewise-linear inner approximations of the nonlinear function of flow +rate. If the pump is inactive, the head gain is restricted to a value of zero. +Here, `wm` is the WaterModels object, `n` is the subnetwork (or time) index +that is considered, `a` is the index of the pump, `node_fr` is the index of the +tail node of the pump, `node_to` is the index of the head node of the pump, and +`q_min_forward` is the _minimum_ (positive) amount of flow when the pump is +active. Head gain is assumed to be nonnegative and is directed from `node_fr` +to `node_to`. """ -function constraint_on_off_pump_power( +function constraint_on_off_pump_head_gain_ne( wm::AbstractPWLRDModel, n::Int, a::Int, + node_fr::Int, + node_to::Int, q_min_forward::Float64 ) - # Gather pump power and status variables. - P, lambda = var(wm, n, :P_pump, a), var(wm, n, :lambda_pump) + # Gather pump flow, head gain, and status variables. + qp, g, z = var(wm, n, :qp_ne_pump, a), var(wm, n, :g_ne_pump, a), var(wm, n, :z_ne_pump, a) + + # Gather convex combination variables. + lambda, x = var(wm, n, :lambda_ne_pump), var(wm, n, :x_con_ne_pump) # Get the corresponding flow partitioning. - @assert haskey(ref(wm, n, :pump, a), "flow_partition") - partition = ref(wm, n, :pump, a, "flow_partition") - bp_range = 1:length(partition) - - # Add a constraint that lower-bounds the power variable. - if minimum(partition) == 0.0 && maximum(partition) == 0.0 - c_1 = JuMP.@constraint(wm.model, P == 0.0) - append!(con(wm, n, :on_off_pump_power)[a], [c_1]) - else - f_ua = _calc_pump_power_ua(wm, n, a, partition) - power_lb_expr = sum(f_ua[k] * lambda[a, k] for k in bp_range) - c_1 = JuMP.@constraint(wm.model, power_lb_expr <= P) + @assert haskey(ref(wm, n, :ne_pump, a), "flow_partition") + partition = ref(wm, n, :ne_pump, a, "flow_partition") + bp_range, bp_range_m1 = 1:length(partition), 1:length(partition)-1 + + # Add the required SOS constraints. + c_1 = JuMP.@constraint(wm.model, sum(lambda[a, k] for k in bp_range) == z) + c_2 = JuMP.@constraint(wm.model, sum(x[a, k] for k in bp_range_m1) == z) + c_3 = JuMP.@constraint(wm.model, lambda[a, 1] <= x[a, 1]) + c_4 = JuMP.@constraint(wm.model, lambda[a, bp_range[end]] <= x[a, bp_range_m1[end]]) + + # Add a constraint for the flow piecewise approximation. + qp_lhs = sum(partition[k] * lambda[a, k] for k in bp_range) + scalar = _get_scaling_factor(vcat(qp_lhs.terms.vals, [1.0])) + c_5 = JuMP.@constraint(wm.model, scalar * qp_lhs == scalar * qp) + + # Get pump head curve function and its derivative. + head_curve_function = ref(wm, n, :ne_pump, a, "head_curve_function") + head_curve_derivative = ref(wm, n, :ne_pump, a, "head_curve_derivative") + + # Add a constraint that lower-bounds the head gain variable. + f_all = map(x -> x = x < 0.0 ? 0.0 : x, head_curve_function.(partition)) + gain_lb_expr = sum(f_all[k] .* lambda[a, k] for k in bp_range) + scalar = _get_scaling_factor(vcat(gain_lb_expr.terms.vals, [1.0])) + c_6 = JuMP.@constraint(wm.model, scalar * gain_lb_expr <= scalar * g) - # Add a constraint that upper-bounds the power variable. - f_oa = _calc_pump_power_oa(wm, n, a, partition) - power_ub_expr = sum(f_oa[k] * lambda[a, k] for k in bp_range) - c_2 = JuMP.@constraint(wm.model, P <= power_ub_expr) + # Append the constraint array. + append!(con(wm, n, :on_off_pump_head_gain_ne, a), [c_1, c_2, c_3, c_4, c_5, c_6]) + + for flow_value in filter(q_val -> q_val > 0.0, partition) + # Add head gain outer (i.e., upper) approximations. + f = head_curve_function(flow_value) + df = head_curve_derivative(flow_value) + + if abs(df) >= 1.0e-4 # Only add an outer-approximation if the derivative isn't too small. + f_rhs = f * z + df * (qp - flow_value * z) + scalar = _get_scaling_factor(vcat(f_rhs.terms.vals, [1.0])) + c = JuMP.@constraint(wm.model, scalar * g <= scalar * f_rhs) + append!(con(wm, n, :on_off_pump_head_gain_ne, a), [c]) + end + end - # Append the :on_off_pump_power constraint array. - append!(con(wm, n, :on_off_pump_power)[a], [c_1, c_2]) + for k in 2:length(partition)-1 + # Add the adjacency constraints for piecewise variables. + adjacency = x[a, k-1] + x[a, k] + c = JuMP.@constraint(wm.model, lambda[a, k] <= adjacency) + append!(con(wm, n, :on_off_pump_head_gain_ne, a), [c]) end end + + +""" + constraint_on_off_pump_power( + wm::AbstractPWLRDModel, + n::Int, + a::Int, + q_min_forward::Float64 + ) + +Adds constraints that model the pump's power consumption, if operating, as +being bounded between piecewise-linear inner and outer approximations of a +(potentially) nonlinear function of a more accurate model. If the pump is +inactive, the power is restricted to a value of zero. Here, `wm` is the +WaterModels object, `n` is the subnetwork (or time) index that is considered, +`a` is the index of the pump, and `q_min_forward` is the _minimum_ (positive) +amount of flow when the pump is active. +""" +# function constraint_on_off_pump_power( +# wm::AbstractPWLRDModel, +# n::Int, +# a::Int, +# q_min_forward::Float64 +# ) +# # Gather pump power and status variables. +# P, lambda = var(wm, n, :P_pump, a), var(wm, n, :lambda_pump) +# +# # Get the corresponding flow partitioning. +# @assert haskey(ref(wm, n, :pump, a), "flow_partition") +# partition = ref(wm, n, :pump, a, "flow_partition") +# bp_range = 1:length(partition) +# +# # Add a constraint that lower-bounds the power variable. +# if minimum(partition) == 0.0 && maximum(partition) == 0.0 +# c_1 = JuMP.@constraint(wm.model, P == 0.0) +# append!(con(wm, n, :on_off_pump_power)[a], [c_1]) +# else +# f_ua = _calc_pump_power_ua(wm, n, a, partition) +# power_lb_expr = sum(f_ua[k] * lambda[a, k] for k in bp_range) +# c_1 = JuMP.@constraint(wm.model, power_lb_expr <= P) +# +# # Add a constraint that upper-bounds the power variable. +# f_oa = _calc_pump_power_oa(wm, n, a, partition) +# power_ub_expr = sum(f_oa[k] * lambda[a, k] for k in bp_range) +# c_2 = JuMP.@constraint(wm.model, P <= power_ub_expr) +# +# # Append the :on_off_pump_power constraint array. +# append!(con(wm, n, :on_off_pump_power)[a], [c_1, c_2]) +# end +# end diff --git a/src/prob/wf.jl b/src/prob/wf.jl index 9b21f6d8..b023fc3a 100644 --- a/src/prob/wf.jl +++ b/src/prob/wf.jl @@ -71,6 +71,7 @@ function build_wf(wm::AbstractWaterModel) constraint_on_off_pump_head(wm, a) constraint_on_off_pump_head_gain(wm, a) constraint_on_off_pump_flow(wm, a) + # println("*******Disabling pump power constraint*****************") constraint_on_off_pump_power(wm, a) end @@ -160,6 +161,13 @@ function build_mn_wf(wm::AbstractWaterModel) variable_reservoir_flow(wm; nw=n) variable_tank_flow(wm; nw=n) + # + if(haskey(wm.ref[:it][wm_it_sym],:sol_from_relaxation)) + sol_from_relaxation = wm.ref[:it][wm_it_sym][:sol_from_relaxation] + println("Fixing Binary Values from relaxed solution") + fix_variables_to_relaxed_solutions(wm, sol_from_relaxation; nw=n, continuous_fixing = false) + end + # Flow conservation at all nodes. for i in ids(wm, :node; nw=n) constraint_flow_conservation(wm, i; nw=n) @@ -227,6 +235,7 @@ function build_mn_wf(wm::AbstractWaterModel) constraint_on_off_valve_head(wm, a; nw=n) constraint_on_off_valve_flow(wm, a; nw=n) end + end # Start with the first network, representing the initial time step. @@ -268,6 +277,7 @@ function build_mn_wf(wm::AbstractWaterModel) end end + # println(wm.model) # Add the objective. objective_wf(wm) end diff --git a/src/util/fixing_variables.jl b/src/util/fixing_variables.jl new file mode 100644 index 00000000..45fc057f --- /dev/null +++ b/src/util/fixing_variables.jl @@ -0,0 +1,154 @@ +function fix_variables_to_relaxed_solutions(wm::AbstractWaterModel, sol_val_dict::Dict{String, Any}; nw::Int=nw_id_default, continuous_fixing::Bool = false) +# println("Using relaxed solutions") +# println(sol_val_dict) + +############################### +# Integer Variables ########### +############################### +fixing_direction = true + if(haskey(ref(wm,nw),:pipe)) + # println("fixing pipe flow direction") + for (a, pipe) in ref(wm, nw, :pipe) + if(fixing_direction == true) + val = sol_val_dict["nw"][string(nw)]["pipe"][string(a)]["y"] + y_pipe = var(wm, nw, :y_pipe, a) + JuMP.unset_binary(y_pipe) + JuMP.fix(y_pipe, val; force = true) + end + # partition = pipe["flow_partition"] + # qp = var(wm, nw, :qp_pipe, a) + # qn = var(wm, nw, :qn_pipe, a) + # qval = sol_val_dict["nw"][string(nw)]["pipe"][string(a)]["q"] + # for i in 1:length(partition)-1 + # if(qval >= partition[i] && qval <= partition[i+1]) + # if(qval >=0) + # JuMP.set_lower_bound(qp,partition[i]) + # JuMP.set_upper_bound(qp,partition[i+1]) + # else + # JuMP.set_lower_bound(qn,partition[i]) + # JuMP.set_upper_bound(qn,partition[i+1]) + # end + # end + # end + end + end + + if(haskey(ref(wm,nw),:short_pipe)) + # println("fixing pipe flow direction") + for (a, short_pipe) in ref(wm, nw, :short_pipe) + if(fixing_direction == true) + val = sol_val_dict["nw"][string(nw)]["short_pipe"][string(a)]["y"] + val = Int(round(val)) + y_short_pipe = var(wm, nw, :y_short_pipe, a) + # println("Unsetting binary for y pipe") + JuMP.unset_binary(y_short_pipe) + JuMP.fix(y_short_pipe, val; force = true) + end + + end + end + if(haskey(ref(wm,nw),:pump)) + println("fixing pump activation status") + # for (a, pump) in ref(wm, nw, :pump) + # val = sol_val_dict["nw"][string(nw)]["pump"][string(a)]["status"] + # val = Int(round(val)) + # z_pump = var(wm, nw, :z_pump, a) + # JuMP.unset_binary(z_pump) + # JuMP.fix(z_pump, val; force = true) + # end + println("fixing pump flow direction") + for (a, pump) in ref(wm, nw, :pump) + # if(fixing_direction == true) + # val = sol_val_dict["nw"][string(nw)]["pump"][string(a)]["y"] + # val = Int(round(val)) + # y_pump = var(wm, nw, :y_pump, a) + # JuMP.unset_binary(y_pump) + # JuMP.fix(y_pump, val; force = true) + # end + + # partition = pump["flow_partition"] + # qp_pump = var(wm, nw, :qp_pump, a) + # # qn_pump = var(wm, nw, :qn_pump, a) + # qpumpval = val = sol_val_dict["nw"][string(nw)]["pump"][string(a)]["q"] + # for i in 1:length(partition)-1 + # if(qpumpval >= partition[i] && qpumpval <= partition[i+1]) + # if(qpumpval >=0) + # JuMP.set_lower_bound(qp_pump,partition[i]) + # JuMP.set_upper_bound(qp_pump,partition[i+1]) + # else + # JuMP.set_lower_bound(qn_pump,partition[i]) + # JuMP.set_upper_bound(qn_pump,partition[i+1]) + # end + # end + # end + end + end + # + # + if(haskey(ref(wm,nw),:ne_pump)) + # println("fixing ne pump build status") + for (a, ne_pump) in ref(wm, nw, :ne_pump) + val = sol_val_dict["nw"][string(nw)]["ne_pump"][string(a)]["build_status"] + val = Int(round(val)) + x_ne_pump = var(wm, nw, :x_ne_pump, a) + JuMP.unset_binary(x_ne_pump) + JuMP.fix(x_ne_pump, val; force = true) + end + # println("fixing ne pump activation status") + for (a, ne_pump) in ref(wm, nw, :ne_pump) + val = sol_val_dict["nw"][string(nw)]["ne_pump"][string(a)]["status"] + val = Int(round(val)) + z_ne_pump = var(wm, nw, :z_ne_pump, a) + JuMP.unset_binary(z_ne_pump) + JuMP.fix(z_ne_pump, val; force = true) + + # # println("fixing ne pump flow direction") + if(fixing_direction == true) + val = sol_val_dict["nw"][string(nw)]["ne_pump"][string(a)]["y"] + val = Int(round(val)) + y_ne_pump = var(wm, nw, :y_ne_pump, a) + JuMP.unset_binary(y_ne_pump) + JuMP.fix(y_ne_pump, val; force = true) + end + end + end + # + if(haskey(ref(wm,nw),:ne_short_pipe)) + # println("fixing ne short pipe flow direction") + for (a, ne_short_pipe) in ref(wm, nw, :ne_short_pipe) + # if(fixing_direction == true) + # val = sol_val_dict["nw"][string(nw)]["ne_short_pipe"][string(a)]["y"] + # val = Int(round(val)) + # y_ne_short_pipe = var(wm, nw, :y_ne_short_pipe, a) + # JuMP.unset_binary(y_ne_short_pipe) + # JuMP.fix(y_ne_short_pipe, val; force = true) + # end + + # println("fixing ne short pipe build status") + val = sol_val_dict["nw"][string(nw)]["ne_short_pipe"][string(a)]["status"] + val = Int(round(val)) + z_ne_short_pipe = var(wm, nw, :z_ne_short_pipe, a) + JuMP.unset_binary(z_ne_short_pipe) + JuMP.fix(z_ne_short_pipe, val; force = true) + + end + end + if(haskey(ref(wm,nw),:valve)) + for (a, valve) in ref(wm, nw, :valve) + val = sol_val_dict["nw"][string(nw)]["valve"][string(a)]["status"] + val = Int(round(val)) + z_valve = var(wm, nw, :z_valve, a) + JuMP.unset_binary(z_valve) + JuMP.fix(z_valve, val; force = true) + end + for (a, valve) in ref(wm, nw, :valve) + if(fixing_direction == true) + yval = sol_val_dict["nw"][string(nw)]["valve"][string(a)]["y"] + yval = Int(round(yval)) + y_valve = var(wm, nw, :y_valve, a) + JuMP.unset_binary(y_valve) + JuMP.fix(y_valve, yval; force = true) + end + end + end +end