问题描述
我正在尝试找到可以最大化我的总和值的最佳组合，但是必须在2个特定约束下，因此我假设线性编程将是最合适的。
问题如下： 一些教育性的世界性事件希望聚集世界上最聪明的青少年学生。 每个州都通过以下考试对10万名学生进行了测试：“ MATH”，“ ENGLISH”，“ COmpuTERS”，“ HISTORY”，“ PHYSICS” .. EACH考试的得分为0-100。
要求每个州从经过测试的10万名学生中向他们发送最好的1万名。
您作为法国代表，被要求从经过测试的您所在国家的10万名学生中选出排名前1万名的学生。为此，您需要从中优化MAX VALUE，以获得最佳的总得分。
但是有 2个主要约束：
1-从所选的1万名学生中，您需要分配特定的学生，仅针对上述5个主题中的1个特定主题对事件进行测试。 所需的分配是：['MATH'：4000，'ENGLISH'：3000，'COmpuTERS'：2000，'HISTORY'：750，'PHYSICS'：250]
2-每个“考试科目”得分的权重必须不同。.例如：97是数学，其历史价值超过97。 高度为：['MATH'：1.9，'ENGLISH'：1.7，'COmpuTERS'：1.5，'HISTORY'：1.3，'PHYSICS'：1.1]
我的解决方案： 我尝试将PULP（python）用作LP库并正确解决了该问题，但运行耗时超过2小时。 您能找到一种更好（更快，更简单..）的解决方法吗？ 有一些NUMPY LP函数可以代替使用，也许会更快一些？ 它应该是一个简单的优化问题，因为我使它变得太慢且太复杂了。 -解决方案只需要使用Python
例如，让我们从小规模来看同样的问题： 有30名学生，您只需要选择15名学生，就以下学科分配需求而言，这将为我们提供最佳组合。 所需的分配是-['MATH'：5，'ENGLISH'：4，'COmpuTERS'：3，'HISTORY'：2，'PHYSICS'：1]
这是所有30名学生及其成绩：
这是我针对原始问题的完整代码（10万名学生）：
import pandas as pd
import numpy as np
import pulp as p
import time
t0=time.time()
demand = [4000,3000,2000,750,250]
weight = [1.9,1.7,1.5,1.3,1.1]
original_data= pd.read_csv('GRADE_100K.csv') #created simple csv file with random scores
data_c=original_data.copy()
data_c.index = np.arange(1,len(data_c)+1)
data_c.columns
data_c=data_c[['STUDENT_ID','MATH','ENGLISH','COmpuTERS','HISTORY','PHYSICS']]
#DataFrame Shape
m=data_c.shape[1]
n=data_c.shape[0]
data=[]
sublist=[]
for j in range(0,n):
for i in range(1,m):
sublist.append(data_c.iloc[j,i])
data.append(sublist)
sublist=[]
def _get_num_students(data):
return len(data)
def _get_num_subjects(data):
return len(data[0])
def _get_weighted_data(data,weight):
return [
[a*b for a,b in zip(row,weight)]
for row in data
]
data = _get_weighted_data(data,weight)
num_students = _get_num_students(data)
num_subjects = _get_num_subjects(data)
# Create a LP Minimization problem
Lp_prob = p.LpProblem('Problem',p.LpMaximize)
# Create problem Variables
variables_matrix = [[0 for i in range(num_subjects)] for j in range(num_students)]
for i in range(0,num_students):
for j in range(0,num_subjects):
variables_matrix[i][j] = p.LpVariable(f"X({i+1},{j+1})",1,cat='Integer')
df_d=pd.DataFrame(data=data)
df_v=pd.DataFrame(data=variables_matrix)
ml=df_d.mul(df_v)
ml['coeff'] = ml.sum(axis=1)
coefficients=ml['coeff'].tolist()
# DEALING WITH TARGET FUNCTION VALUE
suming=0
k=0
sumsum=[]
for z in range(len(coefficients)):
suming +=coefficients[z]
if z % 2000==0:
sumsum.append(suming)
suming=0
if z<2000:
sumsum.append(suming)
sumsuming=0
for s in range(len(sumsum)):
sumsuming=sumsuming+sumsum[s]
Lp_prob += sumsuming
# DEALING WITH the 2 CONSTRAINS
# 1-subject constraints
con1_suming=0
for e in range(num_subjects):
L=df_v.iloc[:,e].to_list()
for t in range(len(L)):
con1_suming +=L[t]
Lp_prob += con1_suming <= demand[e]
con1_suming=0
# 2- students constraints
con2_suming=0
for e in range(num_students):
L=df_v.iloc[e,:].to_list()
for t in range(len(L)):
con2_suming +=L[t]
Lp_prob += con2_suming <= 1
con2_suming=0
print("time taken for TARGET+CONSTRAINS %8.8f seconds" % (time.time()-t0) )
t1=time.time()
status = Lp_prob.solve() # Solver
print("time taken for SOLVER %8.8f seconds" % (time.time()-t1) ) # 632 SECONDS
print(p.LpStatus[status]) # The solution status
print(p.value(Lp_prob.objective))
df_v=pd.DataFrame(data=variables_matrix)
# Printing the final solution
lst=[]
val=[]
for i in range(0,num_students):
lst.append([p.value(variables_matrix[i][j]) for j in range(0,num_subjects)])
val.append([sum([p.value(variables_matrix[i][j]) for j in range(0,num_subjects)]),i])
ones_places=[]
for i in range (0,len(val)):
if val[i][0]==1:
ones_places.append(i+1)
len(ones_places)
data_once=data_c[data_c['STUDENT_ID'].isin(ones_places)]
IDs=[]
for i in range(len(ones_places)):
IDs.append(data_once['STUDENT_ID'].to_list()[i])
course=[]
sub_course=[]
for i in range(len(lst)):
j=0
sub_course='x'
while j<len(lst[i]):
if lst[i][j]==1:
sub_course=j
j=j+1
course.append(sub_course)
coures_ones=[]
for i in range(len(course)):
if course[i]!= 'x':
coures_ones.append(course[i])
# adding the COURSE name to the final table
# NUMBER OF DICTIONARY KEYS based on number of COURSES
col=original_data.columns.values[1:].tolist()
dic = {0:col[0],1:col[1],2:col[2],3:col[3],4:col[4]}
cc_name=[dic.get(n,n) for n in coures_ones]
one_c=[]
if len(IDs)==len(cc_name):
for i in range(len(IDs)):
one_c.append([IDs[i],cc_name[i]])
prob=[]
if len(IDs)==len(cc_name):
for i in range(len(IDs)):
prob.append([IDs[i],cc_name[i],data_once.iloc[i][one_c[i][1]]])
scoring_table = pd.DataFrame(prob,columns=['STUDENT_ID','COURES','score'])
scoring_table.sort_values(by=['COURES','score'],ascending=[False,False],inplace=True)
scoring_table.index = np.arange(1,len(scoring_table)+1)
print(scoring_table)
解决方法
关于我使用最低成本流的想法，还有更多想法。
我们通过采用四层有向图来建模此问题，其中每一层都完全连接到下一层。
节点
-
第一层：单个节点
s将成为我们的来源。 -
第二层：每个学生一个节点。
-
第三层：每个主题一个节点。
-
第四层：OA单节点
t将成为我们的消耗。
边缘容量
-
第一->第二：所有边的容量为1。
-
第二->第三：所有边的容量为1。
-
第三->第四：所有边的容量与必须分配给该学科的学生人数相对应。
边缘成本
-
第一->第二：所有边的成本均为0。
-
第二个->第三点：请记住，该层的边缘将学生与主题联系在一起。这些费用将与学生在该主题上获得的加权分数成反比。
cost = -subject_weight*student_subject_score。 -
第三->第四：所有边的成本均为0。
然后我们要求从s到t的流量等于我们必须选择的学生人数。
为什么这样做？
通过将第三层和第四层之间的所有边作为赋值，最小成本流问题的解决方案将与您的问题的解决方案相对应。
每个学生最多可以选择一个学科，因为相应的节点只有一个入学边缘。
每个学科都有确切的所需学生人数，因为传出容量对应于我们必须为该学科选择的学生人数，我们必须利用这些优势的全部容量，因为我们无法满足流动需求否则。
对MCF问题的最小解决方案与您问题的最大解决方案相对应，因为成本与它们给出的价值相对应。
当您在python中寻求解决方案时，我使用ortools实现了最小成本流问题。在我的colab笔记本中，找到解决方案的时间不到一秒钟。需要“很长时间”是提取溶液。但是，包括设置和解决方案提取在内，对于100000名学生的全部问题，我的运行时间仍不到20秒。
代码
# imports
from ortools.graph import pywrapgraph
import numpy as np
import pandas as pd
import time
t_start = time.time()
# setting given problem parameters
num_students = 100000
subjects = ['MATH','ENGLISH','COMPUTERS','HISTORY','PHYSICS']
num_subjects = len(subjects)
demand = [4000,3000,2000,750,250]
weight = [1.9,1.7,1.5,1.3,1.1]
# generating student scores
student_scores_raw = np.random.randint(101,size=(num_students,num_subjects))
# setting up graph nodes
source_nodes = [0]
student_nodes = list(range(1,num_students+1))
subject_nodes = list(range(num_students+1,num_subjects+num_students+1))
drain_nodes = [num_students+num_subjects+1]
# setting up the min cost flow edges
start_nodes = [int(c) for c in (source_nodes*num_students + [i for i in student_nodes for _ in subject_nodes] + subject_nodes)]
end_nodes = [int(c) for c in (student_nodes + subject_nodes*num_students + drain_nodes*num_subjects)]
capacities = [int(c) for c in ([1]*num_students + [1]*num_students*num_subjects + demand)]
unit_costs = [int(c) for c in ([0.]*num_students + list((-student_scores_raw*np.array(weight)*10).flatten()) + [0.]*num_subjects)]
assert len(start_nodes) == len(end_nodes) == len(capacities) == len(unit_costs)
# setting up the min cost flow demands
supplies = [sum(demand)] + [0]*(num_students + num_subjects) + [-sum(demand)]
# initialize the min cost flow problem instance
min_cost_flow = pywrapgraph.SimpleMinCostFlow()
for i in range(0,len(start_nodes)):
min_cost_flow.AddArcWithCapacityAndUnitCost(start_nodes[i],end_nodes[i],capacities[i],unit_costs[i])
for i in range(0,len(supplies)):
min_cost_flow.SetNodeSupply(i,supplies[i])
# solve the problem
t_solver_start = time.time()
if min_cost_flow.Solve() == min_cost_flow.OPTIMAL:
print('Best Value:',-min_cost_flow.OptimalCost()/10)
print('Solver time:',str(time.time()-t_solver_start)+'s')
print('Total Runtime until solution:',str(time.time()-t_start)+'s')
#extracting the solution
solution = []
for i in range(min_cost_flow.NumArcs()):
if min_cost_flow.Flow(i) > 0 and min_cost_flow.Tail(i) in student_nodes:
student_id = min_cost_flow.Tail(i)-1
subject_id = min_cost_flow.Head(i)-1-num_students
solution.append([student_id,subjects[subject_id],student_scores_raw[student_id,subject_id]])
assert(len(solution) == sum(demand))
solution = pd.DataFrame(solution,columns = ['student_id','subject','score'])
print(solution.head())
else:
print('There was an issue with the min cost flow input.')
print('Total Runtime:',str(time.time()-t_start)+'s')
通过以下列表理解（不是每次迭代也不使用列表查找）来替换上面代码中用于解决方案提取的for循环，可以显着改善运行时间。但是出于可读性原因，我也将旧解决方案留在这里。这是新的：
solution = [[min_cost_flow.Tail(i)-1,subjects[min_cost_flow.Head(i)-1-num_students],student_scores_raw[min_cost_flow.Tail(i)-1,min_cost_flow.Head(i)-1-num_students]
]
for i in range(min_cost_flow.NumArcs())
if (min_cost_flow.Flow(i) > 0 and
min_cost_flow.Tail(i) <= num_students and
min_cost_flow.Tail(i)>0)
]
以下输出为新的更快的实现提供了运行时。
输出
Best Value: 1675250.7
Solver time: 0.542395830154419s
Total Runtime until solution: 1.4248979091644287s
student_id subject score
0 3 ENGLISH 99
1 5 MATH 98
2 17 COMPUTERS 100
3 22 COMPUTERS 100
4 33 ENGLISH 100
Total Runtime: 1.752336025238037s
请指出我可能犯的任何错误。
我希望这会有所帮助。 ;）
,我认为您对此很接近。这是一个相当标准的整数线性程序（ILP）分配问题。由于问题的结构，它会有点慢。
您没有在帖子中说设置和解决时间的细目分类。我看到您正在读取文件并正在使用熊猫。我认为熊猫会因优化问题而很快变得笨拙，但这只是个人喜好。
我使用pyomo解算器在cbc中编写了您的问题，我敢肯定，pulp使用的解算器与您的问题相同。 （见下文）。我认为您有两个约束和一个双索引的二进制决策变量就可以了。
如果我将其削减为10,000名学生（不懈怠，只需一对一配对），它就能在14秒内解决比较问题。我的设置是一台5岁的iMac，带有很多ram。
在池中有10万名学生的情况下运行，在调用求解器之前，它会在大约25分钟的时间内完成10秒的“设置”时间。所以我不太确定为什么您的编码要花2个小时。如果您可以减少求解器的时间，那将会有所帮助。其余的应该是微不足道的。我并没有在输出中花太多时间，但是980K的OBJ函数值似乎是合理的。
其他想法：
如果您可以使求解器选项正确配置并且将mip间隙设置为0.05左右，则可以接受，如果您可以接受稍微不理想的解决方案，则可以加快速度。我对像Gurobi这样的付费求解器只有很不错的求解器选项。我使用免费赠品解算器YMMV并没有那么幸运。
import pyomo.environ as pyo
from random import randint
from time import time
# start setup clock
tic = time()
# exam types
subjects = ['Math','English','Computers','History','Physics']
# make set of students...
num_students = 100_000
students = [f'student_{s}' for s in range(num_students)]
# make 100K fake scores in "flat" format
student_scores = { (student,subj) : randint(0,100)
for student in students
for subj in subjects}
assignments = { 'Math': 4000,'English': 3000,'Computers': 2000,'History': 750,'Physics': 250}
weights = {'Math': 1.9,'English': 1.7,'Computers': 1.5,'History': 1.3,'Physics': 1.1}
# Set up model
m = pyo.ConcreteModel('exam assignments')
# Sets
m.subjects = pyo.Set(initialize=subjects)
m.students = pyo.Set(initialize=students)
# Parameters
m.assignments = pyo.Param(m.subjects,initialize=assignments)
m.weights = pyo.Param(m.subjects,initialize=weights)
m.scores = pyo.Param(m.students,m.subjects,initialize=student_scores)
# Variables
m.x = pyo.Var(m.students,within=pyo.Binary) # binary selection of pairing student to test
# Objective
m.OBJ = pyo.Objective(expr=sum(m.scores[student,subject] * m.x[student,subject]
for student in m.students
for subject in m.subjects),sense=pyo.maximize)
### Constraints ###
# fill all assignments
def fill_assignments(m,subject):
return sum(m.x[student,subject] for student in m.students) == assignments[subject]
m.C1 = pyo.Constraint(m.subjects,rule=fill_assignments)
# use each student at most 1 time
def limit_student(m,student):
return sum(m.x[student,subject] for subject in m.subjects) <= 1
m.C2 = pyo.Constraint(m.students,rule=limit_student)
toc = time()
print (f'setup time: {toc-tic:0.3f}')
tic = toc
# solve it..
solver = pyo.SolverFactory('cbc')
solution = solver.solve(m)
print(solution)
toc = time()
print (f'solve time: {toc-tic:0.3f}')
输出
setup time: 10.835
Problem:
- Name: unknown
Lower bound: -989790.0
Upper bound: -989790.0
Number of objectives: 1
Number of constraints: 100005
Number of variables: 500000
Number of binary variables: 500000
Number of integer variables: 500000
Number of nonzeros: 495094
Sense: maximize
Solver:
- Status: ok
User time: -1.0
System time: 1521.55
Wallclock time: 1533.36
Termination condition: optimal
Termination message: Model was solved to optimality (subject to tolerances),and an optimal solution is available.
Statistics:
Branch and bound:
Number of bounded subproblems: 0
Number of created subproblems: 0
Black box:
Number of iterations: 0
Error rc: 0
Time: 1533.8383190631866
Solution:
- number of solutions: 0
number of solutions displayed: 0
solve time: 1550.528