Constraints Consensus

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

examples/setups/constraints_consensus/launcher.py

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

# get MPI info
NN = MPI.COMM_WORLD.Get_size()
agent_id = MPI.COMM_WORLD.Get_rank()

# Generate a common graph (everyone uses the same seed)
Adj = binomial_random_graph(NN, p=0.03, seed=1)

#####################
# Problem parameters
#####################
np.random.seed(10*agent_id)

# linear objective function
dim = 2
z = Variable(dim)
c = np.ones([dim,1])
obj_func = c @ z

# constraints are circles of the form (z-p)^\top (z-p) <= 1
# equivalently z^\top z - 2(A p)^\top z + p^\top A p <= 1
I = np.eye(dim)
p = np.random.rand(dim,1)
r = 1 # unitary radius

constr = []
ff = QuadraticForm(z,I,- 2*(I @ p),(p.transpose() @ I @ p) - r**2)
constr.append(ff<= 0)

#####################
# Distributed algorithms
#####################

# local agent and problem
agent = Agent(
    in_neighbors=np.nonzero(Adj[agent_id, :])[0].tolist(),
    out_neighbors=np.nonzero(Adj[:, agent_id])[0].tolist())
pb = Problem(obj_func, constr)
agent.set_problem(pb)

# instantiate the algorithm
algorithm = ConstraintsConsensus(agent=agent,
                                 enable_log=True)

# run the algorithm
n_iter = NN*3
x_sequence = algorithm.run(iterations=n_iter, verbose=True)

# print results
x_final = algorithm.get_result()
print("Agent {}: {}".format(agent_id, x_final.flatten()))

# save results to file
if agent_id == 0:
    with open('info.pkl', 'wb') as output:
        pickle.dump({'N': NN, 'size': dim, 'iterations': n_iter}, output, pickle.HIGHEST_PROTOCOL)
    with open('objective_function.pkl', 'wb') as output:
        pickle.dump(obj_func, output, pickle.HIGHEST_PROTOCOL)

with open('agent_{}_constr.pkl'.format(agent_id), 'wb') as output:
    pickle.dump(constr, output, pickle.HIGHEST_PROTOCOL)
np.save("agent_{}_seq.npy".format(agent_id), x_sequence)

examples/setups/constraints_consensus/results.py

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

# initialize
with open('info.pkl', 'rb') as inp:
    info = pickle.load(inp)
NN = info['N']
iters = info['iterations']
size = info['size']
# load agent data
sequence = np.zeros((NN, iters, size))
local_constr = {}
for i in range(NN):
    sequence[i, :, :] = np.load("agent_{}_seq.npy".format(i), allow_pickle=True).reshape((iters, size))
    with open('agent_{}_constr.pkl'.format(i), 'rb') as inp:
        local_constr[i] = pickle.load(inp)
with open('objective_function.pkl', 'rb') as inp:
    obj_func = pickle.load(inp)

# solve centralized problem
global_constr = []
for i in range(NN):
    global_constr.extend(local_constr[i])
global_pb = Problem(obj_func, global_constr)
x_centr = global_pb.solve()
cost_centr = obj_func.eval(x_centr)

# compute cost errors
cost_err = np.zeros((NN, iters))

for i in range(NN):
    for t in range(iters):
        cost_err[i, t] = abs(obj_func.eval(sequence[i, t, :].reshape((size, 1))) - cost_centr)

# compute max violation
vio_err = np.zeros((NN, iters))
for i in range(NN):
    for t in range(iters):
        xt = sequence[i, t, :].reshape((size, 1))
        max_err = np.zeros((len(global_constr), 1))
        for c in range(len(global_constr)):
            max_err[c] = global_constr[c].function.eval(xt)
        vio_err[i, t] = np.max(max_err)

# Plot the evolution of the local estimates
# generate N colors
colors = {}
for i in range(NN):
    colors[i] = np.random.rand(3, 1).flatten()

# plot cost error
plt.figure()
plt.title('Evolution of cost error')
plt.xlabel(r"iteration $k$")
plt.ylabel(r"$|f(x_i^k) - f^\star|$")

for i in range(NN):
    plt.plot(np.arange(iters), cost_err[i, :], color=colors[i])

# plot violation error
plt.figure()
plt.title('Maximum constraint violation')
plt.xlabel(r"iteration $k$")
plt.ylabel(r"$max_j g(x_i^k)$")

for i in range(NN):
    plt.plot(np.arange(iters), vio_err[i, :], color=colors[i])

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