简单线性规划问题的求解#

在看算法设计与分析的书中，练习题里面有个线性规划的问题。

利用坐标和画出线性图像，可以解决。但是这样效率不高。

经过搜索，可以使用python的scipy这个库进行解决。

很早前安装这个库是为了求多元多次方程的。

过程#

问题

当x,y满足下面条件：

x+y<=4
x+3y<=6
x>=0,y>=0 求 3x+5y的最大值。

解决

要说明的是求解的时候要给各个方程乘以-1，确保左边的值>=右边的值，而不是小于。故有： -x-y>=-4 -1-3y>=-6

1
import numpy as np
2
from scipy import optimize as op
3
x=(0,None)
4
y=(0,None)
5
a_ub=np.array([[-1,-1],[-1,-3]])
6
b_ub=np.array([-4,-6])
7
c=np.array([3,5])
8
res=op.linprog(c,a_ub,b_ub,bounds(x,y))
9
res
10
     con: array([], dtype=float64)
11
     fun: 13.999999999646015
12
 message: 'Optimization terminated successfully.'
13
     nit: 4
14
   slack: array([-9.11426490e-11, -1.71698211e-10])
15
  status: 0
16
 success: True
17
       x: array([3., 1.])

上面的scipy求出来的是float类型的，如果想更精确，可以这么做： res=op.linprog(c,a_ub,b_ub,bounds(x,y),method=“revised simplex”)

即当x=3,y=1的时候，3x+5y可以拿到最大值，14.

更进一步#

PuLP is an LP modeler written in python. PuLP can generate MPS or LP files and call GLPK, COIN CLP/CBC, CPLEX, and GUROBI to solve linear problems.

这个库更强大一些，专门为各种优化问题提供的，最大，最小，最优等等。。。

更详细的信息可以看这个： https://coin-or.github.io/pulp/main/the_optimisation_process.html

case study开始：https://coin-or.github.io/pulp/CaseStudies/a_blending_problem.html

1
import pulp as p
2
# 求最大值为LpMaximize,求最小值为LpMinimize
3
lp1=p.LpProblem('problem',p.LpMaximize)
4
# x,y的边界3
5
x=p.LpVariable('x',lowBound=0)
6
y=p.LpVariable('y',lowBound=0)
7
# 输入结构方程
8
lp1+=3*x+5*y
9
lp1+=x+y<=4
10
lp1+=x+3*y<=6
11
# print(lp1)
12
# status=lp1.solve()
13
# 求解
14
print(p.value(x),p.value(y),p.value(lp1.objective))

输出为：3.0 1.0 14.0

可以看到pulp的方式比scipy的更加适合处理规划问题。简单而且不含糊。

再进一步#

在进行多元多项式的处理的时候，scipy会存在只考虑正整数，忽略负数的情况，而且不能用无理数替换float的数，导致求出来的值无法在理论意义上相等。而sympy这个模块可以实现。

1
from sympy import *
2
init_printing(use_unicode=True)
3
sqrt(8)
4
>>> 2⋅√2
5

6
# 矩阵计算
7

8
M=Matrix([[1,-1,2],[3,4,5]])
9
N=Matrix([0,1,1])
10
M*N
11

12
⎡1⎤
13
⎢ ⎥
14
⎣9⎦

微分方程#

我們看範例，解微分方程式 y’-3y=2

1
t = symbols('t')
2

3
y = Function('y')
4
dsolve(Eq(Derivative(y(t), t) -3*y(t), 2), y(t))

1
           3⋅t
2
       C₁⋅ℯ      2
3
y(t) = ─────── - ─
4
          3      3

Derivative(y(t), t)就是微一階y’(t)；若微二階 y”(t) 就是Derivative(y(t), t, t).

解微分方程 y’=(y+1)/(t-3)#

转化为 y’-(y+1)/(t-3)=0

1
t=symbols('t')
2
y=Function('y')
3
dsolve(Eq(Derivative(y(t),t)-(y(t)+1)/(t-3),0),y(t))

1
y(t) = C₁⋅t - 3⋅C₁ - 1

解二阶微分方程#

y‘‘+9y=1

1
t=symbols('t')
2
y=Function('y')
3
dsolve(Eq(Derivative(y(t),t,t)+9*y(t),1),y(t))

1
y(t) = C₁⋅sin(3⋅t) + C₂⋅cos(3⋅t) + 1/9

求解不定积分#

积分2x/(x*x-1)

1
import math
2
from sympy import *
3
x=symbols('x')
4
integrate(2*x/(x**2-1),x)

求解定积分#

更好的资料参看这里：

http://web.ntnu.edu.tw/~tsungwu/Python_DevOps/Part1_Basics&Math/section5_integration.html