欢迎关注我的知乎账号和过程系统工程专栏，专注于化工、过控、制药、优化的学习心得体会分享。如需转载，请私信。收藏前，欢迎点赞！

声明：代码仅供参考，作者不对代码产生的后果负法律责任。

上文写到，用Matlab与Aspen Plus（AP）通信，上文链接。本文继续前文，用同样的例子三效蒸发器说明如何用Python与AP通信。Matlab虽好，但毕竟是商业软件有运用成本，而且现在很多学校也被美国禁用了，Python作为一种广泛使用的胶水语言、开源语言，有很多独到的优势。

用Python与AP通信的详细步骤和思路，完全和Matlab一模一样，只是在具体的代码上有不同。所以详细细节请参看上文，这里只把代码放上来，代码有详细注释。

# -*- coding: utf-8 -*-
"""
Created on Tue Nov 10 11:36:45 2020

@author: Devin @Zhihu https://www.zhihu.com/people/devins/posts
"""

# 导入需要的库
import os
import win32com.client as win32
import time
import numpy as np

# 创建AP的本地服务器

# AP8.8版替换下面的36.0为34.0; 9.0替换为35.0; 10.0替换为36.0； 11.0替换为37.0
Application = win32.Dispatch('Apwn.Document.36.0')

# 获取当前文件夹的地址，三效蒸发器文件和本程序文件需放置在同一个文件夹
address = os.getcwd()

# AP的bkp文件的文件名，即三效蒸发器bkp文件的文件名
SimulationName = 'Triple-Effect Evaporator'

# 打开三效蒸发器文件
Application.InitFromArchive2(os.path.abspath(SimulationName + '.bkp'))

# 设置AP用户界面的可见性，1为可见，0为不可见
Application.Visible = 1

# 压制对话框的弹出，1为压制；0为不压制
Application.SuppressDialogs = 1

# 试运行三效蒸发器模拟
Application.Engine.Run2(1)

# 因为AP运行比较慢，所以用一个循环语句，每两秒钟检查一次是否运行完毕
while Application.Engine.IsRunning == 1:
    time.sleep(2)

# 定义输入、输出变量、收敛标签
FeedFlow = [1200,1400]
HeatDuty = []
SimulationConvergency = []

 
# 程序核心部分

# 用for循环依次改变流量为1200、1400，并运行模拟，然后存储数据，检查历史文件
for ii in range(len(FeedFlow)):
    # 设置输入流量
    Application.Tree.FindNode("\Data\Streams\BRINE\Input\TOTFLOW\MIXED").Value = FeedFlow[ii]
    
    # 初始化，用于清除之前的数据
    Application.Reinit
    
    # 运行程序
    Application.Engine.Run2(1)

    # 定义本次的输出数据
    heatduty = []

    # 每两秒钟检查一次是否运行完毕
    while Application.Engine.IsRunning == 1:
        time.sleep(2)

    # 获取三个蒸发器的热负荷
    heatduty.append(Application.Tree.FindNode("\Data\Blocks\STAGE-1\Output\QCALC").Value)
    heatduty.append(Application.Tree.FindNode("\Data\Blocks\STAGE-2\Output\QCALC").Value)
    heatduty.append(Application.Tree.FindNode("\Data\Blocks\STAGE-3\Output\QCALC").Value)

    # 存储到HeatDuty中
    HeatDuty.append(heatduty)

    # 获取运行ID，并以文本格式打开相应的历史文件，搜索文件中是否有关键词
    name = Application.Tree.FindNode("\Data\Results Summary\Run-Status\Output\RUNID").Value
    Filename = address + '\\' + name + '.his'

    # 如果关键词出现，则标记这一次结果为没有收敛
    with open(Filename,'r') as f: 
        isError = np.any(np.array([line.find('SEVERE ERROR') for line in f.readlines()])>=0)
        SimulationConvergency.append(not isError)

# 关闭AP        
Application.Close
Application.Quit

参考文献

Aspen Plus User Guide

培养基里有各种各样的物质，细胞需要摄取这些物质，然后通过氧化过程，为生长和生存提供能量和各种物质需求。这些转换过程都可以表示成化学反应，如果绘制起来就是代谢网络图，如下图所示。

代谢网络图（图片来自网络）

这些网络中的化学反应，都可以表示成我们高中学的化学反应的形式：

aA+bB\rightarrow cC+dD

这里的a, b, c, d就是化学计量数。

同样的细胞生长本身也可以看作一个化学反应。细胞干重中，主要是蛋白质，蛋白质是由各种氨基酸组成的，合成细胞的过程中，需要消耗一定量的能量（ATP）：

aAla+bMet + cPhe... +eATP \rightarrow c细胞

这些系数，是有办法通过实验测到的。而且非常有意思的是，这些系数对于很多实验室常见的细菌而言，有差别，但差别不是十分的大，所以很多文献经常用大肠杆菌的来做，一些参数通过拟合，就能得到不错的效果。

例如下图，有这样一个细胞，里面还有一个细胞器。外部的A （即途中EX_A[e]）可以流入细胞变成A[e]，再留入细胞器变成A[c]。这些不同空间的物质交换，也可以表示成化学反应的形式。

图片来自书Optimization methods in metabolic networks

涉及的反应如下：

图片来自书Optimization methods in metabolic networks

其中Sink（阱），可以被理解成物质的最终出路，比如合成细胞（biomass，生物质）本身。

所有的反应都可以打包起来，写成一个矩阵，即化学计量数矩阵 \bm{S} 。需要注意的是，化工里，把反应物的系数记为负值，产物的系数写成正值，因为负表示体系中在消耗。

图片来自书Optimization methods in metabolic networks

如果反应速率写成向量 \bm{v} ，化学计量数矩阵是 \bm{S} 。那么 \bm{Sv}= \bm{0} 就表示，代谢流守恒，也就是流入的物质，流出细胞的物质，还有相互之间的转化是遵循物质守恒的。这个式子成立需要假设细胞短时间内，处于拟稳态（quasi-steady state）。如果不假设拟稳态，某些代谢物会积累。可以这样假设的另一个原因是，通常这些代谢物反应，相对宏观的培养基物质消耗快速。在上文提到的物质守恒和拟稳态的基础上，细胞的目的利用现有的资源，去生长，去合成新的细胞。所以这就变成一个经典的线性规划问题。目标是生长，约束是物质守恒和其他生长条件的约束。

拟稳态代谢流模型（Flux balance analysis，FBA），通常情况下表示为线性规划问题，但不限于线性规划。形式如下：

\max_\bm{v} ~~\bm{c}^T\bm{v}\\ \bm{Sv}= \bm{0}\\ \bm{l}\le \bm{Sv}\le \bm{u} (1)

其中 \bm{v} 是代谢流（也就是反应速率），\bm{S}是化学计量数矩阵。 \bm{l} 和 \bm{u} ，表示代谢流的速率限制条件，比如，氧气、葡萄糖的摄入是不可以无限快的。最常见的目标函数 \bm{c}^T\bm{v} 就是最大化细胞的生长。解这个线性规划问题，就可以得到细胞内的代谢流的分布。很多时候解不唯一。所以有的时候会用多个线性规划，去找理想的解。原理是，进化导致了代谢流分布倾向于高效，所以高效的解更符合实际情况。

上面的例子的公式最后可以写成

图片来自书Optimization methods in metabolic networks

代谢流中的逆优化问题

为了能够更好的得到合适的目标函数，不同的人提出了不同的目标函数。比如，最大生长速率、最大熵增、ATP最大合成速率、CO2最大合成速率等。这些目标函数用在代谢流里面，某种程度上都可以得到符合实验的数据。

• 但是如何得到一个最适的目标函数呢？

• 如何能使得目标函数，这个LP得到唯一解呢？

• 如何得到合适的反应速率约束条件呢？

这个时候，逆优化就能起作用了。一个优化方面的大牛教授，就建议我用逆优化去试试这些问题。在有代谢流的数据（ \bm{v} ）下，能不能反过来去找合适的 \bm{c} ，还有反应速率上下限的函数，\bm{l} 和 \bm{u}。

上文的是最基本的代谢流模型，真正在工业中有应用价值的是动态代谢流守恒模型（Dynamic FBM），可以利用这样的模型，去分析优化培养基、控制监控发酵过程。如果又引入一些，热力学或者其他复杂的约束条件或者目标函数，这个问题会变成高度非线性的逆优化问题。

参考文献

Maranas C D, Zomorrodi A R. Optimization methods in metabolic networks[M]. John Wiley & Sons, 2016.

化工、制药、过控、优化的学习心得体会