CVXPY无法获取向量表达式的元素并将其用作数组索引

Question

问题陈述：我正在尝试建模和优化一个放置问题，即将无向图的节点放置在单元格网格中，以使加权欧几里德长度最小化，但受制于每个网格单元格根据其加权容量只能包含一定数量的节点的约束。我试图使用CVXPY将模型框架化为加权边长度最小化和优化的凸问题。网格组织为2X3矩阵。 当我尝试使用表达式将节点位置转换为网格编号，并尝试使用表达式向量的元素作为数组索引（以便我可以汇总该网格单元中的所有节点权重）时，CVXPY似乎不起作用

这是我试过的一些代码。带有加权边矩阵的10个节点的示例图。第97行（gn=gridNum[i]）和第99行似乎有问题（由问题注释突出显示）。最后的打印行输出应显示所有节点仅位于网格单元格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)

来自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

Answer 1

嗯。我找到了一种让CVXPY比较值然后生成索引的方法

for i in range(num_nodes.value):
    idx = 0
    for j in range(num_grid_cells.value):
        if(gridNum[i].value.__eq__(j)):
            idx = j
            break
    gridW[idx] += nodeWts[i]

虽然这会给出正确的网格单元权重-gridW，但当添加权重约束时，优化器无法最小化（最小化距离为inf），如：

constraints = [locX >= 0.1, locY >= 0.1, locX <= max_X.value, locY <= max_Y.value, gridW <= gridCapacity]

---

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.]]
MINWL inf
X:  None
Y:  None
GN:  None
GW: [130   0   0   0   0   0]