Using an abstract layer for admm algorithms now. Implemented sharing as variant besides general consensus.

Author: Rico Schrage

Project.toml

@@ -8,6 +8,7 @@ version = "0.1.0"
 AutoHashEquals = "15f4f7f2-30c1-5605-9d31-71845cf9641f"
 Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
 JuMP = "4076af6c-e467-56ae-b986-b466b2749572"
+LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Mango = "5e49fdec-d473-4d14-b295-7bff2fcf1925"
 OSQP = "ab2f91bb-94b4-55e3-9ba0-7f65df51de79"
 Revise = "295af30f-e4ad-537b-8983-00126c2a3abe"
@@ -16,6 +17,7 @@ Revise = "295af30f-e4ad-537b-8983-00126c2a3abe"
 AutoHashEquals = "2.2.0"
 Distributions = "0.25.115"
 JuMP = "1.25.0"
+LinearAlgebra = "1.11.0"
 Mango = "^0.4"
 OSQP = "0.8.1"
 Revise = "3.6.4"

README.md

@@ -1,7 +1,7 @@
 
 # Distributed Optimization for Julia
 
-The package DistributedOptimization.jl (DO) aims to provide a collection of distributed algorithms. The algorithms are implemented without considering one special communication technique or package. DO provides abstract types and function interfaces to implement so-called carriers, which are able to execute the distributed algorithms asynchronous. All algorithms can also be used without carrier using fitting @spawn or @async statements.
+The package DistributedOptimization.jl (DO) aims to provide a collection of distributed optimization algorithms. The algorithms are implemented without considering one special communication technique or package. DO provides abstract types and function interfaces to implement so-called carriers, which are able to execute the distributed algorithms asynchronous. All algorithms can also be used without carrier using fitting @spawn or @async statements.
 
 Currently there are two tested algorithms:
 * ADMM multi-value consensus such that the sum of all resp values equals a target vector
@@ -10,3 +10,4 @@ Currently there are two tested algorithms:
 There is one carrier implemented:
 * Mango.jl, agent framework for the simulation of distributed systems, DO provides roles to which the specific algorithms can be assigned to
 
+Note that the package is highly work in progress.

src/DistributedOptimization.jl

@@ -5,7 +5,12 @@ include("carrier/core.jl")
 include("algorithm/core.jl")
 include("algorithm/heuristic/cohda/core.jl")
 include("algorithm/heuristic/cohda/decider.jl")
-include("algorithm/admm/flex_admm.jl")
+
+include("algorithm/admm/core.jl")
+include("algorithm/admm/flex_actor.jl")
+include("algorithm/admm/consensus_admm.jl")
+include("algorithm/admm/sharing_admm.jl")
+
 include("carrier/mango.jl")
 
 end

src/algorithm/admm/consensus_admm.jl

@@ -0,0 +1,47 @@
+export create_consensus_target_reach_admm_coordinator, ADMMConsensusGlobalActor, create_admm_start_with_target
+
+
+@kwdef struct ADMMConsensusGlobalActor <: ADMMGlobalActor
+    α::Int = 100
+end
+
+function z_update(actor::ADMMConsensusGlobalActor, input::Vector{<:Real}, x, u, z, ρ, N)
+    m = length(z)
+    S = zeros(m)
+    for i in 1:N
+        S .+= x[i] .+ u[i]
+    end
+    δ = (input .- S) ./ (N + actor.α/ρ)
+    for i in 1:N
+        z[i] = x[i] .+ u[i] .+ δ
+    end
+    return z
+end
+
+function u_update(actor::ADMMConsensusGlobalActor, x, u, z, ρ, N)
+    return u + x - z
+end
+
+function init_z(actor::ADMMConsensusGlobalActor, n::Int, m::Int)
+    return [ones(m) for _ in 1:n]
+end
+
+function init_u(actor::ADMMConsensusGlobalActor, n::Int, m::Int)
+    return [zeros(m) for _ in 1:n]
+end
+
+function actor_correction(actor::ADMMConsensusGlobalActor, x, z, u, i)
+    return - z[i] + u[i]
+end
+
+function primal_residual(actor::ADMMConsensusGlobalActor, x, z)
+    return maximum(norm.(x .- z))
+end
+
+function create_consensus_target_reach_admm_coordinator()
+    return ADMMGenericCoordinator(global_actor=ADMMConsensusGlobalActor())
+end
+
+function create_admm_start_with_target(target::Vector{<:Real})
+    return ADMMStart(target, length(target))
+end

src/algorithm/admm/core.jl

@@ -0,0 +1,125 @@
+export ADMMStart, ADMMAnswer, ADMMAnswer, ADMMGlobalActor, ADMMGlobalObjective, ADMMGenericCoordinator
+
+struct ADMMStart
+    data::Any
+    solution_length::Int
+end
+
+struct ADMMMessage 
+    v::Vector{Float64}
+    ρ::Float64
+end
+
+struct ADMMAnswer
+    x::Vector{Float64}
+end
+
+abstract type ADMMGlobalActor end
+
+function z_update(actor::ADMMGlobalActor, x, u, z, ρ, N)
+    throw("NotImplemented")
+end
+
+function u_update(actor::ADMMGlobalActor, x, u, z, ρ, N)
+    throw("NotImplemented")
+end
+
+function init_z(actor::ADMMGlobalActor, n::Int, m::Int)
+    throw("NotImplemented")
+end
+
+function init_u(actor::ADMMGlobalActor, n::Int, m::Int)
+    throw("NotImplemented")
+end
+
+function actor_correction(actor::ADMMGlobalActor, x, z, u, i)
+    throw("NotImplemented")
+end
+
+function primal_residual(actor::ADMMGlobalActor, x, z)
+    throw("NotImplemented")
+end
+
+abstract type ADMMGlobalObjective end
+
+function objective(global_objective::ADMMGlobalObjective, x, u, z, N) 
+    throw("NotImplemented")
+end
+
+@kwdef struct ADMMGenericCoordinator <: Coordinator
+    global_actor::ADMMGlobalActor
+    ρ::Float64 = 1.0 
+    max_iters::Int64 = 100
+    slack_penalty::Int64 = 100
+    abs_tol::Float64 = 1e-4
+    rel_tol::Float64 = 1e-3
+    μ::Real = 10
+    τ::Real = 2
+end
+
+# ADMM solver
+function _start_coordinator(admm::ADMMGenericCoordinator, carrier::Carrier, input::Any, m::Int)
+    ρ = admm.ρ
+    max_iters = admm.max_iters
+    abs_tol = admm.abs_tol
+    rel_tol = admm.rel_tol
+    n = length(others(carrier, "coordinator"))
+    μ = admm.μ
+    τ = admm.τ
+
+    # Initialize
+    x = [zeros(m) for i in 1:n]
+    z = init_z(admm.global_actor, n, m)
+    u = init_u(admm.global_actor, n, m)
+
+    for k in 1:max_iters
+        # 1. Local x-updates (in parallel)
+        awaitables = []
+        # send all async and get awaitable
+        for (i,addr) in enumerate(others(carrier, "c"))
+            push!(awaitables, send_awaitable(carrier, ADMMMessage(actor_correction(admm.global_actor, x, z, u, i), ρ), addr))
+        end
+        # await all awaitables and update x
+        for (i,awaitable) in enumerate(awaitables)
+            x[i] = wait(carrier, awaitable).x
+        end
+
+        # 2. Global z-update
+        z_old = deepcopy(z)
+
+        z = z_update(admm.global_actor, input, x, u, z, ρ, n)
+        u = u_update(admm.global_actor, x, u, z, ρ, n)
+
+        # 4. Check convergence
+        # primal residual: max over i of ||x_i - z_i||
+        r_norm = primal_residual(admm.global_actor, x, z)
+        # dual residual: ρ * max over i of ||z_i - z_old_i||
+        s_norm = ρ * maximum(norm.(z .- z_old))
+        # tolerances
+        ϵ_pri = sqrt(m*n)*abs_tol + rel_tol*max(maximum(norm.(x)), maximum(norm.(z)))
+        ϵ_dual = sqrt(m*n)*abs_tol + rel_tol*maximum(norm.(u))
+        if r_norm < ϵ_pri && s_norm < ϵ_dual
+            @warn "Converged in $k iterations."
+            break
+        end
+        
+        # Varying penalty paramter according to B. S. He, H. Yang, and S. L. Wang, “Alternating direction method with self
+        # adaptive penalty parameters for monotone variational inequalities,”
+        if r_norm > μ * s_norm
+            ρ = ρ * τ
+        elseif s_norm > μ * r_norm
+            ρ = ρ / τ
+        end
+        
+        if k == max_iters
+            @warn "Reached max iterations ($max_iters) without full convergence."
+            throw("ADMM not converged $x, $u")
+        end
+    end
+    return x, z, u
+end
+
+function start_optimization(coordinator::ADMMGenericCoordinator, carrier::Carrier, start_data::ADMMStart, meta::Any)
+    x,_,_ = _start_coordinator(coordinator, carrier, start_data.data, start_data.solution_length)
+    return x
+end

src/algorithm/admm/flex_actor.jl

@@ -0,0 +1,98 @@
+export ADMMFlexActor, create_admm_flex_actor_one_to_many, result
+
+using JuMP
+using OSQP
+using LinearAlgebra
+
+mutable struct ADMMFlexActor <: DistributedAlgorithm
+    l::Vector{Float64} # lower bounds
+    u::Vector{Float64} # upper bounds
+    C::Matrix{Float64} # coupling matrix
+    d::Vector{Float64} # coupling RHS
+    S::Vector{Float64} # prio penalty per sector
+    x::Vector{Float64} # intermediate result/result
+    ADMMFlexActor(l::Vector{Float64},
+        u::Vector{Float64},
+        C::Matrix{Float64},
+        d::Vector{Float64},
+        S::Vector{Float64}) = new(l, u, C, d, S, Vector{Float64}())
+end
+
+function _create_C_and_d(u::Vector{<:Real})
+    m = length(u)
+    R = 1 + 2*(m-1)
+    C = zeros(eltype(u), R, m)
+    d = zeros(eltype(u), R)
+
+    for i in 1:length(u)
+        C[1, i] = u[i] < 0 ? -1 : 1
+    end
+    d[1] = sum(abs.(u))
+
+    for j in 1:(m-1)
+        
+        r1 = 1 + 2*(j-1) + 1
+        r2 = r1 + 1
+
+        C[r1, j] = u[j] == 0 || u[m] == 0 ? 0 :  1/u[j]
+        C[r1, m] = u[j] == 0 || u[m] == 0 ? 0 : -1/u[m]
+
+        C[r2, j] = u[j] == 0 || u[m] == 0 ? 0 : -1/u[j]
+        C[r2, m] = u[j] == 0 || u[m] == 0 ? 0 :  1/u[m]
+    end
+
+    return C, d
+end
+
+function create_admm_flex_actor_one_to_many(in_capacity::Real, η::Vector{Float64}, S::Union{Nothing,Vector{<:Real}}=nothing)
+    tech_capacity = in_capacity .* η
+
+    if isnothing(S)
+        S = zeros(length(η))
+    end
+    l = min.(zeros(length(tech_capacity)), tech_capacity)
+    u = max.(tech_capacity, zeros(length(tech_capacity)))
+    C, d = _create_C_and_d(tech_capacity)
+
+    return ADMMFlexActor(l, u, C, d, S)
+end
+
+function result(actor::ADMMFlexActor)
+    return actor.x
+end
+
+# Solve the projection / local update via QP: minimize 1/2||x - v||^2 s.t. l <= x <= u, Cx <= d
+function _local_update(actor::ADMMFlexActor, v::Vector{Float64}, ρ::Float64)
+    m = length(v)
+    model = Model(OSQP.Optimizer)
+    set_silent(model)
+    
+    @variable(model, x[1:m])
+    
+    # admm objective: (ρ/2)||x - v||^2 + S_i·x
+    # Note that v = - (correction) in the implementation as a result the equation is based on  (ρ/2)||x + v||^2 + S_i·x
+    # transform to standard form Hx^2 - hx
+    # in some cases about 100-1000x faster
+    H = Diagonal(ones(m) * ρ)
+    h = ρ*v
+
+    # + priority cost term: S_i·x
+    @objective(model, Min, 0.5 * sum(H[i,i]*x[i]^2 for i=1:m) + sum(h[i]*x[i] for i=1:m) + sum(actor.S[i] * x[i] for i in 1:m))
+    #@objective(model, Min, (ρ/2)*sum((x[i] + v[i])^2 + x[i]*actor.S[i] for i in 1:m))
+    
+    # box constraints
+    @constraint(model, [i=1:m], actor.l[i] <= x[i] <= actor.u[i])
+    
+    # coupling constraints
+    @constraint(model, actor.C * x .<= actor.d)
+
+    optimize!(model)
+
+    return value.(x)
+end
+
+function on_exchange_message(actor::ADMMFlexActor, carrier::Carrier, message_data::ADMMMessage, meta::Any)
+    actor.x = _local_update(actor, message_data.v, message_data.ρ)
+
+    reply_to_other(carrier, ADMMAnswer(actor.x), meta)
+end