python求解运输问题-【数学建模】线性规划各种问题的Python调包方法

python求解运输问题_【数学建模】线性规划各种问题的Python调包⽅法关键词：Python、调包、线性规划、指派问题、运输问题、pulp、混合整数线性规划(MILP)注：此⽂章是线性规划的调包实现，具体步骤原理请搜索具体解法。本⽂章的各个问题可能会采⽤多种调⽤⽅法，为什么？因为这些包各有特点，有些语法特别像matlab，只要稍稍改变即可达成代码交换；⽽有些包利⽤了python本⾝的特性，在灵活度与代码的可读性上更⾼。我认为这些包各有优劣，各位各持所需吧。看了本⽂章能做到什么？你可以在本⽂章内学到线性规划的⼏个问题的求解⽅式，并学会如何⽤pulp包解决线性规划问题。⽆论是整数规划(IntegerProgram)、01规划(BinaryProgram)还是混合整数线性规划(MILP)，你都可以得到很好的解题⽅法。⼀、线性规划该问题引⽤⾃《数学建模算法与应⽤-司守奎》第⼀章线性规划3.线性规划包的具体使⽤可参考scipy官⽹求解最普通的线性规划问题：scipy调包代码importnumpyasnpz=np.array([2,3,1])a=np.array([[1,4,2],[3,2,0]])b=np.array([8,6])x1_bound=x2_bound=x3_bound=(0,None)fromscipyimportoptimizeres=optimize.linprog(z,A_ub=-a,b_ub=-b,bounds=(x1_bound,x2_bound,x3_bound))print(res)#output:#fun:7.0#message:'Optimizationterminatedsuccessfully.'#nit:2#slack:array([0.,0.])#status:0#success:True#x:array([0.8,1.8,0.])注意，此函数和matlab的linprog有很多相似之处。⾸先默认就是求解最⼩值，如果想要求最⼤值，需要对⽬标函数的各参数取反(既对z取反)，并在得出的结果(func)前取反。并且所有约束条件都为≤，因此例⼦中传⼊参数是前⾯都加了负号。bound为边界的⼆元元组，None时为⽆边界。如果存在类似

这种情况，可以：A_eq=[[1,2,4]]b_eq=[101]并在linprog中传⼊。得出的结果为scipy.optimize.optimize.OptimizeResult(优化结果)类型，是类似字典的结构，例如想要优化结果代⼊⽬标函数的值，可以res.fun或res['fun']这样取值。pulp调包代码importpulp#⽬标函数的系数z=[2,3,1]#约束a=[[1,4,2],[3,2,0]]b=[8,6]#确定最⼤化最⼩化问题，最⼤化只要把Min改成Max即可m=pulp.LpProblem(sense=pulp.LpMinimize)#定义三个变量放到列表中x=[pulp.LpVariable(f'x{i}',lowBound=0)foriin[1,2,3]]#定义⽬标函数，lpDot可以将两个列表的对应位相乘再加和#相当于z[0]*x[0]+z[1]*x[0]+z[2]*x[2]m+=pulp.lpDot(z,x)#设置约束条件foriinrange(len(a)):m+=(pulp.lpDot(a[i],x)>=b[i])#求解m.solve()#输出结果print(f'优化结果：{pulp.value(m.objective)}')print(f'参数取值：{[pulp.value(var)forvarinx]}')#output:#优化结果：7.0#参数取值：[2.0,0.0,3.0]每⼀步的说明已经注释在代码中，可以看到输出结果，两者的变量取值并不⼀致，但代⼊⽬标函数的结果都是⼀样的。同样的，如果存在类似这种情况，可以：

A_eq=[1,2,4]b_eq=101m+=(pulp.lpDot(A_eq,x)==b_eq)⼆、运输问题某商品有个产地、个销地，各产地的产量分别为，各销地的需求量分别为。若该商品由产地运到销地的单位运价为，问应该如何调运才能使总运费最省？引⼊变量，其取值为由产地运往销地的该商品数量，数学模型为:例题：pulp调包代码importpulpimportnumpyasnpfrompprintimportpprintdeftransportation_problem(costs,x_max,y_max):row=len(costs)col=len(costs[0])prob=pulp.LpProblem('TransportationProblem',sense=pulp.LpMaximize)var=[[pulp.LpVariable(f'x{i}{j}',lowBound=0,cat=pulp.LpInteger)forjinrange(col)]foriinrange(row)]flatten=lambdax:[yforlinxforyinflatten(l)]iftype(x)islistelse[x]prob+=pulp.lpDot(flatten(var),costs.flatten())foriinrange(row):prob+=(pulp.lpSum(var[i])<=x_max[i])forjinrange(col):prob+=(pulp.lpSum([var[i][j]foriinrange(row)])<=y_max[j])prob.solve()

return{'objective':pulp.value(prob.objective),'var':[[pulp.value(var[i][j])forjinrange(col)]foriinrange(row)]}然后构造参数传递进去：if__name__=='__main__':costs=np.array([[500,550,630,1000,800,700],[800,700,600,950,900,930],[1000,960,840,650,600,700],[1200,1040,980,860,880,780]])max_plant=[76,88,96,40]max_cultivation=[42,56,44,39,60,59]res=transportation_problem(costs,max_plant,max_cultivation)print(f'最⼤值为{res["objective"]}')print('各变量的取值为：')pprint(res['var'])#output:#最⼤值为284230.0#各变量的取值为：#[[0.0,0.0,6.0,39.0,31.0,0.0],#[0.0,0.0,0.0,0.0,29.0,59.0],#[2.0,56.0,38.0,0.0,0.0,0.0],#[40.0,0.0,0.0,0.0,0.0,0.0]]三、指派问题拟分配n⼈去⼲项⼯作，没⼈⼲且仅⼲⼀项⼯作，若分配第⼈去⼲第项⼯作，需花费单位时间，问应如何分配⼯作才能使公认花费的总时间最少？假设指派问题的系数矩阵为：  引⼊变量，若分配⼲⼯作，则取，否则取，上述指派问题的数学模型为  指派问题的可⾏解⽤矩阵表⽰，每⾏每列有且只有⼀个元素为1，其余元素均为0。

调⽤scipy解决原书使⽤匈⽛利算法解决的，在这⾥我们⽤scipy的优化模块解决importnumpyasnpfromscipy.optimizeimportlinear_sum_assignment引⼊包，linear_sum_assignment是专门⽤于解决指派问题的模块。efficiency_matrix=np.array([[12,7,9,7,9],[8,9,6,6,6],[7,17,12,14,12],[15,14,6,6,10],[4,10,7,10,6]])row_index,col_index=linear_sum_assignment(efficiency_matrix)print(row_index+1)print(col_index+1)print(efficiency_matrix[row_index,col_index])#output：#[12345]#[23145]#[76766]定义了开销矩阵(指派问题的系数矩阵)efficiency_matrix，传⼊linear_sum_assignment，结果返回的是最优指派的⾏和列，例如第⼀⾏选择第⼆列，意为：将第⼀个⼈派往第⼆个⼯作。⽽根据numpy.array的性质，传⼊⾏和列就会返回⾏列所对应的值，即为输出的第三列print(efficiency_matrix[row_index,col_index].sum())#output：#32对其求和，即可得到指派问题的最⼩时间。调⽤pulp解决先定义通⽤解决⽅法，其中的flatten是递归展开列表⽤的。defassignment_problem(efficiency_matrix):row=len(efficiency_matrix)col=len(efficiency_matrix[0])flatten=lambdax:[yforlinxforyinflatten(l)]iftype(x)islistelse[x]m=pulp.LpProblem('assignment',sense=pulp.LpMinimize)var_x=[[pulp.LpVariable(f'x{i}{j}',cat=pulp.LpBinary)forjinrange(col)]foriinrange(row)]

m+=pulp.lpDot(efficiency_matrix.flatten(),flatten(var_x))foriinrange(row):m+=(pulp.lpDot(var_x[i],[1]*col)==1)forjinrange(col):m+=(pulp.lpDot([var_x[i][j]foriinrange(row)],[1]*row)==1)m.solve()print(m)return{'objective':pulp.value(m.objective),'var':[[pulp.value(var_x[i][j])forjinrange(col)]foriinrange(row)]}然后定义矩阵，输⼊求解efficiency_matrix=np.array([[12,7,9,7,9],[8,9,6,6,6],[7,17,12,14,9],[15,14,6,6,10],[4,10,7,10,9]])res=assignment_problem(efficiency_matrix)print(f'最⼩值{res["objective"]}')print(res['var'])#output#最⼩值32.0#[[0.0,1.0,0.0,0.0,0.0],[0.0,0.0,0.0,1.0,0.0],[0.0,0.0,0.0,0.0,1.0],[0.0,0.0,1.0,0.0,0.0],[1.0,0.0,0.0,0.0,0.0]]四、pulp的使⽤⽅式基本使⽤可以看出，pulp在线性规划这⼀部分，有更多的通⽤性，编写程序更⾃由。前⾯的例⼦已经⾜够丰富了，但是如果读到这⾥还不去清楚pulp的使⽤⽅式的话……可以去百度⼀下，我这⾥也简单讲⼀讲。⾸先是定义⼀个问题，第⼀个参数为问题的名称，不过并⾮必要的参数，⽽通过sense参数可决定优化的是最⼤值还是最⼩值prob=pulp.LpProblem('problemname',sense=pulp.LpMinimize)然后是定义变量：x0=pulp.LpVariable('x0',lowBound=0,upBound=None,cat=pulp.LpInteger)x1=pulp.LpVariable('x1',lowBound=0,upBound=None,cat=pulp.LpInteger)x2=pulp.LpVariable('x2',lowBound=0,upBound=None,cat=pulp.LpInteger)这⾥定义了三个变量，第⼀个参数为变量名，lowBound、upBound为上下限，cat为变量类型，默认为连续变量，还可以设为离散变量或⼆值变量，具体怎么设置？

如上述代码所⽰，pulp.LpInteger为离散变量，pulp.LpBinary为⼆值变量，同时也可以传⼊'Integer'字符串的⽅式指明变量类型。从上⾯⼏个问题的代码可以看出，我⼏乎没有单个定义变量，⽽是批量定义的。然后是添加⽬标函数：prob+=2*x0-5*x1+4*x2只要让问题(prob)直接加变量的表达式即可添加⽬标函数。再之后是添加约束：prob+=(x0+x1-6*x2<=120)只要让问题(prob)直接加变量的判断式即可添加约束prob.solve()调⽤solve⽅法解出答案，如果省去这⼀步，后⾯的变量和结果值都会显⽰None。print(pulp.value(prob.objective))print(pulp.value(x0))打印优化结果，并显⽰当优化达到结果时x0的取值。思考程序本质problem对象是如何通过不断加来获得⽬标函数和约束的？熟悉python或者c++的朋友可能会想到⼀个词：操作符重载。没错，就是这么实现的，上述的对象⼏乎都