CVXPY/Scipy/Octave线性规划问题求解

相关软件/工具包：
（1）Python3.6
（2）cvxpy1.0.21
（3）scipy1.2.1
（4） Octave 5.1

博主相关博文：
（1）python线性规划（linear programming）与分配问题（assignment problem）
（2）Scipy/CVXPY/Octave非线性约束优化

线性规划问题，在计算求解过程主要确定下面4点：
（1）目标函数
（2）等式约束
（3）不等式约束
（4）变量约束

CVXPY是一个基于python环境的功能强大的凸优化计算软件包，使用其求解一个简单的线性规划问题 ，

min f=-4a+b+7c (1)目标函数
s.t.
a+b-c=5 （2）等式约束
3a-b+c<=4 （3）不等式约束
a+b-4c<=-7
a,b>=0 （4）变量约束

import cvxpy as cp
import numpy as np

'''
min f=-4a+b+7c

s.t.
a+b-c=5 

3a-b+c<=4
a+b-4c<=-7 

a,b>=0
'''
c = np.array([-4, 1, 7])
A = np.array([[3, -1, 1], [1,1,-4]])
b = np.array([4, -7])

x = cp.Variable(3)
prob=cp.Problem(cp.Minimize(c.T@x),
            [x[0]+x[1]-x[2]==5,A@x<=b,x[0]>=0,x[1]>=0])
prob.solve()

print(prob.value)
print(x.value)

结果为：

25.75
[2.25 6.75 4.  ]

从结果可以看出，当a,b,c取值分别为2.25,6.75,4时，目标函数取值最新，最小值为25.75，和matlab计算结果一致，博主在之前一篇介绍scipy.optimize.linprog的博客中的线性规划问题没有涉及到等式约束，这里同样以该函数实现本文中的线性规划问题。

scipy.optimize.linprog的函数原型如下：

scipy.optimize.linprog(c, A_ub=None, b_ub=None, A_eq=None, b_eq=None, bounds=None, method='simplex', callback=None, options=None)[source]¶

求解代码如下：

import cvxpy as cp
import numpy as np
from scipy.optimize import linprog
'''
min f=-4a+b+7c

s.t.
a+b-c=5 

3a-b+c<=4
a+b-4c<=-7 

a,b>=0
'''

c = np.array([-4, 1, 7])
A = np.array([[3, -1, 1], [1,1,-4]])
b = np.array([4, -7])
A_eq = np.array([[1,1,-1]]) # 注意这里的维度是(1,3)
b_eq = np.array([5.])

x0_bounds = (0, None)
x1_bounds = (0, None)
x2_bounds = (None, None)

res = linprog(c, 
              A_ub=A,       # 不等式约束
              b_ub=b, 
              A_eq = A_eq,  # 等式约束
              b_eq = b_eq,
              bounds=(x0_bounds, x1_bounds,x2_bounds), # 变量边界约束
              options={"disp": True})

print(res)

结果如下：

Optimization terminated successfully.
         Current function value: 25.750000   
         Iterations: 3
     con: array([0.])
     fun: 25.75
 message: 'Optimization terminated successfully.'
     nit: 3
   slack: array([0., 0.])
  status: 0
 success: True
       x: array([2.25, 6.75, 4.  ])

和cvxpy计算结果一致。

同样，我们使用Octave提供的glpk函数计算线性优化问题，glpk函数原型如下：

[xopt, fmin, errnum, extra] = glpk (c, A, b, lb, ub, ctype, vartype, sense, param)

参数说明（重要）：
（1）c:现象规划问题的系数，注意转置操作
（2）A：系数矩阵，这个不同于cvxpy或scipy.optimize.linprog将等式约束和不等式约束分隔开来，在Octave中glpk是将全部约束条件（等式约束和不等式约束）写在同一系数矩阵中，那么如何区分是等式约束还是不等式约束呢？下面一个参数会介绍
（3）b，约束条件的右边值
（4）lb：变量约束左边界
（5）ub：upper boundary，变量右边界约束
（6）ctype：constraint type约束类型，有下面几种类型和意义

"F"
A free (unbounded) constraint (the constraint is ignored). 
"U"
An inequality constraint with an upper bound (A(i,:)*x <= b(i)). 
"S"
An equality constraint (A(i,:)*x = b(i)). 
"L"
An inequality with a lower bound (A(i,:)*x >= b(i)). 
"D"
An inequality constraint with both upper and lower bounds (A(i,:)*x >= -b(i)) and (A(i,:)*x <= b(i)). 

前面说到，就是该参数来控制约束类型
（7）vartype：变量类型，'C’表示连续性变量，'I’表示整数变量
（8）sense：1为minimize问题，-1为maximize问题，这个一定不能搞错
（9）param：其他的一些参数，详细可以参照Octave文档

Octave代码：

c = [-4, 1, 7]';
A = [ 1, 1, -1;
      3, -1, 1;
      1,1,-4];
b = [5, 4, -7]';
lb = [0, 0, 0]';
ub = [];
ctype = "SUU"; % 表示系数矩阵第一行代表'S'，即等式约束，第二三行代表'U'，即Ax<=b形式的约束
vartype = "CCC"; %表示三个求解变量都是连续型
%% If sense is 1, the problem is a minimization. 
%% If sense is -1, the problem is a maximization. The default value is 1.
s = 1; % minimize问题

param.msglev = 1;
param.itlim = 10000;

[xmin, fmin, status, extra] = ...
   glpk (c, A, b, lb, ub, ctype, vartype, s, param)

计算结果如下：

xmin =

   2.2500
   6.7500
   4.0000

fmin =  25.750
status = 0
extra =

  scalar structure containing the fields:

    lambda =

       2.4167
      -1.2500
      -2.6667

    redcosts =

       0
       0
       0

    time = 0
    status =  5

由结果可见和前面两种方法求解结果相同。

和scipy.optimize.linprog不同的是，其约束条件可以灵活控制，不一定非得转换为Ax>=b的形式。

Octave的glpk有个优点是变量可以控制为整数，这个比较符合我们实际问题中遇到的情况，比如货物运输问题，需要找到货车的最优配置方案，这就需要约束为整数规划问题，因为货车不可能为半辆（前面的问题中将vartype修改为’III’便能求出整数解）。cvxpy也可以进行整数规划，但笔者还没学习到。