根据给定的约束创建矩阵
Create a matrix based on given constraints
我正在尝试创建具有以下约束的矩阵。
- 列总和应介于 300 和 390 之间，包括这两个值。
- 行总和应等于每行用户指定的值。
- 矩阵中的非零值不得小于 10。
- 给定列中非零值的数量不应超过 4。
- 列应按对角线顺序排列。
如果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
输出
如果您已经知道哪些近对角元素是非零的，它是线性方程组（对于列总和 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 技巧让求解器对该列尝试更努力一些。您可以尝试一个一个地禁用小元素（例如 <10）。您也可以尝试使用 linear programming optimizer，它更适合同时具有硬相等要求和边界的问题。
我认为您将其设置为 LP 的做法是正确的。可以表述为MIP。
我在这里没有修改任何类型的输入，我不确定在您的约束条件下，您是否能保证所有输入的结果都是可行的。
我惩罚了 off-diagonal 选择以鼓励对角线上的东西，并设置了一些“选择完整性”约束来强制执行 block-selection。
在大约 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
我正在尝试创建具有以下约束的矩阵。
- 列总和应介于 300 和 390 之间，包括这两个值。
- 行总和应等于每行用户指定的值。
- 矩阵中的非零值不得小于 10。
- 给定列中非零值的数量不应超过 4。
- 列应按对角线顺序排列。
如果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
输出
如果您已经知道哪些近对角元素是非零的，它是线性方程组（对于列总和 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 技巧让求解器对该列尝试更努力一些。您可以尝试一个一个地禁用小元素（例如 <10）。您也可以尝试使用 linear programming optimizer，它更适合同时具有硬相等要求和边界的问题。
我认为您将其设置为 LP 的做法是正确的。可以表述为MIP。
我在这里没有修改任何类型的输入，我不确定在您的约束条件下，您是否能保证所有输入的结果都是可行的。
我惩罚了 off-diagonal 选择以鼓励对角线上的东西，并设置了一些“选择完整性”约束来强制执行 block-selection。
在大约 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