编写：@谭胖子重新启航

修订：@_散沙_Python玩家_

Concrete Math with Python IV

Python带你游览具体世界的数学基础

一般来说，高深的数学总会让人有种抗拒心理，但跟具体世界紧密联系的简单数学知识会则不容易让人有如下的疑问：

"我们知道这些知识到底有什么用处？"

这些数学知识第一次出现在高中、大学课本时可能会让人摸不到头脑，而他们往往是基础的计算机科学、统计学、经济学知识的基石。

这部分内容包括：

数列与金融基础
排列组合与概率基础
方程组/矩阵
线性规划
质数与加密
随机数与概率分布

注意：具备扎实的高中数学、理工科本科数学基础的朋友可以简单阅览或者跳过这一部分。



In [22]:

    
import numpy as np
import scipy as sp
import scipy.signal as ssig
import bokeh as bk
import bokeh.plotting as bplt
import bokeh.palettes as bpal



In [23]:

    
# 工具配置
TOOLS = "pan,wheel_zoom,box_zoom,box_select"
# 作图配置
plot_config = dict(plot_width=800, plot_height=500, tools=TOOLS)
bplt.output_notebook()









    




    
        
        
        
    
        
        BokehJS successfully loaded.

1 运筹学的分支: 线性规划

考虑一个家庭要买一定量的食材来做饭，而且今晚决定合理的搭配牛肉和羊肉的量。

首先是热量：牛肉1000千卡/kg,羊肉2000/kg，总热量不要超过4000千卡
首先是价格：牛肉40元/kg，羊肉20元/kg，总价不要超过120元
最后是孩子的喜好：家里的孩子喜欢吃牛肉胜过吃羊肉，尽管羊肉便宜，但是需要保证羊肉不要比牛肉多出1kg以上。

全家的目标：出门买一次菜不容易，所以要追求一次买的牛肉与羊肉的总重量达到最大。

$$Max:{x_1} + {x_2}$$$$s.t.$$$$\eqalign{ & {x_1} + 2{x_2} \leqslant 4 \cr & 4{x_1} + 2{x_2} \leqslant 12 \cr & - {x_1} + {x_2} \leqslant 1 \cr} $$

变量的线性不等式，称为约束（热量、价格、孩子的喜好），因为可以不取等号，所以矩阵法求解就无能为力了
在满足所有约束前提下，满足所有不等式条件的变量取值成为可行域
在可行域中，可以搜寻使得目标函数（总重量）取得最大值的点



In [46]:

    
xlim,ylim = np.array([0,3]),np.array([0,3])
res1y = (4-xlim)/2.0
res2x = (12-2*ylim)/4.0
res3x = -1+ylim
cross1,cross2,cross3 = np.array([8.0/3,2.0/3]),np.array([5.0/3,8.0/3]),np.array([2.0/3,5.0/3])
tarfuncoff,refp1 = np.array([1,-1]),np.array([0,1])

p1 = bplt.figure(title='线性规划：约束的可行域与优化目标',background_fill='#FFFFFF',**plot_config)
p1.set(x_range=bk.models.Range1d(-2,4))

p1.line(xlim,res1y,line_width=3)
p1.line(res2x,ylim,line_width=3)
p1.line(res3x,ylim,line_width=3)
p1.line(*zip(refp1+tarfuncoff,refp1-tarfuncoff),line_width=3,line_color='darkred',line_dash='8 8')
p1.line(*zip(cross1+tarfuncoff,cross1-tarfuncoff),line_width=3,line_color='darkred',line_dash='8 8')
p1.line(*zip(cross3+tarfuncoff,cross3-tarfuncoff),line_width=3,line_color='darkred',line_dash='8 8')
p1.patch(*zip(cross1,cross3,np.array([-1,0]),np.array([3,0])),alpha=0.3)
bplt.show(p1)

通过代数运算，我们可以做到：

先画出约束表示的三条线，这三条线在坐标系中围出的范围就是蓝色区域的可行域
再水平画出我们的最大化目标函数的等高线，用红色虚线来表示。

在同一条虚线上的不同点，目标函数取值一样。随着红色虚线向右移动，目标函数取值越来越大，直到触碰到$(\frac{8}{3},\frac{2}{3})$这个点。

这样如果我们还想让目标函数取到更大的值也不可能了，因为已经到达了蓝色区域的角点。我们用不严密的数学语言来说明这个现象：

优化函数在约束条件围成闭区域的时候，最优值一定出现在角点
如果约束条件围成开区域（如上例）的时候，最优值有两种可能：
- 最优值是无界的：它能取到正无穷大或者负无穷大
- 最优值是有界的：它是有限值，并且出现在角点

2 线性规划问题的对偶性

刚才提到，所谓最优值出现在角点的意思，就是构成那个角点的约束都取到了等号

标准最大化问题

$$\eqalign{ & Max:{c^T}x \cr & Ax \leqslant b \cr} $$

其中小写字母为列向量，大写字母为矩阵。

如果我们能够找到约束取到等号的那些行，并且让行向量通过系数$y$进行线性组合成为$C^T$（没有形成约束的第$i$行，对应的$y_i$取0即可）：

$$y^TA = c^T$$

这样我们就知道$$c^Tx = y^TAx = y^Tb$$

这样我们就把刚才"寻找合适的$x$"的问题改变成了"寻找合适的$y$"，也就找到了$c^Tx$能达到的最大值$y^Tb$。

这个问题跟如下的问题写法不一致，但是内涵是一致的，因此跟如下的问题被称作线性规划的对偶问题：

标准最小化问题

$$\eqalign{ & Min:{y^T}b \cr & {y^T}A \geqslant c \cr} $$

为什么找到对偶问题很重要

有时，直接求解现有问题不一定比求解对偶问题更直观。

对于偷懒的我们来说，哪个问题方便，就求解哪个，反正最后的最优化函数值都会殊途同归。

3 线性规划问题的定义：pulp的AMPL

也许你在学校里已经学习过大名鼎鼎的LINGO/LINDO软件来求解线性规划问题了。当然没接触过这类软件也没关系。

下面我们需要用Python中的pulp来解决这类问题。要想解决问题，先要用求解器所能读懂的语言定义这些问题。这个语言就是我们先来学习一个好玩的东西，AMPL(A Math Programming Language),AMPL支持集合(Sets)与参数(Parameters)两种类型，参数和Python内置的list表现很类似。



In [64]:

    
import pulp
data = pulp.Amply("set CITIES := Auckland Wellington Christchurch;")
print data.CITIES
print data['CITIES']
for c in data.CITIES:
    print(c)









    



<SetObject: ['Auckland', 'Wellington', 'Christchurch']>
<SetObject: ['Auckland', 'Wellington', 'Christchurch']>
Auckland
Wellington
Christchurch

Sets可以有不同类型的参数



In [53]:

    
data = pulp.Amply("""
set BitsNPieces := 0 3.2 -6e4 Hello "Hello, World!";
""")
print data.BitsNPieces









    



<SetObject: [0.0, 3.2, -60000.0, 'Hello', 'Hello, World!']>

Sets可以包含多维的数据



In [54]:

    
data = pulp.Amply("""
set pairs dimen 2;
set pairs := (1,2) (2,3) (4,2) (3,1) (2,2) (4,4) (3,4);
""")
print data.pairs









    



<SetObject: [(1.0, 2.0), (2.0, 3.0), (4.0, 2.0), (3.0, 1.0), (2.0, 2.0), (4.0, 4.0), (3.0, 4.0)]>

Sets本身也可以是多维的



In [56]:

    
data = pulp.Amply("""
... set CITIES{COUNTRIES};
... set CITIES[Australia] := Adelaide Melbourne Sydney;
... set CITIES[Italy] := Florence Milan Rome;
... """)
print data.CITIES['Italy']









    



['Florence', 'Milan', 'Rome']

Sets也可以通过矩阵来指定，矩阵的 + 代表这个pair存在，- 则代表不存在



In [55]:

    
data = pulp.Amply("""
set ROUTES dimen 2;
set ROUTES : A B C D :=
            E + - - +
            F + + - -;
""")
print data.ROUTES









    



<SetObject: [('E', 'A'), ('E', 'D'), ('F', 'A'), ('F', 'B')]>

AMPL的参数(Parameters)也是得到支持的



In [58]:

    
data = pulp.Amply("""
param T := 30;
param n := 5;
""")
print data.T
print data.n



In [60]:

    
data = pulp.Amply("""
param COSTS{PRODUCTS};
param COSTS :=
FISH 8.5
CARROTS 2.4
POTATOES 1.6
;
""")
print data.COSTS









    



<ParamObject: {'POTATOES': 1.6, 'FISH': 8.5, 'CARROTS': 2.4}>

参数本身也是支持多维的



In [62]:

    
data = pulp.Amply("""
param COSTS{CITIES, PRODUCTS};
param COSTS :=
    Auckland FISH 5
    Auckland CHIPS 3
    Wellington FISH 4
    Wellington CHIPS 1
;
""")
print data.COSTS['Auckland']['CHIPS']
print data.COSTS['Auckland','FISH']

二维参数可以通过矩阵的方式来指定，高维的参数可以通过切片来实现：



In [63]:

    
data = pulp.Amply("""
param COSTS{CITIES, PRODUCTS, SIZE};
param COSTS :=
[Auckland, *, *] :   SMALL LARGE :=
                FISH  5     9
                CHIPS 3     5
[Wellington, *, *] : SMALL LARGE :=
                FISH  4     7
                CHIPS 1     2
;
""")
print(data.COSTS)









    



<ParamObject: {'Wellington': {'FISH': {'SMALL': 4.0, 'LARGE': 7.0}, 'CHIPS': {'SMALL': 1.0, 'LARGE': 2.0}}, 'Auckland': {'FISH': {'SMALL': 5.0, 'LARGE': 9.0}, 'CHIPS': {'SMALL': 3.0, 'LARGE': 5.0}}}>

4.Show me the code?

定义好了问题之后，我们就要大干一场了。

The Diet Problem V1

问题描述

现在有5中食物${F_1}, \cdots ,{F_5}$，每种食物都含义不同数量的4种营养${N_1}, \cdots ,{N_4}$，每种营养都很重要，人体对营养$j$的需求量为${c_j}$，每种食物的价格是${b_i}$，第$i$种食物中含有$j$中营养的量为${a_{ij}}$

数学抽象 $$\min :0.013{x_1} + 0.008{x_2} + 0.01{x_3} + 0.005{x_4} + 0.015{x_5}$$

$$\eqalign{ & 0.100{x_1} + 0.200{x_2} + 0.150{x_3} + 0.000{x_4} + 0.040{x_5} \ge 8.0 \cr & 0.080{x_1} + 0.100{x_2} + 0.110{x_3} + 0.010{x_4} + 0.010{x_5} \ge 6.0 \cr & 0.001{x_1} + 0.005{x_2} + 0.003{x_3} + 0.100{x_4} + 0.150{x_5} \ge 2.0 \cr & 0.002{x_1} + 0.005{x_2} + 0.007{x_3} + 0.002{x_4} + 0.008{x_5} \ge 0.4 \cr} $$$${x_1} \ge 0;{x_2} \ge 0;{x_3} \ge 0;{x_4} \ge 0;{x_5} \ge 0$$

解数学问题



In [67]:

    
Ingredients = ['F1', 'F2', 'F3', 'F4', 'F5']
#每单位食物的成本
costs = {
    "F1":0.013,
    "F2":0.008,
    "F3":0.010,
    "F4":0.005,
    "F5":0.015
}
#每单位食物第一种营养的量
N1={
    "F1":0.100,
    "F2":0.200,
    "F3":0.150,
    "F4":0.000,
    "F5":0.040
}
#每单位食物第二种营养的量
N2={
    "F1":0.080,
    "F2":0.100,
    "F3":0.110,
    "F4":0.010,
    "F5":0.010
}
#每单位食物第三种营养的量
N3={
    "F1":0.001,
    "F2":0.005,
    "F3":0.003,
    "F4":0.100,
    "F5":0.150
}
#每单位食物第四种营养的量
N4={
    "F1":0.002,
    "F2":0.005,
    "F3":0.007,
    "F4":0.002,
    "F5":0.008
}
#选择最小化优化问题
prob = pulp.LpProblem("The Diet Problem",LpMinimize)
#我们可以通过两种方式指定变量
# 方法1：逐一指定
# x1 = LpVariable('F1',0,None)
# 方法2：通过dict直接指定
vars = pulp.LpVariable.dicts("Food",Ingredients,0)
print vars









    



{'F1': Food_F1, 'F2': Food_F2, 'F3': Food_F3, 'F4': Food_F4, 'F5': Food_F5}

将目标函数加入到prob中：



In [68]:

    
prob += pulp.lpSum([costs[i]*vars[i] for i in Ingredients])

将四个预算约束加入到prob中：



In [69]:

    
prob += pulp.lpSum(N1[i]*vars[i] for i in Ingredients)>=8.0
prob += pulp.lpSum(N2[i]*vars[i] for i in Ingredients)>=6.0
prob += pulp.lpSum(N3[i]*vars[i] for i in Ingredients)>=2.0
prob += pulp.lpSum(N4[i]*vars[i] for i in Ingredients)>=0.4
#将线性规划问题保持到文件
prob.writeLP("WhiskasModel.lp")
#求解该问题
prob.solve();
print u"状态:", LpStatus[prob.status]
print "="*32
for v in prob.variables():
    print v.name, "=", v.varValue
print "="*32
print u"所有食物的总成本：", value(prob.objective)









    



状态: Optimal
================================
Food_F1 = 0.0
Food_F2 = 4.9673203
Food_F3 = 48.366013
Food_F4 = 18.300654
Food_F5 = 0.0
================================
所有食物的总成本： 0.6149019624

The Diet Problem V2

当然这个问题可以用上面的AMPL更方便的表示出来：



In [70]:

    
data = pulp.Amply("""
set Food := F1 F2 F3 F4 F5;
set Ntr  := N1 N2 N3 N4; 
param Costs{Food};
param Costs :=
        F1 0.013
        F2 0.008
        F3 0.010
        F4 0.005
        F5 0.015;
param Need{Ntr};
param Need :=
        N1 8.0
        N2 6.0
        N3 2.0
        N4 0.4;
param Content{Food, Ntr};
param Content :  F1    F2    F3    F4    F5 :=
              N1 0.100 0.200 0.150 0.000 0.040
              N2 0.080 0.100 0.110 0.010 0.010
              N3 0.001 0.005 0.003 0.100 0.150
              N4 0.002 0.005 0.007 0.002 0.008;
""")

现在用data中包含的变量与数据构造线性规划问题：



In [71]:

    
prob2 = pulp.LpProblem("The Diet Problem2",LpMinimize)
vars = pulp.LpVariable.dicts("Food",data.Food,0)
prob2 += pulp.lpSum([data.Costs[i]*vars[i] for i in data.Food])
for j in data.Ntr:
    prob2 += pulp.lpSum(data.Content[j,i]*vars[i] for i in data.Food)>=data.Need[j]
#将线性规划问题保持到文件
prob2.writeLP("DietModel2.lp")
#求解该问题
prob2.solve()

输出优化的结果



In [72]:

    
print u"状态:", LpStatus[prob2.status]
print "="*32
for v in prob2.variables():
    print v.name, "=", v.varValue
print "="*32
print u"所有食物的总成本：", value(prob2.objective)









    



状态: Optimal
================================
Food_F1 = 0.0
Food_F2 = 4.9673203
Food_F3 = 48.366013
Food_F4 = 18.300654
Food_F5 = 0.0
================================
所有食物的总成本： 0.6149019624

The Diet Problem V3

我们换个思考问题的方式：还是刚才的营养配餐问题，原有问题中4个约束被写成了4个菜单。

我希望能够通过将4个菜单进行线性组合的方式得到一个不等式，而组合后不等式中$F1,F2,...F5$的系数刚好是它们的价格。

只要能够组合出来，因为这是一个左边大于右边的不等式，食品的最小可能成本也就根据我们的菜单重新组合出来了。



In [74]:

    
data = pulp.Amply("""
set Food := F1 F2 F3 F4 F5;
set Ntr  := N1 N2 N3 N4; 
param Need{Ntr};
param Need :=
        N1 8.0
        N2 6.0
        N3 2.0
        N4 0.4;
param Costs{Food};
param Costs :=
        F1 0.013
        F2 0.008
        F3 0.010
        F4 0.005
        F5 0.015;
param Content{Ntr, Food};
param Content :  N1    N2    N3    N4 :=
              F1 0.100 0.080 0.001 0.002
              F2 0.200 0.100 0.005 0.005
              F3 0.150 0.110 0.003 0.007
              F4 0.000 0.010 0.100 0.002
              F5 0.040 0.010 0.150 0.008;
""")



In [75]:

    
prob3 = pulp.LpProblem("The Diet Problem3",LpMaximize)
vars = pulp.LpVariable.dicts("Ntr",data.Ntr,0)
prob3 += pulp.lpSum([data.Need[i]*vars[i] for i in data.Ntr])
for j in data.Food:
    prob3 += pulp.lpSum(data.Content[j,i]*vars[i] for i in data.Ntr)<=data.Costs[j]
#将线性规划问题保持到文件
prob3.writeLP("DietModel3.lp")
#求解该问题
prob3.solve()









    Out[75]:





1



In [78]:

    
print u"状态:", LpStatus[prob3.status]
print "="*32
for v in prob3.variables():
    print v.name, "=", v.varValue
print "="*32
print u"所有营养菜单的总成本：", value(prob2.objective)









    



状态: Optimal
================================
Ntr_N1 = 0.0
Ntr_N2 = 0.036078431
Ntr_N3 = 0.029411765
Ntr_N4 = 0.84901961
================================
所有营养菜单的总成本： 0.6149019624

对偶问题总结

求解了这两个问题，我们发现：

求解原问题时，发现最后不购买$F1,F5$，对偶问题的等式矩阵只有3维。
求解对偶问题，发现最后不使用菜单$N1$，原问题的矩阵也只有3维。
不管按照分别买食物的方式还是按照菜单成本购买的方式，核算出来的成本都是同一个值。

5 总结

这一节讲了线性规划的概念，图形求解法与对偶问题。
并且提到一种全新的数学问题描述语言AMPL
讲解了利用Python求解线性规划方法的方法
并且给出了一个简单的实例。



In [ ]: