Distributed Simplex

This is an example on how to use the DistributedSimplex class. See also the reference [BuNo11].

examples/setups/distributed_simplex/launcher.py

import dill as pickle
import numpy as np
from mpi4py import MPI
from disropt.agents import Agent
from disropt.algorithms.constraintexchange import DistributedSimplex
from disropt.functions import Variable
from disropt.utils.graph_constructor import ring_graph
from disropt.utils.LP_utils import generate_LP
from disropt.problems import LinearProblem

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

# Generate a ring graph (for which the diameter is nproc-1)
Adj = ring_graph(nproc)
graph_diam = nproc-1

# reset local seed
np.random.seed(1)

# number of constraints
n_constr = 3

# number of columns for each processor
k = 2

# generate a feasible optimization problem of size k * nproc
c_glob, A_glob, b_glob, x_glob = generate_LP(k * nproc, n_constr, 50, constr_form='eq')

# extract the columns assigned to this agent
local_indices = list(np.arange(k*local_rank, k*(local_rank+1)))

c_loc = c_glob[local_indices, :]
A_loc = A_glob[:, local_indices]
b_loc = b_glob

# define the local problem data
x = Variable(k)
obj = c_loc @ x
constr = A_loc.transpose() @ x == b_loc
problem = LinearProblem(objective_function=obj, constraints=constr)

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

# instantiate the algorithm
algorithm = DistributedSimplex(agent, enable_log=True, problem_size=nproc*k,
    local_indices=local_indices, stop_iterations=2*graph_diam+1)

# run the algorithm
x_sequence, J_sequence = algorithm.run(iterations=100, verbose=True)

# print results
_, _, _, J_final = algorithm.get_result()
print("Agent {} - {} iterations - final cost {}".format(agent.id, len(J_sequence), J_final))

# save results to file
if agent.id == 0:
    with open('info.pkl', 'wb') as output:
        pickle.dump({'N': nproc, 'n_vars': k * nproc, 'n_constr': n_constr, 'c': c_glob,
            'A': A_glob, 'b': b_glob, 'opt_sol': x_glob}, output, pickle.HIGHEST_PROTOCOL)

np.save("agent_{}_x_seq.npy".format(agent.id), x_sequence)
np.save("agent_{}_J_seq.npy".format(agent.id), J_sequence)

examples/setups/distributed_simplex/results.py

import numpy as np
import matplotlib.pyplot as plt
import dill as pickle
from disropt.functions import Variable
from disropt.problems import LinearProblem

# initialize
with open('info.pkl', 'rb') as inp:
    info = pickle.load(inp)
NN        = info['N']
c         = info['c']
opt_sol   = info['opt_sol']

# load agent data
sequence_J = {}
for i in range(NN):
    sequence_J[i] = np.load("agent_{}_J_seq.npy".format(i))

# compute optimal cost
opt_cost = (c.transpose() @ opt_sol).flatten()

# plot cost evolution
plt.figure()
plt.title("Cost evolution")
colors = np.random.rand(NN+1, 3)
max_iters = 0

for i in range(NN):
    seq_J_i = sequence_J[i]
    n_iters_i = len(seq_J_i)
    max_iters = max(max_iters, n_iters_i)
    plt.plot(np.arange(n_iters_i), seq_J_i, color=colors[i])

# plot optimal cost
plt.plot(np.arange(max_iters), np.ones(max_iters)*opt_cost, '--', color=colors[NN])
plt.show()

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

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

Distributed Simplex (dual problem)

This is an example on how to use the DualDistributedSimplex class, which forms the dual problem of the given optimization problem and solves it with the Distributed Simplex algorithm.

examples/setups/distributed_simplex/launcher_dual.py

import dill as pickle
import numpy as np
from mpi4py import MPI
from disropt.agents import Agent
from disropt.algorithms.constraintexchange import DualDistributedSimplex
from disropt.functions import Variable
from disropt.utils.graph_constructor import ring_graph
from disropt.utils.LP_utils import generate_LP
from disropt.problems import LinearProblem

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

# Generate a ring graph (for which the diameter is nproc-1)
Adj = ring_graph(nproc)
graph_diam = nproc-1

# reset local seed
np.random.seed(1)

# number of variables
n_vars = 3

# number of constraints for each processor
k = 2

# generate a feasible optimization problem with k * nproc constraints
c_glob, A_glob, b_glob, x_glob = generate_LP(n_vars, k * nproc, 50, direction='max')

# extract the constraints assigned to this agent
local_indices = list(np.arange(k*local_rank, k*(local_rank+1)))

c_loc = c_glob
A_loc = A_glob[local_indices, :]
b_loc = b_glob[local_indices, :]

# define the local problem data
x = Variable(n_vars)
obj = -c_loc @ x # minus sign for maximization
constr = A_loc.transpose() @ x <= b_loc
problem = LinearProblem(objective_function=obj, constraints=constr)

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

# instantiate the algorithm
algorithm = DualDistributedSimplex(agent, enable_log=True, num_constraints=nproc*k,
    local_indices=local_indices, stop_iterations=2*graph_diam+1)

# run the algorithm
x_sequence, J_sequence = algorithm.run(iterations=100, verbose=True)

# print results
_, _, _, J_final = algorithm.get_result()
print("Agent {} - {} iterations - final cost {}".format(agent.id, len(J_sequence), J_final))

# save results to file
if agent.id == 0:
    with open('info.pkl', 'wb') as output:
        pickle.dump({'N': nproc, 'n_constr': k * nproc, 'n_vars': n_vars, 'c': c_glob,
            'A': A_glob, 'b': b_glob, 'opt_sol': x_glob}, output, pickle.HIGHEST_PROTOCOL)

np.save("agent_{}_x_seq.npy".format(agent.id), x_sequence)
np.save("agent_{}_J_seq.npy".format(agent.id), J_sequence)

examples/setups/distributed_simplex/results_dual.py

import numpy as np
import matplotlib.pyplot as plt
import dill as pickle
from disropt.functions import Variable
from disropt.problems import LinearProblem

# initialize
with open('info.pkl', 'rb') as inp:
    info = pickle.load(inp)
NN        = info['N']
c         = info['c']
opt_sol   = info['opt_sol']

# load agent data
sequence_J = {}
for i in range(NN):
    sequence_J[i] = np.load("agent_{}_J_seq.npy".format(i))

# compute optimal cost
opt_cost = (c.transpose() @ opt_sol).flatten()

# plot cost evolution
plt.figure()
plt.title("Cost evolution")
colors = np.random.rand(NN+1, 3)
max_iters = 0

for i in range(NN):
    seq_J_i = sequence_J[i]
    n_iters_i = len(seq_J_i)
    max_iters = max(max_iters, n_iters_i)
    plt.plot(np.arange(n_iters_i), seq_J_i, color=colors[i])

# plot optimal cost
plt.plot(np.arange(max_iters), np.ones(max_iters)*opt_cost, '--', color=colors[NN])
plt.show()

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

> mpirun -np 10 --oversubscribe python launcher_dual.py
> python results_dual.py