diff --git a/src/solution/methods/ProgressiveHedging/optimize.jl b/src/solution/methods/ProgressiveHedging/optimize.jl index c08e391..03831d8 100644 --- a/src/solution/methods/ProgressiveHedging/optimize.jl +++ b/src/solution/methods/ProgressiveHedging/optimize.jl @@ -6,106 +6,95 @@ using TimerOutputs import JuMP const to = TimerOutput() -function optimize!(model::JuMP.Model, method::ProgressiveHedging)::FinalResult +function optimize!(model::JuMP.Model, method::ProgressiveHedging)::Nothing mpi = MpiInfo(MPI.COMM_WORLD) - iterations = Array{IterationInfo,1}(undef, 0) - if method.consensus_vars === nothing - method.consensus_vars = - [var for var in all_variables(model) if is_binary(var)] + iterations = PHIterationInfo[] + consensus_vars = [var for var in all_variables(model) if is_binary(var)] + nvars = length(consensus_vars) + weights = ones(nvars) + if method.initial_weights !== nothing + weights = copy(method.initial_weights) end - nvars = length(method.consensus_vars) - if method.weights === nothing - method.weights = [1.0 for _ in 1:nvars] + target = zeros(nvars) + if method.initial_target !== nothing + target = copy(method.initial_target) end - if method.initial_global_consensus_vals === nothing - method.initial_global_consensus_vals = [0.0 for _ in 1:nvars] - end - - ph_sp_params = SpParams( + params = PHSubProblemParams( ρ = method.ρ, - λ = [method.λ_default for _ in 1:nvars], - global_consensus_vals = method.initial_global_consensus_vals, + λ = [method.λ for _ in 1:nvars], + target = target, ) - ph_subproblem = - SubProblem(model, model[:obj], method.consensus_vars, method.weights) - set_optimizer_attribute(model, "Threads", method.num_of_threads) + sp = PHSubProblem(model, model[:obj], consensus_vars, weights) while true - it_time = @elapsed begin - solution = solve_subproblem(ph_subproblem, ph_sp_params) + iteration_time = @elapsed begin + solution = solve_subproblem(sp, params, method.inner_method) MPI.Barrier(mpi.comm) global_obj = compute_global_objective(mpi, solution) - global_consensus_vals = compute_global_consensus(mpi, solution) - update_λ_and_residuals!( - solution, - ph_sp_params, - global_consensus_vals, - ) + target = compute_target(mpi, solution) + update_λ_and_residuals!(solution, params, target) global_infeas = compute_global_infeasibility(solution, mpi) global_residual = compute_global_residual(mpi, solution) - if has_numerical_issues(global_consensus_vals) + if has_numerical_issues(target) break end end - total_elapsed_time = compute_total_elapsed_time(it_time, iterations) - it = IterationInfo( - it_num = length(iterations) + 1, - sp_consensus_vals = solution.consensus_vals, - global_consensus_vals = global_consensus_vals, - sp_obj = solution.obj, + total_elapsed_time = + compute_total_elapsed_time(iteration_time, iterations) + current_iteration = PHIterationInfo( + global_infeas = global_infeas, global_obj = global_obj, - it_time = it_time, - total_elapsed_time = total_elapsed_time, global_residual = global_residual, - global_infeas = global_infeas, + iteration_number = length(iterations) + 1, + iteration_time = iteration_time, + sp_vals = solution.vals, + sp_obj = solution.obj, + target = target, + total_elapsed_time = total_elapsed_time, ) - iterations = [iterations; it] - print_progress(mpi, it, method.print_interval) - if should_stop(mpi, iterations, method.termination_criteria) + push!(iterations, current_iteration) + print_progress(mpi, current_iteration, method.print_interval) + if should_stop(mpi, iterations, method.termination) break end end - - return FinalResult( - last(iterations).global_obj, - last(iterations).sp_consensus_vals, - last(iterations).global_infeas, - last(iterations).it_num, - last(iterations).total_elapsed_time, - ) + return end function compute_total_elapsed_time( - it_time::Float64, - iterations::Array{IterationInfo,1}, + iteration_time::Float64, + iterations::Array{PHIterationInfo,1}, )::Float64 length(iterations) > 0 ? current_total_time = last(iterations).total_elapsed_time : current_total_time = 0 - return current_total_time + it_time + return current_total_time + iteration_time end -function compute_global_objective(mpi::MpiInfo, s::SpSolution)::Float64 +function compute_global_objective( + mpi::MpiInfo, + s::PhSubProblemSolution, +)::Float64 global_obj = MPI.Allreduce(s.obj, MPI.SUM, mpi.comm) global_obj /= mpi.nprocs return global_obj end -function compute_global_consensus(mpi::MpiInfo, s::SpSolution)::Array{Float64,1} - sp_consensus_vals = s.consensus_vals - global_consensus_vals = MPI.Allreduce(sp_consensus_vals, MPI.SUM, mpi.comm) - global_consensus_vals = global_consensus_vals / mpi.nprocs - return global_consensus_vals +function compute_target(mpi::MpiInfo, s::PhSubProblemSolution)::Array{Float64,1} + sp_vals = s.vals + target = MPI.Allreduce(sp_vals, MPI.SUM, mpi.comm) + target = target / mpi.nprocs + return target end -function compute_global_residual(mpi::MpiInfo, s::SpSolution)::Float64 - n_vars = length(s.consensus_vals) +function compute_global_residual(mpi::MpiInfo, s::PhSubProblemSolution)::Float64 + n_vars = length(s.vals) local_residual_sum = abs.(s.residuals) global_residual_sum = MPI.Allreduce(local_residual_sum, MPI.SUM, mpi.comm) return sum(global_residual_sum) / n_vars end function compute_global_infeasibility( - solution::SpSolution, + solution::PhSubProblemSolution, mpi::MpiInfo, )::Float64 local_infeasibility = norm(solution.residuals) @@ -113,9 +102,13 @@ function compute_global_infeasibility( return global_infeas end -function solve_subproblem(sp::SubProblem, ph_sp_params::SpParams)::SpSolution +function solve_subproblem( + sp::PHSubProblem, + params::PHSubProblemParams, + method::SolutionMethod, +)::PhSubProblemSolution G = length(sp.consensus_vars) - if norm(ph_sp_params.λ) < 1e-3 + if norm(params.λ) < 1e-3 @objective(sp.mip, Min, sp.obj) else @objective( @@ -124,40 +117,31 @@ function solve_subproblem(sp::SubProblem, ph_sp_params::SpParams)::SpSolution sp.obj + sum( sp.weights[g] * - ph_sp_params.λ[g] * - (sp.consensus_vars[g] - ph_sp_params.global_consensus_vals[g]) - for g in 1:G + params.λ[g] * + (sp.consensus_vars[g] - params.target[g]) for g in 1:G ) + - (ph_sp_params.ρ / 2) * sum( - sp.weights[g] * - ( - sp.consensus_vars[g] - - ph_sp_params.global_consensus_vals[g] - )^2 for g in 1:G + (params.ρ / 2) * sum( + sp.weights[g] * (sp.consensus_vars[g] - params.target[g])^2 for + g in 1:G ) ) end - optimize!(sp.mip, XavQiuWanThi2019.Method()) + optimize!(sp.mip, method) obj = objective_value(sp.mip) - sp_consensus_vals = value.(sp.consensus_vars) - return SpSolution( - obj = obj, - consensus_vals = sp_consensus_vals, - residuals = zeros(G), - ) + sp_vals = value.(sp.consensus_vars) + return PhSubProblemSolution(obj = obj, vals = sp_vals, residuals = zeros(G)) end function update_λ_and_residuals!( - solution::SpSolution, - ph_sp_params::SpParams, - global_consensus_vals::Array{Float64,1}, + solution::PhSubProblemSolution, + params::PHSubProblemParams, + target::Array{Float64,1}, )::Nothing - n_vars = length(solution.consensus_vals) - ph_sp_params.global_consensus_vals = global_consensus_vals + n_vars = length(solution.vals) + params.target = target for n in 1:n_vars - solution.residuals[n] = - solution.consensus_vals[n] - ph_sp_params.global_consensus_vals[n] - ph_sp_params.λ[n] += ph_sp_params.ρ * solution.residuals[n] + solution.residuals[n] = solution.vals[n] - params.target[n] + params.λ[n] += params.ρ * solution.residuals[n] end end @@ -179,22 +163,22 @@ end function print_progress( mpi::MpiInfo, - iteration::IterationInfo, + iteration::PHIterationInfo, print_interval, )::Nothing if !mpi.root return end - if iteration.it_num % print_interval != 0 + if iteration.iteration_number % print_interval != 0 return end @info @sprintf( - "Current iteration %8d %20.6e %20.6e %12.2f %% %8.2f %8.2f", - iteration.it_num, + "%8d %20.6e %20.6e %12.2f %% %8.2f %8.2f", + iteration.iteration_number, iteration.global_obj, iteration.global_infeas, iteration.global_residual * 100, - iteration.it_time, + iteration.iteration_time, iteration.total_elapsed_time ) end @@ -209,21 +193,21 @@ end function should_stop( mpi::MpiInfo, - iterations::Array{IterationInfo,1}, - criteria::TerminationCriteria, + iterations::Array{PHIterationInfo,1}, + termination::PHTermination, )::Bool - if length(iterations) >= criteria.max_iterations + if length(iterations) >= termination.max_iterations if mpi.root @info "Iteration limit reached. Stopping." end return true end - if length(iterations) < criteria.min_iterations + if length(iterations) < termination.min_iterations return false end - if last(iterations).total_elapsed_time > criteria.max_time + if last(iterations).total_elapsed_time > termination.max_time if mpi.root @info "Time limit reached. Stopping." end @@ -233,9 +217,9 @@ function should_stop( curr_it = last(iterations) prev_it = iterations[length(iterations)-1] - if curr_it.global_infeas < criteria.min_feasibility + if curr_it.global_infeas < termination.min_feasibility obj_change = abs(prev_it.global_obj - curr_it.global_obj) - if obj_change < criteria.min_improvement + if obj_change < termination.min_improvement if mpi.root @info "Feasibility limit reached. Stopping." end diff --git a/src/solution/methods/ProgressiveHedging/structs.jl b/src/solution/methods/ProgressiveHedging/structs.jl index 39d68a4..46f68f2 100644 --- a/src/solution/methods/ProgressiveHedging/structs.jl +++ b/src/solution/methods/ProgressiveHedging/structs.jl @@ -4,81 +4,34 @@ using JuMP, MPI, TimerOutputs -mutable struct TerminationCriteria - max_iterations::Int - max_time::Float64 - min_feasibility::Float64 - min_improvement::Float64 - min_iterations::Int - - function TerminationCriteria(; - max_iterations::Int = 1000, - max_time::Float64 = 14400.0, - min_feasibility::Float64 = 1e-3, - min_improvement::Float64 = 1e-3, - min_iterations::Int = 2, - ) - return new( - max_iterations, - max_time, - min_feasibility, - min_improvement, - min_iterations, - ) - end +Base.@kwdef mutable struct PHTermination + max_iterations::Int = 1000 + max_time::Float64 = 14400.0 + min_feasibility::Float64 = 1e-3 + min_improvement::Float64 = 1e-3 + min_iterations::Int = 2 end -Base.@kwdef mutable struct IterationInfo - it_num::Int - sp_consensus_vals::Array{Float64,1} - global_consensus_vals::Array{Float64,1} - sp_obj::Float64 +Base.@kwdef mutable struct PHIterationInfo + global_infeas::Float64 global_obj::Float64 - it_time::Float64 - total_elapsed_time::Float64 global_residual::Float64 - global_infeas::Float64 -end - -mutable struct ProgressiveHedging <: SolutionMethod - consensus_vars::Union{Array{VariableRef,1},Nothing} - weights::Union{Array{Float64,1},Nothing} - initial_global_consensus_vals::Union{Array{Float64,1},Nothing} - num_of_threads::Int - ρ::Float64 - λ_default::Float64 - print_interval::Int - termination_criteria::TerminationCriteria - - function ProgressiveHedging(; - consensus_vars::Union{Array{VariableRef,1},Nothing} = nothing, - weights::Union{Array{Float64,1},Nothing} = nothing, - initial_global_consensus_vals::Union{Array{Float64,1},Nothing} = nothing, - num_of_threads::Int = 1, - ρ::Float64 = 1.0, - λ_default::Float64 = 0.0, - print_interval::Int = 1, - termination_criteria::TerminationCriteria = TerminationCriteria(), - ) - return new( - consensus_vars, - weights, - initial_global_consensus_vals, - num_of_threads, - ρ, - λ_default, - print_interval, - termination_criteria, - ) - end + iteration_number::Int + iteration_time::Float64 + sp_vals::Array{Float64,1} + sp_obj::Float64 + target::Array{Float64,1} + total_elapsed_time::Float64 end -struct FinalResult - obj::Float64 - vals::Any - infeasibility::Float64 - total_iteration_num::Int - wallclock_time::Float64 +Base.@kwdef mutable struct ProgressiveHedging <: SolutionMethod + initial_weights::Union{Vector{Float64},Nothing} = nothing + initial_target::Union{Vector{Float64},Nothing} = nothing + ρ::Float64 = 1.0 + λ::Float64 = 0.0 + print_interval::Int = 1 + termination::PHTermination = PHTermination() + inner_method::SolutionMethod = XavQiuWanThi2019.Method() end struct SpResult @@ -86,23 +39,23 @@ struct SpResult vals::Array{Float64,1} end -Base.@kwdef mutable struct SubProblem +Base.@kwdef mutable struct PHSubProblem mip::JuMP.Model obj::AffExpr consensus_vars::Array{VariableRef,1} weights::Array{Float64,1} end -Base.@kwdef struct SpSolution +Base.@kwdef struct PhSubProblemSolution obj::Float64 - consensus_vals::Array{Float64,1} + vals::Array{Float64,1} residuals::Array{Float64,1} end -Base.@kwdef mutable struct SpParams +Base.@kwdef mutable struct PHSubProblemParams ρ::Float64 λ::Array{Float64,1} - global_consensus_vals::Array{Float64,1} + target::Array{Float64,1} end struct MpiInfo @@ -118,9 +71,3 @@ struct MpiInfo return new(comm, rank, is_root, nprocs) end end - -Base.@kwdef struct Callbacks - before_solve_subproblem::Any - after_solve_subproblem::Any - after_iteration::Any -end diff --git a/test/src/solution/methods/ProgressiveHedging/ph.jl b/test/src/solution/methods/ProgressiveHedging/ph.jl index 1e6acd0..0cc1f49 100644 --- a/test/src/solution/methods/ProgressiveHedging/ph.jl +++ b/test/src/solution/methods/ProgressiveHedging/ph.jl @@ -1,4 +1,4 @@ -using Cbc +using HiGHS using MPI using JuMP using UnitCommitment @@ -9,29 +9,32 @@ function fixture(path::String)::String return "$basedir/../../../../fixtures/$path" end -# 1. Initialize MPI +# Initialize MPI MPI.Init() -# 2. Configure progressive hedging method +# Configure progressive hedging method ph = UnitCommitment.ProgressiveHedging() -# 3. Read problem instance +# Read problem instance instance = UnitCommitment.read( [fixture("case14.json.gz"), fixture("case14.json.gz")], ph, ) -# 4. Build JuMP model +# Build JuMP model model = UnitCommitment.build_model( instance = instance, - optimizer = optimizer_with_attributes(Cbc.Optimizer, "LogLevel" => 0), + optimizer = optimizer_with_attributes( + HiGHS.Optimizer, + MOI.Silent() => true, + ), ) -# 5. Run the decentralized optimization algorithm +# Run the decentralized optimization algorithm UnitCommitment.optimize!(model, ph) -# 6. Fetch the solution +# Fetch the solution solution = UnitCommitment.solution(model, ph) -# 7. Close MPI +# Close MPI MPI.Finalize()