0


分段模型线性化(PWL)【Python|Gurobi实现】

1 概述

在许多问题中,可能既包含着纯线性函数,也包含着分段线性函数。在这一节,我们使用运输问题解释处理分段线性模型的各种方法。

从几何的观点,下图为一个经典的分段线性函数。这个函数由四条线段构成,分段点为4、5 和 7。这三个分段点,将整个函数分为(-∞, 4) 、 [4, 5) 、 [5, 7) 和 [7, ∞) 四个区间。

每个区间内线段的斜率是一个常数,但不同区间之间的线段斜率是不同的。斜率的变化的点就是分段点。

分段线性函数常用于表示或逼近非线性一元函数。

2 算例及Python|Gurobi实现

2.1 算例

考虑经典的运输问题,即:有一些供应商(节点i)和一些需求(节点j),确定每个供应商的产量(P_{i}) 、用户的需求量(D_{i}) ,以及在给定每条路径运输成本(C_{ij})的情况下,最小化运输成本。

\begin{array}{l} \min _{x_{i j}} \mathrm{OF}=\sum_{i j} C_{i j} x_{i j}^{2}\\ \left\{\begin{array}{ll} P_{i}^{\min } \leq P_{i} \leq P_{i}^{\max }, & \text { if unit } i \text { is on, } \\ P_{i}=0 & \text { if unit } i \text { is } \end{array}\right.\\ \sum_{i} x_{i j} \geq D_{j}\\ \sum_{j} x_{i j}=P_{i}\\ 0 \leq x_{i j} \leq x_{i j}^{\max } \end{array}

2.2 Python|Gurobi实现

  1. # -*- coding: utf-8 -*-
  2. from gurobipy import *
  3. import numpy as np
  4. N_i=3
  5. N_j=4
  6. k=1.8
  7. C=np.array([
  8. [0.0755,0.0655,0.0498,0.0585],
  9. [0.0276,0.0163,0.0960,0.0224],
  10. [0.0680,0.0119,0.0340,0.0751]
  11. ])
  12. Pmin=np.array([100,50,30])
  13. Pmax=np.array([450,350,500])
  14. demand=np.array([217,150,145,244])
  15. # 创建模型
  16. M_PWL=Model("PWL-1")
  17. # 变量声明
  18. x =M_PWL.addVars(N_i,N_j,lb=0,ub=100, name="x")
  19. P =M_PWL.addVars(N_i,lb=0,ub=Pmax, name="P")
  20. U =M_PWL.addVars(N_i,vtype=GRB.BINARY,name="U") #采用0-1变量对分段约束进行线性化处理。即:通过定义0-1变量“U”,
  21. # 设置目标函数
  22. M_PWL.setObjective(quicksum(C[i,j]*x[i,j]*x[i,j] for i in range(N_i) for j in range(N_j)),GRB.MINIMIZE)
  23. # 添加约束
  24. M_PWL.addConstrs((P[i]<=Pmax[i]*U[i] for i in range(N_i)),"Con_P1")
  25. M_PWL.addConstrs((P[i]>=Pmin[i]*U[i] for i in range(N_i)),"Con_P2")
  26. M_PWL.addConstrs((sum(x[i,j] for i in range(N_i))>=demand[j] for j in range(N_j)),"Con_x1")
  27. M_PWL.addConstrs((sum(x[i,j] for j in range(N_j))==P[i] for i in range(N_i)),"Con_x2")
  28. # Optimize model
  29. M_PWL.optimize()
  30. M_PWL.write("PWL1.lp")
  31. P_c=np.zeros(N_i)
  32. x_c=np.zeros((N_i,N_j))
  33. for i in range(N_i):
  34. P_c[i]=P[i].x
  35. for j in range(N_j):
  36. x_c[i,j]=x[i,j].x
  37. print('P_c is',P_c)
  38. print('x_c is',x_c)

2.3 运行结果

  1. Gurobi Optimizer version 9.5.2 build v9.5.2rc0 (win64)
  2. Thread count: 8 physical cores, 16 logical processors, using up to 16 threads
  3. Optimize a model with 13 rows, 18 columns and 39 nonzeros
  4. Model fingerprint: 0xc29ac6cf
  5. Model has 12 quadratic objective terms
  6. Variable types: 15 continuous, 3 integer (3 binary)
  7. Coefficient statistics:
  8. Matrix range [1e+00, 5e+02]
  9. Objective range [0e+00, 0e+00]
  10. QObjective range [2e-02, 2e-01]
  11. Bounds range [1e+00, 5e+02]
  12. RHS range [1e+02, 2e+02]
  13. Presolve removed 7 rows and 6 columns
  14. Presolve time: 0.00s
  15. Presolved: 6 rows, 12 columns, 20 nonzeros
  16. Presolved model has 12 quadratic objective terms
  17. Variable types: 12 continuous, 0 integer (0 binary)
  18. Root relaxation presolve time: 0.00s
  19. Root relaxation presolved: 6 rows, 12 columns, 20 nonzeros
  20. Root relaxation presolved model has 12 quadratic objective terms
  21. Root barrier log...
  22. Ordering time: 0.00s
  23. Barrier statistics:
  24. AA' NZ : 8.000e+00
  25. Factor NZ : 2.100e+01
  26. Factor Ops : 9.100e+01 (less than 1 second per iteration)
  27. Threads : 1
  28. Objective Residual
  29. Iter Primal Dual Primal Dual Compl Time
  30. 0 3.49914954e+06 -3.87161263e+06 3.49e+03 8.83e+02 1.01e+06 0s
  31. 1 1.27698128e+05 -5.62657044e+05 3.70e+02 1.04e+02 1.31e+05 0s
  32. 2 4.07264925e+03 -3.53732018e+05 0.00e+00 1.04e-04 1.19e+04 0s
  33. 3 4.02551091e+03 -4.00359989e+03 0.00e+00 2.04e-06 2.68e+02 0s
  34. 4 2.47855475e+03 -8.87063350e+02 0.00e+00 2.03e-12 1.12e+02 0s
  35. 5 2.20928849e+03 1.86235277e+03 0.00e+00 1.07e-14 1.16e+01 0s
  36. 6 2.16290222e+03 2.16128475e+03 0.00e+00 3.55e-15 5.39e-02 0s
  37. 7 2.16264766e+03 2.16264604e+03 0.00e+00 1.78e-15 5.40e-05 0s
  38. 8 2.16264740e+03 2.16264740e+03 0.00e+00 1.78e-15 5.40e-08 0s
  39. 9 2.16264740e+03 2.16264740e+03 0.00e+00 7.11e-15 5.40e-11 0s
  40. Barrier solved model in 9 iterations and 0.00 seconds (0.00 work units)
  41. Optimal objective 2.16264740e+03
  42. Root relaxation: objective 2.162647e+03, 0 iterations, 0.00 seconds (0.00 work units)
  43. Nodes | Current Node | Objective Bounds | Work
  44. Expl Unexpl | Obj Depth IntInf | Incumbent BestBd Gap | It/Node Time
  45. * 0 0 0 2162.6474028 2162.64740 0.00% - 0s
  46. Explored 1 nodes (0 simplex iterations) in 0.00 seconds (0.00 work units)
  47. Thread count was 16 (of 16 available processors)
  48. Solution count 1: 2162.65
  49. Optimal solution found (tolerance 1.00e-04)
  50. Best objective 2.162647402807e+03, best bound 2.162647402807e+03, gap 0.0000%
  51. P_c is [199.24499511 282.49440796 274.26059693]
  52. x_c is [[ 55.44250871 14.2550231 48.60135551 80.94610778]
  53. [100. 57.28245479 25.21195317 100. ]
  54. [ 61.55749129 78.46252211 71.18669131 63.05389222]]
  55. Process finished with exit code 0

Gurobi Optimizer version 9.5.2 build v9.5.2rc0 (win64)
Thread count: 8 physical cores, 16 logical processors, using up to 16 threads
Optimize a model with 13 rows, 18 columns and 39 nonzeros
Model fingerprint: 0xc29ac6cf
Model has 12 quadratic objective terms
Variable types: 15 continuous, 3 integer (3 binary)
Coefficient statistics:
Matrix range [1e+00, 5e+02]
Objective range [0e+00, 0e+00]
QObjective range [2e-02, 2e-01]
Bounds range [1e+00, 5e+02]
RHS range [1e+02, 2e+02]
Presolve removed 7 rows and 6 columns
Presolve time: 0.00s
Presolved: 6 rows, 12 columns, 20 nonzeros
Presolved model has 12 quadratic objective terms
Variable types: 12 continuous, 0 integer (0 binary)
Root relaxation presolve time: 0.00s
Root relaxation presolved: 6 rows, 12 columns, 20 nonzeros
Root relaxation presolved model has 12 quadratic objective terms
Root barrier log...

Ordering time: 0.00s

Barrier statistics:
AA' NZ : 8.000e+00
Factor NZ : 2.100e+01
Factor Ops : 9.100e+01 (less than 1 second per iteration)
Threads : 1

  1. Objective Residual

Iter Primal Dual Primal Dual Compl Time
0 3.49914954e+06 -3.87161263e+06 3.49e+03 8.83e+02 1.01e+06 0s
1 1.27698128e+05 -5.62657044e+05 3.70e+02 1.04e+02 1.31e+05 0s
2 4.07264925e+03 -3.53732018e+05 0.00e+00 1.04e-04 1.19e+04 0s
3 4.02551091e+03 -4.00359989e+03 0.00e+00 2.04e-06 2.68e+02 0s
4 2.47855475e+03 -8.87063350e+02 0.00e+00 2.03e-12 1.12e+02 0s
5 2.20928849e+03 1.86235277e+03 0.00e+00 1.07e-14 1.16e+01 0s
6 2.16290222e+03 2.16128475e+03 0.00e+00 3.55e-15 5.39e-02 0s
7 2.16264766e+03 2.16264604e+03 0.00e+00 1.78e-15 5.40e-05 0s
8 2.16264740e+03 2.16264740e+03 0.00e+00 1.78e-15 5.40e-08 0s
9 2.16264740e+03 2.16264740e+03 0.00e+00 7.11e-15 5.40e-11 0s

Barrier solved model in 9 iterations and 0.00 seconds (0.00 work units)
Optimal objective 2.16264740e+03

Root relaxation: objective 2.162647e+03, 0 iterations, 0.00 seconds (0.00 work units)

  1. Nodes | Current Node | Objective Bounds | Work

Expl Unexpl | Obj Depth IntInf | Incumbent BestBd Gap | It/Node Time

  • 0 0 0 2162.6474028 2162.64740 0.00% - 0s

Explored 1 nodes (0 simplex iterations) in 0.00 seconds (0.00 work units)
Thread count was 16 (of 16 available processors)

Solution count 1: 2162.65

Optimal solution found (tolerance 1.00e-04)
Best objective 2.162647402807e+03, best bound 2.162647402807e+03, gap 0.0000%
P_c is [199.24499511 282.49440796 274.26059693]
x_c is [[ 55.44250871 14.2550231 48.60135551 80.94610778]
[100. 57.28245479 25.21195317 100. ]
[ 61.55749129 78.46252211 71.18669131 63.05389222]]

Process finished with exit code 0

3 运输问题——包含分段线性约束的优化问题(Python+Cplex实现)

这个博主写得很好:

包含分段线性约束的优化问题(Python+Cplex实现)

  1. import sys
  2. import cplex
  3. from cplex.exceptions import CplexError
  4. def transport(convex):
  5. # 定义工厂供给量
  6. supply = [1000.0, 850.0, 1250.0]
  7. nbSupply = len(supply)
  8. # 定义展厅的需求量向量
  9. demand = [900.0, 1200.0, 600.0, 400.0]
  10. nbDemand = len(demand)
  11. #
  12. n = nbSupply * nbDemand
  13. # 判断成本斜率的凹凸性
  14. if convex:
  15. pwl_slope = [120.0, 80.0, 50.0]
  16. else:
  17. pwl_slope = [30.0, 80.0, 130.0]
  18. def varindex(m, n):
  19. return m * nbDemand + n
  20. # 设置分段线性函数的断点 x 坐标
  21. k = 0
  22. pwl_x = [[0.0] * 4] * n
  23. pwl_y = [[0.0] * 4] * n
  24. for i in range(nbSupply):
  25. for j in range(nbDemand):
  26. if supply[i] < demand[j]:
  27. midval = supply[i]
  28. else:
  29. midval = demand[j]
  30. pwl_x[k][1] = 200.0
  31. pwl_x[k][2] = 400.0
  32. pwl_x[k][3] = midval
  33. pwl_y[k][1] = pwl_x[k][1] * pwl_slope[0]
  34. pwl_y[k][2] = pwl_y[k][1] + \
  35. pwl_slope[1] * (pwl_x[k][2] - pwl_x[k][1])
  36. pwl_y[k][3] = pwl_y[k][2] + \
  37. pwl_slope[2] * (pwl_x[k][3] - pwl_x[k][2])
  38. k = k + 1
  39. # 建立数学模型
  40. model = cplex.Cplex()
  41. model.set_problem_name("transport_py")
  42. model.objective.set_sense(model.objective.sense.minimize)
  43. # 定义决策变量 x_{ij} 表示工厂 i 到展厅 j 的运输量
  44. colname_x = ["x{0}".format(i + 1) for i in range(n)]
  45. model.variables.add(obj=[0.0] * n, lb=[0.0] * n,
  46. ub=[cplex.infinity] * n, names=colname_x)
  47. # y(varindex(i, j)) is used to model the PWL cost associated with
  48. # this shipment.
  49. colname_y = ["y{0}".format(j + 1) for j in range(n)]
  50. model.variables.add(obj=[1.0] * n, lb=[0.0] * n,
  51. ub=[cplex.infinity] * n, names=colname_y)
  52. # 供给必须满足需求约束
  53. for i in range(nbSupply):
  54. ind = [varindex(i, j) for j in range(nbDemand)]
  55. val = [1.0] * nbDemand
  56. row = [[ind, val]]
  57. model.linear_constraints.add(lin_expr=row,
  58. senses="E", rhs=[supply[i]])
  59. # 需求必须符合供给约束
  60. for j in range(nbDemand):
  61. ind = [varindex(i, j) for i in range(nbSupply)]
  62. val = [1.0] * nbSupply
  63. row = [[ind, val]]
  64. model.linear_constraints.add(lin_expr=row,
  65. senses="E", rhs=[demand[j]])
  66. # 添加分段线性约束
  67. for i in range(n):
  68. # preslope is the slope before the first breakpoint. Since the
  69. # first breakpoint is (0, 0) and the lower bound of y is 0, it is
  70. # not meaningful here. To keep things simple, we re-use the
  71. # first item in pwl_slope.
  72. # Similarly, postslope is the slope after the last breakpoint.
  73. # We just use the same slope as in the last segment; we re-use
  74. # the last item in pwl_slope.
  75. model.pwl_constraints.add(vary=n + i,
  76. varx=i,
  77. preslope=pwl_slope[0],
  78. postslope=pwl_slope[-1],
  79. breakx=pwl_x[i],
  80. breaky=pwl_y[i],
  81. name="p{0}".format(i + 1))
  82. # solve model
  83. model.solve()
  84. model.write('transport_py.lp')
  85. # Display solution
  86. print()
  87. print("Solution status :", model.solution.get_status())
  88. print("Cost : {0:.2f}".format(
  89. model.solution.get_objective_value()))
  90. print()
  91. print("Solution values:")
  92. for i in range(nbSupply):
  93. print(" {0}: ".format(i), end='')
  94. for j in range(nbDemand):
  95. print("{0:.2f}\t".format(
  96. model.solution.get_values(varindex(i, j))),
  97. end='')
  98. print()
  99. if __name__ == "__main__":
  100. transport(True)

本文转载自: https://blog.csdn.net/weixin_66436111/article/details/127429729
版权归原作者 wlz249 所有, 如有侵权,请联系我们删除。

“分段模型线性化(PWL)【Python|Gurobi实现】”的评论:

还没有评论