根据给定的约束创建矩阵

Question

我正在尝试创建一个具有以下约束的矩阵。

1. 列总和应介于 300 和 390 之间，包括两个值。
2. 行总和应等于用户为每行指定的值。
3. 矩阵中的非零值不应小于 10。
4. 给定列中非零值的计数不应超过 4。
5. 列应按对角线顺序排列。

如果UserInput = [427.7, 12.2, 352.7, 58.3, 22.7, 31.9, 396.4, 29.4, 171.5, 474.5, 27.9, 200]

我想要这样的输出矩阵，

编辑 1

我使用 Pyomo 尝试了以下方法，但是，我遇到了第 5 个约束，即 列值应在矩阵中对角对齐

import sys
import math
import numpy as np
import pandas as pd

from pyomo.environ import *

solverpath_exe= 'glpk-4.65\\w64\\glpsol.exe'
solver=SolverFactory('glpk',executable=solverpath_exe)

# Minimize the following:
# Remaining pieces to be zero for all et values
# The number of cells containg non-zero values

# Constraints
# 1) Column sum, CS, is: 300 <= CS <= 390
# 2) Row sum, RS, is equal to user-specified values, which are present in the E&T ticket column of the file
# 3) Number of non-zero values, NZV, in each column, should be: 0 < NZV <= 4
# 4) The NZV in the matrix should be: NZV >= 10
# 5) The pieces are stacked on top of each other. So, a the cell under a non-zero value cell is zero, than all cells underneath should have zeros.

maxlen = 390
minlen = 300
npiece = 4
piecelen = 10

# Input data: E&T Ticket values
etinput = [427.7, 12.2, 352.7, 58.3, 22.7, 31.9,
           396.4, 29.4, 171.5, 474.5, 27.9, 200]


# Create data structures to store values
etnames  = [f'et{i}' for i in range(1,len(etinput) + 1)]
colnames = [f'col{i}' for i in range(1, math.ceil(sum(etinput)/minlen))] #+1 as needed

et_val = dict(zip(etnames, etinput))

# Instantiate Concrete Model
model2 = ConcreteModel()

# define variables and set upper bound to 390 
model2.vals = Var(etnames, colnames, domain=NonNegativeReals,bounds = (0, maxlen), initialize=0)

# Create Boolean variables
bigM = 10000
model2.y = Var(colnames, domain= Boolean)
model2.z = Var(etnames, colnames, domain= Boolean)


# Minimizing the sum of difference between the E&T Ticket values and rows 
model2.minimizer = Objective(expr= sum(et_val[r] - model2.vals[r, c]
                                      for r in etnames for c in colnames),
                             sense=minimize)

model2.reelconstraint = ConstraintList()
for c in colnames:
    model2.reelconstraint.add(sum(model2.vals[r,c] for r in etnames) <= bigM * model2.y[c])
    

# Set constraints for row sum equal to ET values
model2.rowconstraint = ConstraintList()
for r in etnames:
    model2.rowconstraint.add(sum(model2.vals[r, c] for c in colnames) <= et_val[r])

    
# Set contraints for upper bound of column sums
model2.colconstraint_upper = ConstraintList()
for c in colnames:
    model2.colconstraint_upper.add(sum(model2.vals[r, c] for r in etnames) <= maxlen)
    

# Set contraints for lower bound of column sums
model2.colconstraint_lower = ConstraintList()
for c in colnames:
    model2.colconstraint_lower.add(sum(model2.vals[r, c] for r in etnames) + bigM * (1-model2.y[c]) >= minlen)
    

model2.bool = ConstraintList()
for c in colnames:
    for r in etnames:
        model2.bool.add(model2.vals[r,c] <= bigM * model2.z[r,c])
    

model2.npienceconstraint = ConstraintList()
for c in colnames:
    model2.npienceconstraint.add(sum(model2.z[r, c] for r in etnames) <= npiece)

# Call solver for model
solver.solve(model2);

# Create dataframe of output
pdtest = pd.DataFrame([[model2.vals[r, c].value for c in colnames] for r in etnames],
                        index=etnames,
                        columns=colnames)

pdtest

输出

Answer 1

我认为您将其设置为 LP 是正确的。它可以表述为 MIP。

我在这里没有修改任何种类的输入，我不确定在你有限制的情况下，所有输入都能保证可行的结果。

我惩罚非对角线选择以鼓励对角线选择，并设置一些“选择完整性”约束来强制块选择。

在大约 1/10 秒内解决...

# magic matrix

# Constraints
# 1) Column sum, CS, is: 300 <= CS <= 390
# 2) Row sum, RS, is equal to user-specified values, which are present in the E&T ticket column of the file
# 3) Number of non-zero values, NZV, in each column, should be: 0 < NZV <= 4
# 4) The NZV in the matrix should be: NZV >= 10
# 5) The pieces are stacked on top of each other. So, a the cell under a non-zero value cell is zero, than all cells underneath should have zeros.

import pyomo.environ as pyo

# user input
row_tots = [427.7, 12.2, 352.7, 58.3, 22.7, 31.9, 396.4, 29.4, 171.5, 474.5, 27.9, 200]
min_col_sum = 300
max_col_sum = 390
max_non_zero = 4
min_size = 10
bigM = max(row_tots)

m = pyo.ConcreteModel()

# SETS
m.I = pyo.Set(initialize=range(len(row_tots)))
m.I_not_first = pyo.Set(within=m.I, initialize=range(1, len(row_tots)))
m.J = pyo.Set(initialize=range(int(sum(row_tots)/min_col_sum)))

# PARAMS
m.row_tots = pyo.Param(m.I, initialize={k:v for k,v in enumerate(row_tots)})

# set up weights (penalties) based on distance from diagonal line
# between corners using indices as points and using distance-to-line formula
weights = { (i, j) : abs((len(m.I)-1)/(len(m.J)-1)*j - i) for i in m.I for j in m.J}
m.weight  = pyo.Param(m.I * m.J, initialize=weights)

# VARS
m.X = pyo.Var(m.I, m.J, domain=pyo.NonNegativeReals)
m.Y = pyo.Var(m.I, m.J, domain=pyo.Binary)          # selection indicator
m.UT = pyo.Var(m.I, m.J, domain=pyo.Binary)         # upper triangle of non-selects

# C1: col min sum
def col_sum_min(m, j):
    return sum(m.X[i, j] for i in m.I) >= min_col_sum
m.C1 = pyo.Constraint(m.J, rule=col_sum_min)

# C2: col max sum
def col_sum_max(m, j):
    return sum(m.X[i, j] for i in m.I) <= max_col_sum
m.C2 = pyo.Constraint(m.J, rule=col_sum_max)

# C3: row sum 
def row_sum(m, i):
    return sum(m.X[i, j] for j in m.J) == m.row_tots[i]
m.C3 = pyo.Constraint(m.I, rule=row_sum)

# C4: max nonzeros
def max_nz(m, j):
    return sum(m.Y[i, j] for i in m.I) <= max_non_zero
m.C4 = pyo.Constraint(m.J, rule=max_nz)


# selection variable enforcement
def selection_low(m, i, j):
    return min_size*m.Y[i, j] <= m.X[i, j]
m.C10 = pyo.Constraint(m.I, m.J, rule=selection_low)
def selection_high(m, i, j):
    return m.X[i, j] <= bigM*m.Y[i, j]
m.C11 = pyo.Constraint(m.I, m.J, rule=selection_high)

# continuously select blocks in columns.  Use markers for "upper triangle" to omit them

# a square may be selected if previous was, or if previous is in upper triangle
def continuous_selection(m, i, j):
    return m.Y[i, j] <= m.Y[i-1, j] + m.UT[i-1, j]
m.C13 = pyo.Constraint(m.I_not_first, m.J, rule=continuous_selection)
# enforce row-continuity in upper triangle
def upper_triangle_continuous_selection(m, i, j):
    return m.UT[i, j] <= m.UT[i-1, j]
m.C14 = pyo.Constraint(m.I_not_first, m.J, rule=upper_triangle_continuous_selection)
# enforce either-or for selection or membership in upper triangle
def either(m, i, j):
    return m.UT[i, j] + m.Y[i, j] <= 1
m.C15 = pyo.Constraint(m.I, m.J, rule=either)

# OBJ:  Minimze number of selected cells, penalize for off-diagonal selection
def objective(m):
    return sum(m.Y[i, j]*m.weight[i, j] for i in m.I for j in m.J)
#   return sum(sum(m.X[i,j] for j in m.J) - m.row_tots[i] for i in m.I) #+\
#           sum(m.Y[i,j]*m.weight[i,j] for i in m.I for j in m.J)
m.OBJ = pyo.Objective(rule=objective)
    

solver = pyo.SolverFactory('cbc')
results = solver.solve(m)

print(results)
for i in m.I:
    for j in m.J:
        print(f'{m.X[i,j].value : 3.1f}', end='\t')
    print()
print('\npenalty matrix check...')
for i in m.I:
    for j in m.J:
        print(f'{m.weight[i,j] : 3.1f}', end='\t')
    print()

结果

 300.0   127.7   0.0     0.0     0.0     0.0     0.0    
 0.0     12.2    0.0     0.0     0.0     0.0     0.0    
 0.0     165.6   187.1   0.0     0.0     0.0     0.0    
 0.0     0.0     58.3    0.0     0.0     0.0     0.0    
 0.0     0.0     22.7    0.0     0.0     0.0     0.0    
 0.0     0.0     31.9    0.0     0.0     0.0     0.0    
 0.0     0.0     0.0     300.0   96.4    0.0     0.0    
 0.0     0.0     0.0     0.0     29.4    0.0     0.0    
 0.0     0.0     0.0     0.0     171.5   0.0     0.0    
 0.0     0.0     0.0     0.0     10.0    390.0   74.5   
 0.0     0.0     0.0     0.0     0.0     0.0     27.9   
 0.0     0.0     0.0     0.0     0.0     0.0     200.0

Answer 2

如果您已经知道哪些近对角线元素不为零，则它是线性方程组（对于列总和 345 和指定的行总和），但您必须迭代组合。您有 19 个方程，其中有 10 个未知数（非零项的数量），这通常是不可解的。它变得更容易了，因为您可以选择 10 个未知数有帮助，并且只需要大致满足其中的 7 个方程，但我认为只有当您幸运时才存在解决方案（或者这是一个旨在有一个解决方案）。

鉴于 12 行中的每一行都必须有正确的总和，因此您至少需要 12 个非零元素。最有可能的是，每行至少需要两个，每列至少需要两个。

找到具有解的最优集合可能是一个 NP 完全问题，这意味着您必须系统地迭代所有组合，直到找到解。

对于您的示例，大约有 m=31 个矩阵元素；不可能遍历所有组合。您需要反复试验。

这是一个示例代码，允许使用 numpy 的最小二乘求解器优化所有 31 个元素。

import numpy as np

rowsums = np.array([427.7, 12.2, 352.7, 58.3, 22.7, 31.9, 396.4, 29.4, 171.5, 474.5, 27.9, 200])
nrows = len(rowsums)
ncols = 7
colsum_target = 345 # fuzzy target
    
mask = np.array([
       [1, 1, 0, 0, 0, 0, 0],
       [1, 1, 0, 0, 0, 0, 0],
       [1, 1, 1, 0, 0, 0, 0],
       [0, 1, 1, 0, 0, 0, 0],
       [0, 1, 1, 1, 0, 0, 0],
       [0, 0, 1, 1, 1, 0, 0],
       [0, 0, 1, 1, 1, 0, 0],
       [0, 0, 0, 1, 1, 1, 0],
       [0, 0, 0, 1, 1, 1, 0],
       [0, 0, 0, 0, 1, 1, 1],
       [0, 0, 0, 0, 0, 1, 1],
       [0, 0, 0, 0, 0, 1, 1]]).astype(bool)
assert mask.shape == (nrows, ncols)

m = mask.sum() # number of elements to fit

# idx is the index matrix, referring to the element in the x-vector
idx = np.full(mask.shape, -1, dtype=int)
k = 0
for i in range(nrows):
    for j in range(ncols):
        if mask[i, j]:
            idx[i, j] = k
            k += 1
print(f'Index matrix:\n{idx}')

# We're going to solve A @ x = b, where x are the near-diagonal elements
# Shapes: A (nrows+ncols, m); b (nrows+ncols,); x: (m,)
# and b are the ocnditions on the row and column sums.
# Rows A[:nrows] represent the conditions on row sums.
# Rows A[-ncols:] represent the conditions on the column sums.
A = np.zeros((ncol + nrow, m))
for i in range(nrows):
    for j in range(ncols):
        if mask[i, j]:
            A[i, idx[i, j]] = 1
            A[nrows+j, idx[i, j]] = 1
            
b = np.concatenate((rowsums, np.full(ncols, colsum_target, dtype=np.float64)))

# Force priority on row sums (>>1 to match row sums, <<1 to match column sums)
priority = 1000
A[:nrows, :] *= priority
b[:nrows] *= priority

# Get the solution vector x
x, _, _, _ = np.linalg.lstsq(A, b, rcond=None)

# map the elements of x into the matrix template
mat = np.concatenate((x, [0]))[idx] # extra [0] is for the -1 indices
round_mat = np.around(mat, 1)

row_sum_errors = np.around(mat.sum(axis=1)-rowsums, 6)
col_sums = np.around(mat.sum(axis=0), 2)

print(f'mat:\n{round_mat}\nrow_sums error:\n{row_sum_errors}')
print(f'column sums:\n{col_sums}')

这会产生输出：

Index matrix:
[[ 0  1 -1 -1 -1 -1 -1]
 [ 2  3 -1 -1 -1 -1 -1]
 [ 4  5  6 -1 -1 -1 -1]
 [-1  7  8 -1 -1 -1 -1]
 [-1  9 10 11 -1 -1 -1]
 [-1 -1 12 13 14 -1 -1]
 [-1 -1 15 16 17 -1 -1]
 [-1 -1 -1 18 19 20 -1]
 [-1 -1 -1 21 22 23 -1]
 [-1 -1 -1 -1 24 25 26]
 [-1 -1 -1 -1 -1 27 28]
 [-1 -1 -1 -1 -1 29 30]]
mat:
[[210.8 216.9   0.    0.    0.    0.    0. ]
 [  3.1   9.1   0.    0.    0.    0.    0. ]
 [101.1 107.1 144.4   0.    0.    0.    0. ]
 [  0.   10.5  47.8   0.    0.    0.    0. ]
 [  0.  -28.6   8.7  42.6   0.    0.    0. ]
 [  0.    0.   -3.7  30.1   5.5   0.    0. ]
 [  0.    0.  117.8 151.6 127.    0.    0. ]
 [  0.    0.    0.   21.6  -3.   10.8   0. ]
 [  0.    0.    0.   69.   44.3  58.2   0. ]
 [  0.    0.    0.    0.  141.3 155.1 178.1]
 [  0.    0.    0.    0.    0.    2.5  25.4]
 [  0.    0.    0.    0.    0.   88.5 111.5]]
row_sums error:
[-0. -0. -0. -0. -0. -0. -0. -0. -0. -0. -0. -0.]
column sums:
[315.03 315.03 315.03 315.03 315.03 315.03 315.03]

最小二乘求解器无法处理硬约束；如果您发现某一列只是有点超出范围（例如 299），您可以使用相同的 priority 技巧让求解器对该列更加努力。您可以尝试逐个禁用较小的元素（例如 linear programming optimizer，它更适合同时具有硬相等要求和边界的问题。