16 混合

在之前的章节，我们学习了如何模拟咖啡冷却的过程 ，只要我们知道咖啡初始温度，空气温度，和冷却速率（r）我们就能预测咖啡在规定时间冷却了多少，但是，一般而言，我们是不知道冷却速率（r）的，我们可以通过数据进行估测它，只要给出初始温度和最终温度，以及时间，我们就可以通过测试和纠错来算出r。在这一章，我们就通过 二分法来测出 r 。

16.1 计算根

ModSim 里提供了一个方法叫 root_bisect 来求非线性方程的根，比如你想算：

那么就可以算出x=1，x=2，x=3 这些事方程的根（默认计算f(x) = 0 的值）

下面我就用一个例子来解释一下如何用这个函数：

#首先自己定义一函数，待会算这个函数的根
def func(x):
    return (x-1) * (x-2) * (x-3)

#下面就是使用求根函数，两个参数，一个参数是我们的要求根的方程，
#第二是限制根的范围，指定我们要算哪个范围的根
#此函数是用二分法来算，用1.5 和2.5 的中间值 x2 作为func的参数
#得到一个数，如果大于0就把x2-1.5作为新的范围进行迭代（反之则x2-2.5）直到算出值
res = root_bisect(func, [1.5, 2.5])
print(res.root)

下面我们如何用 root_bisect() 函数类估测我们的 r 呢？

首先我们得定义一个函数，这个函数的输入是我们要算的r，我们再测出一个粗略的r所在的范围，然后让root_bisect() 来迭代计算

#这个函数通过给定的r 计算预测温度，和实际温度做比较（做差）作为输出 
#然后让root_bisect()计算r
def error_func1(r):
    system = make_system(T_init=90, r=r, volume=300, t_end=30)
    results = run_simulation(system, update_func)
    T_final = get_last_value(results.T)
    return T_final - 70

定义好我们的方程（函数） 就是求方程的跟啦（r） 下面就是用root_bisect() 来求根

#首先我们得粗略算出r所在区间，作者 把0.01带入 得到值是2.3℃（正）所以需要把 r调大
error_func1(r=0.01)

#把0.02带入得 -10.9℃ （负）  所以区间就定好了，r真实值在0.01到0.02之间 
#剩下的就可以交给root_bisect()啦
error_func1(r=0.02)

#通过函数求得我们的r_coffe 的值为0.012
res = root_bisect(error_func1, [0.01, 0.02])
r_coffee = res.root

通过同样的方法计算 r_milk 的值

16_1.png

图16.1： 咖啡与牛奶的随时间变化的温度

16.2 混合液体（牛奶和咖啡）

下面我们混合两个液体（牛奶和咖啡），混合液体的温度取决于混合的材料，并且我们假设没有化学反应或者其他因素来影响温度，我们的热量是守恒的。

如果第一个液体的温度是 T₁ 第二个液体的温度是 T₂ ，混合后的温度是 T ，那么与第一个液体进行的热量交换是 C₁(T - T₁) ,同理，与第二个液体进行的热量交换是 C₂(T - T₂) , C₁ 和 C₂ 是两个液体的比热容(热质量)，根据热量守恒，两个交换的热量相加应该为0，如下：

下面我们把 T 提取出来：

如果我们知道液体的密度和体积以及比热容cp，那么我们就可以算出来热质量C：

如果液体有相同的比热容和密度，那么我们就可以得到以下式子(v1 v2 是两个液体的体积)：

下面就是对混合系统的模拟了:

def mix(system1, system2):
    #判断两个系统的结束时间，确认是否相同，如果不相同就返回错误并且终止运行其他代码
    assert system1.t_end == system2.t_end
    #分别赋值两个系统的体积
    V1, V2 = system1.volume, system2.volume
    #分别赋值两个系统的温度
    T1, T2 = system1.temp, system2.temp
    #混合后的体积和温度
    V_mix = V1 + V2
    T_mix = (V1 * T1 + V2 * T2) / V_mix
    #返回一个新的系统，温度为混合后的温度，体积为混合后的体积
    return make_system( T_init=T_mix,
                        r=system1.r,
                        volume=V_mix,
                        t_end=30)

主义我们的mix（）函数要求我们的两个系统的温度有 temp变量，并且是已经更新过的温度，我用了这个函数：

def run_and_set(system):
    #更新温度
    results = run_simulation(system, update_func)
    system.temp = get_last_value(results.T)
    #返回的是对象
    return results

现在我们有了所有我们需要的东西！

16.3 先混合还是后混合？

首先我先创建两个对象来代表咖啡和奶油

#注意两者的冷却时间都是一样的，是30（分钟）
#咖啡
coffee = make_system(T_init=90, r=r_coffee,volume=300, t_end=30)
#奶油
milk = make_system(T_init=5, r=r_milk, volume=50, t_end=30)

先混合：

mix_first = mix(coffee, milk)
mix_results = run_and_set(mix_first)
#先混合 得到混合后冷却得到的温度是61.4℃

后混合：

coffee_results = run_and_set(coffee)
milk_results = run_and_set(milk)
mix_last = mix(coffee, milk)
#先各自放置30分钟再混合得到的温度是 63.1℃

显然 后混合效果好一些，咖啡温度也高1.7℃！

16.4 优化

过30分钟再添加牛奶比立即添加牛奶要好一些，但是可能有更好的选择，下面我用新的函数来调节添加牛奶的时间，让T_add 作为参数

def run_and_mix(t_add, t_total):
    #咖啡  
    coffee = make_system(T_init=90, r=r_coffee,volume=300, t_end=t_add)
    coffee_results = run_and_set(coffee)
    #牛奶
    milk = make_system(T_init=5, r=r_milk,volume=50, t_end=t_add)
    milk_results = run_and_set(milk)
    #混合
    mixture = mix(coffee, milk)
    #可以看到 混合液体冷却的温度是： 总温度-我们尝试加牛奶的温度
    mixture.t_end = t_total - t_add
    results = run_and_set(mixture)
    return mixture.temp

t_add =0 时 表示立即添加牛奶（也就是混合液最后得到温度：61.4℃），t_end = 30表示最后添加牛奶（也就是混合液体最终温度为：63.1℃），下面我们尝试不同的添加牛奶时间，来看看什么时间添加会好一些：

sweep = SweepSeries()
    #表示0到30之间选择了11次（numpy.linespace()函数中表示11个）
    for t_add in linspace(0, 30, 11):
        sweep[t_add] = run_and_mix(t_add, 30)

得到的结果作图得到如下图：

16_2.png

图16.2：以牛奶加入的时间为变量最终的温度曲线

图片横纵表是我们添加牛奶的时机，可见在最后添加牛奶是最好的！

16.5 分析

用牛顿冷却定律来求解是不太合理的，我们可以用微分方程来求解：

一般求解如下：

特解：T₍₀₎ = T_init ：

现在我们可以用观测的数据来计算估测我们的参数 r， 如果T(tend) = Tend 我们可以把tend 和Tend带入到特解了 来求出我们的 r ， 公式如下：

我们让 tend = 30 ， Tend = 70 ，Tinit = 90 ，Tenv = 22 我们就可以计算出估测 r 为 0.0116

接下来我们可以用下面的方程来计算出 我们的时间序列：

def run_analysis(system):
    T_env, r = system.T_env, system.r
    T_init = system.init.T
    #numpy.linspace(num = 50)缺省为50
    ts = linspace(0, system.t_end)
    T_array = T_env + (T_init - T_env) * exp(-r * ts)
    return TimeFrame(T_array, index=ts, columns=['T'])

这个函数类似于run_simulation 它把System当做参数，然后返回TimeFrame作为结果

因为linerange返回的是Numpy 矩阵， T_array也是Numpy矩阵，为了让run_simulation()里的数据格式一致，我们把它放到 TimeFrame , 下面是代码：

#之前计算得到r
r_coffee2 = 0.0116
coffee2 = make_system(T_init=90, r=r_coffee2,volume=300, t_end=30)
results = run_analysis(coffee2)

最终温度是70℃，我建议您先阅读chap16.ipynb 然后自己下载数据，进行实践一下，这有助于您的学习和更深刻的了解！

本书的中文翻译由南开大学医学院智能医学工程专业2018级、2019级的师生完成，方便后续学生学习《Python仿真建模》课程。翻译人员（排名不分前后）：薛淏源、金钰、张雯、张莹睿、赵子雨、李翀、慕振墺、许靖云、李文硕、尹瀛寰、沈纪辰、迪力木拉、樊旭波、商嘉文、赵旭、连煦、杨永新、樊一诺、刘志鑫、彭子豪、马碧婷、吴晓玲、常智星、陈俊帆、高胜寒、韩志恒、刘天翔、张艺潇、刘畅。

整理校订由刘畅完成，如果您发现有翻译不当或者错误，请邮件联系changliu@nankai.edu.cn。

本书中文版不用于商业用途，供大家自由使用。

未经允许，请勿转载。