模型
代码
import re
from pyscipopt import Model, quicksumclass Constrain:def __init__(self, n, rowConstrain, columnConstrain):self.n = nself.rowConstrain = rowConstrainself.columnConstrain = columnConstrain#X2constrain
def X2constrain(X):n = len(X)rowConstrain = []for i in range(n):row = []num = None#若遇到第一个1，则开始计数，若遇到0，则停止计数for j in range(n):if X[i][j] == 1:if num is None:num = 1else:num += 1else:if num is None:continueelse:row.append(num)num = None#把最后一个数加进去if num is not None:row.append(num)rowConstrain.append(row)columnConstrain = []for j in range(n):column = []num = None#若遇到第一个1，则开始计数，若遇到0，则停止计数for i in range(n):if X[i][j] == 1:if num is None:num = 1else:num += 1else:if num is None:continueelse:column.append(num)num = None#把最后一个数加进去if num is not None:column.append(num)columnConstrain.append(column)constrain = Constrain(n, rowConstrain, columnConstrain)return constraindef case():#1n = 5rowConstrain = [[2],[5],[2,1],[2],[3]]columnConstrain = [[1,2],[5],[3,1],[1],[2]]constrain = Constrain(n, rowConstrain,columnConstrain)X0 = [[0,1,1,0,0],[1,1,1,1,1],[0,1,1,0,1],[1,1,0,0,0],[1,1,1,0,0]]constrain.X0 = X0#2n = 15rowConstrain = [[4,2,3],[4,6],[5,1,2],[1,4,1],[6,1,1],[6,1],[4,4],[6,1,1],[5,1,2],[1,1,2],[2,1,1],[1,2,1],[3,2,1],[7,1,4],[7,6]]columnConstrain = [[1,1,5],[1,1,3],[2,2,3],[4,2,2],[6,4,2],[9,2],[7,2],[1,1,5],[2,3,1],[1,3,2],[3,3,1],[2,1,4],[3,4],[2,2,2],[2,8]]constrain = Constrain(n, rowConstrain,columnConstrain)return constrain#打印X
def XPrint(X):for i in range(len(X)):print("[" + " ".join([str(c) for c in X[i]])+"]")print("")def IP(constrain, timeLimit=100, isPrint=True):n = constrain.nR = constrain.rowConstrainC = constrain.columnConstrainM = 1000# 建模求解model = Model("IP")X = {(i, j): model.addVar(vtype="B", name="X[%s, %s]" % (i, j)) for i in range(n) for j in range(n)}Y = {(i, j, k): model.addVar(vtype="B", name="Y[%s, %s, %s]" % (i, j, k)) for i in range(n) for j in range(n) for k in range(len(R[i]))}Z = {(i, j, l): model.addVar(vtype="B", name="Z[%s, %s, %s]" % (i, j, l)) for i in range(n) for j in range(n) for l in range(len(C[j]))}S = {(i, k): model.addVar(vtype="I", name="S[%s, %s]" % (i, k)) for i in range(n) for k in range(len(R[i]))}T = {(j, l): model.addVar(vtype="I", name="T[%s, %s]" % (j, l)) for j in range(n) for l in range(len(C[j]))}#随便设置个目标model.setObjective(quicksum(X[i, j] for i in range(n) for j in range(n)), "maximize")#设置可行范围for i in range(n):for k in range(len(R[i])):model.addCons(S[i, k] >= 0)model.addCons(S[i, k] <= n-1)for j in range(n):for l in range(len(C[j])):model.addCons(T[j, l] >= 0)model.addCons(T[j, l] <= n-1)#每行在某一个区间范围内时，X为1for i in range(n):for k in range(len(R[i])):for j in range(n):model.addCons(S[i, k] <= j + M * (1 - Y[i, j, k]))model.addCons(S[i, k]+R[i][k]-1 >= j - M * (1 - Y[i, j, k]))#Y转Xfor i in range(n):for j in range(n):model.addCons(quicksum(Y[i, j, k] for k in range(len(R[i]))) >= X[i, j])model.addCons(quicksum(Y[i, j, k] for k in range(len(R[i]))) <= M*X[i, j])#每行的某个起点的取值与前面一个的尾巴中间至少隔一个for i in range(n):for k in range(len(R[i])-1):model.addCons(S[i, k] + R[i][k] + 1 <= S[i,k+1])#每列在某一个区间范围内时，X为1for j in range(n):for l in range(len(C[j])):for i in range(n):model.addCons(T[j, l] <= i + M * (1 - Z[i, j, l]))model.addCons(T[j, l]+C[j][l]-1 >= i - M * (1 - Z[i, j, l]))#Z转Xfor i in range(n):for j in range(n):model.addCons(quicksum(Z[i, j, l] for l in range(len(C[j]))) >= X[i, j])model.addCons(quicksum(Z[i, j, l] for l in range(len(C[j]))) <= M*X[i, j])#每行的某个起点的取值与前面一个的尾巴中间至少隔一个for j in range(n):for l in range(len(C[j])-1):model.addCons(T[j, l] + C[j][l] + 1 <= T[j,l+1])# 设置求解时间model.setRealParam("limits/time", timeLimit)# model.hideOutput()model.optimize()if isPrint:print("\ngap:", model.getGap())X1 = [[round(model.getVal(X[i, j])) for j in range(n)] for i in range(n)]return X1if __name__ == '__main__':#数据录入constrain = case()#IP求解X = IP(constrain)XPrint(X)
