精华内容
下载资源
问答
  • 文章目录初识随机规划 (2):两阶段随机规划模型两阶段随机规划模型生产计划的例子参数的不确定性随机规划模型(Stochastic Programming)Python调用Gurobi求解随机规划模型参考文献OlittleRer 初识随机规划 (2):两...

    初识随机规划 (2):两阶段随机规划(一个详细的例子)

    本文中的确定性问题的例子来源于《运筹学》第四版,运筹学编写组编著。

    两阶段随机规划模型

    上篇文章介绍了随机规划的一个小例子。

    文章链接: https://mp.weixin.qq.com/s/9RsWE6E00W_dz7bpThHInA 或者 https://blog.csdn.net/HsinglukLiu/article/details/112150067

    为了方便阅读,我们把那个例子拿过来,再复习一遍。

    一个工厂生产2种产品 A A A B B B A A A产品需要1个单位人工工时和4个单位的设备1的工时, B B B产品需要2个单位的人工工时和4个单位的设备2的工时。且该厂一共可提供人工工时、设备1工时和设备2的工时分别为8,16和12。且生产单位产品 A A A B B B的净利润分别为2和3。假设该工厂生产的产品是供不应求的 ,问,该工厂应该如何生产(可以是连续的量),使得净利润最大化。
    根据上述描述,我们引入下面的决策变量:
    x 1 , x 2 x_1, x_2 x1,x2分别代表生产产品 A 、 B A、B AB的量,则该问题的数学模型为
    max ⁡    2 x 1 + 3 x 2 s . t . x 1 + 2 x 2 ⩽ 8 4 x 1 ⩽ 16 4 x 2 ⩽ 12 x 1 , x 2 ⩾ 0 \begin{aligned} \max \,\,&2x_1 + 3x_2 \\ s.t. \quad &x_1 + 2x_2 \leqslant 8 \\ &4x_1 \qquad \leqslant 16 \\ &\qquad 4x_2 \leqslant 12 \\ &x_1, x_2 \geqslant 0 \end{aligned} maxs.t.2x1+3x2x1+2x284x1164x212x1,x20

    但是上面的文章只是让读者直观的感受一下随机规划究竟是什么,文中的随机规划模型如下:

    max ⁡    ∑ k ∈ K p k ( c 1 x 1 k + c 2 x 2 k ) s . t . a 11 k x 1 k + a 12 k x 2 k ⩽ 8 , ∀ k ∈ K a 21 k x 1 k         ⩽ 16 , ∀ k ∈ K            a 32 k x 2 ⩽ 12 , ∀ k ∈ K x 1 k , x 2 k ⩾ 0 , ∀ k ∈ K \begin{aligned} \max \,\,& \sum_{k \in K}p_k (c_1 x_1^k + c_2 x_2^k ) \\ s.t. \quad & a_{11}^k x_1^k + a_{12}^kx_2^k \leqslant 8, \quad \forall k \in K \\ &a_{21}^kx_1^k \qquad\,\,\,\,\,\,\, \leqslant 16, \quad \forall k \in K \\ &\,\,\,\,\,\,\,\,\,\,\qquad a_{32}^kx_2 \leqslant 12, \quad \forall k \in K \\ &x_1^k, x_2^k \geqslant 0, \quad \forall k \in K \end{aligned} maxs.t.kKpk(c1x1k+c2x2k)a11kx1k+a12kx2k8,kKa21kx1k16,kKa32kx212,kKx1k,x2k0,kK
    其紧凑形式为
    max ⁡    ∑ k ∈ K p k ∑ i ∈ I c i x i k s . t . ∑ i ∈ I a j i k x i k ⩽ b j , ∀ j ∈ J , ∀ k ∈ K x i k ⩾ 0 , ∀ i ∈ I , ∀ k ∈ K \begin{aligned} \max \,\,& \sum_{k \in K}p_k \sum_{i \in I}c_i x_i^k \\ s.t. &\quad \sum_{i \in I}{a_{ji}^k x_i^k} \leqslant b_j, \qquad \forall j \in J, \forall k \in K \\ &x_i^k \geqslant 0, \qquad \qquad \forall i \in I, \forall k \in K \end{aligned} maxs.t.kKpkiIcixikiIajikxikbj,jJ,kKxik0,iI,kK
    其中 I = { 1 , 2 } I = \{1, 2\} I={1,2}表示产品的集合, J = { 1 , 2 , 3 } J = \{1, 2, 3\} J={1,2,3}表示资源的集合。 a j i k a_{ji}^k ajik表示第 j ∈ J j \in J jJ种资源被第 i ∈ I i \in I iI种商品在第 k ∈ K k \in K kK个场景下需要的量。

    细心的读者会发现,上面的模型,其实可以等价为:分别求解 K K K个线性规划,然后将 K K K个单独的线性规划的目标值加权求和一下。原因是这 K K K个场景下的模型没有任何联系,是相互独立的。

    难道说随机规划就是很多个确定性线性规划简单的分别求解一下,然后把目标函数加和一下这么简单吗?当然不是。

    3层决策:战略层(strategic)战术层(tactical)作业层(operational)

    一般情况下,随机规划都是有多个阶段的决策的。一般来讲,决策有三个层次:
    战略层(strategic)战术层(tactical)作业层(operational)。从strategictactical再到operational,决策的周期越来越短,决策的事情越来越具体化。随机规划中,比较常见的是两个阶段的决策two-stage decision,因此,读论文的时候,经常看到关于随机规划的文章是把问题建模为一个two-stage mixed integer programming的模型求解的。

    所谓的two-stage,一般可以认为包含两种决策,一种是比较长期的 ,上层的决策,比如服务网络怎么搭建,物流网络怎么搭建,基础设施怎么搭建等,这种决策会影响企业较长时间的发展,我们称其为战略层决策(strategic level),有的决策也许没有到战略层的级别,可以称为战术层决策(tactical level)。第二种决策就是非常具体化的,比如生产计划,采购计划,VRP的路径规划等,都是非常具体、细节的决策,这一类决策对企业的发展影响周期较短,可以称为作业层决策(operational level)。为了方便统一,本文中我们统一使用tactical leveloperational level来表示上层决策和下层决策。

    Two-stage Stochastic Programming

    Two-stage Stochastic Programming的直观写法

    还是回到随机规划上来。我们回去看之前的例子,例子中说,我们有生产产品 A A A和产品 B B B的设备1和设备2,有了设备,我们可以决定生产 A A A B B B。但是在实际中,也许我们可以再往前座一层决策,那就是,我们决策一下,要不要接产品 A A A和产品 B B B的订单,也就是要不要生产产品 A A A和产品 B B B。例如,如果要生产 A A A,我们就需要准备生产 A A A所需的设备等,这些活动是需要花费一些成本的,我们可以认为这些成本是固定的,即固定成本或者启动成本,固定成本可以认为就是需要去购置机器,雇佣人员等。这个决策,可以看做是上层决策(tactical level),很明显,这个决策并没有涉及具体生产多少产品的决策。

    在决定了要生产 A A A以后,也准备好了所有的所需材料和人工,我们基于这些材料,再去制定具体要生产多少商品的生产计划。这个生产计划,就是非常具体的决策,可以看做是作业层决策(operational level), 即下层决策。

    假设决定生产产品 A A A和产品 B B B的固定成本分别为 f A f_A fA f B f_B fB。我们引入上层决策(tactical level)的决策变量:
    z i ∈ { 0 , 1 } , i = A  or  B \begin{aligned} z_i \in \{ 0, 1\}, i = A \text{ or } B \end{aligned} zi{0,1},i=A or B
    表示我们决定是否生产产品 A A A或者 B B B

    之后,具体操作层的决策,就是每种产品生产多少,即为 x i k x_i^k xik,在第 k k k个场景下,生产产品 i i i的数量。

    因此,Two-stage Stochastic Programming的模型可以写成

    max ⁡    ∑ k ∈ K p k ( ∑ i ∈ I c i x i k ) − ∑ i ∈ I f i z i s . t . ∑ i ∈ I a j i k x i k ⩽ b j , ∀ j ∈ J , ∀ k ∈ K ∑ k ∈ K x i k − M z i ⩽ 0 , ∀ i ∈ I z i ∈ { 0 , 1 } , x i k ⩾ 0 , ∀ i ∈ I , ∀ k ∈ K \begin{aligned} \max \,\,&\sum_{k\in K}{p_k}\left( \sum_{i\in I}{c_i}x_{i}^{k} \right) -\sum_{i\in I}{f_i}z_i \\ s.t.\quad &\sum_{i\in I}{a_{ji}^{k}x_{i}^{k}}\leqslant b_j, \quad \forall j\in J,\forall k\in K \\ &\sum_{k\in K}{x_{i}^{k}}-Mz_i\leqslant 0, \quad \forall i\in I \\ &z_i\in \left\{ 0,1 \right\} ,x_{i}^{k}\geqslant 0, \quad \forall i\in I,\forall k\in K\end{aligned} maxs.t.kKpk(iIcixik)iIfiziiIajikxikbj,jJ,kKkKxikMzi0,iIzi{0,1},xik0,iI,kK
    第2条约束表示,只要有一个场景下, x i k > 0 x_i^k >0 xik>0,则 z i = 1 z_i=1 zi=1。也就是,只要有一种情况下,生产了产品 i ∈ { A , B } i \in \{A, B\} i{A,B},则我们就需要在一开始就要决定生产产品 i ∈ { A , B } i \in \{A, B\} i{A,B},且准备好生产 i ∈ { A , B } i \in \{A, B\} i{A,B}的所有预备材料。

    上面的模型虽然包含了两阶段决策:上层决策(tactical level)作业层决策(operational level)。但是模型看上去看是单阶段的模型。不过上面的写法是没问题的,只是没有体现出两阶段的思想。

    Two-stage Stochastic Programming的更学术的写法

    为了体现出两阶段的思想,我们换一种更专业,更学术的写法,学术论文里大家一般也都是按照下面的方法写的。

    我们按照事情发生的先后顺序的逻辑来写。

    首先,我们需要做一个上层决策(tactical level),这个决策即为:要不要生产产品 A A A和产品 B B B。做出的决策必须要使得总的期望净利润最大化。即
    max ⁡ E p [ h ( z , K ) ] − ∑ i ∈ I f i z i z i ∈ { 0 , 1 } , ∀ i ∈ I \begin{aligned} \max \quad &\mathbb{E}_p\left[ h\left( \mathbf{z},K \right) \right] -\sum_{i\in I}{f_i}z_i \\ &z_i\in \left\{ 0,1 \right\} ,\qquad \forall i\in I \end{aligned} maxEp[h(z,K)]iIfizizi{0,1},iI
    这里,符号上可能有些难理解,我在这里详细的解释一下。 E p [ h ( z , K ) ] \mathbb{E}_p\left[ h\left( \mathbf{z},K \right) \right] Ep[h(z,K)],就是在各个场景出现的概率为 p p p时(例如 p = [ 0.2 , 0.2 , 0.2 , 0.2 , 0.2 ] p = [0.2, 0.2, 0.2, 0.2, 0.2] p=[0.2,0.2,0.2,0.2,0.2]),决策变量 z = [ z A , z B ] T \mathbf{z}=[z_A, z_B]^T z=[zA,zB]T的取值确定以后,在所有场景 K K K下的收入的期望值。实际上,就是第二阶段决策产生的收入的期望值,但是现在我们还没有决策第二阶段的生产计划,因此我们先用一个符号去表示第二阶段的期望收入值

    当我们确定了是否要生产产品 A A A B B B的决策 z = [ z A , z B ] T \mathbf{z}=[z_A, z_B]^T z=[zA,zB]T以后,第二阶段,也就是制定生产计划的阶段,我们会面临不确定的情况,也就是面临多种场景 ∀ k ∈ K \forall k \in K kK。对于每一种场景 k ∈ K k \in K kK当这种场景 k k k真的发生了的时候,我们需要做出在这种场景下的最优的生产计划。最优的生产计划可以用下面的模型规划出来

    注意,之前是 h ( z , K ) h\left( \mathbf{z},K \right) h(z,K), 是指在所有的场景 K K K下的所有决策,这里我们固定了某一种场景 k k k因此变成了 h ( z , k ) h\left( \mathbf{z},k \right) h(z,k)。此时 z \mathbf{z} z已经可以看做是一个给定的值,因为在做第二阶段生产计划的时候,决定生产还是不生产 A A A B B B是已经知道的事实。

    h ( z , k ) = max ⁡ x ( k )    ∑ i ∈ I c i x i k s . t . ∑ i ∈ I a j i k x i k ⩽ b j , ∀ j ∈ J , ∀ k ∈ K ∑ k ∈ K ∑ i ∈ I x i k − M z i ⩽ 0 , ∀ i ∈ I x i k ⩾ 0 , ∀ i ∈ I , ∀ k ∈ K \begin{aligned} h\left( \mathbf{z},k \right) &=\underset{\mathbf{x}\left( k \right)}{\max}\,\,\sum_{i\in I}{c_i}x_{i}^{k} \\ s.t.&\quad \sum_{i\in I}{a_{ji}^{k}x_{i}^{k}}\leqslant b_j, \qquad \quad \forall j\in J,\forall k\in K \\ &\quad \sum_{k\in K}{\sum_{i\in I}{x_{i}^{k}}}-Mz_i\leqslant 0, \quad \forall i\in I \\ &\quad x_{i}^{k}\geqslant 0,\quad \qquad \qquad \qquad \forall i\in I,\forall k\in K \end{aligned} h(z,k)s.t.=x(k)maxiIcixikiIajikxikbj,jJ,kKkKiIxikMzi0,iIxik0,iI,kK

    怎么样,这样写,是不是就比较好理解两个阶段的决策上层决策(tactical level)作业层决策(operational level)之间是什么关系了。

    我们来汇总一下上面的两阶段模型。

    首先,第一阶段,作为企业的决策者,我们需要作出是否生产产品 A A A B B B的上层决策(tactical level decisions),使得总的期望净利润最大化。这一阶段的决策,是和随机场景无关的,无论场景怎么变,这一阶段的决策产生的固定成本都是相同的。模型为
    max ⁡ E p [ h ( z , K ) ] − ∑ i ∈ I f i z i z i ∈ { 0 , 1 } , ∀ i ∈ I \begin{aligned} \max \quad &\mathbb{E}_p\left[ h\left( \mathbf{z},K \right) \right] -\sum_{i\in I}{f_i}z_i \\ &z_i\in \left\{ 0,1 \right\} ,\qquad \forall i\in I \end{aligned} maxEp[h(z,K)]iIfizizi{0,1},iI
    制定好了上层决策之后,真正到了要生产的时候,我们观察到某种可能的场景 k ∈ K k \in K kK发生了,我们需要作出该场景下的最优生产决策,即:
    第二阶段的决策:亦即下层决策(operational level decisions),对于场景 k ∈ K k \in K kK,我们有
    h ( z , k ) = max ⁡ x ( k )    ∑ i ∈ I c i x i k s . t . ∑ i ∈ I a j i k x i k ⩽ b j , ∀ j ∈ J , ∀ k ∈ K ∑ k ∈ K ∑ i ∈ I x i k − M z i ⩽ 0 , ∀ i ∈ I x i k ⩾ 0 , ∀ i ∈ I , ∀ k ∈ K \begin{aligned} h\left( \mathbf{z},k \right) &=\underset{\mathbf{x}\left( k \right)}{\max}\,\,\sum_{i\in I}{c_i}x_{i}^{k} \\ s.t.&\quad \sum_{i\in I}{a_{ji}^{k}x_{i}^{k}}\leqslant b_j, \quad \forall j\in J,\forall k\in K \\ &\quad \sum_{k\in K}{\sum_{i\in I}{x_{i}^{k}}}-Mz_i\leqslant 0, \qquad \forall i\in I \\ &\quad x_{i}^{k}\geqslant 0,\quad \quad \quad \forall i\in I,\forall k\in K \end{aligned} h(z,k)s.t.=x(k)maxiIcixikiIajikxikbj,jJ,kKkKiIxikMzi0,iIxik0,iI,kK
    怎么样,这么写是不是高大上一些。以后小伙伴们写论文,大概都需要这么写的。这样比较专业一些。如果按照之前的那种写法,可能会被喷不专业。哈哈哈。

    两阶段的写法比较容易理解两阶段决策的逻辑,但在真正求解的时候,是和前面那种写法是完全一样的。为了方便大家查看,我再抄一遍过来。

    合并书写的方式 max ⁡    ∑ k ∈ K p k ( ∑ i ∈ I c i x i k ) − ∑ i ∈ I f i z i s . t . ∑ i ∈ I a j i k x i k ⩽ b j , ∀ j ∈ J , ∀ k ∈ K ∑ k ∈ K x i k − M z i ⩽ 0 , ∀ i ∈ I z i ∈ { 0 , 1 } , x i k ⩾ 0 , ∀ i ∈ I , ∀ k ∈ K \begin{aligned} \max \,\,&\sum_{k\in K}{p_k}\left( \sum_{i\in I}{c_i}x_{i}^{k} \right) -\sum_{i\in I}{f_i}z_i \\ s.t.\quad &\sum_{i\in I}{a_{ji}^{k}x_{i}^{k}}\leqslant b_j, \quad \forall j\in J,\forall k\in K \\ &\sum_{k\in K}{x_{i}^{k}}-Mz_i\leqslant 0, \quad \forall i\in I \\ &z_i\in \left\{ 0,1 \right\} ,x_{i}^{k}\geqslant 0, \quad \forall i\in I,\forall k\in K\end{aligned} maxs.t.kKpk(iIcixik)iIfiziiIajikxikbj,jJ,kKkKxikMzi0,iIzi{0,1},xik0,iI,kK
    我们真正进行求解的时候,还是用这种写法去求解就可以。还是那句话,求解器根本不知道你的决策变量哪个是上层决策,哪个是下层决策。当你把模型丢给求解器,人家直接就求解了,管你什么一阶段二阶段呢,哈哈哈。除非你自己根据一阶段二阶段设计了专门的算法,先固定哪个,后做哪个。比如说,你用Benders Decomposition,把第一阶段的z作为复杂变量,把第二阶段的x作为简单变量,那就另当别论。

    至于随机规划更为复杂的例子,大家可以去读我的母系,工业工程系毕业的黄一潇师兄2017年发在TRB上的论文Time-dependent vehicle routing problem with path flexibility

    下面我们继续用Python调用Gurobi来求解上述Two-stage Stochastic Programming

    Python调用Gurobi求解Two-stage Stochastic Programming

    我们沿用前一篇文章的所有关于均值和方差的假设,即
    E ( a ~ 11 ) = μ 11 = 1 , σ 11 = 0.2 E ( a ~ 12 ) = μ 12 = 2 , σ 12 = 0.2 E ( a ~ 21 ) = μ 21 = 4 , σ 21 = 0.2 E ( a ~ 32 ) = μ 32 = 4 , σ 32 = 0.2 \begin{aligned} \mathbb{E}(\tilde{a}_{11}) = \mu_{11} = 1, \quad \sigma_{11} = 0.2 \\ \mathbb{E}(\tilde{a}_{12}) = \mu_{12} = 2, \quad \sigma_{12} = 0.2 \\ \mathbb{E}(\tilde{a}_{21}) = \mu_{21} = 4, \quad \sigma_{21} = 0.2 \\ \mathbb{E}(\tilde{a}_{32}) = \mu_{32} = 4, \quad \sigma_{32} = 0.2 \end{aligned} E(a~11)=μ11=1,σ11=0.2E(a~12)=μ12=2,σ12=0.2E(a~21)=μ21=4,σ21=0.2E(a~32)=μ32=4,σ32=0.2
    我们设置 f A = 6 , f B = 4 f_A = 6, f_B=4 fA=6,fB=4
    我们用Python调用Gurobi对其进行求解。

    为了方便与只有操作层决策的版本进行对比,我们首先将不确定参数固定住,即保证二者使用同一组约束系数。

    scenario_num = 5  
    # generate parameters under scenarios 
    a_11 = np.random.normal(1, 0.2, scenario_num)  # mu, sigma, sample_number 
    a_12 = np.random.normal(2, 0.2, scenario_num)
    a_21 = np.random.normal(4, 0.2, scenario_num)
    a_32 = np.random.normal(4, 0.2, scenario_num) 
    print('a_11:', a_11)
    print('a_12:', a_12)
    print('a_21:', a_21)
    print('a_32:', a_32) 
    

    产生的随机输入参数为

    a_11:  [0.90039351 1.38590641 1.18988416 1.01751025 0.7549129 ]
    a_12:  [2.1688726  1.79995693 1.69104578 2.23760596 2.06338852]
    a_21:  [4.18417176 4.06374553 4.17136612 3.86979488 3.79315143]
    a_32:  [4.1363189  3.83931807 3.86209004 3.9088935  4.00349583]
    

    Two-stage Stochastic Programming

    from gurobipy import * 
    
    
    sto_model = Model('Stochastic Programming')
    prob = [0.2, 0.2, 0.2, 0.2, 0.2] 
    # prob = [0.1, 0.2, 0.3, 0.15, 0.25]
    profit = [2, 3]
    fix_cost = [6, 4] 
    b = [8, 16, 12] 
    big_M = 1000 
    
    
    x_sto = {}
    z = {} 
    for i in range(2): 
        # creat variables  
        z[i] = sto_model.addVar(lb = 0, ub = 1, vtype = GRB.BINARY, name = 'z_' + str(i) + '_' + str(i)) 
        for s in range(scenario_num):
            x_sto[i, s] = sto_model.addVar(lb = 0, ub = GRB.INFINITY, vtype = GRB.CONTINUOUS, name = 'x_' + str(i) + '_' + str(s))
        
    # set objective 
    obj = LinExpr(0)
    for key in x_sto.keys():
        product_ID = key[0]
        scenario_ID = key[1] 
        obj.addTerms(prob[scenario_ID] * profit[product_ID], x_sto[key])
    for key in z.keys():
        obj.addTerms(-fix_cost[key], z[key])  
        
    sto_model.setObjective(obj, GRB.MAXIMIZE)
    
    # add constraints:
    # constraints 1-1
    for s in range(scenario_num):
        lhs = LinExpr(0)
        lhs.addTerms(a_11[s], x_sto[0, s]) 
        lhs.addTerms(a_12[s], x_sto[1, s])  
        sto_model.addConstr(lhs <= b[0], name = 'cons_' + str(0))  
        
    # constraints 1-2      
    for s in range(scenario_num):
        lhs = LinExpr(0)
        lhs.addTerms(a_21[s], x_sto[0, s]) 
        sto_model.addConstr(lhs <= b[1], name = 'cons_' + str(1))   
        
    # constraints 1-3      
    for s in range(scenario_num):
        lhs = LinExpr(0)
        lhs.addTerms(a_32[s], x_sto[1, s])  
        sto_model.addConstr(lhs <= b[2], name = 'cons_' + str(2))     
    
    # constraints 2 
    for i in range(2):
        lhs = LinExpr(0)
        for s in range(scenario_num):
            lhs.addTerms(1, x_sto[i, s])
        sto_model.addConstr(lhs <= big_M * z[i]) 
    
    # solve 
    sto_model.optimize()
    
    print('Obj = ', sto_model.ObjVal)
    print('\n\n  ------  tactical level decision  ------- ') 
    for key in z.keys():
        print(z[key].varName, ' = ', z[key].x)   
    print('\n\n ------  operational level decision  ------- ')
    for s in range(scenario_num): 
        print(' ------  scenario   ', s, '  ------- ') 
        for i in range(2):
            print(x_sto[i,s].varName, ' = ', x_sto[i,s].x)  
        print('\n\n') 
    

    求解结果为

     -- Obj =  5.120668471999894
    
    
      ------  tactical level decision  ------- 
    z_0_0  =  0.0
    z_1_1  =  1.0
    
    
     ------  operational level decision  ------- 
     ------  scenario    0   ------- 
    x_0_0  =  0.0
    x_1_0  =  2.901130275374204
    
    
    
     ------  scenario    1   ------- 
    x_0_1  =  0.0
    x_1_1  =  3.1255550569323436
    
    
    
     ------  scenario    2   ------- 
    x_0_2  =  0.0
    x_1_2  =  3.107125898642543
    
    
    
     ------  scenario    3   ------- 
    x_0_3  =  0.0
    x_1_3  =  3.069922473497831
    
    
    
     ------  scenario    4   ------- 
    x_0_4  =  0.0
    x_1_4  =  2.997380415552898
    

    可以看到,我们最终是决定,不生产产品 A A A,只生产产品 B B B,因为产品 A A A的固定成本高( f A = 6 f_A=6 fA=6),但是利润比较低(为2)。而产品 B B B的固定成本低( f A = 4 f_A=4 fA=4),而且利润高(为3)。

    接下来,我们就对比一下之前那种只有操作层的版本。

    Single stage version

    scenario_num = 5 
    sto_model = Model('Stochastic Programming')
    prob = [0.2, 0.2, 0.2, 0.2, 0.2] 
    # prob = [0.1, 0.2, 0.3, 0.15, 0.25]
    profit = [2, 3] 
    b = [8, 16, 12] 
    
    # generate parameters under scenarios 
    # a_11 = np.random.normal(1, 0.2, scenario_num)  # mu, sigma, sample_number 
    # a_12 = np.random.normal(2, 0.2, scenario_num)
    # a_21 = np.random.normal(4, 0.2, scenario_num)
    # a_32 = np.random.normal(4, 0.2, scenario_num) 
        
    x_sto = {}
    for i in range(2): 
        # creat variables 
        for s in range(scenario_num):
            x_sto[i, s] = sto_model.addVar(lb = 0, ub = GRB.INFINITY, vtype = GRB.CONTINUOUS, name = 'x_' + str(i) + '_' + str(s))
        
    # set objective 
    obj = LinExpr(0)
    for key in x_sto.keys():
        product_ID = key[0]
        scenario_ID = key[1] 
        obj.addTerms(prob[scenario_ID] * profit[product_ID], x_sto[key])
    sto_model.setObjective(obj, GRB.MAXIMIZE)
    
    # add constraints:
    # constraints 1 
    for s in range(scenario_num):
        lhs = LinExpr(0)
        lhs.addTerms(a_11[s], x_sto[0, s]) 
        lhs.addTerms(a_12[s], x_sto[1, s])  
        sto_model.addConstr(lhs <= b[0], name = 'cons_' + str(0))  
        
    # constraints 2      
    for s in range(scenario_num):
        lhs = LinExpr(0)
        lhs.addTerms(a_21[s], x_sto[0, s]) 
        sto_model.addConstr(lhs <= b[1], name = 'cons_' + str(1))   
        
    # constraints 3      
    for s in range(scenario_num):
        lhs = LinExpr(0)
        lhs.addTerms(a_32[s], x_sto[1, s])  
        sto_model.addConstr(lhs <= b[2], name = 'cons_' + str(2))     
    
    # solve 
    sto_model.optimize()
    
    print('Obj = ', sto_model.ObjVal)
    # for key in x_sto.keys():
    #     print(x_sto[key].varName, ' = ', x_sto[key].x)   
    for s in range(scenario_num): 
        print(' ------  scenario   ', s, '  ------- ') 
        for i in range(2):
            print(x_sto[i,s].varName, ' = ', x_sto[i,s].x)  
        print('\n\n') 
    

    求解结果如下

     --- Obj =  13.896545328523917
     ------  scenario    0   ------- 
    x_0_0  =  3.8239347951176703
    x_1_0  =  2.101070361744041
    
    
    
     ------  scenario    1   ------- 
    x_0_1  =  1.71305615957035
    x_1_1  =  3.1255550569323436
    
    
    
     ------  scenario    2   ------- 
    x_0_2  =  2.307542152652112
    x_1_2  =  3.107125898642543
    
    
    
     ------  scenario    3   ------- 
    x_0_3  =  4.134586067378273
    x_1_3  =  1.6951225436724826
    
    
    
     ------  scenario    4   ------- 
    x_0_4  =  4.218128458180001
    x_1_4  =  2.333869931282845
    

    是吧,这两个的结果还是完全不一样的吧。

    此时,由于上层决策变量 z i , ∀ i ∈ I z_i, \forall i \in I zi,iI和下层决策变量 x i k , ∀ i ∈ I , k ∈ K x_i^k, \forall i \in I, k \in K xik,iI,kK是有耦合关系的,是一揽子决策。场景和场景之间也是有一些关系的,比如在场景1下,我们生产了产品 A A A,即 x A 1 > 0 x_A^1 >0 xA1>0,那么 z A z_A zA就一定为1,因此在场景2下,也许就会生产 x A 2 x_A^2 xA2也许就大于0了。但是,也许你单独求解多个模型的时候, x A 2 = 0 x_A^2 = 0 xA2=0。这种情形是有可能会出现的。在这种情况下,单独求解每个场景下的模型,然后加权求和一下,就和随机规划的模型的结果是不等价的。

    所以说,随机规划肯定不是单独求解多个线性规划,然后加权求和一下。它是有自己独特的强大功能在的。

    好啦,这篇文章就是完整的给大家展示了两阶段随机规划的一个小例子,从构建决策,到建模,到代码实现。相信大家读完这篇文章,应该对随机规划有了一个非常清楚的初步认识。之后如果想要继续深入,欢迎阅读下面的两本书。或者持续关注我们的公众号。

    参考文献

    [1]: Shapiro, A and Ruszczynski, A, Handbooks in operations research and management science: Vol. 10. Stochastic programming, 2003
    [2]: Birge, John R., and Francois Louveaux. Introduction to stochastic programming. Springer Science & Business Media, 2011.
    [3]: Huang, Y., Zhao, L., Van Woensel, T., & Gross, J. P. (2017). Time-dependent vehicle routing problem with path flexibility. Transportation Research Part B: Methodological, 95, 169-195.


    作者:刘兴禄,清华大学,清华伯克利深圳学院,博士在读
    联系方式:hsinglul@163.com
    博客主页:https://blog.csdn.net/HsinglukLiu/article/details/107848461

    OlittleRer

    运小筹公众号是致力于分享运筹优化(LP、MIP、NLP、随机规划、鲁棒优化)、凸优化、强化学习等研究领域的内容以及涉及到的算法的代码实现。编程语言和工具包括Java、Python、Matlab、CPLEX、Gurobi、SCIP 等。

    关注我们: 运筹小公众号

    在这里插入图片描述

    也欢迎广大同行投稿!欢迎大家关注我们的公众号“运小筹”以及粉丝群!

    微信群:
    在这里插入图片描述

    QQ群:

    QQ群里有我们共享的一些书籍、论文、算例等资料,欢迎大家加入。

    在这里插入图片描述

    展开全文
  • 数组的完全随机排列算法

    千次阅读 2016-12-21 16:46:04
    数组的完全随机排列算法 Array.prototype.sort 方法被许多 JavaScript 程序员误用来随机排列数组。最近做的前端星计划挑战项目中,一道实现 blackjack 游戏的问题,就发现很多同学使用了 Array....

    数组的完全随机排列算法

    Array.prototype.sort 方法被许多 JavaScript 程序员误用来随机排列数组。最近做的前端星计划挑战项目中,一道实现 blackjack 游戏的问题,就发现很多同学使用了 Array.prototype.sort 来洗牌。就连最近一期 JavaScript Weekly上推荐的一篇文章也犯了同样的错误。

    以下就是常见的完全错误的随机排列算法:

    function shuffle(arr){
        return arr.sort(function(){ return Math.random() - 0.5; }); } 

    以上代码看似巧妙利用了 Array.prototype.sort 实现随机,但是,却有非常严重的问题,甚至是完全错误。

    证明 Array.prototype.sort 随机算法的错误

    为了证明这个算法的错误,我们设计一个测试的方法。假定这个排序算法是正确的,那么,将这个算法用于随机数组 [0, 1, 2, 3, 4, 5, 6, 7, 8, 9],如果算法正确,那么每个数字在每一位出现的概率均等。因此,将数组重复洗牌足够多次,然后将每次的结果在每一位相加,最后对每一位的结果取平均值,这个平均值应该约等于 (0 + 9) / 2 = 4.5,测试次数越多次,每一位上的平均值就都应该越接近于 4.5。所以我们简单实现测试代码如下:

    var arr = [0,1,2,3,4,5,6,7,8,9]; var res = [0,0,0,0,0,0,0,0,0,0]; var t = 10000; for(var i = 0; i < t; i++){ var sorted = shuffle(arr.slice(0)); sorted.forEach(function(o,i){ res[i] += o; }); } res = res.map(function(o){ return o / t; }); console.log(res); 

    将上面的 shuffle 方法用这段测试代码在 chrome 浏览器中测试一下,可以得出结果,发现结果并不随机分布,各个位置的平均值越往后越大,这意味着这种随机算法越大的数字出现在越后面的概率越大。

    为什么会产生这个结果呢?我们需要了解 Array.prototype.sort 究竟是怎么作用的。

    首先我们知道排序算法有很多种,而 ECMAScript 并没有规定 Array.prototype.sort 必须使用何种排序算法。在这里,有兴趣的同学不妨看一下 JavaScriptCore 的源码实现:

    排序不是我们今天讨论的主题,但是不论用何种排序算法,都是需要进行两个数之间的比较和交换,排序算法的效率和两个数之间比较和交换的次数有关系。

    最基础的排序有冒泡排序和插入排序,原版的冒泡或者插入排序都比较了 n(n-1)/2 次,也就是说任意两个位置的元素都进行了一次比较。那么在这种情况下,如果采用前面的 sort 随机算法,由于每次比较都有 50% 的几率交换和不交换,这样的结果是随机均匀的吗?我们可以看一下例子

    function bubbleSort(arr, compare){
      var len = arr.length; for(var i = 0; i < len - 1; i++){ for(var j = 0; j < len - 1 - i; j++){ var k = j + 1; if(compare(arr[j], arr[k]) > 0){ var tmp = arr[j]; arr[j] = arr[k]; arr[k] = tmp; } } } return arr; } function shuffle(arr){ return bubbleSort(arr, function(){ return Math.random() - 0.5; }); } var arr = [0,1,2,3,4,5,6,7,8,9]; var res = [0,0,0,0,0,0,0,0,0,0]; var t = 10000; for(var i = 0; i < t; i++){ var sorted = shuffle(arr.slice(0)); sorted.forEach(function(o,i){ res[i] += o; }); } res = res.map(function(o){ return o / t; }); console.log(res); 

    上面的代码的随机结果也是不均匀的,测试平均值的结果越往后的越大。(笔者之前没有复制原数组所以错误得出均匀的结论,已更正于 2016-05-10)

    冒泡排序总是将比较结果较小的元素与它的前一个元素交换,我们可以大约思考一下,这个算法越后面的元素,交换到越前的位置的概率越小(因为每次只有50%几率“冒泡”),原始数组是顺序从小到大排序的,因此测试平均值的结果自然就是越往后的越大(因为越靠后的大数出现在前面的概率越小)。

    我们再换一种算法,我们这一次用插入排序

    function insertionSort(arr, compare){
      var len = arr.length; for(var i = 0; i < len; i++){ for(var j = i + 1; j < len; j++){ if(compare(arr[i], arr[j]) > 0){ var tmp = arr[i]; arr[i] = arr[j]; arr[j] = tmp; } } } return arr; } function shuffle(arr){ return insertionSort(arr, function(){ return Math.random() - 0.5; }); } var arr = [0,1,2,3,4,5,6,7,8,9]; var res = [0,0,0,0,0,0,0,0,0,0]; var t = 10000; for(var i = 0; i < t; i++){ var sorted = shuffle(arr.slice(0)); sorted.forEach(function(o,i){ res[i] += o; }); } res = res.map(function(o){ return o / t; }); console.log(res); 

    由于插入排序找后面的大数与前面的数进行交换,这一次的结果和冒泡排序相反,测试平均值的结果自然就是越往后越小。原因也和上面类似,对于插入排序,越往后的数字越容易随机交换到前面。

    所以我们看到即使是两两交换的排序算法,随机分布差别也是比较大。除了每个位置两两都比较一次的这种排序算法外,大多数排序算法的时间复杂度介于 O(n) 到 O(n2) 之间,元素之间的比较次数通常情况下要远小于 n(n-1)/2,也就意味着有一些元素之间根本就没机会相比较(也就没有了随机交换的可能),这些 sort 随机排序的算法自然也不能真正随机。

    我们将上面的代码改一下,采用快速排序

    function quickSort(arr, compare){
      arr = arr.slice(0); if(arr.length <= 1) return arr; var mid = arr[0], rest = arr.slice(1); var left = [], right = []; for(var i = 0; i < rest.length; i++){ if(compare(rest[i], mid) > 0){ right.push(rest[i]); }else{ left.push(rest[i]); } } return quickSort(left, compare).concat([mid]) .concat(quickSort(right, compare)); } function shuffle(arr){ return quickSort(arr, function(){ return Math.random() - 0.5; }); } var arr = [0,1,2,3,4,5,6,7,8,9]; var res = [0,0,0,0,0,0,0,0,0,0]; var t = 10000; for(var i = 0; i < t; i++){ var sorted = shuffle(arr.slice(0)); sorted.forEach(function(o,i){ res[i] += o; }); } res = res.map(function(o){ return o / t; }); console.log(res); 

    快速排序并没有两两元素进行比较,它的概率分布也不随机。

    所以我们可以得出结论,用 Array.prototype.sort 随机交换的方式来随机排列数组,得到的结果并不一定随机,而是取决于排序算法是如何实现的,用 JavaScript 内置的排序算法这么排序,通常肯定是不完全随机的。

    经典的随机排列

    所有空间复杂度 O(1) 的排序算法的时间复杂度都介于 O(nlogn) 到 O(n2) 之间,因此在不考虑算法结果错误的前提下,使用排序来随机交换也是慢的。事实上,随机排列数组元素有经典的 O(n) 复杂度的算法:

    function shuffle(arr){
      var len = arr.length; for(var i = 0; i < len - 1; i++){ var idx = Math.floor(Math.random() * (len - i)); var temp = arr[idx]; arr[idx] = arr[len - i - 1]; arr[len - i -1] = temp; } return arr; } 

    在上面的算法里,我们每一次循环从前 len - i 个元素里随机一个位置,将这个元素和第 len - i 个元素进行交换,迭代直到 i = len - 1 为止。

    我们同样可以检验一下这个算法的随机性

    function shuffle(arr){
      var len = arr.length; for(var i = 0; i < len - 1; i++){ var idx = Math.floor(Math.random() * (len - i)); var temp = arr[idx]; arr[idx] = arr[len - i - 1]; arr[len - i -1] = temp; } return arr; } var arr = [0,1,2,3,4,5,6,7,8,9]; var res = [0,0,0,0,0,0,0,0,0,0]; var t = 10000; for(var i = 0; i < t; i++){ var sorted = shuffle(arr.slice(0)); sorted.forEach(function(o,i){ res[i] += o; }); } res = res.map(function(o){ return o / t; }); console.log(res); 

    从结果可以看出这个算法的随机结果应该是均匀的。不过我们的测试方法其实有个小小的问题,我们只测试了平均值,实际上平均值接近只是均匀分布的必要而非充分条件,平均值接近不一定就是均匀分布。不过别担心,事实上我们可以简单从数学上证明这个算法的随机性。

    随机性的数学归纳法证明

    对 n 个数进行随机:

    1. 首先我们考虑 n = 2 的情况,根据算法,显然有 1/2 的概率两个数交换,有 1/2 的概率两个数不交换,因此对 n = 2 的情况,元素出现在每个位置的概率都是 1/2,满足随机性要求。

    2. 假设有 i 个数, i >= 2 时,算法随机性符合要求,即每个数出现在 i 个位置上每个位置的概率都是 1/i。

    3. 对于 i + 1 个数,按照我们的算法,在第一次循环时,每个数都有 1/(i+1) 的概率被交换到最末尾,所以每个元素出现在最末一位的概率都是 1/(i+1) 。而每个数也都有 i/(i+1) 的概率不被交换到最末尾,如果不被交换,从第二次循环开始还原成 i 个数随机,根据 2. 的假设,它们出现在 i 个位置的概率是 1/i。因此每个数出现在前 i 位任意一位的概率是 (i/(i+1)) * (1/i) = 1/(i+1),也是 1/(i+1)。

    4. 综合 1. 2. 3. 得出,对于任意 n >= 2,经过这个算法,每个元素出现在 n 个位置任意一个位置的概率都是 1/n。

    总结

    一个优秀的算法要同时满足结果正确和高效率。很不幸使用 Array.prototype.sort 方法这两个条件都不满足。因此,当我们需要实现类似洗牌的功能的时候,还是应该采用巧妙的经典洗牌算法,它不仅仅具有完全随机性还有很高的效率。

    除了收获这样的算法之外,我们还应该认真对待这种动手分析和解决问题的思路,并且捡起我们曾经学过而被大多数人遗忘的数学(比如数学归纳法这种经典的证明方法)。

    有任何问题欢迎与作者探讨~

    展开全文
  • 目录 ...生活中有许多例子是牛顿力学基本能完全描述的:苹果落地,炮弹轨迹等。这些简单的过程可以由微分方程和初始条件唯一确定。比如知道了炮弹的初始速度和角度,那么后续的轨迹就能用牛顿第...

    目录

    一杯咖啡是随机的吗?

    生活中有许多例子是牛顿力学基本能完全描述的:苹果落地,炮弹轨迹等。这些简单的过程可以由微分方程和初始条件唯一确定。比如知道了炮弹的初始速度和角度,那么后续的轨迹就能用牛顿第二定律唯一解出来。但是生活中还有一些其他的例子,虽然也由牛顿力学控制,但是由于自由度太大,需要描述的参数实在太多,以致理论上的描述和预测根本行不通。自由度很大的系统有哪些?一小勺咖啡里的分子数量就足以说明问题:约有 N = 1 0 23 N=10^{23} N=1023(阿伏伽德罗常数)个分子,每个分子有6个自由度(x,y,z方向的位置和速度),每个自由度由一个微分方程控制,用目前的计算机解这么多个方程(6N个)组成的微分方程组是不可能完成的任务。

    图1. 一杯咖啡里的分子(包括水分子和咖啡分子)个数约为阿伏伽德罗常数。图片摘自http://toddyoungblood.com/。

    怎么解决这个问题呢?仔细考察这样的系统会发现,“随机性"承担了很重要的角色。这里的“随机性"可以从不同角度来理解:首先,在这么多粒子的系统中,我们不可能精确知道每个粒子的信息(位置和速度),这带来了一定的随机性;其次,这个系统可能和外界热源接触或有复杂的相互作用,用统计力学的概率分布和系综平均的方法来描述更合理。因此,即使这一小勺咖啡里的分子都受牛顿力学控制,对于信息的缺失迫使我们不得不使用统计力学描述(“随机性"描述)。实际上,这样的描述是现实世界很好的近似,而且与实验结果能很好符合。

    我们就拿咖啡在水中的扩散作为例子。从宏观上看,当我们把咖啡磨成的粉稳稳地倒进一杯热水中时(不人为搅拌),可以看到咖啡粉不断地向外扩散,直到整杯水都变成褐色。从微观上来看,当我们把一堆咖啡分子(严格来说应该是块状小碎屑,我们不可能把咖啡磨到分子层面,但为了说明方便,暂且用这个叫法)倒进水中后,如果我们能站在其中某个分子上跟着它运动,就会看到如下一幕:这个分子在水中走的轨迹很杂乱,一会儿撞到其他水分子或咖啡分子改变了运动方向和速度,走了一段直线距离(牛顿第一定律)后又撞到另一个分子,又改变了运动方向和速度……前面已经解释过,如果用牛顿运动方程组来求解每个粒子的运动信息,目前的计算能力根本不可能完成;即使能够求解,我们得到的也只是每个粒子的位置和速度随时间的变化而已,对于这个扩散过程本身并不能带来深刻的理解。什么意思呢?举个简单的例子,假设我们求解出了咖啡粉在100ml水中均匀扩散开所需的时间,那么如果想知道在200ml的水中完全扩散的时间,是不是还得重新计算一遍?多麻烦呀!没关系,这个咖啡杯问题我们马上就能解答。

    图2. 一团咖啡分子(红色)在水中的溶解,注意到红色粒子之间会相互碰撞改变速度方向和大小。不过这个动画描述的过于简单,没有画出水分子,也没有展现水分子和咖啡分子碰撞的效应。图片摘自http://cronodon.com/BioTech/Diffusion.html。

    换个角度看待上述不可能求解的方程:受牛顿力学控制的方程组自由度那么大,其实里面包含的信息是冗余的;我们没必要知道那么多信息,只要用一种统计的描述就够了,也就是对所有分子轨迹的平均行为进行描述即可。因为粒子数足够多,这个求平均的操作也就能很好定义。

    现在我们到了最关键的一步:既然我们只想知道平均意义下的系统的行为,我们就可以对真实系统建模,用一个随机过程来近似一个(受牛顿力学控制的)确定的过程,建模的要求是这两个过程具有相同的平均行为。这里的“确定的过程"指的是上述庞大的牛顿方程组的确定解。具体来说,我们把每个粒子不断受撞击而产生的复杂轨迹建模成一个“随机游走"过程,就像深夜一个醉汉在马路上晃悠一样,没有别人推搡,自己在那乱走;而且我们假设这种粒子的随机“乱走"在行为表现上等效于真实的不断受到撞击的走法。这种近似在多大程度上是合理的,我们在后面会提到。

    图3. 随机行走(random walk),或醉汉行走。图片摘自George Gamov, 1961。
    图4. 800个粒子构成的热力学系统的模拟,蓝线为其中5个粒子的轨迹,红色为其中一个粒子的速度向量。当系统里的粒子数足够大时,单个粒子的行为可以用随机游走模型来近似。图片摘自Wikipedia:Brownian motion。

    从上面的讨论中,我们已经差不多达成了一致:可以用随机游走模型来等效描述咖啡在水中的扩散行为。**一杯咖啡是确定的,但是为了理解它的行为,我们通过建立模型引入了随机性。**现在我们需要讲一点硬核的物理了。

    随机游走与扩散

    限于篇幅,我们在这里只讨论简单的情况:一维固定步长的无漂移随机游走(取那么长的名字是为了更严谨)。假设每个粒子只能沿着一条直线运动,每次只能向左或向右随机走一步,那么在走了第N步后,我们可以对这些粒子的位置 x N = ∑ n = 1 N s n x_N=\sum_{n=1}^Ns_n xN=n=1Nsn有哪些认知呢?上式中 s n = k n L s_n=k_nL sn=knL L L L表示每一步的固定步长, k N k_N kN相同的可能性取 ± 1 \pm1 ±1

    首先,从对称性的角度来讲,每一步位移向左或向右是等概率的,均值为零( &lt; s n &gt; = 0 &lt;s_n&gt;=0 <sn>=0),因此对所有粒子随机游走的轨迹求平均不会在任何方向上有系统性的偏移(你会看到一杯咖啡自己沿着桌子跑吗?)。于是, &lt; x N &gt; = ∑ n = 1 N &lt; s n &gt; = 0 &lt;x_N&gt;=\sum_{n=1}^N&lt;s_n&gt;=0 <xN>=n=1N<sn>=0这里的平均严格来说指的是统计力学中的系综平均。

    其次,我们来看位移的方均值(在平均值附近的涨落;或二阶矩,因为均值为零)。 &lt; x N 2 &gt; = &lt; ( x N − 1 + k N L ) 2 &gt; = &lt; x N − 1 2 &gt; + 2 L &lt; x N − 1 k N &gt; + L 2 &lt; k N 2 &gt; &lt;x_N^2&gt;=&lt;(x_{N-1}+k_NL)^2&gt;=&lt;x_{N-1}^2&gt;+2L&lt;x_{N-1}k_N&gt;+L^2&lt;k_N^2&gt; <xN2>=<(xN1+kNL)2>=<xN12>+2L<xN1kN>+L2<kN2>
    因为每一步的 k N k_N kN与前一步的位置 x N − 1 x_{N-1} xN1相互独立,而 &lt; k N &gt; = 0 &lt;k_N&gt;=0 <kN>=0,所以最右边中间一项的贡献为零( &lt; x N − 1 k N &gt; = &lt; x N − 1 &gt; &lt; k N &gt; = 0 &lt;x_{N-1}k_N&gt;=&lt;x_{N-1}&gt;&lt;k_N&gt;=0 <xN1kN>=<xN1><kN>=0)。最后一项正好等于 L 2 L^2 L2,因为 ( ± 1 ) 2 = 1 (\pm1)^2=1 (±1)2=1

    上式表明第N步的 &lt; x N 2 &gt; &lt;x_N^2&gt; <xN2>比上一步的 &lt; x N − 1 2 &gt; &lt;x_{N-1}^2&gt; <xN12> L 2 L^2 L2,以此类推,得到 &lt; x N 2 &gt; = N L 2 &lt;x_N^2&gt;=NL^2 <xN2>=NL2
    把上式改写一下,如果我们等待的时间为 t t t,那么粒子走了 N = t / Δ t N=t/\Delta t N=t/Δt步,其中 Δ t = 1 s \Delta t=1s Δt=1s。于是,一维固定步长的无漂移随机游走的方均位移随着时间线性增加: &lt; x N 2 &gt; = 2 D t &lt;x_N^2&gt;=2Dt <xN2>=2Dt 其中扩散常数 D = L 2 / ( 2 Δ t ) D=L^2/(2\Delta t) D=L2/(2Δt),这个常数也可通过实验来追踪大量咖啡颗粒,对他们的位移方均值随时间的变化曲线取平均得到(基本上是一条过原点的直线)。

    图5. 一维固定步长随机游走的方均位移随步数(等效于时间)变化关系,左边是两个独立的粒子位移随时间变化关系,右边是对500个这样的粒子位移的方均值随时间的变化曲线求平均,基本上可以用一条过原点的直线拟合。图片摘自Nicholas J. Giordano & Hisao Nakanishi[^Nicholas]。
    [^Nicholas]: Giordano, Nicholas J. Computational physics. Pearson Education India, 2006.

    这里的求平均符号需要仔细对待——任意某一个粒子的随机游走并不满足上式的扩散定律,我们在讨论的是足够多个粒子的平均行为。此外要说明的是,一个匀速直线运动的物体位移随时间的变化是线性的 x = v t x=vt x=vt,而一个随机游走的粒子不一样,位移的方均根随时间变化关系在平均行为上表现为 &lt; x 2 &gt; ∼ t 1 / 2 \sqrt{&lt;x^2&gt;}\sim t^{1/2} <x2> t1/2,也就是随时间的增长率要慢一些(指数比1小)。现在我们可以回答前面的咖啡杯问题了,我们假设当 &lt; x 2 &gt; \sqrt{&lt;x^2&gt;} <x2> 增长到与杯子尺寸相当时,咖啡粉在水中完全扩散开,那么当杯子尺寸加倍时, t t t要变为 4 t 4t 4t,即需要等待原来四倍的时间他们才能充分混合。

    从更深层的视角来看,上述扩散定律连接了宏观和微观世界。我们将咖啡粉的布朗运动作为微观世界(主要是水分子热运动造成的效应)和宏观世界(可用肉眼测量的量)之间联系的桥梁。具体来说,我们已经得到了咖啡粉布朗运动的微观参数(步长 L L L和间隔时间 Δ t \Delta t Δt)和宏观可测量量(扩散常数D)之间的关系。

    这个扩散定律还可推广到更高维的情形。在二维棋盘上,每步的位移 r r r x x x, y y y两个独立方向,因此有 &lt; r N 2 &gt; = &lt; x N 2 &gt; + &lt; y N 2 &gt; = 4 D t &lt;r_N^2&gt;=&lt;x_N^2&gt;+&lt;y_N^2&gt;=4Dt <rN2>=<xN2>+<yN2>=4Dt 类似的,三维情形得到 &lt; r N 2 &gt; = 6 D t &lt;r_N^2&gt;=6Dt <rN2>=6Dt

    注:维度越高,时间t前面的系数越大,意味着越容易扩散。这是合理的:更高的维度意味着粒子可以从其他维度的路径走到目的地。

    最后,我们来看一下 x N 2 x_N^2 xN2的涨落。这个量说的是粒子个体之间的涨落差异大小。这么说有点不好理解,具体一点,这个量可以理解为两个从同一地点出发,独立做随机游走的粒子的位移方均值的差异。我们想要知道的是,这个差异随着步数(或时间)是越来越大(发散)还是越来越小(收敛)?在继续阅读前,我建议大家花几秒时间思考一下这个问题。
    ……1s
    ……2s
    ……3s
    好了,我们直接上公式:
    &lt; ( x N 2 − &lt; x N 2 &gt; ) 2 &gt; = &lt; x N 4 &gt; − &lt; x N 2 &gt; 2 \sqrt{&lt;(x_N^2-&lt;x_N^2&gt;)^2&gt;}=\sqrt{&lt;x_N^4&gt;-&lt;x_N^2&gt;^2} <(xN2<xN2>)2> =<xN4><xN2>2
    上式中 &lt; x N 2 &gt; &lt;x_N^2&gt; <xN2>是前面已经算过的,只需要计算四阶矩 &lt; x N 4 &gt; &lt;x_N^4&gt; <xN4>
    &lt; x N 4 &gt; = ∑ i , j , k , l = 1 N &lt; s i s j s k s l &gt; &lt;x_N^4&gt;=\sum_{i,j,k,l=1}^N&lt;s_is_js_ks_l&gt; <xN4>=i,j,k,l=1N<sisjsksl>
    这个求和号一共包含 N 4 N^4 N4项,但是由于每一步移动相互独立,导致了不同下标 i i i s i s_i si相互独立。因此只要有至少一个下标“落单",这一项就为零。非零项只有四个下指标均相同,或下指标能两两配对(但不相同,与前一种情况区别开)这两种情况。所以只剩下
    &lt; x N 4 &gt; = ∑ i = 1 N &lt; s i 4 &gt; + 3 ∑ p = 1 N [ s p 2 ∑ q ≠ p s q 2 ] = 3 N 2 L 4 − 2 N L 4 &lt;x_N^4&gt;=\sum_{i=1}^N&lt;s_i^4&gt;+3\sum_{p=1}^N\left[s_p^2\sum_{q\ne p}s_q^2 \right]=3N^2L^4-2NL^4 <xN4>=i=1N<si4>+3p=1Nsp2q̸=psq2=3N2L42NL4
    代入前面 x N 2 x_N^2 xN2的涨落公式,得
    &lt; ( x N 2 ) 2 − &lt; x N 2 &gt; 2 = 2 N 2 L 4 − 2 N L 4 ∼ 2 N L 2 \sqrt{&lt;(x_N^2)^2-&lt;x_N^2&gt;^2}=\sqrt{2N^2L^4-2NL^4}\sim\sqrt{2}NL^2 <(xN2)2<xN2>2 =2N2L42NL4 2 NL2
    最后一个近似当 N N N足够大时成立。
    这个结果有没有出乎你的意料?它说明 x N 2 x_N^2 xN2的涨落竟然也随着 N N N的增大而变大!也就是两个从同一地点出发独立地做随机游走的粒子越往下走就会离的越远!这和我们之前的对于位置的方均值的讨论是类似的:单个粒子的随机游走,与多个粒子随机游走的平均行为,并不相同

    扩散定律是普适的

    至此我们都在讨论最简单的模型:一维固定步长的无漂移随机游走。当然还有其他更接近真实情况也更复杂的模型,比如每一步的步长在某个连续或离散区间里取值,或者每一步的步长可以无限制地取任意长度( L e v y   f l i g h t Levy\ flight Levy flight)。但是不管模型如何建立,只要给定独立随机移动的步长分布,扩散定律依旧成立,只是系数会有所不同。

    我们举步长取某个离散分布的情况作简要说明,假设 P k P_k Pk定义为步长为 k L kL kL的概率,其中 k k k是一个整数,可正可负,分别代表前进后退。这个模型中,位移的方均值为
    v a r i a n c e ( x N ) = 2 D t ,   其 中 D = L 2 2 Δ t × v a r i a n c e ( k ) variance(x_N)=2Dt,\ 其中D=\frac{L^2}{2\Delta t}\times variance(k) variance(xN)=2Dt, D=2ΔtL2×variance(k)
    依旧满足扩散定律。其中 v a r i e n c e varience varience表示方差。具体推导可以参考Philip Nelson书1的第四章。

    这个结论说明,扩散定律是普适的,定律本身并不依赖于模型,只有系数才依赖于具体的细节。换句话说,对模型的假设做修正只是改变数学的复杂程度,并不改变其物理实质。

    ##由扩散支配的亚细胞世界
    前面说了这么多枯燥的数学,这些结论有什么用?其实,在远比我们小的尺度上,扩散是非常重要的一个过程。没有它,生命可能无法存活。

    图6. 不同尺度下的世界。 图片摘自https://microbiologyinfo.com/different-size-shape-and-arrangement-of-bacterial-cells/。
    我们从上图的中间可以看到,细菌的尺寸在$1\mu m$左右(约头发丝的1/50倍)。假设一个细菌可以看成半径为$1\mu m$的小水球,并假设在水球中心放入一些糖分子,大概需要多长时间,这些糖分子才能均匀扩散到整个细菌?换做是细胞大小的水球($半径10\mu m$)呢?

    对于室温下水中的小分子,扩散系数的大小为 D ≈ 1 0 − 9 m 2 / s − 1 D\approx 10^{-9}m^2/s^{-1} D109m2/s1,或者写为 D ≈ 1 μ m 2 m s − 1 D\approx 1\mu m^2ms^{-1} D1μm2ms1。于是,在细菌中的糖分子需要时间 ( 1 μ m ) 2 / ( 6 D ) ≈ 0.2 m s (1\mu m)^2/(6D)\approx 0.2ms (1μm)2/(6D)0.2ms就能扩散开来。而在比它尺寸大一个数量级的细胞中,要花100倍于此的时间。这说明,如果想要在微米级别传递物质或信息,仅仅扩散就够了,不需要花力气设计复杂的结构,这或许也是细菌结构下相对简单的原因之一吧。

    另一方面,如果想要在长距离传输物质,比如从你看到远处的狮子后大脑发出信号,再通过脊髓传递到腿部的肌肉神经(约1m长的距离)撒腿逃跑,这个信号传递过程如果仅凭扩散来完成的话你早就被吃掉了,于是神经元细胞自己发展出了一套高速传递信息的方式:电信号。当然,在不同神经元之间的突触间隙(约20-30nm)中,神经递质也是通过扩散来传递的。

    多说一句,这个扩散的粒子流(前一个神经元流到后一个神经元)并不是由于浓度大的地方粒子之间相互推挤而向把他们浓度低的地方驱赶,因为我们假设所有递质粒子相互独立,也就是没有任何相互作用(我们认为神经递质的数量远小于周围水分子的数量)。导致这个粒子流的合理解释是:由于每个粒子向各个方向运动的概率是均等的,那么粒子多的地方会有更多粒子跑到粒子少的地方。**这是一个由“概率"驱动的流动。**这个概念与“熵"以及“熵力"紧密相连,有机会再细说。

    图7. 神经突触示意图。(右下角)神经递质被释放后,向另一边扩散。 图片摘自https://upload.wikimedia.org/wikipedia/commons/3/30/Chemical_synapse_schema_cropped.jpg。

    摩擦与扩散是一对双胞胎——爱因斯坦关系

    最开始我们对咖啡分子在水中溶解的过程有过唯象的描述,现在我们再来审视一下这个过程。依旧做一些比较强的假设来简化计算。考察一维运动,假设每次碰撞间隔都是 Δ t \Delta t Δt,在两次碰撞之间颗粒的运动遵循牛顿定律 Δ x = v 0 Δ t + 1 2 f m ( Δ t ) 2 \Delta x=v_0\Delta t+\frac{1}{2}\frac{f}{m}(\Delta t)^2 Δx=v0Δt+21mf(Δt)2其中施加在粒子上的外力 f f f可以是电场对该带电粒子的力,也可以是重力等。

    另假设每次碰撞消去了碰撞之前的所有记忆,因此 v 0 v_0 v0随机地指向左边或右边,均值为零。对上式求平均得到 &lt; Δ x &gt; = f 2 m ( Δ t ) 2 &lt;\Delta x&gt;=\frac{f}{2m}(\Delta t)^2 <Δx>=2mf(Δt)2这个式子表明,在施加了力场后,粒子沿着力的方向有个恒定速度的漂移,平均漂移速度为 v d r i f t = &lt; Δ x &gt; Δ t = f / ζ v_{drift}=\frac{&lt;\Delta x&gt;}{\Delta t}=f/\zeta vdrift=Δt<Δx>=f/ζ上式中粘性摩擦系数 ζ = 2 m / Δ t \zeta=2m/\Delta t ζ=2m/Δt可以看到这个漂移速度正比于外力的大小。

    ζ \zeta ζ可以通过实验测量。球体的粘性摩擦系数和半径服从简单的关系,即斯托克斯公式 ζ = 6 π η R \zeta=6\pi\eta R ζ=6πηR。室温下水的粘度约为 η = 1 0 − 3 k g   m − 1 s − 1 \eta=10^{-3}kg\ m^{-1}s^{-1} η=103kg m1s1。因此只要测得了颗粒的尺寸,就能得到 ζ \zeta ζ,进一步测量密度可以间接测量颗粒的质量 m m m。有了实验能够测量的 ζ \zeta ζ m m m,那么碰撞的时间间隔 Δ t \Delta t Δt就可以由 ζ = 2 m / Δ t \zeta=2m/\Delta t ζ=2m/Δt计算出来。进一步,从扩散定律 D = L 2 / ( 2 Δ t ) D=L^2/(2\Delta t) D=L2/(2Δt)就可以得到分子尺度的特征步长 L L L

    目前为止,上述理论有一个致命的缺陷:我们不可能在实验中观测到 Δ t \Delta t Δt L L L,因此这个理论无法证伪

    所幸的是,爱因斯坦注意到在 Δ t \Delta t Δt L L L之间还有另一个关系。记特征速度为 v = L / Δ t v=L/\Delta t v=L/Δt,统计力学告诉我们,(一维)理想气体分子热运动满足 m &lt; v 2 &gt; = k B T m&lt;v^2&gt;=k_BT m<v2>=kBT

    现在我们未知的只是 Δ t \Delta t Δt L L L两个变量,但是有粘性摩擦定律,扩散定律,还有热运动速度的关系式三个方程构成了超定方程,唯一的可能是参数 ζ \zeta ζ D D D并非独立,而是有个特定关系,这个关系就是爱因斯坦关系 ζ D = k B T   ( 爱 因 斯 坦 关 系 ) \zeta D=k_BT\ (爱因斯坦关系) ζD=kBT 

    这个关系由爱因斯坦在1905年得到,它表明颗粒扩散系数(描述位置的涨落)和粘性摩擦系数(描述颗粒受到的耗散)相互联系。

    爱因斯坦关系具有深刻的含义。通过它我们可以从宏观测量推测出分子的大小,而不需要用高端的显微镜来完成。我们知道1mol理想气体的体积为22.4L(0℃,1标准大气压),通过理想气体状态方程 P V = n R T PV=nRT PV=nRT可以算出气体常数 R = k B N a R=k_BN_a R=kBNa。而实验参数能够给出 k B k_B kB,于是我们知道了阿伏伽德罗常数 N a N_a Na,从而可以估算分子的尺寸。

    此外,爱因斯坦关系等号右边没有质量 m m m,这说明该关系不依赖于颗粒的质量。颗粒质量越大,扩散起来越难(D越小),且摩擦阻力越大( ζ \zeta ζ越大)。这种普适性也提供了可证伪的预言:我们可以检查不同颗粒质量、温度、尺寸条件下,是否都能给出相同的 k B k_B kB值。另外,尽管 ζ \zeta ζ D D D与温度的关系很复杂,但他们的乘积与温度的关系却非常简单。值得称道的是,这个关系已经被实验完美证实。

    爱因斯坦关系溯源——涨落耗散定理

    一个优美的物理关系背后一定有客观存在的真理。爱因斯坦关系其实是一个更普适的定理在随机游走问题中的体现——涨落耗散定理(Fluctuation-dissipation theorem)。涨落耗散定理是统计物理中非常重要而且有意思的话题,一方面在于深入地理解近平衡体系中响应与关联的关系;另一方面在于非平衡态体系中的涨落耗散关系。关于这些内容在其他地方讨论。


    1. Nelson, Philip. Biological physics. New York: WH Freeman, 2004. ↩︎

    展开全文
  • 巧用随机区组设计让你实验更轻松

    千次阅读 2020-02-26 13:50:13
    值得注意的是,试验者得到多个批次的试验数据后,可能未意识到这是随机区组设计,统计分析时没有加上批次这一因素,导致误差变大:①当批次间差异较小时,单因素与单因素随机区组的统计结果是相差不多的,下次进行...

    1、单因素随机区组方差分析  VS  t检验

    当你想要比较A、B两种处理的效果时,第一次实验数据如下: 

    A1 <- c(2.68, 2.47, 2.62 )
    B1 <- c(2.69, 2.78, 2.66 )
    t.test(A1, B1)

    p-value = 0.189 > 0.05

    第二次实验数据如下:

    A2 <- c(3.07, 3.14, 3.23)
    B2 <- c(3.26, 3.26, 3.29)
    t.test(A2, B2)

    p-value = 0.1108 > 0.05

    第三次实验数据如下:

    A3 <- c(2.84, 2.85,2.92)
    B3 <- c(2.92, 3.02,2.99)
    t.test(A3, B3)

    p-value = 0.05318 > 0.05

    实验做到这一步,你可能要疯了,毕竟做了三次,p值都大于0.05。当然,你可能会想到把三次实验的数据放一起分析。

    A <- c(2.68, 2.47, 2.62, 3.07, 3.14, 3.23, 2.84, 2.85, 2.92)
    B <- c(2.69, 2.78, 2.66, 3.26, 3.26, 3.29, 2.92, 3.02, 2.99)
    t.test(A, B)

    p-value = 0.3344 > 0.05,终于,你忍不住要发pyq了。

    事实上,如果你认真观察会发现,每次实验得到的数据存在些许差异,你看看下面均值的变化:

    A1:  2.59    B1: 2.71

    A2:  3.15    B2:   3.27

    A3:  2.87    B3:   2.98

    不同“批次”实验间的差异,可能是由于每次实验的材料、测量仪器的工作状态或者实验员实验熟练程度等非实验因素条件造成的,你可以把每一次实验都看做一个“区组”,区组内非试验因素差异最小而区组间内非试验因素差异最大,当你直接用t检验的时候,会将区组间非试验因素的差异并入误差项,使得原本显著差异的结果变得不显著。

    如果你用单因素随机区组的方法分析,便可以将批次间的差异剔除出。

    res <- c(2.68, 2.47, 2.62, 3.07, 3.14, 3.23, 2.84, 2.85, 2.92, 
             2.69, 2.78, 2.66,3.26, 3.26, 3.29, 2.92, 3.02, 2.99)
    trt <- factor(c(rep('A', 9), rep('B', 9)))  # 处理因素
    bk <- factor(c(1,1,1,2,2,2,3,3,3,1,1,1,2,2,2,3,3,3))  #  区组因素
    data <- data.frame(trt, bk, res)
    aov <- aov(res ~ trt + bk, data = data)
    summary(aov)

                        Df         Sum Sq         Mean Sq      F value         Pr(>F)    
    trt                 1            0.0613           0.0613         15.89         0.00135     ** 
    bk                2            0.9353           0.4677         121.32       1.44e-09    ***
    Residuals   14            0.0540           0.0039                     
    ---
    Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

    结果表明,处理因素间存在极显著差异(p = 0.00135 < 0.01),同时,结果也表明区组间也存在极显著差异(p = 1.44e-09 < 0.01),说明采用这种方法进行分析是一种d明智且正确的选择。

    2、随机区组设计助你合理规划实验

    随机区组设计(randomized blocks design):指根据局部控制和随机原则进行的,将试验单位按性质不同分成与重复数一样多的区组(窝组),使区组内非试验因素差异最小而区组间内非试验因素差异最大,每个区组均包括全部的处理。区组内各处理随机排列,各区组独立随机排列。

    随机区组设计通过划分区组降低对试验单位的整体要求,只要求各区组内非试验因素条件尽量一致,能有效减少非试验因素的单向差异,降低试验误差,广泛应用于生物实验设计。但也存在不足之处,即不允许处理数太多。当处理数太多,或需要在区组内设置重复时,必然导致区组规模增大,组内误差相应增大,就达不到局部控制的效果。——常规应用场景

    把随机区组的观念灵活应用在实验中,可以是这样的:试验单位满足试验的整体要求,即各试验材料差异很小,且材料充足。但试验的重复数较多,在一次试验中难以完成。试验者可分批次进行试验,每一批次即为一个区组,根据自身情况合理分配每一批次试验的重复数,从而将繁重的工作平均分配到各批次试验。

    在之前的一篇文章《原创 R语言告诉你该做几个实验重复》中,给出了一种估计实验重复数的方法,通过功效分析,你想要有90%的把握得到p<0.05的数据需要做9个重复,如果实验不复杂,不耗时耗力,也不缺材料,你当然可以在一次实验中做9个重复,如果实验一次做完9个重复,或者你不想一次做那么多,想要做半天玩半天,完全可以分3次,每次3个重复来开展实验。

    以下提供一个轻松做实验的“套路”:
    1)当你做完一次实验后,p>0.05,先别急着继续干,可根据《原创 R语言告诉你该做几个实验重复》中的方法做一个功效分析,看看自己该做几个重复比较合适;

    2)当你知道该做几个重复数之后,可以这样安排实验:

    N=4~6:一次做4~6个重复

    N=7~123~4批次,每次3~4个重复

    N=13~203~4批次,每次4~5个重复

    如果重复数很大,你可能需要衡量下这个实验是否有继续的必要性,或许你的处理本来就是没有效果的,就算你做再多也做不出来。每个人可以根据实验的复杂性,个人精力,时间的紧迫性合理安排自己的实验,以上仅是一些建议。

    3、小结

    本文的例子是单因素两水平的,你还可以将其扩展到单因素多水平,双因素甚至三因素的应用上。

    值得注意的是,试验者得到多个批次的试验数据后,可能未意识到这是随机区组设计,统计分析时没有加上批次这一因素,导致误差变大:①当批次间差异较小时,单因素与单因素随机区组的统计结果是相差不多的,下次进行类似试验的设计时可不考虑批次的影响;②当批次间差异较大时,可以通过统计分析剔除批次间的差异,说明在重复数多的情况下,采用随机区组设计,进行多批次试验做法是可行的;③多批次试验数据的分析,倘若忽略“批次”这个因素,可能使原本显著差异的结果变得不显著。

    展开全文
  • 随机产生20个单词

    千次阅读 2019-10-14 14:30:49
    随机产生20个单词 一、问题来源: 老师给了一份专业单词word,说第二天要全背下来。错了就五十遍啊五十遍。 然后,有人提出要做一个产生随机单词的Demo,来测试自己。 老师表示呵呵,做出...
  • 数值随机化算法和舍伍德随机算法

    千次阅读 2020-06-01 19:25:04
    随机化算法的一个基本特征就是对所求解的问题的同一实例用同一随机化算法求解两次可能得到完全不同的结果。 随机化算法可以分为四类:数值随机化、蒙特卡洛算法、拉斯维加斯算法和舍伍德算法。这次我们主要介绍数值...
  • 一文通透优化算法:从随机梯度、随机梯度下降法到牛顿法、共轭梯度       1 什么是梯度下降法 经常在机器学习中的优化问题中看到一个算法,即梯度下降法,那到底什么是梯度下降法呢? 维基百科给出的定义...
  • 随机数C语言编程

    千次阅读 2018-06-10 22:38:03
    在程序设计中,有时会用到随机数。本文介绍在 Linux 编程环境下,如何生成伪随机数。 什么是伪随机数 伪随机数是通过一个确定性的算法计算出来的“似乎”是随机的数序,因此伪随机数实际上并不随机。在计算伪...
  • 随机效应与固定效应&面板数据回归

    万次阅读 多人点赞 2014-11-13 21:07:45
    随机效应与固定效应 方差分析主要有三种模型:即固定效应模型(fixed effects model),随机效应模型(random effects model),混合效应模型(mixed effects model)。 所谓的固定、随机、混合,主要是针对...
  • 随机森林

    千次阅读 2015-10-28 17:12:14
    Random Forest(s),随机森林,又叫Random Trees[2][3],是一种由多棵决策树组合而成的联合预测模型,天然可以作为快速且有效的多类分类模型。如下图所示,RF中的每一棵决策树由众多split和node组成:split通过输入的...
  • 加权随机抽样算法

    千次阅读 2019-02-18 16:17:24
    1. 基于均匀分布概率的算法   例如,3等奖抽中的概率是70%,2等奖是20%,1等奖是10%,这样,大部分人都只能中3等奖,小...另一个例子:按权重均一化后,编号3被抽中的概率要求是70%,5出现的概率为25%,0出现的...
  • 随机模型,估计与控制 ——介绍

    千次阅读 2019-04-19 10:27:51
    在考虑系统分析或控制器设计时,工程师可以运用确定性系统和控制理论的丰富知识。人们自然会问,为什么我们必须超越这些结果,提出随机系统模型,并基于这些随机模型提出估计和控制概念?为了回答这个问题,让我们...
  • 全网目前最全python例子(附源码)

    万次阅读 多人点赞 2019-12-30 14:55:05
    告别枯燥,60秒学会一个小例子,系统学习Python,从入门到大师。Python之路已有190个例子: 第零章:感受Python之美 第一章:Python基础 第二章:Python之坑 第三章:Python字符串和正则 第四章:Python文件 第五章...
  • Java生成随机数、随机种子

    万次阅读 2018-06-01 10:12:27
    Java里面有一个随机函数——Random,刚开始只是知道这个函数具有随机取值的作用,于是上网搜索了资料一番,做了一下一些关于Random函数的总结: Java中存在着两种Random函数:一、java.lang.Math.Random;...
  • 5.8 决策树和随机森林 原文:In-Depth: Decision Trees and Random Forests 译者:飞龙 协议:CC BY-NC-SA 4.0 译文没有得到原作者授权,不保证与原文的意思严格一致。 之前,我们深入研究了简单的生成...
  • 研究一类带乘性噪声的离散时间随机Markov跳跃系统的有限时间控制问题. 首先,定义系统的有限时间.稳定和有限时间有界,通过逐次迭代和条件期望给出系统有限时间稳定的充分必要条件;其次,针对含干扰的系.统,利用...
  • 三天研读《中兴电路设计规范》精华总结

    万次阅读 多人点赞 2020-05-16 18:25:52
    本博客将简述中兴通讯股份有限公司在原理图设计中需要注意的一些事项,其中包含了中兴设计开发部积累的大量硬件开发知识和经验,可以作为学习使用。硬件工程师可以学习并掌握检查条目的内容以及对条目的详细说明,...
  • 入门python,看完这个300行代码的例子,足矣~

    万次阅读 多人点赞 2020-04-28 21:35:26
    一个例子全搞定! 一个300行的代码,竟然包含了138个知识点。列表,元组,字典,集合,字符串,也有他们的基本操作,有面向对象的类,循环语句,选择语句,函数的创建,包的导入,文件的读取,切片,表达式推导。 还...
  • 条件随机场(Conditional random fields),是一种判别式图模型,因为其强大的表达能力和出色的性能,得到了广泛的应用。从最通用角度来看,CRF本质上是给定了观察值集合 (observations)的马尔可夫随机场(MRF)。在...
  • 系统设计:关于高可用系统的一些技术方案

    万次阅读 多人点赞 2017-09-17 09:22:32
    系统设计关于高可用系统的一些技术方案 高可用方法论 扩展 隔离 解耦 限流 分类 漏桶算法 令牌桶算法 滑动窗口计数法 动态限流 降级 熔断 发布相关 模块级自动化测试 灰度发布 回滚 其他 总结 参考资料 ...
  • 一 arcgis10.6种有这两种方法深度学习与图像分类图像分类是基于不同地物光谱、形状等特征差异将影像...ArcGIS Pro集成的机器学习图像分类方法包括:最大似然分类,随机树,支持向量机等。在Pro 2.1中,ArcGIS提供2个...
  • 源码的下载地址时http://yunpan.cn/cjwwij3FcBtZV 访问密码3579 本列表源码永久免费下载地址...卷 yunpan 的文件夹 PATH 列表 卷序列号为 0000-73EC E:. ...│ 例子大全说明.txt │ 本例子永久更新地址~.url │
  • 标准模型与随机预言模型的比较

    千次阅读 2016-11-08 13:25:59
    随机预言机模型:在安全证明中,随机预言机模型通常是现实中哈希函数的理想化的替身。哈希函数是一个输入为任意长度,输出为固定长度的函数,...在随机预言机模型下,通常设计一个方案并证明是安全的;而在方案的实际执
  • 随机数加密思路

    千次阅读 2006-06-14 15:23:00
    随机数加密思路 伪随机数加密是伪随机数组作为密钥加密的例子。 它工作时需要用户输入密码。密码是干什么的?密码就像一个种子,有了它才能完成加密或解密。密码的数值渗透到整个过程,它决定,密钥的发生和使用...
  • 随机访问的队列

    千次阅读 2017-02-09 17:10:10
    最近在做毕设,关于手势识别的,自己yy了一种算法,实现里面为了方便,需要一种可供我随机访问的大小固定的队列
  • 随机场(Random field)

    万次阅读 2012-07-25 20:38:21
    一、随机场定义 http://zh.wikipedia.org/zh-cn/随机场  随机场(Random field)定义如下: 在概率论中, 由样本空间Ω = {0, 1, …, G − 1}n取样构成的随机变量Xi所组成的S = {X1, …, Xn}。若对所有的ω...
  • 这个挑战的目标是设计和实现一个系统,通过概括在输入文本中发现的模式随机生成合理的输出。 该系统应该能够根据我使用输入文本设计的模型生成独特的句子。 这样的系统可用于随机生成诗歌、歌词或神秘的学术论文。 ...
  • 常用概率分布函数及随机特征

    万次阅读 2018-02-08 20:28:52
    常见分布的随机特征离散随机变量分布伯努利分布(二点分布)伯努利分布亦称“零一分布”、“两点分布”。称随机变量X有伯努利分布, 参数为p(0&lt;p&lt;1),如果它分别以概率p和1-p取1和0为值。EX= p,DX=p(1-p...
  • 在很多时候,我们会使用rand()函数与srand()配合来达到产生随机数的效果,srand初始化随机种子,rand产生随机数,下面进行展开的分析(当然我们在这里先不考虑某些游戏引擎会另外设计自己的随机数产生机制。...
  • Spark设计理念与基本架构

    万次阅读 热门讨论 2016-01-22 13:52:31
    MLlib目前已经提供了基础统计、分类、回归、决策树、随机森林、朴素贝叶斯、保序回归、协同过滤、聚类、维数缩减、特征提取与转型、频繁模式挖掘、预言模型标记语言、管道等多种数理统计、概率论、数据挖掘方面的...

空空如也

空空如也

1 2 3 4 5 ... 20
收藏数 98,441
精华内容 39,376
关键字:

完全随机设计的例子