2021 January 11 书籍, 阅读

简单线性规划问题的求解

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

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

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

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

过程

问题

当x,y满足下面条件：

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

解决

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

import numpy as np
from scipy import optimize as op
x=(0,None)
y=(0,None)
a_ub=np.array([[-1,-1],[-1,-3]])
b_ub=np.array([-4,-6])
c=np.array([3,5])
res=op.linprog(c,a_ub,b_ub,bounds(x,y))
res
     con: array([], dtype=float64)
     fun: 13.999999999646015
 message: 'Optimization terminated successfully.'
     nit: 4
   slack: array([-9.11426490e-11, -1.71698211e-10])
  status: 0
 success: True
       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

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

输出为：3.0 1.0 14.0

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

再进一步

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

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

# 矩阵计算

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

⎡1⎤
⎢ ⎥
⎣9⎦

微分方程

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

t = symbols('t')

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

           3⋅t    
       C₁⋅ℯ      2
y(t) = ─────── - ─
          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

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

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

解二阶微分方程

y‘‘+9y=1

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

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

求解不定积分

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

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

求解定积分

更好的资料参看这里：

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