矩阵建模功能

杉数求解器COPT的Python API提供矩阵建模方式,支持 NumPy 的多维数组,二维的 NumPy 矩阵,以及 SciPy 列压缩矩阵( csc_matrix ) 和行压缩矩阵( csr_matrix )的运算,并且可以和常规(标量)变量和常规约束结合。 NumPy 版本要求需要为1.23及以上, 并且不高于2.0版本( NumPy 2.0 引入了不兼容的改动)。 Python 最低版本要求为3.8。主要提供如下功能:

  1. 添加多维变量( MVar 类对象)等相关操作;

  2. 构建多维线性表达式( MLinExpr 类对象),添加多维线性约束( MConstr 类对象)等相关操作;

  3. 构建多维二次表达式( MQuadExpr 类对象),添加多维凸二次约束( QConstraint 类对象)等相关操作。

两种矩阵建模模式

COPT提供的矩阵建模功能目前支持两种模式:原始模式和内测模式。其中,原始模式(为默认模式)依赖于Python的 NumPy 库。 7.1.4新增内测模式,基于COPT C++接口的原生矩阵建模类实现,建模速度相较前者有一定程度提升。

在Python接口中,两种模式通过 Model.matrixmodelmode 进行控制和切换,默认是原始模式。请根据使用需要设置:

  • Model.matrixmodelmode = "legacy" (默认模式),依赖于 NumPy 库;

  • Model.matrixmodelmode = "experimental" ,基于COPT的C++原生矩阵类( NdArray ),无外部依赖。

两种模式支持的函数以及操作类型略有差异。

1.MConstr.reshape()/MVar.reshape() 的 order(存储)方式,即函数参量 order 的取值:

  • 原始模式:同时支持 'C' (按行)和 'F' (按列)

  • 内测模式:只支持 'C' (按行)

2.高级索引(advanced indexing)支持情况:

以nqueen问题中,对角线上最多只能有一个queen对应的约束为例:

  • 原始模式:能够兼容 NumPy 的高级索引,可具体实现如下:

    import coptpy as cp
    import numpy as np
    
    env = cp.Envr()
    model = env.createModel()
    model.matrixmodelmode = "legacy"
    for i in range(1, 2*n):
        # At most one queen per diagonal
        diagn = (range(max(0, i-n), min(n, i)), range(min(n, i)-1, max(0, i-n)-1, -1))
        model.addConstrs(x[diagn].sum() <= 1, nameprefix="diag"+str(i))
        # At most one queen per anti-diagonal
        adiagn = (range(max(0, i-n), min(n, i)), range(max(0, n-i), min(n, 2*n-i)))
        model.addConstrs(x[adiagn].sum() <= 1, nameprefix="adiag"+str(i))
    
  • 内测模式:不支持上述高级索引方式,需要通过Python原生的索引方式,可具体实现如下:

    import coptpy as cp
    import numpy as np
    
    env = cp.Envr()
    model = env.createModel()
    model.matrixmodelmode = "experimental"
    for i i in range(-n+1, n):
        # At most one queen per diagonal
        model.addConstrs(x.diagonal(i).sum() <= 1, nameprefix="diag"+str(i))
        # At most one queen per anti-diagonal
        model.addConstrs(x[:, ::-1].diagonal(i).sum() <= 1, nameprefix="adiag"+str(i))
    

3.获取矩阵类对象( MConstrMVar 等对象)相关信息,返回值类型不同:

  • 原始模式: numpy.ndarray

  • 内测模式: coptcore.NdArray

4.支持的维度情况有所不同:

  • 原始模式:对维度无限制,可支持更高维矩阵。

  • 内测模式:支持的维度最高为三维。

5.添加多维二次约束时,对参与构成约束的变量维度支持情况,以及返回值的类型有所不同:

  • 原始模式:参与构成二次约束的变量维度需要为1(即向量),不支持更高的维度,返回值类型为 QConstraint类 对象。

    import coptpy as cp
    import numpy as np
    
    env = cp.Envr()
    model = env.createModel()
    model.matrixmodelmode = "legacy"
    Q = np.full((3, 3), 1)
    mx = model.addMVar(3, nameprefix="mx")
    # mqc <coptpy.QConstraint: >
    mqc = model.addQConstr(mx@Q@mx<=1)
    
  • 内测模式:参与构成二次约束的变量维度可以为一维或者二维,返回值类型为 MQConstr 类 对象。

    import coptpy as cp
    import numpy as np
    
    env = cp.Envr()
    model = env.createModel()
    model.matrixmodelmode = "experimental"
    Q = np.full((3, 3), 1)
    mx = model.addMVar((3, 3), nameprefix="mx")
    # mqc <coptpy.MQConstr: shape=(3, 3)>
    mqc = model.addQConstr(mx@Q@mx <= 1)
    

多维变量

  1. 添加多维变量 MVar

MVar 是多维变量相关操作的封装,用户可以调用 Model.addMVar() 向模型中添加任意维度和形状的多维变量 MVar 。除了需要指定参数 shape 形状之外,其余参数与添加普通变量时需指定的一致,包括: lb , ub , vtype , nameprefix

  • 添加一维连续型多维变量: x = Model.addMVar(3)

  • 添加二维 3x3 二进制型多维变量: y = Model.addMVar((3,3), vtype=COPT.BINARY)

此外,还可以对 MVar 多维变量进行切片操作,如: y1 = y[:,0:2]

  1. 获取多维变量相关属性

    • 维度数量:MVar.ndim

    • 多维变量的形状:MVar.shape

    • 元素数量:MVar.size

多维数组运算及表达式

多维线性表达式

多维变量与其系数(可以为 ndarray )构成多维线性表达式( MLinExpr ),支持的运算主要有:

  1. 矩阵乘法运算:A @ x

x = model.addMVar(3)
A = np.array([[1, 0, 1],[0, 0, 1]])
expr1 = A @ x
  1. 向量内积运算

x = model.addMVar(3)
c = np.array([1, 2, 3])
expr2 = c @ x

多维二次表达式

常见的多维二次表达式及对应的数学形式如下:

  • x @ Q @ x:\(x^TQx\)

  • x @ x:\(x^Tx\)

  • x @ Q @ x + c @ x + b:\(x^TQx+c^Tx+b\)

其他多维数组运算

  1. 与常规变量、常规线性表达式,以及常数结合:

x = model.addMVar(3)
y = model.addVar()
c = np.array([1, 2, 3])
Q = np.full((3, 3), 1)
expr3 = 2 * x @ Q @ x + c @ x + 2 * y + 1
  1. 自加/自减/自乘运算:

mx = m.addMVar((3, 3))
B = np.array([[1, 0, 1], [0, 1, 1]])
expr_add = B @ mx
expr_add += 1
expr_add *= 2

注意

  • 我们直接打印多维表达式 print(MLinExpr)/print(MQuadExpr) 时,同时会输出表达式的形状 shape 。当 shape=() 时,表示该表达式为标量(单一线性/二次表达式),对应的 ndim=0size=1 ,对于多维变量 MVarshape 也是如此;

  • 当进行矩阵乘法运算(A@x)时, 需要满足矩阵乘法运算法则,A的列数和X的行数需要相同;

  • COPT支持 MLinExpr 类对象和 LinExpr 结合使用,但需注意此时的 MLinExpr 类对象需要 shape=() ,最终返回的表达式为 MLinExpr 类对象且 shape=()

矩阵约束

矩阵线性约束

COPT支持添加矩阵线性约束的两种方式,函数参数提供的格式有所不同:

  1. 专门添加矩阵线性约束的 Model.addMConstr() ,可以指定的参数有:

    • A :线性约束的系数矩阵

    • x :决策变量( MVar 类对象)

    • sense :线性约束的类型,可取值有: 'L' (小于等于)、 'G' (大于等于)、 'E' (等于)等

    • b :线性约束右端项(向量,维度与矩阵 A 的行数一致)

    • name :约束名称

x = model.addMVar(shape=3, vtype=COPT.BINARY, nameprefix='x')
A = np.array([[1, 2, 3], [3, 2, 1]])
b = np.array([2, 5])
mconstrs = model.addMConstr(A, x, 'L', b, nameprefix='c')
obj = np.array([1, 2, 1])
model.setObjective(obj @ x, COPT.MINIMIZE)
  1. 多维线性约束可视为一组线性约束,故 Model.addConstrs() 也可以添加多维线性约束:

x = model.addMVar(shape=3, vtype=COPT.BINARY, nameprefix='x')
A = np.array([[1, 2, 3], [3, 2, 1]])
b = np.array([2, 5])
mconstrs = model.addConstrs(A @ x <= b, nameprefix='c')
obj = np.array([1, 2, 1])
model.setObjective(obj @ x, COPT.MINIMIZE)

二次约束

COPT支持添加矩阵二次约束的两种方式,函数参数提供的格式有所不同:

  1. 专门添加矩阵二次约束的 Model.addMQConstr() ,可以指定的参数有:

    • Q :二次项系数矩阵

    • c :线性项系数向量,若无线性项,则为 None

    • sense :二次约束的类型,可取值有: 'L' (小于等于)、 'G' (大于等于)、 'E' (等于)等。

    • rhs :二次约束右端项

    • xQ_L :二次项系数矩阵 Q 左端变量(向量,长度与矩阵 Q 的行数一致)

    • xQ_R :二次项系数矩阵 Q 右端变量(向量,长度与矩阵 Q 的列数一致)

    • xc :线性项的变量,若无线性项,则为 None

    • name : 约束名称

Q = np.diag([3, 2, 1])
x = model.addMVar(3)
c1 = model.addMQConstr(Q, None, 'L', 1.0, x, x)
  1. Model.addQConstr() ,直接给出多维二次表达式

  • lhs : 多维二次表达式

  • sense :约束类型

  • rhs :约束右端项

Q = np.diag([3, 2, 1])
x = model.addMVar(3)
c2 = model.addQConstr(x@Q@x<=1.0)

由多维变量构成的目标函数

COPT支持设置线性和二次目标函数,并提供设置由多维变量构成目标函数的两种方式,函数参数的格式有所不同:

  1. 专门设置由多维变量构成的目标函数的 Model.setMObjective() ,可以指定的参数有:

    • Q :二次项系数矩阵,若目标函数是线性的,则为 None

    • c :线性项系数向量,若无线性项,则为 None

    • constant :目标函数的常数项

    • xQ_L :二次项系数矩阵 Q 左端变量(向量,长度与矩阵 Q 的行数一致),若目标函数是线性的,则为 None

    • xQ_R :二次项系数矩阵 Q 右端变量(向量,长度与矩阵 Q 的列数一致),若目标函数是线性的,则为 None

    • xc :线性项的变量,若无线性项,则为 None

    • sense :优化方向,可取值为: COPT.MINIMIZECOPT.MAXIMIZE

  2. 直接给出目标函数表达式的 Model.setObjective()

    • expr :目标函数表达式,可以为线性或者二次表达式

    • sense :优化方向,可取值为: COPT.MINIMIZECOPT.MINIMIZE

x = model.addMVar(shape=3, vtype=COPT.BINARY, nameprefix="x")
obj = np.array([1, 2, 1])
model.setObjective(obj @ x, COPT.MINIMIZE)

关于矩阵建模方式,COPT Python接口分别提供多维变量及其(线性和凸二次)表达式、矩阵约束的类,对相关操作进行封装,其中包含的成员方法及具体介绍请参考 Python API 对应部分: