pythonnumpyoptimizationpyomoconstraint-programming
Create a matrix based on given constraints
I am trying to create a matrix with the following constraints.
- The column sum should be between 300 and 390, both values inclusive.
- Row sum should be equal to user-specified values per row.
- No non-zero value in the matrix should be less than 10.
- The count of non-zero values in a given column should not exceed 4.
- The columns should be arranged in a diagonal order.
if UserInput = [427.7, 12.2, 352.7, 58.3, 22.7, 31.9, 396.4, 29.4, 171.5, 474.5, 27.9, 200]
I want output matrix something like this,
Edit 1
I have tried the following approach using Pyomo, however, I got stuck on 5th constraint that column values should be diagonally aligned in the matrix
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
Output
Solution
I think you were on the right track with setting this up as an LP. It can be formulated as a MIP.
I haven't tinkered with any variety of inputs here, and I'm not sure you are guaranteed feasible results for all inputs with the constraints you have.
I penalized off-diagonal selection to encourage things on diagonal, and set up some "selection integrality" constraints to enforce block-selection.
Solves in about 1/10 of second...
# 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()
Result
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