Distributed ADMM

This is an example on how to use the ADMM class.

examples/setups/distributed_ADMM/launcher.py

import dill as pickle
import numpy as np
from mpi4py import MPI
from disropt.agents import Agent
from disropt.algorithms.admm import ADMM
from disropt.functions import QuadraticForm, Variable
from disropt.utils.graph_constructor import binomial_random_graph
from disropt.problems import Problem

# get MPI info
comm = MPI.COMM_WORLD
nproc = comm.Get_size()
local_rank = comm.Get_rank()

# Generate a common graph (everyone use the same seed)
Adj = binomial_random_graph(nproc, p=0.3, seed=1)

# reset local seed
np.random.seed(10*local_rank)

agent = Agent(
    in_neighbors=np.nonzero(Adj[local_rank, :])[0].tolist(),
    out_neighbors=np.nonzero(Adj[:, local_rank])[0].tolist())

# variable dimension
n = 2

# generate a positive definite matrix
P = np.random.randn(n, n)
P = P.transpose() @ P
bias = 3*np.random.randn(n, 1)

# declare a variable
x = Variable(n)

# define the local objective function
fn = QuadraticForm(x - bias, P)

# define a (common) constraint set
constr = [x <= 1, x >= -1]

# local problem
pb = Problem(fn, constr)
agent.set_problem(pb)

# instantiate the algorithms
initial_z = np.ones((n, 1))
# initial_lambda = {local_rank: 10*np.random.rand(n, 1)}
initial_lambda = {local_rank: np.ones((n, 1))}

for j in agent.in_neighbors:
    # initial_lambda[j] = 10*np.random.rand(n, 1)
    initial_lambda[j] = np.ones((n, 1))

algorithm = ADMM(agent=agent,
                 initial_lambda=initial_lambda,
                 initial_z=initial_z,
                 enable_log=True)

# run the algorithm
x_sequence, lambda_sequence, z_sequence = algorithm.run(iterations=100, penalty=0.1, verbose=True)
x_t, lambda_t, z_t = algorithm.get_result()
print("Agent {}: primal {} dual {} auxiliary primal {}".format(agent.id, x_t.flatten(), lambda_t, z_t.flatten()))

np.save("agents.npy", nproc)

# save agent and sequence
if local_rank == 0:
    with open('constraints.pkl', 'wb') as output:
        pickle.dump(constr, output, pickle.HIGHEST_PROTOCOL)
with open('agent_{}_function.pkl'.format(agent.id), 'wb') as output:
    pickle.dump(agent.problem.objective_function, output, pickle.HIGHEST_PROTOCOL)
with open('agent_{}_dual_sequence.pkl'.format(agent.id), 'wb') as output:
    pickle.dump(lambda_sequence, output, pickle.HIGHEST_PROTOCOL)
np.save("agent_{}_primal_sequence.npy".format(agent.id), x_sequence)
np.save("agent_{}_auxiliary_primal_sequence.npy".format(agent.id), z_sequence)

examples/setups/distributed_ADMM/results.py

import dill as pickle
import numpy as np
import matplotlib.pyplot as plt
from disropt.problems import Problem

N = np.load("agents.npy")
n = 2

lambda_sequence = {}
x_sequence = {}
z_sequence = {}
local_obj_function = {}
for i in range(N):
    with open('agent_{}_dual_sequence.pkl'.format(i), 'rb') as input:
        lambda_sequence[i] = pickle.load(input)
    x_sequence[i] = np.load("agent_{}_primal_sequence.npy".format(i))
    z_sequence[i] = np.load("agent_{}_auxiliary_primal_sequence.npy".format(i))
    with open('agent_{}_function.pkl'.format(i), 'rb') as input:
        local_obj_function[i] = pickle.load(input)
with open('constraints.pkl', 'rb') as input:
    constr = pickle.load(input)
iters = x_sequence[0].shape[0]

# solve centralized problem
global_obj_func = 0
for i in range(N):
    global_obj_func += local_obj_function[i]

global_pb = Problem(global_obj_func, constr)
x_centr = global_pb.solve()
cost_centr = global_obj_func.eval(x_centr)
x_centr = x_centr.flatten()

# compute cost errors
cost_err = np.zeros((N, iters)) - cost_centr

for i in range(N):
    for t in range(iters):
        # add i-th cost
        cost_err[i, t] += local_obj_function[i].eval(x_sequence[i][t, :])

# plot maximum cost error
plt.figure()
plt.title('Maximum cost error (among agents)')
plt.xlabel(r"iteration $k$")
plt.ylabel(r"$\max_{i} \: \left|\sum_{j=1}^N f_j(x_i^k) - f^\star \right|$")
plt.semilogy(np.arange(iters), np.amax(np.abs(cost_err), axis=0))

# plot maximum solution error
sol_err = np.zeros((N, iters))
for i in range(N):
    sol_err[i] = np.linalg.norm(np.squeeze(x_sequence[i]) - x_centr[None, :], axis=1)

plt.figure()
plt.title('Maximum solution error (among agents)')
plt.xlabel(r"iteration $k$")
plt.ylabel(r"$\max_{i} \: \|x_i^k - x^\star\|$")
plt.semilogy(np.arange(iters), np.amax(sol_err, axis=0))

plt.show()

The two files can be executed by issuing the following commands in the example folder:

> mpirun -np 30 --oversubscribe python launcher.py
> python results.py