Problem Statement: I am trying to model and optimize a placement problem of placing the nodes of an undirected graph in a grid of cells such that the weighted Euclidean length is minimized subject to the constraint each grid cell can contain only a certain number of nodes based on its weighted-capacity. I am trying to use CVXPY to frame the model as a convex problem of weighted edge length minimization and optimize. The grid is organized as a 2X3 matrix. When I try to convert the node locations to a grid number with an expression and try to use the elements of the expression vector as an array index ( so that I can sum up all node weights in that grid cell), CVXPY does not seem to work.
Here is some code I tried. Sample graph of 10 nodes with weighted edge matrix. Line 97 (gn = gridNum[i]) and 99 seem to have an issue( highlighted by PROBLEM comment. The last print line output is expected to show all the nodes are in only grid cell 0.
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
@author: rajkumar
"""
import numpy as np
import cvxpy as cvxp
import matplotlib.pyplot as plt
# Problem data creation
# Given a set of nodes and weighted edges,
# place them on a grid of cells which have a link capacity,
# such that the total euclidean distance is minimized
# subject to the constraint that number of nodes in a grid cell
# dont exceed a certain capacity
# Will deal with maximizing grid cell link capacity objective function later
# once I get the grid cell number issue fixed
# E.g graph - 10 nodes with connectivity matrix
# 6 grid cells arranged as 2x3
num_nodes = cvxp.Parameter(nonneg=True,value=10)
num_grid_cells = cvxp.Parameter(nonneg=True,value=6)
# Max X and Y grid coordinate values
max_X = cvxp.Parameter(nonneg=True, value=(num_grid_cells.value/2))
max_Y = cvxp.Parameter(nonneg=True,value=(num_grid_cells.value/3))
#Weight of each node
nodeWts = np.array([10,20,15,12,19,11,14,9,12,8])
#Capacity of each grid cell
gridCellCapacities = np.array([30,30,30,30,30,30])
# created a Variable intialized with above
gridCapacity = cvxp.Parameter(shape=num_grid_cells.value, value=gridCellCapacities, nonneg=True)
# Node adjacency matrix
cellConnectivity = np.matrix([[0.,5.,2.,3.,4.,5.,0.,0.,0.,0.],
[5.,0.,0.,2.,3.,4.,5.,6.,7.,8.],
[2.,0.,0.,2.,3.,4.,5.,6.,7.,8.],
[3.,2.,2.,0.,2.,3.,4.,5.,6.,7.],
[4.,3.,3.,2.,0.,2.,3.,4.,5.,6.],
[5.,4.,4.,3.,2.,0.,2.,3.,4.,5.],
[0.,5.,5.,4.,3.,2.,0.,2.,3.,4.],
[0.,6.,6.,5.,4.,3.,2.,0.,2.,3.],
[0.,7.,7.,6.,5.,4.,3.,2.,0.,8.],
[0.,8.,8.,7.,6.,5.,4.,3.,8.,0.]])
cellWeightedDeg = np.matrix([[19,0,0,0,0,0,0,0,0,0],
[0,40,0,0,0,0,0,0,0,0],
[0,0,37,0,0,0,0,0,0,0],
[0,0,0,34,0,0,0,0,0,0],
[0,0,0,0,32,0,0,0,0,0],
[0,0,0,0,0,32,0,0,0,0],
[0,0,0,0,0,0,28,0,0,0],
[0,0,0,0,0,0,0,31,0,0],
[0,0,0,0,0,0,0,0,42,0],
[0,0,0,0,0,0,0,0,0,49]])
# Positive semi definite laplacian
cellLP = cellWeightedDeg - cellConnectivity
locX = cvxp.Variable(num_nodes.value,pos=True)
locY = cvxp.Variable(num_nodes.value,pos=True)
# Variable to store sum of node weights in a grid cell
gridWt = cvxp.Variable(num_grid_cells.value, pos=True)
#Number of nodes in each grid cell - num_nodes x num_grid_cells matrix
nodesInGridCells = cvxp.Variable(shape=(num_nodes.value,num_grid_cells.value))
#### Some workaround if not positive semi definite
cellLP = 0.5*(cellLP+cellLP.T) # make Q symmetric
w, v = np.linalg.eig(cellLP) # eigen decomposition
print("Eigenvalues\n",w)
w1 = min(w) # first eigen value
print("Smallest eigenvalue",w1)
tol = 1.0e-10 # tolerance
f = 0 # factor for diagonal perturbation
if w1<tol:
f = -w1 + tol
cellLP += f*np.eye(10)
print('cellLP - sh', cellLP)
# Get grid number from X,Y locations of nodes
# Grid arrangement
# 3,4,5
# 0,1,2
gridW = np.array([0,0,0,0,0,0])
gridNum = cvxp.ceil(locX) + cvxp.abs(cvxp.square(cvxp.ceil(locY))) - 2
# The following loop is to get the gridNum from the above expression vector
# so that it can be used to index into gridW - weight array
# Elementwise extraction from gridNum expression is not working
# PROBLEM IN THIS LOOP
for i in range(num_nodes.value):
gn = gridNum[i]
print('ISGV', gn.is_vector())
gridW[gn.value] += nodeWts[i]
print("GN_IS_VEC: ", gridNum.is_vector())
constraints = [locX >= 0.1, locY >= 0.1, locX <= max_X.value, locY <= max_Y.value]
objectiveX = (1/2)*cvxp.quad_form(locX,cellLP)
objectiveY = (1/2)*cvxp.quad_form(locY,cellLP)
# PSD based objective
#objectiveX = (1/2)*cvxp.quad_form(locX, cvxp.Parameter(shape=cellLP.shape, value=cellLP, PSD=True))
#objectiveY = (1/2)*cvxp.quad_form(locY, cvxp.Parameter(shape=cellLP.shape, value=cellLP, PSD=True))
prob = cvxp.Problem(cvxp.Minimize(objectiveX+objectiveY), constraints)
prob.solve();
print('MINWL',prob.value)
print('X: ', locX.value)
print('Y: ', locY.value)
print('GN: ', gridNum.value)
# Wrong values printed - Only grid number 0 should have all the nodes
print('GW: ', gridW)
OUTPUT FROM CVXPY
Eigenvalues
[-7.10542736e-15 5.82050153e+01 1.91148884e+01 5.18253316e+01
3.09023706e+01 4.05651162e+01 3.85097439e+01 3.38040401e+01
3.48251251e+01 3.62483688e+01]
Smallest eigenvalue -7.105427357601002e-15
cellLP - sh [[19. -5. -2. -3. -4. -5. 0. 0. 0. 0.]
[-5. 40. 0. -2. -3. -4. -5. -6. -7. -8.]
[-2. 0. 37. -2. -3. -4. -5. -6. -7. -8.]
[-3. -2. -2. 34. -2. -3. -4. -5. -6. -7.]
[-4. -3. -3. -2. 32. -2. -3. -4. -5. -6.]
[-5. -4. -4. -3. -2. 32. -2. -3. -4. -5.]
[ 0. -5. -5. -4. -3. -2. 28. -2. -3. -4.]
[ 0. -6. -6. -5. -4. -3. -2. 31. -2. -3.]
[ 0. -7. -7. -6. -5. -4. -3. -2. 42. -8.]
[ 0. -8. -8. -7. -6. -5. -4. -3. -8. 49.]]
ISGV True
ISGV True
ISGV True
ISGV True
ISGV True
ISGV True
ISGV True
ISGV True
ISGV True
ISGV True
GN_IS_VEC: True
MINWL 2.872635462836115e-11
X: [0.16948144 0.16948144 0.16948144 0.16948144 0.16948144 0.16948144
0.16948144 0.16948144 0.16948144 0.16948144]
Y: [0.16948144 0.16948144 0.16948144 0.16948144 0.16948144 0.16948144
0.16948144 0.16948144 0.16948144 0.16948144]
GN: [0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]
GW: [130 130 130 130 130 130] <---- Problem here
