<p>我认为你把它设置为LP是正确的。它可以表示为MIP</p>
<p>我没有在这里修改过任何种类的输入,我也不确定您是否能保证所有具有约束的输入都有可行的结果</p>
<p>我惩罚非对角选择以鼓励对角选择,并设置一些“选择完整性”约束以强制块选择</p>
<p>在大约1/10秒内解决</p>
<pre><code># 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()
</code></pre>
<h3>结果</h3>
<pre><code> 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
</code></pre>