python3实现线性规划求解

Background

对于数学规划问题，有很多的实现。Matlab+YALMIP+CPLEX这个组合应该是比较主流的，尤其是在电力相关系统中占据着比较重要的地位。
MATLAB是一个强大的数值计算工具，用于数学建模、算法开发和数据分析。
Yalmip是一个MATLAB工具箱，用于建模和解决凸优化问题。它提供了一个简单的语法，使用户能够轻松地定义优化问题，并使用各种内置求解器求解这些问题。
Cplex是一个商业优化求解器，由IBM公司开发。它可以用于解决各种优化问题，包括线性规划、混合整数线性规划和二次规划等。在MATLAB中，用户可以使用Yalmip接口轻松地与Cplex集成。
目前 cplex 对 python 的支持目前还不是太全，相关的学习资料比较少，ibm 自己出的资料对 python 包的介绍也很简略，例子及相关类方法的介绍也不详细，这一点远没有对 java 或 c++ 支持地好。
关于一些求解简单线性规划问题，python中也有一些库可以实现：
z3-solver是由Microsoft Research(微软)开发的SMT求解器，它用于检查逻辑表达式的可满足性，可以找到一组约束中的其中一个可行解，缺点是无法找出所有的可行解（对于规划求解问题可以是scipy）。z3-solver可应用于软/硬件的验证与测试、约束求解、混合系统的分析、安全、生物，以及几何求解等问题。Z3 主要由 C++ 开发，提供了 .NET、C、C++、Java、Python 等语言调用接口。
scipy库中的函数scipy.optimize.linprog也可以进行线性规划求解，但不支持整数约束，只能求解出实数。
pulp库是一个专门进行规划求解的库。pulp库也不是万能的，虽然可以解决线性规划问题，但不能进行非线性的规划求解。当然对于规划求解，95%以上的场景都是线性规划求解，pulp就足够应对我们需要应对的场景。pulp库它将优化问题描述为数学模型，生成MPS或者LP文件，然后调用LP求解器，如CBC、GLPK、CPLEX、Gurobi等来进行求解。
安装完pulp库默认就拥有了CBC求解器，其他求解器需要额外安装才能使用，这里我没有测试这个。
cvxpy库是一个用于凸优化的Python库，它提供了一个简洁的、符合数学约束的方式来定义和求解各种凸优化问题。cvxpy的目标是提供一个易于使用的界面，使用户能够以简洁的方式表达凸优化问题，并通过优化求解器快速求解。内建的凸优化问题集类：cvxpy包含了一系列常见的凸优化问题模型，如线性规划、二次规划、半正定规划等。支持多种求解器：cvxpy可以与多个优化求解器集成，包括open-source的solver（如SCS、ECOS等）以及商业求解器（如Gurobi、CPLEX等）。

1、单个决策变量-`scipy`

线性规划案例

1.1、`matlab+yalmip+cplex`实现

%% yalmip+cplex
%（1）设定决策变量X（1）、X（2）
%（2）sdpvar：实数变量；binvar：0—1变量；intvar：整型变量
%（3）Yalmip默认是对称的，要求非对称用full

% 清除工作区
clear;clc;close all;
% function [m1, m2]=test_cplex()
% 创建决策变量
x=sdpvar(2,1,'full'); 
% 添加约束条件
st=[];
st=[st,-3*x(1)+x(2)<=6];
st=[st,x(1)+2*x(2)>=4];
st=[st,x(1)+3*x(2)==4];
st=[st,x(2)>=-3];
% 配置求解器
ops = sdpsettings('solver','cplex','verbose',0);
% 目标函数，默认最小
z=-2*x(1)+4*x(2); 
% 求解，如果最大值用-z
reuslt = optimize(st,z,ops);
% 求解结果x的值
x=value(x);
% 目标函数的最优解，即最小值
z=value(z);
% 打印结果
fprintf('x:%s\n', join(string(x),', '))
fprintf('z:%d\n', z)

运行结果

x:13, -3
z:-38

在这里插入图片描述

1.2、`python3+scipy`实现

scipy.optimize.linprog

这里只介绍后面要用到的几个参数。

在这里插入图片描述

代码实现

import numpy as np
from scipy.optimize import linprog


def main():
    """主函数"""
    # 要最小化的线性目标函数的系数
    c = np.array([-2, 4])
    # 不等式约束矩阵
    A_ub = np.array([[-3, 1], [-1, -2]])
    B_ub = np.array([6, -4])
    # 等式约束矩阵
    A_ed = np.array([[1, 3]])
    B_ed = np.array([4])
    # 决策变量的最小值和最大值对序列
    bounds = [None, None], [-3, None]
    # 求解
    res = linprog(c, A_ub, B_ub, A_ed, B_ed, bounds=bounds)
    # 求解结果x的值
    x = list(res.x.round())
    # 目标函数的最优解，即最小值
    z = round(res.fun)
    print(f'x: {x}')
    print(f'z: {z}')
    print(f'msg: {res.message}')


if __name__ == '__main__':
    main()

运行结果

在这里插入图片描述

2、多个决策变量-`cvxpy`

cvxpy几个决策变量都行哈，单个也行。我这里只是用不同的方式实现记录下。

2.1、`matlab+yalmip+cplex`实现

% 创建决策变量
x=sdpvar(2,1,'full'); 
y=sdpvar(2,1,'full'); 
% 添加约束条件
st=[];
st=[st,-3*x(1)+x(2)<=6];
st=[st,x(1)+2*x(2)>=4];
st=[st,x(1)+2*x(2)<=2*y(1)];
st=[st,x(1)+3*x(2)==4];
st=[st,x(1)+4*x(2)==y(2)];
st=[st,x(2)>=-3];
% 配置求解器
ops = sdpsettings('solver','cplex','verbose',0);
% 目标函数，默认最小
z=-2*x(1)+4*y(2); 
% 求解，如果最大值用-z
reuslt = optimize(st,z,ops);
% 求解结果x的值
x=value(x);
y=value(y);
% 目标函数的最优解，即最小值
z=value(z);
% 打印结果
fprintf('x:%s\n', join(string(x),', '))
fprintf('y:%s\n', join(string(y),', '))
fprintf('z:%d\n', z)

运行结果

x:13, -3
y:3.5, 1
z:-22

在这里插入图片描述

2.2、`python3+cvxpy+cplex`实现

需要先安装cvxpy和cplex库

pip3 install cvxpy
pip3 install cplex

实现代码

import cvxpy as cp

# 创建决策变量
x = cp.Variable(2)
y = cp.Variable(2)

# 添加约束条件
constraints = []
constraints += [-3 * x[0] + x[1] <= 6]
constraints += [x[0] + 2 * x[1] >= 4]
constraints += [x[0] + 2 * x[1] <= 2 * y[0]]
constraints += [x[0] + 3 * x[1] == 4]
constraints += [x[0] + 4 * x[1] == y[1]]
constraints += [x[1] >= -3]

# 定义目标函数
objective = -2 * x[0] + 4 * y[1]

# 创建模型求解
problem = cp.Problem(cp.Minimize(objective), constraints)
problem.solve(solver=cp.CPLEX, verbose=False)

# 打印结果
x = list(x.value)
y = list(y.value)
# 目标函数的最优解，即最小值
z = round(problem.value)
print(f'x: {x}')
print(f'y: {y}')
print(f'z: {z}')