【数学建模学习手册】第三章：规划模型（二）

本专栏内容为：数学建模原理 记录学习数学建模

💓博主csdn个人主页：小小unicorn
⏩专栏分类：数学建模
🚚代码仓库：小小unicorn的代码仓库🚚
🌹🌹🌹关注我带你学习编程知识

非线性规划

如果目标函数或约束条件中至少有一个是非线性函数时的最优化问题就叫做非线性规划问题。

二次规划

二次规划是非线性规划中的一类特殊规划问题，如果目标函数是二次函数，约束函数是线性函数时，这就是一个二次规划问题。一个有n个变量与m个限制的二次规划问题可以用以下的形式描述。首先给定：


那么，当维数增大时，是否具有类似的性质呢？再此不做证明，直接给出结论：

• 如果𝐻是半正定矩阵，那么𝑓(𝑥)是一个凸函数。相应的二次规划为凸二次规划问题；此时如果有局部最优解，那这个局部最优解就是全局最优解。但这个全局最小值可能是不唯一的，如下图所示。
  
• 若𝐻为非正定矩阵，则目标函数是有多个驻点和局部极小点的NP难问题。
• 如果𝐻=0，二次规划问题就变成线性规划问题。

下面我们想一下为什么要引入二次规划，以及为什么二次规划问题被写成了上面那种标准形式。首先要明确两个概念，维度是自变量的个数，次数是自变量最高的阶数，两个不能弄混。我们解决问题的思路是将不熟悉的问题，转换成熟悉的问题。其中优化轨迹点是不熟悉的问题，但是求解多项式的极值是熟悉的问题。因此我们希望将复杂问题转换为一个与多项式求极值类似的问题。而求解多项式极值也要看多项式的次数，如果是一次，那就是直线，没什么好讨论的。如果是二次，那函数将有唯一的极值，也是最值，这是我们最想看到的。如果超过2次，那问题就复杂了。

多项式求极值是二维空间的问题，只有𝑥和𝑦两个变量。而运动规划的优化问题是高维的。高维空间的优化问题是十分复杂的，我们如果将其转换为一个类似于我们熟悉的多项式求极值问题，我们当然希望选择对优化来说最理想的二次。如果问题的维度也高，次数也高，就很难求解了。于是我们希望将问题转化为这样一个高维空间中，类似于二次函数求极值的问题，所以目标函数这样定义：

非线性规划的求解

类比线性规划的基本形式，我们把所有的线性方程约束中系数做成系数矩阵𝐴_eq，等号右边的常数作为列向量𝑏_eq；线性不等式约束中的系数矩阵𝐴和不等号右边的常数b，为了方便起见通常将不等式统一为小于等于；变量𝑥在向量𝑙_𝑏到𝑢_𝑏之间取值；非线性不等式和非线性方程分别用𝐶(𝑥)和𝐶_eq(𝑥)，则线性规划的标准形式如下所示：

我们仍然将问题统一为函数极小值问题，不等约束统一为小于等于。如果原问题是最大值或者有大于等于，那就乘−1进行取反即可。

注意：非线性规划中只要破坏了等式约束、不等约束和目标函数当中任何一个的线性就可以说是一个非线性规划。仅仅破坏一条问题的求解难度会上来很多。

从线性到非线性，多元函数可能初来乍到的同学并不一定理解。在高中我们说，函数是从一个集合（定义域）到另一个集合（值域）的一一对应的映射，那现在多元函数？不就变成多个集合到一个集合的映射了嘛？这，怎么也能叫一一对应呢？诶，这个时候你就注意了，多元函数仍然是一个集合到另一个集合的映射，只不过自变量的集合不是数集，而是点集，或者说是𝑛维空间里面向量的集合。

如果高中学过导数的同学可能知道，一元函数求极值的一个方法就是求导导数为0。对于多元函数也类似，如果是无约束就只是给了一个多元函数求极值，只需要针对每个变量分别求导数（我们把这个过程叫做求偏导），所有偏导为0解方程组得到的就是极值点的候选解。当然，我们有可能解出来的是𝑥₃=0的这种极值解，这种情况下我们会通过二阶导数的符号去进一步判定。

注意：在求对某一个变量的偏导数的时候通常把其他变量视作常数，我们把这种策略叫主元策略。

当我们碰上了有约束条件下的非线性函数求极值的时候，我们通常使用 Lagrange 法。例如，对于广义的含等式条件（先暂时只考虑等式问题）的极值问题：

我们通过引入𝑛个称之为Lagrange 乘子的常数把原问题改写为新的函数：

接下来就像解无约束极值一样对每个𝑥和乘子求偏导即可解决问题。而当我们在问题中考虑不等条件，那么这个问题的求解策略其实类似，我们称其为KKT条件。对于问题

我们分别引入两个不同乘子，函数𝐿将在𝑥取得极值当且仅当它在下面的问题中取得极值：

当然，如果想求的不是一个精确解而是一个近似的数值解，那么我们的方法同样有很多。除了在下一章中会介绍的一些数值方法外，Monte Carlo模拟几乎是用的最广泛的一种。

整数规划与指派问题

整数规划概念

离散和连续是一对重要的概念。这一对概念其实非常好理解：比如说，连续的变量取值是连续的实数，取值可以是任意的浮点数（小数），但离散变量取值是不连续的，是有间隔的。离散量的取值往往是整数，也可能是有限的取值。

离散优化的例子其实也很简单，如果把前面的线性规划或者非线性规划当中加上一条约束：自变量取整数，这个问题就开始有意思了。可能最优解是一个全浮点数，但加上整数约束以后究竟在哪个整数点上取到最优解那还真说不好，你如果用枚举的方式去解那复杂度是成倍上涨。我们想，一定有更快的求解策略。

另一类离散优化的典型例子是匹配问题和组合优化问题，比如有100个人匹配100项任务，百配百就有五千种匹配模式。而这么多的匹配模式中究竟哪一个是最优解呢？那就得在五千种匹配模式中搜索，变量就是某一个人是否匹配某一项任务，取值只能是{0,1}。所以这就是一种离散优化。下面我们也会对这类问题进行介绍。

分支界定法

同单纯形法与Monte Carlo法之于线性规划，解整数规划的基础原理其实更多。最典型的两种算法就是分支定界法和割平面法。

分支定界法是一种经典的搜索算法，这里把它用在规划当中主要是为了对上下界进行搜索。

分支定界本质上是构造一棵搜索树进行上下界搜索，它会把问题的搜索空间组织成一棵树，而分支就是从根出发将原始问题按照整数约束去分支为左子树和右子树，通过不断检查子树的上下界去搜索最优解的过程。

我们举一个例子：
例题：求规划的最优解

首先我们忽略取值{0,1}这个条件，把它就当做一个取值范围[0,1]之间的线性规划去做，它肯定是有最优解的，这毋庸置疑。现在，我们从x₃开始分支，分别按取0和取1去对原始问题进行划分。

如果x_3​取值为1，那么原始问题变成了 7𝑥₁+8𝑥₂⩽7 的条件下最大𝑥₁+ 𝑥₂+1；如果取值为0那么原始问题变成了 7𝑥₁+8𝑥₂⩽14 的条件下最大化 𝑥₁+𝑥₂ 。然后我们分别计算两边的最优解，两边都可能存在最优整数解于是对这两种情况再进行划分。每划分一次我们就会对同一层的子问题求解对应的线性规划（把整数条件换成区间条件）观察谁最小谁可分，到最后遍历完成就得到了最优解。如图所示：

注意：我们每经过一层就会更新一个线性规划的最优解（也可以叫松弛解），在根节点我们将它设置为负无穷，而在子节点中只要能够比上一次的最优解更优我们就会更新这个最优解。比如说在第三层，节点4的最优解2就是最优，而节点6的最优解1和节点7的最优解$1.8571$再怎么分不会比节点4更优，所以我们下一步只对节点4分支并按照这一分支为子问题定界。

import math
from scipy.optimize import linprog
import sys

def integerPro(c, A, b, Aeq, beq, t=1.0E-8):
    res = linprog(c, A_ub=A, b_ub=b, A_eq=Aeq, b_eq=beq)
    bestVal = sys.maxsize  # 很大一个数
    bestX = res.x
    if not (type(res.x) is float or res.status != 0):
        bestVal = sum([x * y for x, y in zip(c, bestX)])
    if all(((x - math.floor(x)) <= t or (math.ceil(x) - x) <= t) for x in bestX):
        return bestVal, bestX
    else:
        ind = [i for i, x in enumerate(bestX) if (x - math.floor(x)) > t and (math.ceil(x) - x) > t][0]
        newCon1 = [0] * len(A[0])
        newCon2 = [0] * len(A[0])
        newCon1[ind] = -1
        newCon2[ind] = 1
        newA1 = A.copy()
        newA2 = A.copy()
        newA1.append(newCon1)
        newA2.append(newCon2)
        newB1 = b.copy()
        newB2 = b.copy()
        newB1.append(-math.ceil(bestX[ind]))
        newB2.append(math.floor(bestX[ind]))
        r1 = integerPro(c, newA1, newB1, Aeq, beq)
        r2 = integerPro(c, newA2, newB2, Aeq, beq)
        if r1[0] < r2[0]:
            return r1
        else:
            return r2
if __name__ == '__main__':
    c = [3, 4, 1]
    A = [[-1, -6, -2], [-2, 0, 0]]
    b = [-5, -3]
    Aeq = [[0, 0, 0]]
    beq = [0]
    print(integerPro(c, A, b, Aeq, beq))

指派问题与匈牙利法

0-1规划是整数规划中最特殊的一种。它的限制不仅仅是要求变量是整数，而且只能是0或1，故而名曰0-1规划。事实上，在数学建模竞赛当中，0-1规划可以说是最常见的整数规划，是学习的重点。指派问题又是怎么一回事呢？

我们可以看到这样一个例子：
例题：现在有4个人{A,B,C,D}可以做四项工作{1,2,3,4}，他们每个人只能做一项工作，所需要的时间按照表2.3给出：

我们把四个人与4项工作对应的4⋅4=16项安排作为决策变量，变量取值为{0,1}，表示某个人是否执行某项工作。一个人只能做一项工作，一项工作也只能一个人来做，比如对𝐴而言，𝐴如果做了2就不能做1,3,4，所以𝐴匹配的四个变量只能有一个是1其余三个都是0。同样地，对任务2而言如果它被安排给了𝐴那么B,C,D也都不能去做，所以任务2匹配的四个变量也只能一个是1其余三个是0。由此，我们给出定义：

指派问题和0-1规划的解法可以从matlab的整数规划函数中进行约束，但还有一种经典的算法可以解0-1规划。这个方法被称作匈牙利法。匈牙利法的操作比较有趣，对于时间排布表，首先将每行减去当前行的最小值，然后将每一列的值减去当前列的最小值。接下来一步比较有趣，需要用最少的水平线和竖直线覆盖所有的0项。如果线条总数为4那么算法停止，给出指派方案；如果少于4条那么则计算没有被覆盖的最小值，将没有被覆盖的每行减去最小值，被覆盖的每列加上最小值，然后重新进行覆盖。整个过程如图所示：

from scipy.optimize import linear_sum_assignment
import numpy as np
T=np.array([[25,29,31,42],[39,38,26,20],[34,27,28,40],[24,42,36,23]])
row_ind,col_ind=linear_sum_assignment(T)
print(row_ind)
print(col_ind)
print(T[row_ind,col_ind])
print(T[row_ind,col_ind].sum())

使用Scipy和Cvxpy解决规划问题

使用scipy求解函数优化问题

在 scipy.optimize 模块中，提供了多种用于非线性规划问题的方法，适用于不同类型的问题：

• brent(): 适用于单变量无约束优化问题，结合了Newton 法和二分法。

• fmin(): 适用于多变量无约束优化问题，采用单纯形法，只需利用函数值，无需函数的导数或二阶导数。

• leastsq(): 用于解决非线性最小二乘问题，用于求解非线性最小二乘拟合问题。

• minimize(): 适用于约束优化问题，利用 Lagrange 乘子法将约束优化问题转化为无约束优化问题。

• minimize(fun, x_0[, args, method, jac, hess, ...])：用于对一个或多个变量的标量函数进行最小化。

对于无约束问题优化算法：

• method='CG': 非线性共轭梯度算法，只能处理无约束优化问题，需要使用一阶导数函数。
• method='BFGS': 拟 Newton 法，只能处理无约束优化问题，需要使用一阶导数函数。BFGS 算法性能良好，是无约束优化问题的默认算法。
• method='Newton-CG': 截断 Newton 法，只能处理无约束优化问题，需要使用一阶导数函数，适合处理大规模问题。
• method='dogleg': Dog-leg 信赖域算法，需要使用梯度和 Hessian 矩阵（必须正定），只能处理无约束优化问题。
• method='trust-ncg': 采用 Newton 共轭梯度信赖域算法，需要使用梯度和 Hessian 矩阵（必须正定），只能处理无约束优化问题，适合大规模问题。
• method='trust-exact': 求解无约束极小化问题的信赖域方法，需要梯度和 Hessian 矩阵（不需要正定）。
• method='trust-krylov': 使用 Newton-GLTR 信赖域算法，需要使用梯度和 Hessian 矩阵（必须正定），只能处理无约束优化问题，适合中大规模问题。

对于边界约束条件问题优化算法：

• method='Nelder-Mead': 下山单纯形法，可以处理边界约束条件（决策变量的上下限），只使用目标函数，不使用导数函数或二阶导数，具有较强的鲁棒性。

• method='L-BFGS-B': 改进的 BFGS 拟Newton 法，“L” 指有限内存，“B” 指边界约束，可以处理边界约束条件，需要使用一阶导数函数。L-BFGS-B 算法性能良好，内存消耗小，适合处理大规模问题，是边界约束优化问题的默认算法。

• method='Powell': 改进的共轭方向法，可以处理边界约束条件（决策变量的上下限）。

• method='TNC': 截断Newton 法，可以处理边界约束条件。
  对于带有约束条件问题优化算法：

• method='COBYLA': 线性近似约束优化方法，通过线性逼近目标函数和约束条件来处理非线性问题。只使用目标函数，不需要导数或二阶导数值，可以处理约束条件。

• method='SLSQP': 序贯最小二乘规划算法，可以处理边界约束、等式约束和不等式约束条件。SLSQP 算法性能良好，是带有约束条件优化问题的默认算法。

• method='trust-constr': 信赖域算法，通用的约束最优化方法，适合处理大规模问题。

下面给出了一些利用 scipy.optimize 模块进行求解的代码示例：

from scipy.optimize import brent, fmin, fmin_ncg, minimize
import numpy as np

# 1. Demo1：单变量无约束优化问题(Scipy.optimize.brent)

def objf(x):  # 目标函数
    fx = x**2 - 8*np.sin(2*x+np.pi)
    return fx
xIni = -5.0
xOpt= brent(objf, brack=(xIni,2))
print("xIni={:.4f}\tfxIni={:.4f}".format(xIni,objf(xIni)))
print("xOpt={:.4f}\tfxOpt={:.4f}".format(xOpt,objf(xOpt)))

def objf2(x):  # Rosenbrock benchmark function
    fx = 100.0 * (x[0] - x[1] ** 2.0) ** 2.0 + (1 - x[1]) ** 2.0
    return fx

xIni = np.array([-2, -2])
xOpt = fmin(objf2, xIni)
print("xIni={:.4f},{:.4f}\tfxIni={:.4f}".format(xIni[0],xIni[1],objf2(xIni)))
print("xOpt={:.4f},{:.4f}\tfxOpt={:.4f}".format(xOpt[0],xOpt[1],objf2(xOpt)))

#目标函数：
def func(args):
    fun = lambda x: 60 - 10*x[0] - 4*x[1] + x[0]**2 + x[1]**2 - x[0]*x[1] 
    return fun

#约束条件，包括等式约束和不等式约束
def con(args):
    cons = ({'type': 'eq', 'fun': lambda x: x[0]+x[1]-8})
    return cons
# 解决问题是：当x_1+x_2=8时，求解函数60-10x_1-4x_2+x_1^2+x_2^2-x_1x_2的极小值
if __name__ == "__main__":
    args = ()
    args1 = ()
    cons = con(args1)
    x_0 = np.array((2.0, 1.0)) #设置初始值，初始值的设置很重要，很容易收敛到另外的极值点中，建议多试几个值
    #求解#
    res = minimize(func(args), x_0, method='SLSQP', constraints=cons)
    print(res)
def objF4(x):  # 定义目标函数
    a, b, c, d = 1, 2, 3, 8
    fx = a*x[0]**2 + b*x[1]**2 + c*x[2]**2 + d
    return fx

# 定义约束条件函数
def constraint1(x):  # 不等式约束 f(x)>=0
    return x[0]** 2 - x[1] + x[2]**2
def constraint2(x):  # 不等式约束 转换为标准形式
    return -(x[0] + x[1]**2 + x[2]**3 - 20)
def constraint3(x):  # 等式约束
    return -x[0] - x[1]**2 + 2
def constraint4(x):  # 等式约束
    return x[1] + 2*x[2]**2 -3
# 定义边界约束
b = (0.0, None)
bnds = (b, b, b)
# 定义约束条件
con1 = {'type': 'ineq', 'fun': constraint1}
con2 = {'type': 'ineq', 'fun': constraint2}
con3 = {'type': 'eq', 'fun': constraint3}
con4 = {'type': 'eq', 'fun': constraint4}
cons = ([con1, con2, con3,con4])  # 3个约束条件
# 求解优化问题
x_0 = np.array([1., 2., 3.])  # 定义搜索的初值
res = minimize(objF4, x_0, method='SLSQP', bounds=bnds, constraints=cons)
print("Optimization problem (res):\t{}".format(res.message))  # 优化是否成功
print("xOpt = {}".format(res.x))  # 自变量的优化值
print("min f(x) = {:.4f}".format(res.fun))  # 目标函数的优化值
def objF5(x,args):  # 定义目标函数
    a,b,c,d = args
    fx = lambda x: a*x[0]**2 + b*x[1]**2 + c*x[2]**2 + d
    return a*x[0]**2 + b*x[1]**2 + c*x[2]**2 + d

def constraint1():  # 定义约束条件函数
    cons = ({'type': 'ineq', 'fun': lambda x: (x[0]**2 - x[1] + x[2]**2)},  # 不等式约束 f(x)>=0
            {'type': 'ineq', 'fun': lambda x: -(x[0] + x[1]**2 + x[2]**3 - 20)},  # 不等式约束 转换为标准形式
            {'type': 'eq', 'fun': lambda x: (-x[0] - x[1]**2 + 2)},  # 等式约束
            {'type': 'eq', 'fun': lambda x: (x[1] + 2*x[2]**2 - 3)})  # 等式约束
    return cons

# 定义边界约束
b = (0.0, None)
bnds = (b, b, b)
# 定义约束条件
cons = constraint1()
args1 = (1,2,3,8)  # 定义目标函数中的参数
# 求解优化问题
x_0 = np.array([1., 2., 3.])  # 定义搜索的初值
res1 = minimize(objF5, x_0=x_0, args=[1,2,3,8], method='SLSQP', bounds=bnds, constraints=cons)
print("Optimization problem (res1):\t{}".format(res1.message))  # 优化是否成功
print("xOpt = {}".format(res1.x))  # 自变量的优化值
print("min f(x) = {:.4f}".format(res1.fun))  # 目标函数的优化值

使用cvxpy求解函数优化问题

CVXPY 是一种用于凸优化问题的 Python 嵌入式建模语言。它允许您以遵循数学的自然方式表达您的问题，而不是求解器要求的限制性标准形式

CVXPY 的变量类型

变量可以是标量、向量以及矩阵 cvxpy 中可以做常数使用的有：

• NumPy ndarrays
• NumPy matrices
• SciPy sparse matrices

CVXPY 的约束可以使用==, <=，>= ，不能使用< ，>。也不能使用0 <= x <= 1 or x == y == 2。 parameters可以理解为参数求解问题里的一个常数，可以是标量、向量、矩阵。在没有求解问题前（调用类xxx.solve()），其允许你改变其值。

例1

解决这个问题的代码如下:

#例题1
import cvxpy as cp
import numpy as np

coef = np.array([72, 64])#输入目标函数系数
left = np.array([[1, 1], [12, 8], [3, 0]])#输入约束条件系数
right = np.array([50, 480, 100])#输入约束条件上限值
x = cp.Variable(2)#构造决策变量
obj = cp.Maximize(coef @ x)#构造目标函数
cons = [x >= 0, left @ x <= right]#构造约束条件
prob = cp.Problem(obj, cons)#构建模型
prob.solve(solver='ECOS')#模型求解
print("最优值:", prob.value)
print("最优解:", x.value)
print("剩余牛奶：", right[0] - sum(left[0] * x.value))
print("剩余劳动时间：", right[1] - sum(left[1] * x.value))
print("A1剩余加工能力：", right[2] - sum(left[2] * x.value))

例2

import cvxpy as cp
import numpy as np

coef = np.array([2, 3, 4])#输入目标函数系数
left = np.array([[1.5, 3, 5], [280, 250, 400]])#输入约束条件系数
right = np.array([600, 60000])#输入输入约束条件上限值
x = cp.Variable(3, integer=True)#创建决策变量，并且为整数
obj = cp.Maximize(coef @ x)#构造目标函数
cons = [x >= 0, left @ x <= right]#构造约束条件
prob = cp.Problem(obj, cons)#构建模型
prob.solve(solver="GUROBI")#模型求解
print("最优值：", prob.value)
print("最优解：", x.value)
print("钢材剩余量：", right[0] - sum(left[0] * x.value))
print("劳动时间剩余量：", right[1] - sum(left[1] * x.value))

例题3

import cvxpy as cp
import numpy as np

coef = np.array([2, 3, 4])
left = np.array([[1.5, 3, 5], [280, 250, 400]])
right = np.array([600, 60000])
x = cp.Variable(3, integer=True)
y = cp.Variable(3, integer=True)
obj = cp.Maximize(coef @ x)
cons = [x >= 0, left @ x <= right,
        y >= 0, y <= 1,
        x[0] >= 80 * y[0], x[0] <= 1000 * y[0],
        x[1] >= 80 * y[1], x[1] <= 1000 * y[1],
        x[2] >= 80 * y[2], x[2] <= 1000 * y[2], ]
prob = cp.Problem(obj, cons)
prob.solve(solver="GUROBI")
print("最优值：", prob.value)
print("最优解：", x.value)
print("钢材剩余量：", right[0] - sum(left[0] * x.value))
print("劳动时间剩余量：", right[1] - sum(left[1] * x.value))

例题4

import cvxpy as cp
import numpy as np

coef_x = np.array([4.8, 5.6, 4.8, 5.6])#输入目标函数x对应系数
coef_cx = np.array([0, 5000, 9000, 12000])#输入用z表示cx的系数
coef_buy_x = np.array([0, 500, 1000, 1500])#输入用z表示x的系数
left = np.array([[0, 0, 1, 1], [-1, 0, 1, 0], [0, -2, 0, 3]])#输入约束条件系数
right = np.array([1000, 0, 0])#输入约束条件上限值
x = cp.Variable(4)#创建决策变量x
y = cp.Variable(3, integer=True)#创建0-1变量y
z = cp.Variable(4)#创建变量z
obj = cp.Maximize(coef_x @ x - coef_cx @ z)#构造目标函数
#构造约束条件
cons = [cp.sum(x[0:2]) <= 500 + cp.sum(coef_buy_x @ z),
        left @ x <= right,
        sum(coef_buy_x @ z) <= 1500,
        x >= 0,
        z[0] <= y[0], z[1] <= y[0] + y[1], z[2] <= y[1] + y[2], z[3] <= y[2],
        cp.sum(z[:]) == 1, z >= 0,
        cp.sum(y[:]) == 1,
        y >= 0, y <= 1]
prob = cp.Problem(obj, cons)#构造模型
prob.solve(solver="GUROBI")#求解模型
print("最优值：", prob.value)
print("最优解：", x.value)
print("购买原油A：", sum(coef_buy_x * z.value), "t")

例题5

import cvxpy as cp
import numpy as np
#输入目标函数系数
coef = np.array([66.8, 75.6, 87, 58.6,
                 57.2, 66, 66.4, 53,
                 78, 67.8, 84.6, 59.4,
                 70, 74.2, 69.6, 57.2,
                 67.4, 71, 83.8, 62.4])
x = cp.Variable(20, integer=True)#构造决策变量
#构造目标函数
obj = cp.Minimize(coef @ x)
#输入约束条件
cons = [x >= 0, x <= 1,
        cp.sum(x[0:4]) <= 1,
        cp.sum(x[4:8]) <= 1,
        cp.sum(x[8:12]) <= 1,
        cp.sum(x[12:16]) <= 1,
        cp.sum(x[16:20]) <= 1,
        cp.sum(x[0:20:4]) == 1,
        cp.sum(x[1:20:4]) == 1,
        cp.sum(x[2:20:4]) == 1,
        cp.sum(x[3:20:4]) == 1]
prob = cp.Problem(obj, cons)#构造模型
prob.solve(solver="GUROBI")#模型求解
print("最优值：", prob.value)
print("最优解：", x.value)

例题6

import cvxpy as cp
import numpy as np
#输入目标函数的系数
coef_obj = np.array([-0.8, -0.5, -0.5, -0.2, -0.5, -0.2, 0.1, 0.1, -0.2])
coef_credits = np.array([5, 4, 4, 3, 4, 3, 2, 2, 3])#输入课程学分系数
x = cp.Variable(9, integer=True)#构造决策变量
obj = cp.Minimize(coef_obj @ x)#构造目标函数
#输入约束条件
cons = [cp.sum(x[0:5]) >= 2,
        x[2] + [4] + x[5] + x[7] + x[8] >= 3,
        x[3] + x[5] + x[6] + x[8] >= 2,
        2 * x[2] - x[0] - x[1] <= 0,
        x[3] - x[6] <= 0,
        2 * x[4] - x[0] - x[1] <= 0,
        x[5] - x[6] <= 0,
        x[7] - x[4] <= 0,
        2 * x[8] - x[0] - x[2] <= 0,
        x >= 0, x <= 1]
prob = cp.Problem(obj, cons)#模型构建
prob.solve(solver="GUROBI")#模型求解
print("选课结果：", x.value)
print("学分总和：", sum(coef_credits * x.value))