使用minimize、linprog、GA（遗传算法）、cvxpy（凸优化）实现线性规划

星辰大海936

已于 2024-08-08 08:54:06 修改

阅读量373

点赞数 3

分类专栏： python 文章标签： python numpy 开发语言

于 2024-08-06 00:21:17 首次发布

本文链接：https://blog.csdn.net/weixin_74188799/article/details/140939362

版权

python 专栏收录该内容

10 篇文章 0 订阅

订阅专栏

题目

结合标准形式

( 注意minimize和 $\le$ 符号):

把题目转化为标准形式如下:

$minimize\quad f=c^Tx \\\\ subject\ to\ \begin{matrix}Ax \le b\\ A_{eq} = b_{eq} \\ l_b \le x \le u_b\end{matrix} \\$ $\\ \\c=\begin{bmatrix} 3\\ 2\\ -1 \end{bmatrix} \\ \\A=\begin{bmatrix} 3 & -1 &1 \\ 1& 2& 0\\ -1 & -1& -1 \end{bmatrix} , b=\begin{bmatrix} 18\\ 16\\ -2 \end{bmatrix}\\ \\ A_{eq}=\begin{bmatrix} 1&2&3\\ 1&0&-1\end{bmatrix}, b_{eq}=\begin{bmatrix}15\\ 4\end{bmatrix}\\ \\ l_b=\begin{bmatrix}0\\ 0\\0\end{bmatrix}, u_b=\begin{bmatrix}16\\16\\16\end{bmatrix}$

代码1：minimize

import numpy as np
from scipy.optimize import minimize

# 定义系数
c = np.array([3, 2, -1])
A = np.array([
    [3, -1, 1],
    [1, 2, 0],
    [-1, -1, -1]
])
b = np.array([18, 16, -2])
A_eq = np.array([
    [1, 2, 3],
    [1, 0, -1]
])
b_eq = np.array([15, 4])
l_b = np.array([0, 0, 0])
u_b = np.array([16, 16, 16])

# 定义目标函数
def objective(x):
    return np.dot(c, x)

# 定义约束
def inequality_constraints(x):
    return np.dot(A, x) - b
def equality_constraints(x):
    return np.dot(A_eq, x) - b_eq

# 初始猜测值
x0 = np.array([0, 0, 0])

# 定义约束条件
constraints = [
    {'type': 'ineq', 'fun': inequality_constraints}, # 按照字典的形式分别指定了约束的类型( 等式约束和不等式约束 )
    {'type': 'eq', 'fun': equality_constraints}
]

# 定义变量的上下界
bounds = list(zip(l_b, u_b))
# 例如，zip([0, 0, 0], [16, 16, 16]) 会生成 [(0, 16), (0, 16), (0, 16)]这样一个元组迭代器
# list(zip(l_b, u_b))将 zip 生成的迭代器转换为列表。结果是一个包含每个变量上下界的元组列表，例如 [(0, 16), (0, 16), (0, 16)]。

# 求解优化问题
result = minimize(objective, x0, method='SLSQP', bounds=bounds, constraints=constraints)

# 输出结果
if result.success:
    print(f"x = {result.x}")
    print(f"Objective = {result.fun}")
else:
    print("No solution found")

结果:（结果出错了，见最后）

No solution found

进程已结束，退出代码为 0

代码2：linprog线性规划

from scipy.optimize import linprog
c=[[3], [2], [-1]]
A=[
    [3, -1, 1],
    [1, 2, 0],
    [-1, -1, -1]
]
b=[[18], [16], [-2]]
A_eq=[[1,2,3],[1,0,-1]]
b_eq=[[15],[4]]
bounds=((0,None),(0,None),(0,None))
res=linprog(c=c, A_ub=A, b_ub=b, A_eq=A_eq, b_eq=b_eq, bounds=bounds)
print (res)

结果:

最小值: fun: 19.166666666666668, 此时决策变量x: [ 5.917,1.667,1.917]

E:\anaconda\pythonw.exe E:\pycharm\python_doc\FromBook\0_summer_training\whale\3\线性规划.py 
        message: Optimization terminated successfully. (HiGHS Status 7: Optimal)
        success: True
         status: 0
            fun: 19.166666666666668
              x: [ 5.917e+00  1.667e+00  1.917e+00]
            nit: 0
          lower:  residual: [ 5.917e+00  1.667e+00  1.917e+00]
                 marginals: [ 0.000e+00  0.000e+00  0.000e+00]
          upper:  residual: [       inf        inf        inf]
                 marginals: [ 0.000e+00  0.000e+00  0.000e+00]
          eqlin:  residual: [ 0.000e+00  0.000e+00]
                 marginals: [ 8.333e-01  3.167e+00]
        ineqlin:  residual: [ 0.000e+00  6.750e+00  7.500e+00]
                 marginals: [-3.333e-01 -0.000e+00 -0.000e+00]
 mip_node_count: 0
 mip_dual_bound: 0.0
        mip_gap: 0.0

进程已结束，退出代码为 0

代码3：GA遗传算法

from sko.GA import GA
def func(x):
    return 3*x[0]+2*x[1]-x[2]
cons1=lambda x: x[0]+2*x[1]+3*x[2]-15
cons2=lambda x: x[0]-x[2]-4
con1=lambda x: 3*x[0]-x[1]+x[2]-18
con2=lambda x:x[0]+2*x[1]-16
con3=lambda x: -x[0]-x[1]-x[2]-2
ga=GA(func=func,n_dim=3,size_pop=500,max_iter=500,constraint_eq=[cons1,cons2],constraint_ueq=[con1,con2,con3],lb=[0,0,0],ub=[16,16,16])
# 参数：
    # func:目标函数
    # n_dim:目标函数的维度
    # size_pop:种群的大小，这里指定了种群大小为500，即每次迭代时有500个个体。
    # max_iter:最大迭代次数
    # constraint_eq:等式约束条件，是一个列表，每个元素是一个函数，输入是个体的解向量，输出是该个体是否满足约束条件。
    # constraint_ueq:不等式约束条件，是一个列表，每个元素是一个函数，输入是个体的解向量，输出是该个体是否满足约束条件。
    # lb:每个维度的下界，是一个列表。
    # ub:每个维度的上界，是一个列表。
best_x,best_y=ga.run()
print("best x:\n",best_x,"\nbest_y:\n",best_y)

结果：

E:\anaconda\pythonw.exe E:\pycharm\python_doc\FromBook\0_summer_training\test.py 
best x:
 [5.15894552 3.18210889 1.15894556] 
best_y:
 [20.68210878]

进程已结束，退出代码为 0

结果中可以发现，代码3的结果比代码2要大。这是因为遗传算法作为一种智能算法，不能保证结果是最优的。可以通过提高迭代次数改进结果：

#最大迭代次数：1000
E:\anaconda\pythonw.exe E:\pycharm\python_doc\FromBook\0_summer_training\test.py 
best x:
 [5.76398792 1.97202409 1.76398797] 
best_y:
 [19.47202398]

进程已结束，退出代码为 0

代码4：cvxpy凸优化

import cvxpy as cp
import numpy as np

coef = np.array([3, 2, -1])
left = np.array([[3, -1, 1],[1, 2, 0],[-1, -1, -1]])
right = np.array([18, 16, -2])
x = cp.Variable(3)# 可以指定 integer=True，表示变量只能取整数值；或者boolen=True，表示变量只能取0或1值
obj = cp.Minimize(coef @ x)# 也可以选择 Maximize
cons = [x >= 0, left @ x <= right, # 约束条件可以选择 >=, <=, ==, !=
        np.array([[1,2,3],[1,0,-1]])@x==np.array([15,4])]
prob = cp.Problem(obj, cons)
prob.solve(solver=cp.CLARABEL)
print("最优值：", prob.value)
print("最优解：", x.value)

cvxpy的好处在于可以在定义变量的时候指定变量类型，同时对于约束条件的限制比较少。也就是可以不写成规范形式。

结果：

E:\anaconda\pythonw.exe E:\pycharm\python_doc\FromBook\0_summer_training\test.py 
最优值： 19.16666666699132
最优解： [5.91666667 1.66666667 1.91666667]

进程已结束，退出代码为 0

代码1问题分析

1. result是这样的. 说明结果错误了.

 message: Positive directional derivative for linesearch
 success: False
  status: 8
     fun: 0.0
       x: [ 0.000e+00  0.000e+00  0.000e+00]
     nit: 5
     jac: [ 3.000e+00  2.000e+00 -1.000e+00]
    nfev: 4
    njev: 1
x = [0. 0. 0.]
Objective = 0.0

进程已结束，退出代码为 0

2. 添加模块分析是不是初始值的问题:

for i in tqdm(range(17)):
    for j in range(17):
        for k in range(17):
            x0 = np.array([i, j, k])
            result = minimize(objective, x0, method='SLSQP', bounds=bounds, constraints=constraints)
            if result.success:
                print(result)
                exit()
print("No solution found")

结果发现不是初始值的问题: 仍然找不到函数值的最小值 ... ...

E:\anaconda\pythonw.exe E:\pycharm\python_doc\FromBook\0_summer_training\whale\3\线性规划.py 
100%|██████████| 17/17 [00:54<00:00,  3.22s/it]
No solution found

3. 去问chatgpt：

问题可能出在 minimize 函数使用的约束定义上。在 minimize 函数中，不等式约束默认是要求函数返回的值大于等于0（即 g(x) >= 0），而你的 inequality_constraints 函数返回的是 np.dot(A, x) - b，这实际上是形成了 Ax <= b 的约束。为了符合 minimize 的要求，你需要将不等式约束改为 b - np.dot(A, x)，这样才能确保 g(x) >= 0 的形式。

结果：

E:\anaconda\pythonw.exe E:\pycharm\python_doc\FromBook\0_summer_training\test.py 
x = [5.91666667 1.66666667 1.91666667]
Objective = 19.166666659241173

进程已结束，退出代码为 0

结论：minimize的不等式约束默认是要求函数返回的值大于等于0（即 g(x) >= 0）。与linprog不同。