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

# -*- coding: utf-8 -*-
from gurobipy import *
import numpy as np
N_i=3
N_j=4
k=1.8
C=np.array([
    [0.0755,0.0655,0.0498,0.0585],
    [0.0276,0.0163,0.0960,0.0224],
    [0.0680,0.0119,0.0340,0.0751]
    ])
Pmin=np.array([100,50,30])
Pmax=np.array([450,350,500])
demand=np.array([217,150,145,244])
# 创建模型
M_PWL=Model("PWL-1")
# 变量声明
x  =M_PWL.addVars(N_i,N_j,lb=0,ub=100, name="x")
P  =M_PWL.addVars(N_i,lb=0,ub=Pmax, name="P")
U  =M_PWL.addVars(N_i,vtype=GRB.BINARY,name="U")  #采用0-1变量对分段约束进行线性化处理。即：通过定义0-1变量“U”，
# 设置目标函数
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)
# 添加约束
M_PWL.addConstrs((P[i]<=Pmax[i]*U[i] for i in range(N_i)),"Con_P1")
M_PWL.addConstrs((P[i]>=Pmin[i]*U[i] for i in range(N_i)),"Con_P2")
M_PWL.addConstrs((sum(x[i,j] for i in range(N_i))>=demand[j] for j in range(N_j)),"Con_x1")
M_PWL.addConstrs((sum(x[i,j] for j in range(N_j))==P[i] for i in range(N_i)),"Con_x2")
# Optimize model
M_PWL.optimize()
M_PWL.write("PWL1.lp")
P_c=np.zeros(N_i)
x_c=np.zeros((N_i,N_j))
for i in range(N_i):
    P_c[i]=P[i].x
    for j in range(N_j):
        x_c[i,j]=x[i,j].x
print('P_c is',P_c)
print('x_c is',x_c)

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
                  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)
    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 

import sys
import cplex
from cplex.exceptions import CplexError
def transport(convex):
    # 定义工厂供给量
    supply = [1000.0, 850.0, 1250.0]
    nbSupply = len(supply)
    # 定义展厅的需求量向量
    demand = [900.0, 1200.0, 600.0, 400.0]
    nbDemand = len(demand)
    # 
    n = nbSupply * nbDemand
    # 判断成本斜率的凹凸性
    if convex:
        pwl_slope = [120.0, 80.0, 50.0]
    else:
        pwl_slope = [30.0, 80.0, 130.0]
    def varindex(m, n):
        return m * nbDemand + n
    # 设置分段线性函数的断点 x 坐标
    k = 0
    pwl_x = [[0.0] * 4] * n
    pwl_y = [[0.0] * 4] * n
    for i in range(nbSupply):
        for j in range(nbDemand):
            if supply[i] < demand[j]:
                midval = supply[i]
            else:
                midval = demand[j]
            pwl_x[k][1] = 200.0
            pwl_x[k][2] = 400.0
            pwl_x[k][3] = midval
            pwl_y[k][1] = pwl_x[k][1] * pwl_slope[0]
            pwl_y[k][2] = pwl_y[k][1] + \
                pwl_slope[1] * (pwl_x[k][2] - pwl_x[k][1])
            pwl_y[k][3] = pwl_y[k][2] + \
                pwl_slope[2] * (pwl_x[k][3] - pwl_x[k][2])
            k = k + 1
    # 建立数学模型
    model = cplex.Cplex()
    model.set_problem_name("transport_py")
    model.objective.set_sense(model.objective.sense.minimize)
    # 定义决策变量 x_{ij} 表示工厂 i 到展厅 j 的运输量
    colname_x = ["x{0}".format(i + 1) for i in range(n)]
    model.variables.add(obj=[0.0] * n, lb=[0.0] * n,
                        ub=[cplex.infinity] * n, names=colname_x)
    # y(varindex(i, j)) is used to model the PWL cost associated with
    # this shipment.
    colname_y = ["y{0}".format(j + 1) for j in range(n)]
    model.variables.add(obj=[1.0] * n, lb=[0.0] * n,
                        ub=[cplex.infinity] * n, names=colname_y)
    # 供给必须满足需求约束
    for i in range(nbSupply):
        ind = [varindex(i, j) for j in range(nbDemand)]
        val = [1.0] * nbDemand
        row = [[ind, val]]
        model.linear_constraints.add(lin_expr=row,
                                     senses="E", rhs=[supply[i]])
    # 需求必须符合供给约束
    for j in range(nbDemand):
        ind = [varindex(i, j) for i in range(nbSupply)]
        val = [1.0] * nbSupply
        row = [[ind, val]]
        model.linear_constraints.add(lin_expr=row,
                                     senses="E", rhs=[demand[j]])
    # 添加分段线性约束
    for i in range(n):
        # preslope is the slope before the first breakpoint.  Since the
        # first breakpoint is (0, 0) and the lower bound of y is 0, it is
        # not meaningful here.  To keep things simple, we re-use the
        # first item in pwl_slope.
        # Similarly, postslope is the slope after the last breakpoint.
        # We just use the same slope as in the last segment; we re-use
        # the last item in pwl_slope.
        model.pwl_constraints.add(vary=n + i,
                                  varx=i,
                                  preslope=pwl_slope[0],
                                  postslope=pwl_slope[-1],
                                  breakx=pwl_x[i],
                                  breaky=pwl_y[i],
                                  name="p{0}".format(i + 1))
    # solve model
    model.solve()
    model.write('transport_py.lp')
    # Display solution
    print()
    print("Solution status :", model.solution.get_status())
    print("Cost            : {0:.2f}".format(
        model.solution.get_objective_value()))
    print()
    print("Solution values:")
    for i in range(nbSupply):
        print("   {0}: ".format(i), end='')
        for j in range(nbDemand):
            print("{0:.2f}\t".format(
                model.solution.get_values(varindex(i, j))),
                end='')
        print()
if __name__ == "__main__":
    transport(True)