CN110109430B - 一种间歇式啤酒发酵装置优化控制系统 - Google Patents

Abstract

本发明公开了一种间歇式啤酒发酵装置优化控制系统，根据进入啤酒发酵装置中物料浓度、温度，以及冷媒温度、流量信息自动的计算出给定目标下的最优目标函数、最优发酵时间、冷媒最优控制轨线值，并通过冷媒流量控制实现生产过程最优，实现乙醇产量最大化与批次时间最小化的综合目标。本发明的优化控制系统利用传感测量模块、人机交互模块得到的信息，调用内部存储的啤酒发酵装置优化模型和参数、执行内部的优化控制程序，计算得到给定目标函数下最优发酵时间和冷媒流量最优控制轨线值。然后利用其流量控制模对冷媒流量进行控制，确保冷媒实际流量与最优控制轨线值一致。

Description

一种间歇式啤酒发酵装置优化控制系统

技术领域

本发明属于化工生产过程节能优化技术领域，涉及一种间歇式啤酒发酵装置优化控制系统。

背景技术

啤酒生产过程主要分为：制麦、糖化、发酵、罐装四个部分。其中啤酒生产过程的发酵过程属于化工反应过程，是啤酒生产过程的重要环节，由于发酵过程需要的时间比较长，并且对乙醇和副产品的影响非常大，所以需要很好的进行控制。以往采用人工经验控制，导致啤酒的质量变化比较大，乙醇的产量也无法实现最大化。随着自动化装置的应用，啤酒生产过程自动化程度大幅度升高。但是由于发酵过程受到多种因素影响，发酵过程冷媒流量要根据啤酒发酵罐内部和环境状态变化进行调整，控制不好容易造成发酵过程乙醇产量降低和发酵过程批次时间太长，影响生产效益。而目前更多的是采用固定的控制方式，例如对啤酒发酵罐外侧冷媒的流量进行简单的固定值控制，导致物料浓度等参数变化或者目标函数不同时，乙醇产量不理想或者达到较理想产量目标时物耗多、时间长。

一般的啤酒发酵生产过程是典型批次生产过程，这类过程的动态特性模型可描述为多个微分-代数方程组，其优化控制问题可描述为含有多个微分-代数方程、多个轨线约束方程和多个上下界约束的动态优化问题。求解该类动态优化问题可以得到满足用户目标要求的最优控制轨线和最优发酵时间，但是该问题的求解相对复杂、实现起来也存在困难。

对于啤酒发酵过程的优化控制问题，求解该类优化控制问题采用较多的方法是动态规划算法，但是动态规划算法的求解效率低、精度不够，且对于状态变量含有边界约束和整体含有不等式约束的问题需要处理后才能求解，难以实际应用。本发明可根据进入啤酒发酵罐的物料浓度和冷媒温度等，得到满足用户设定目标的最优控制轨线，并以此为设定值，通过调节冷媒阀控制冷媒流量，实现乙醇产量增加和批次时间缩短的综合目标。

发明内容

本发明的目的就是提供一种间歇式啤酒发酵装置优化控制系统。

本发明包括传感测量模块、人机交互模块、中央控制单元、现场总线网络、A/D和D/A转换模块、流量控制模块，其中传感测量模块用于测量初始时刻进入啤酒发酵罐内葡萄糖浓度、麦芽糖浓度、麦芽三糖浓度、初始时刻啤酒发酵罐内温度，以及用于冷却啤酒发酵罐的冷媒流量和温度。人机交互模块主要用于输入离散网格数、优化求解器容许误差、优化目标权重、进入发酵罐物料体积、酵母添加量。中央控制单元主要用于采集测量数据、存储啤酒发酵装置的优化模型和参数、对模型进行离散化处理、进行优化计算以及发送控制指令。现场总线用于数据通信，A/D和D/A转换模块用于模数/数模转换，流量控制模块接收中央控制单元发送的冷媒最优流量轨线值，并将其作为设定值，将冷媒流量实际值作为反馈值，通过控制冷媒阀开度确保实际流量与给定值一致。

本发明给出的优化控制系统通过以下步骤实现间歇式啤酒发酵装置的最优控制：

步骤A1：用户或者工程师通过人机交互模块输入离散网格数、配置点个数、优化求解器容许误差、优化目标权重、进入发酵罐物料体积、酵母添加量。

步骤A2：传感测量模块测量进入啤酒发酵罐内的麦汁中葡萄糖浓度、麦芽糖浓度、麦芽三糖浓度、啤酒发酵罐内温度和冷媒温度，然后通过A/D转换模块发送给中央控制单元；

步骤A3：中央控制单元从人机交互模块读取用户或者工程师输入的参数、接收从A/D转换模块发送的数据，然后调用内部存储的啤酒发酵装置优化模型和参数、执行内部的优化控制程序，包括将啤酒发酵装置优化模型进行离散化处理，调用非线性求解模块进行优化计算，通过计算得到给定目标函数下最优发酵时间以及在此时间内最优控制轨线值，并将其转换成冷媒流量最优控制轨线值。

步骤A4：中央控制单元通过D/A转换模块将得到的冷媒流量最优控制轨线值发送给流量控制模块，流量控制模块则以该最优控制轨线值作为冷媒流量设定值，以通过传感测量模块得到的冷媒实际流量作为反馈值，然后调用内部的PID控制算法对冷媒阀进行调节，确保冷媒实际流量与设定值一致。

步骤A5：啤酒发酵罐内发酵时间达到最优发酵时间后，系统停止运行。当啤酒发酵罐内物料进入下一环节，并且上一环节得到麦汁进入啤酒发酵罐内后，根据用户指令，转入步骤A1，开始下一批次的最优控制操作。

进一步的，所述的中央控制单元，包括信息采集模块、啤酒发酵装置优化模型模块、模型离散模块、非线性求解模块、控制指令输出模块。其中信息采集模块主要用于接受传感测量模块数据和人机交互模块输入的参数。啤酒发酵装置优化模型模块主要用于存储啤酒发酵过程模型和参数。模型离散模块用于将含有微分-代数方程组的啤酒发酵装置优化模型转化为非线性规划形式的优化模型。非线性求解模块则用于求解得到的非线性规划形式的优化模型。控制指令输出模块则用于将得到的最优控制轨线发送给流量控制模块。

中央控制单元内啤酒发酵装置优化模型如式(1.1)～(1.9)：

u^L≤u(t)≤u^U (1.8)；

temp^L≤temp(t)≤temp^U (1.9)；

其中，t_f表示最优发酵时间，是待优化变量，ω表示平衡两个目标的权重；u(t)表示冷却速率最优控制轨线，u^U和u^L则表示其上下界；e(t_f)表示发酵到t_f时刻时乙醇浓度。e(t)表示乙醇浓度，e(0)表示初始时刻t₀时乙醇浓度，temp(t)表示啤酒发酵罐内温度，temp^U和temp^L则表示其上下界，xb(t)代表酵母浓度，g(t)代表葡萄糖浓度，m(t)代表麦芽糖浓度，n(t)代表麦芽三糖浓度，xb(0)、g(t₀)、m(t₀)、n(t₀)分别表示初始时刻t₀时相应物料的浓度。μ₁(g(t))、μ₂(m(t),g(t))和μ₃(n(t),m(t),g(t))分别为葡糖糖、麦芽糖和麦芽三糖吸收比函数，C_p表示啤酒发酵罐内物料热容，ρ表示啤酒发酵罐内物料密度，TEMP_c表示冷媒温度，

和

表示葡糖糖、麦芽糖和麦芽三糖反应过程的反应热，

Rx_m、Rx_n分别表示葡糖糖、麦芽糖和麦芽三糖对应酵母的化学计量产率，

分别表示葡糖糖、麦芽糖和麦芽三糖对应乙醇的化学计量产率。t为时间，

为关于时间的导数。

糖吸收率见式(1.10)～(1.12)：

其中V_g(t)、V_m(t)和V_n(t)分别表示葡糖糖、麦芽糖和麦芽三糖的最大吸收速率，K_g(t)、K_m(t)和K_n(t)分别表示葡糖糖、麦芽糖和麦芽三糖的米氏常数，K'_g(t)和K'_m(t)分别表示葡糖糖、麦芽糖的反应抑制常数，而且这些常数值均取决于温度。依赖Arrhenius温度有如式(1.13)所示关系：

R表示气体常数，V_i0表示初始时刻t₀所对应的V_i(t₀)值，k_i0 k′_i0和分别表示初始时刻t₀所对应的K_i(t)值和K′_i(t)值，

表示对应的最大速度活化能，

和

则表示对应相应米氏常数的活化能和对应相应抑制常数的活化能。式(1.13)中i表示g，m，n中的一个，当i＝g时，K_i(t)为K_g(t)，K′_i(t)为K'_g(t)，V_i(t)为V_g(t)；当i＝m时，K_i(t)为K_m(t)，K′_i(t)为K'_m(t)，V_i(t)为V_m(t)；当i＝n时，K_i(t)为K_n(t)，K′_i(t)为K'_n(t)，V_i(t)为V_n(t)；其中K'_n(t)为麦芽三糖的反应抑制常数。

中央控制单元内的模型离散模块采用如下方法将啤酒发酵装置优化模型离散为非线性规划问题：

将式(1.1)～(1.13)所示的啤酒发酵装置优化模型转换为式(2.1)～(2.8)所示的动态优化问题：

dz/dt＝f(z(t),y(t),u(t),t,p) (2.2)；

F(z(t),y(t),u(t),t,p)＝0 (2.3)；

z^L≤z(t)≤z^U (2.4)；

u^L≤u(t)≤u^U (2.5)；

y^L≤y(t)≤y^U (2.6)；

t₀≤t≤t_f (2.7)；

z(t₀)＝z⁰ (2.8)；

这里

表示标量目标函数，z(t)、y(t)和u(t)分别表示与时间t相关的微分状态变量、代数状态变量和控制变量值。t₀和t_f表示开始与终端时间，也是啤酒发酵过程开始与结束时间，p表示外界环境参数。z(t_f)、y(t_f)和u(t_f)则分别表示在终端时刻微分状态变量、代数状态变量和控制变量的值。dz/dt表示微分状态变量z(t)对时间t的导数。f(z(t),y(t),u(t),t,p)表示微分方程形式的动态方程，F(z(t),y(t),u(t),t,p)表示代数方程形式的过程轨线束方程，z⁰表示状态变量z(t)在t₀时刻的初值，z^L和z^U表示状态变量z(t)的下界和上界，u^L和u^U分别表示控制变量u(t)的下界和上界，y^L和y^U表示代数状态变量y(t)的下界和上界。

对于式(2.1)-(2.8)所示的动态优化问题，首先将时间区间[t₀,t_f]均匀离散化为ne个网格，每个网格的长度h_i表示为式(2.9)：

h_i＝(t_f-t₀)/ne,i＝1,...,ne (2.9)；

在每个网格内插入K个配置点，配置点的相对位置选择Radau方程的根[ρ₁,ρ₂，…,ρ_K]，在第i个网格内微分状态变量表示为式(2.10)：

代数状态变量表示为式(2.11)：

控制变量表示为：

这里，z_i-1,0表示z(t)在第i个网格内的初始值，h_i是第i个网格的长度，dz/dt_i,q表示在第i个网格第q个配置点处z(t)对时间的导数值，t_i-1表示第i个网格的初始时刻，Ω_q为关于时间的K阶多项式。

y_i,q和u_i,q分别表示在第i个网格第q个配置点处代数变量y(t)和控制变量u(t)的值，ψ_q表示在第i个网格第q个配置点的拉格朗日函数，其形式如式(2.13)：

其中，t_i,j表示在第i个网格第j个配置点处的时间，ρ_q和ρ_j表示第q个和j个Radau方程的根，且满足式(2.14)：

考虑到微分状态变量的连续性，在下一个网格微分状态变量的初值z_i,0等于前一个网格微分状态变量的终值，因此有式(2.15)：

根据以上离散策略，将式(2.1)～(2.8)形式的动态优化问题离散化为式(2.16)：

其中x为nu维变量，ff(x)和c(x)分别表示连续可微的目标函数和mu维约束方程。x中既包含离散化的待优控制变量u_i,q、终端时间t_f，也包含离散的各种状态变量，x^U和x^L分别表示变量的上下界约束。

中央控制单元中的非线性求解模块，采用如下方法对离散化的非线性优化问题进行求解：

步骤B1：对式(2.16)所示的优化问题，首先根据变量上下界约束和经验给x赋予初值。

步骤B2：在x第k次迭代，的迭代点x_k处对式(2.16)进行泰勒展开，k为大于等于零的整数，并忽略高次项和目标函数中的常数项，则原问题求解转为式(3.1)所示的QP子问题的迭代求解：

式中d_k为搜索方向，gg_k ^T和A_k ^T分别表示在x_k处目标函数的导数的转置和约束方程的雅克比矩阵，c_k表示在x_k处c(x_k)的值，W_k为拉格朗日函数的Hessian阵。其中拉格朗日函数为式(3.2)：

L(x,λ,v,π)＝ff(x)+λ^Tc(x)+v^T(x-x^U)-π^T(x-x^L) (3.2)；

λ、v和π分别表示与等式约束、上边界约束和下边界约束相关的拉格朗日乘子，λ^T、v^T和π^T分别表示相应乘子的转置。

步骤B3：在迭代过程中将搜索空间分解为两个子空间Y和Z，将式(3.1)表示的QP子问题转化为低维QP子问题。其中Z∈R^nu×(nu-mu)，由雅克比矩阵的零空间向量组成；Y∈R^nu×mu，由雅克比矩阵的值空间向量组成。在x_k处子空间Y和Z的值表示为Y_k和Z_k，Z_k满足式(3.3)：

A_k ^TZ_k＝0 (3.3)；

搜索方向d_k表示为式(3.4)：

d_k＝Y_kp_y+Z_kp_z (3.4)；

这里p_y和p_z表示值空间和零空间移动的矢量矩阵，且p_y∈R^mu，p_z∈R^nu-mu。将式(3.3)和式(3.4)带入到QP子问题(3.1)中的等式约束，得式(3.5)：

因此根据式(3.5)p_y被唯一确定为式(3.6)：

p_y＝-(A_k ^TY_k)^-1c_k (3.6)；

搜索方向为式(3.7)：

d_k＝-Y_k(A_k ^TY_k)^-1c_k+Z_kp_z (3.7)；

将式(3.7)中的搜索方向d_k代入到以上QP子问题(3.1)中，并去掉与变量p_z无关的常数项，则原QP子问题可表示为以p_z∈R^nu-mu为变量的以下QP子问题形式：

其中w_k为(nu-mu)×1矩阵

B_k为(nu-mu)×(nu-mu)矩阵

p_z为(nu-mu)维变量。

步骤B4：采用积极集方法求解式(3.8)获得p_z，然后根据式(3.4)获得搜索方向矢量的值d_k。

步骤B5：令x_k+1＝x_k+αd_k，得到下一个迭代点。其中α∈(0,1]，通过一维搜索方法求取。

步骤B6：求取搜索方向二范数norm(d_k,2)、拉格朗日函数梯度与搜索方向乘积的值和一阶优化条件值。如果一阶优化条件值小于容许误差ε，或者norm(d_k,2)的值和拉格朗日函数梯度与搜索方向乘积的值同时小于容许误差ε，则停止计算，得到最优目标函数值、最优发酵时间和在此时间段内控制变量值。然后将这些值发送给中央控制单元。否则令x_k＝x_k+1，转步骤B2，继续计算。

中央控制单元中的控制指令输出模块采用如下方法输出冷媒流量最优控制轨线：

步骤C1：中央控制单元利用非线性求解模块计算得到最优值信息，获得最优目标函数值、最优发酵时间和在此时间段内控制变量值，并与模型离散过程对应起来，得到[t₀,t_f]区间内t_i,j时刻对应的控制量u_i,j,u_i,j按照时间连接起来就是最优控制轨线值，其中i＝1,....ne；j＝1,…,K。

步骤C2：中央控制单元根据控制量u_i,j的最优控制轨线值得到冷媒流量的最优控制轨线值。控制量u_i,j表示冷媒对啤酒发酵罐的冷却速率，与冷媒的流量存在式(4.1)所示关系：

Q_i,j＝V_beeru_i,j/(ρ_coolC_pcool) (4.1)；

其中，Q_i,j表示冷媒在t_i,j时刻的流量，V_beer表示啤酒发酵罐内物料体积，ρ_cool表示冷媒密度，C_pcool表示冷媒热容。

中央控制单元根据式(4.1)得到不同时刻冷媒最优流量值Q_i,j，Q_i,j按照时间连接起来就是冷媒流量最优控制轨线值。

步骤C3：中央控制单元利用现场总线，通过A/D和D/A转换模块将得到的冷媒最优流量值Q_i,j和对应的时间发送给流量控制模块。

所述的流量控制模块接受到中央控制单元发送的冷媒流量最优控制轨线值，并以此作为对应时刻的设定值，以通过传感测量模块得到的冷媒实际流量作为反馈值，然后调用内部的PID控制算法对冷媒阀进行调节，使得冷媒流量实际值与设置值保持一致。

本发明根据啤酒发酵罐内初始物料浓度、温度等信息通过最优化计算获得使得目标函数最小的最优冷媒流量控制轨线，并能自动程式化的对冷媒流量进行控制，确保了生产过程最优。即乙醇产量最大化和批次时间最小化的目标综合。另外，用户可以通过改变乙醇产量和批次时间目标之间的权重来优化调整生产过程，使得效益最大化。

附图说明

图1是本发明的功能示意图；

图2是本发明的结构示意图；

图3是本发明的中央控制单元内部结构和流程图；

图4是实施例1获得的冷却速率最优控制轨线图；

图5是实施例1获得的乙醇产品浓度变化情况图。

具体实施方式

下面结合附图说明本发明的具体实施方式。

一种间歇式啤酒发酵装置优化控制系统，对于某间歇式啤酒发酵装置，从上一环节中进入本装置内物料中乙醇浓度为0mol/m³，为确保生产产品质量，啤酒发酵装置内物料温度控制在0℃到12℃之间，冷却速率控制在0kJ h^-1m^-3℃^-1到30kJ h^-1m^-3℃^-1之间。执行如图1～3所示步骤以实现该啤酒发酵装置的最优控制：

步骤A1：用户或者工程师通过人机交互模块输入的离散网格数ne为40、配置点个数K为3、优化求解器容许误差ε为1.0E-7、优化目标权重ω为3、进入发酵罐物料体积为20m³、酵母添加量为3500mol。

步骤A2：传感测量模块测量出进入啤酒发酵罐内的麦汁中葡萄糖浓度、麦芽糖浓度、麦芽三糖浓度分别为70mol/m³、220mol/m³和40mol/m³，啤酒发酵罐内温度和冷媒温度分别为8℃和0℃，然后通过A/D转换模块发送给中央控制单元；

步骤A3：中央控制单元从人机交互模块读取用户或者工程师输入的参数、接收从A/D转换模块发送的数据，然后调用内部存储的啤酒发酵装置优化模型和参数、执行内部的优化控制程序，包括将啤酒发酵装置优化模型进行离散化处理，调用非线性求解模块进行优化计算，通过计算得到给定目标函数下最优发酵时间以及在此时间内最优控制轨线，并将其转换成冷媒流量的最优控制轨线。

步骤A4：中央控制单元通过D/A转换模块将得到的冷媒流量最优控制轨线发送给流量控制模块，流量控制模块则以该最优控制轨线值作为冷媒流量设定值，以通过传感测量模块得到的冷媒实际流量作为反馈值，然后调用内部的PID控制算法对冷媒阀进行调节，确保冷媒实际流量与设定值一致。

步骤A5：啤酒发酵罐内发酵时间达到最优发酵时间后，本系统停止运行。当啤酒发酵罐内物料进入下一环节，并且上一环节得到麦汁进入啤酒发酵罐内后，根据用户指令，转入步骤A1，开始下一批次的最优控制操作。

中央控制单元采集到相关数据后，其存储的啤酒发酵装置优化模型如下：

0≤u(t)≤30 (1.8)；

0≤temp(t)≤12 (1.9)；

其中e(0)、xb(0)、g(t₀)、m(t₀)、n(t₀)的值分别为0、175、70、220和40mol/m³。模型中的物性参数如表1所示。

表1

所述的中央控制单元，包括信息采集模块、啤酒发酵装置优化模型模块、模型离散模块、非线性求解模块、控制指令输出模块。其中信息采集模块主要用于接受传感测量模块数据和人机交互模块输入的参数。啤酒发酵装置优化模型模块主要用于存储啤酒发酵过程模型和参数。模型离散模块用于将含有微分-代数方程组的啤酒发酵装置优化模型转化为非线性规划形式的优化模型。非线性求解模块则用于求解得到的非线性规划形式的优化模型。控制指令输出模块则用于将得到的最优控制轨线发送给流量控制模块。

中央控制单元内啤酒发酵装置优化模型如式(1.1)～(1.9)：

u^L≤u(t)≤u^U (1.8)；

temp^L≤temp(t)≤temp^U (1.9)；

其中，t_f表示最优发酵时间，是待优化变量，ω表示平衡两个目标的权重；u(t)表示冷却速率最优控制轨线，u^U和u^L则表示其上下界；e(t_f)表示发酵到t_f时刻时乙醇浓度。e(t)表示乙醇浓度，e(0)表示初始时刻t₀时乙醇浓度，temp(t)表示啤酒发酵罐内温度，temp^U和temp^L则表示其上下界，xb(t)代表酵母浓度，g(t)代表葡萄糖浓度，m(t)代表麦芽糖浓度，n(t)代表麦芽三糖浓度，xb(0)、g(t₀)、m(t₀)、n(t₀)分别表示初始时刻t₀时相应物料的浓度。μ₁(g(t))、μ₂(m(t),g(t))和μ₃(n(t),m(t),g(t))分别为葡糖糖、麦芽糖和麦芽三糖吸收比函数，C_p表示啤酒发酵罐内物料热容，ρ表示啤酒发酵罐内物料密度，TEMP_c表示冷媒温度，

和

表示葡糖糖、麦芽糖和麦芽三糖反应过程的反应热，

Rx_m、Rx_n分别表示葡糖糖、麦芽糖和麦芽三糖对应酵母的化学计量产率，

分别表示葡糖糖、麦芽糖和麦芽三糖对应乙醇的化学计量产率。t为时间，

为对时间的导数。

糖吸收率见式(1.10)～(1.12)：

其中V_g(t)、V_m(t)和V_n(t)分别表示葡糖糖、麦芽糖和麦芽三糖的最大吸收速率，K_g(t)、K_m(t)和K_n(t)分别表示葡糖糖、麦芽糖和麦芽三糖的米氏常数，K'_g(t)和K'_m(t)分别表示葡糖糖、麦芽糖的反应抑制常数，而且这些常数值均取决于温度。依赖Arrhenius温度有如式(1.13)所示关系：

R表示气体常数，V_i0表示初始时刻t₀所对应的V_i(t₀)值，k_i0和k′_i0分别表示初始时刻t₀所对应的K_i(t)值和K′_i(t)值，

表示对应的最大速度活化能，

和

则表示对应相应米氏常数的活化能和对应相应抑制常数的活化能。因为式(1.13)中i表示g，m，n中的一个，当i＝g时，K_i(t)为K_g(t)，K′_i(t)为K'_g(t)，V_i(t)为V_g(t)；当i＝m时，K_i(t)为K_m(t)，K′_i(t)为K'_m(t)，V_i(t)为V_m(t)；当i＝n时，K_i(t)为K_n(t)，K′_i(t)为K'_n(t)，V_i(t)为V_n(t)；其中K'_n(t)为麦芽三糖的反应抑制常数。

中央控制单元内的模型离散模块采用如下方法将啤酒发酵装置优化模型离散为非线性规划问题：

将式(1.1)～(1.13)所示的啤酒发酵装置优化模型转换为式(2.1)～(2.8)所示的动态优化问题：

dz/dt＝f(z(t),y(t),u(t),t,p) (2.2)；

F(z(t),y(t),u(t),t,p)＝0 (2.3)；

z^L≤z(t)≤z^U (2.4)；

u^L≤u(t)≤u^U (2.5)；

y^L≤y(t)≤y^U (2.6)；

t₀≤t≤t_f (2.7)；

z(t₀)＝z⁰ (2.8)；

这里

表示标量目标函数，z(t)、y(t)和u(t)分别表示与时间t相关的微分状态变量、代数状态变量和控制变量值。t₀和t_f表示开始与终端时间，也是啤酒发酵过程开始与结束时间，p表示外界环境参数。z(t_f)、y(t_f)和u(t_f)则分别表示在终端时刻微分状态变量、代数状态变量和控制变量的值。dz/dt表示微分状态变量z(t)对时间t的导数。f(z(t),y(t),u(t),t,p)表示微分方程形式的动态方程，F(z(t),y(t),u(t),t,p)表示代数方程形式的过程轨线束方程，z⁰表示状态变量z(t)在t₀时刻的初值，z^L和z^U表示状态变量z(t)的下界和上界，u^L和u^U分别表示控制变量u(t)的下界和上界，y^L和y^U表示代数状态变量y(t)的下界和上界。

对于式(2.1)-(2.8)所示的动态优化问题，首先将时间区间[t₀,t_f]均匀离散化为ne个网格(ne＝40)，每个网格的长度h_i表示为式(2.9)：

h_i＝(t_f-t₀)/ne,i＝1,...,ne (2.9)；

在每个网格内插入K个配置点(K＝3)，配置点的相对位置选择Radau方程的根[ρ₁，ρ₂，…，ρ_K]，在第i个网格内微分状态变量表示为式(2.10)：

代数状态变量表示为式(2.11)：

控制变量表示为：

这里，z_i-1,0表示z(t)在第i个网格内的初始值，h_i是第i个网格的长度，dz/dt_i,q表示在第i个网格第q个配置点处z(t)对时间的导数值，t_i-1表示第i个网格的初始时刻，Ω_q为关于时间的K阶多项式。

y_i,q和u_i,q分别表示在第i个网格第q个配置点处代数变量y(t)和控制变量u(t)的值，ψ_q表示在第i个网格第q个配置点的拉格朗日函数，其形式如式(2.13)：

其中，t_i,j表示在第i个网格第j个配置点处的时间，ρ_q和ρ_j表示第q个和j个Radau方程的根，且满足式(2.14)：

考虑到微分状态变量的连续性，在下一个网格微分状态变量的初值z_i,0等于前一个网格微分状态变量的终值，因此有式(2.15)：

根据以上离散策略，将式(2.1)～(2.8)形式的动态优化问题离散化为式(2.16)：

其中x为nu维变量，ff(x)和c(x)分别表示连续可微的目标函数和mu维约束方程。x中既包含离散化的待优控制变量u_i,q、终端时间t_f，也包含离散的各种状态变量，x^U和x^L分别表示变量的上下界约束。

中央控制单元中的非线性求解模块，采用如下方法对离散化的非线性优化问题进行求解：

步骤B1：对式(2.16)所示的优化问题，首先根据变量上下界约束和经验给x赋予初值。

步骤B2：在x第k次迭代，的迭代点x_k处对式(2.16)进行泰勒展开，k为大于等于零的整数，并忽略高次项和目标函数中的常数项，则原问题求解转为式(3.1)所示的QP子问题的迭代求解：

式中d_k为搜索方向，gg_k ^T和A_k ^T分别表示在x_k处目标函数的导数的转置和约束方程的雅克比矩阵，c_k表示在x_k处c(x_k)的值，W_k为拉格朗日函数的Hessian阵。其中拉格朗日函数为式(3.2)：

L(x,λ,v,π)＝ff(x)+λ^Tc(x)+v^T(x-x^U)-π^T(x-x^L) (3.2)；

λ、v和π分别表示与等式约束、上边界约束和下边界约束相关的拉格朗日乘子，λ^T、v^T和π^T分别表示相应乘子的转置。

步骤B3：在迭代过程中将搜索空间分解为两个子空间Y和Z，将式(3.1)表示的QP子问题转化为低维QP子问题。其中Z∈R^nu×(nu-mu)，由雅克比矩阵的零空间向量组成；Y∈R^nu×mu，由雅克比矩阵的值空间向量组成。在x_k处子空间Y和Z的值表示为Y_k和Z_k，Z_k满足式(3.3)：

A_k ^TZ_k＝0 (3.3)；

搜索方向d_k表示为式(3.4)：

d_k＝Y_kp_y+Z_kp_z (3.4)；

这里p_y和p_z表示值空间和零空间移动的矢量矩阵，且p_y∈R^mu，p_z∈R^nu-mu。将式(3.3)和式(3.4)带入到QP子问题(3.1)中的等式约束，得式(3.5)：

因此根据式(3.5)p_y被唯一确定为式(3.6)：

p_y＝-(A_k ^TY_k)^-1c_k (3.6)；

搜索方向为式(3.7)：

d_k＝-Y_k(A_k ^TY_k)^-1c_k+Z_kp_z (3.7)；

将式(3.7)中的搜索方向d_k代入到以上QP子问题(3.1)中，并去掉与变量p_z无关的常数项，则原QP子问题可表示为以p_z∈R^nu-mu为变量的以下QP子问题形式：

其中w_k为(nu-mu)×1矩阵

B_k为(nu-mu)×(nu-mu)矩阵

p_z为(nu-mu)维变量。

步骤B4：采用积极集方法求解式(3.8)获得p_z，然后根据式(3.4)获得搜索方向矢量的值d_k，其中式(3.8)表示的QP子问题的积极集求解方法为成熟技术。

步骤B5：令x_k+1＝x_k+αd_k，得到下一个迭代点。其中α∈(0,1]，通过一维搜索方法求取，该求取方法为成熟技术。

步骤B6：求取搜索方向二范数norm(d_k,2)、拉格朗日函数梯度与搜索方向乘积的值和一阶优化条件值。如果一阶优化条件值小于容许误差ε，或者norm(d_k,2)的值和拉格朗日函数梯度与搜索方向乘积的值同时小于容许误差ε，则停止计算，得到最优目标函数值、最优发酵时间和在此时间段内控制变量值。然后将这些值发送给中央控制单元。否则令x_k＝x_k+1，转步骤B2，继续计算。

所述的中央控制单元通过调用非线性求解模块，得到最优目标函数值、最优发酵时间和此时间段内控制变量值等参数，其中最优发酵时间为167.2小时，在终止时间乙醇产品浓度为877.805mol/m³，最优冷却速率轨线如图4所示，乙醇浓度变化情况如图5所示。

中央控制单元中的控制指令输出模块采用如下方法输出冷媒流量最优控制轨线：

步骤C1：中央控制单元利用非线性求解模块计算得到最优值信息，获得最优目标函数值、最优发酵时间和在此时间段内控制变量值，并与模型离散过程对应起来，得到[0，167.2]区间内t_i,j时刻对应的控制量u_i,j，u_i,j按照时间连接起来就是最优控制轨线值，其中i＝1,....40；j＝1,…,3。

步骤C2：中央控制单元根据控制量u_i,j的最优控制轨线值得到冷媒流量的最优控制轨线值。控制量u_i,j表示冷媒对啤酒发酵罐的冷却速率，与冷媒的流量存在式(4.1)所示关系：

Q_i,j＝V_beeru_i,j/(ρ_coolC_pcool) (4.1)；

其中，Q_i,j表示冷媒在t_i,j时刻的流量，V_beer表示啤酒发酵罐内物料体积，ρ_cool表示冷媒密度，C_pcool表示冷媒热容。

中央控制单元根据式(4.1)得到不同时刻冷媒最优流量值Q_i,j，Q_i,j按照时间连接起来就是冷媒流量最优控制轨线值。

步骤C3：中央控制单元利用现场总线，通过A/D和D/A转换模块将得到的冷媒最优流量值Q_i,j和对应的时间发送给流量控制模块。

所述的流量控制模块接受到中央控制单元发送的冷媒流量最优控制轨线值，并以此作为对应时刻的设定值，以通过传感测量模块得到的冷媒实际流量作为反馈值，然后调用内部的PID控制算法对冷媒阀进行调节，使得冷媒流量实际值与设置值保持一致。

Claims (2)

1.一种间歇式啤酒发酵装置优化控制系统，包括传感测量模块、人机交互模块、中央控制单元、现场总线网络、A/D和D/A转换模块、流量控制模块，其特征在于：所述的传感测量模块用于测量初始时刻进入啤酒发酵罐内葡萄糖浓度、麦芽糖浓度、麦芽三糖浓度、初始时刻啤酒发酵罐内温度，以及用于冷却啤酒发酵罐的冷媒流量和温度；人机交互模块主要用于输入离散网格数、优化求解器容许误差、优化目标权重、进入发酵罐物料体积、酵母添加量；中央控制单元主要用于采集测量数据、存储啤酒发酵装置的优化模型和参数、对模型进行离散化处理、进行优化计算以及发送控制指令；现场总线用于数据通信，A/D和D/A转换模块用于模数/数模转换，流量控制模块接收中央控制单元发送的冷媒最优流量轨线值，并将其作为设定值，将冷媒流量实际值作为反馈值，通过控制冷媒阀开度确保实际流量与给定值一致；

所述的优化控制系统通过以下步骤实现间歇式啤酒发酵装置的最优控制：

步骤A1：用户或者工程师通过人机交互模块输入离散网格数、配置点个数、优化求解器容许误差、优化目标权重、进入发酵罐物料体积、酵母添加量；

步骤A2：传感测量模块测量进入啤酒发酵罐内的麦汁中葡萄糖浓度、麦芽糖浓度、麦芽三糖浓度、啤酒发酵罐内温度和冷媒温度，然后通过A/D转换模块发送给中央控制单元；

步骤A3：中央控制单元从人机交互模块读取用户或者工程师输入的参数、接收从A/D转换模块发送的数据，然后调用内部存储的啤酒发酵装置优化模型和参数、执行内部的优化控制程序，包括将啤酒发酵装置优化模型进行离散化处理，调用非线性求解模块进行优化计算，通过计算得到给定目标函数下最优发酵时间以及在此时间内最优控制轨线值，并将其转换成冷媒流量最优控制轨线值；

步骤A4：中央控制单元通过D/A转换模块将得到的冷媒流量最优控制轨线值发送给流量控制模块，流量控制模块则以该最优控制轨线值作为冷媒流量设定值，以通过传感测量模块得到的冷媒实际流量作为反馈值，然后调用内部的PID控制算法对冷媒阀进行调节，确保冷媒实际流量与设定值一致；

步骤A5：啤酒发酵罐内发酵时间达到最优发酵时间后，系统停止运行；当啤酒发酵罐内物料进入下一环节，并且上一环节得到麦汁进入啤酒发酵罐内后，根据用户指令，转入步骤A1，开始下一批次的最优控制操作。

2.如权利要求1所述的一种间歇式啤酒发酵装置优化控制系统，其特征在于：所述的中央控制单元，包括信息采集模块、啤酒发酵装置优化模型模块、模型离散模块、非线性求解模块、控制指令输出模块；其中信息采集模块主要用于接受传感测量模块数据和人机交互模块输入的参数；啤酒发酵装置优化模型模块主要用于存储啤酒发酵过程模型和参数；模型离散模块用于将含有微分-代数方程组的啤酒发酵装置优化模型转化为非线性规划形式的优化模型；非线性求解模块则用于求解得到的非线性规划形式的优化模型；控制指令输出模块则用于将得到的最优控制轨线发送给流量控制模块；

中央控制单元内啤酒发酵装置优化模型如式(1.1)～(1.9)：

u^L≤u(t)≤u^U (1.8)；

temp^L≤temp(t)≤temp^U (1.9)；

其中，t_f表示最优发酵时间，是待优化变量，ω表示平衡两个目标的权重；u(t)表示冷却速率最优控制轨线，u^U和u^L则表示其上下界；e(t_f)表示发酵到t_f时刻时乙醇浓度；e(t)表示乙醇浓度，e(0)表示初始时刻t₀时乙醇浓度，temp(t)表示啤酒发酵罐内温度，temp^U和temp^L则表示其上下界，xb(t)代表酵母浓度，g(t)代表葡萄糖浓度，m(t)代表麦芽糖浓度，n(t)代表麦芽三糖浓度，xb(0)、g(t₀)、m(t₀)、n(t₀)分别表示初始时刻t₀时相应物料的浓度；μ₁(g(t))、μ₂(m(t),g(t))和μ₃(n(t),m(t),g(t))分别为葡萄糖、麦芽糖和麦芽三糖吸收比函数，C_p表示啤酒发酵罐内物料热容，ρ表示啤酒发酵罐内物料密度，TEMP_c表示冷媒温度，

和

表示葡萄糖、麦芽糖和麦芽三糖反应过程的反应热，

分别表示葡萄糖、麦芽糖和麦芽三糖对应酵母的化学计量产率，

分别表示葡萄糖、麦芽糖和麦芽三糖对应乙醇的化学计量产率；t为时间，

为关于时间的导数；

糖吸收率见式(1.10)～(1.12)：

其中V_g(t)、V_m(t)和V_n(t)分别表示葡萄糖、麦芽糖和麦芽三糖的最大吸收速率，K_g(t)、K_m(t)和K_n(t)分别表示葡萄糖、麦芽糖和麦芽三糖的米氏常数，K'_g(t)和K'_m(t)分别表示葡萄糖、麦芽糖的反应抑制常数，而且这些常数值均取决于温度；依赖Arrhenius温度有如式(1.13)所示关系：

R表示气体常数，V_i0表示初始时刻t₀所对应的V_i(t₀)值，k_i0和k′_i0分别表示初始时刻t₀所对应的K_i(t)值和K′_i(t)值，

表示对应的最大速度活化能，

和

则表示对应相应米氏常数的活化能和对应相应抑制常数的活化能；式(1.13)中i表示g，m，n中的一个，当i＝g时，K_i(t)为K_g(t)，K′_i(t)为K'_g(t)，V_i(t)为V_g(t)；当i＝m时，K_i(t)为K_m(t)，K′_i(t)为K'_m(t)，V_i(t)为V_m(t)；当i＝n时，K_i(t)为K_n(t)，K′_i(t)为K'_n(t)，V_i(t)为V_n(t)；其中K'_n(t)为麦芽三糖的反应抑制常数；

中央控制单元内的模型离散模块采用如下方法将啤酒发酵装置优化模型离散为非线性规划问题：

将式(1.1)～(1.13)所示的啤酒发酵装置优化模型转换为式(2.1)～(2.8)所示的动态优化问题：

dz/dt＝f(z(t),y(t),u(t),t,p) (2.2)；

F(z(t),y(t),u(t),t,p)＝0 (2.3)；

z^L≤z(t)≤z^U (2.4)；

u^L≤u(t)≤u^U (2.5)；

y^L≤y(t)≤y^U (2.6)；

t₀≤t≤t_f (2.7)；

z(t₀)＝z⁰ (2.8)；

这里

表示标量目标函数，z(t)、y(t)和u(t)分别表示与时间t相关的微分状态变量、代数状态变量和控制变量值；t₀和t_f表示开始与终端时间，也是啤酒发酵过程开始与结束时间，p表示外界环境参数；z(t_f)、y(t_f)和u(t_f)则分别表示在终端时刻微分状态变量、代数状态变量和控制变量的值；dz/dt表示微分状态变量z(t)对时间t的导数；f(z(t),y(t),u(t),t,p)表示微分方程形式的动态方程，F(z(t),y(t),u(t),t,p)表示代数方程形式的过程轨线束方程，z⁰表示状态变量z(t)在t₀时刻的初值，z^L和z^U表示状态变量z(t)的下界和上界，u^L和u^U分别表示控制变量u(t)的下界和上界，和y^L y^U表示代数状态变量y(t)的下界和上界；

对于式(2.1)-(2.8)所示的动态优化问题，首先将时间区间[t₀,t_f]均匀离散化为ne个网格，每个网格的长度h_i表示为式(2.9)：

h_i＝(t_f-t₀)/ne,i＝1,...,ne (2.9)；

在每个网格内插入K个配置点，配置点的相对位置选择Radau方程的根[ρ₁，ρ₂，…,ρ_K]，在第i个网格内微分状态变量表示为式(2.10)：

代数状态变量表示为式(2.11)：

控制变量表示为：

这里，z_i-1,0表示z(t)在第i个网格内的初始值，h_i是第i个网格的长度，dz/dt_i,q表示在第i个网格第q个配置点处z(t)对时间的导数值，t_i-1表示第i个网格的初始时刻，Ω_q为关于时间的K阶多项式；

y_i,q和u_i,q分别表示在第i个网格第q个配置点处代数变量y(t)和控制变量u(t)的值，ψ_q表示在第i个网格第q个配置点的拉格朗日函数，其形式如式(2.13)：

其中，t_i,j表示在第i个网格第j个配置点处的时间，ρ_q和ρ_j表示第q个和j个Radau方程的根，且满足式(2.14)：

考虑到微分状态变量的连续性，在下一个网格微分状态变量的初值z_i,0等于前一个网格微分状态变量的终值，因此有式(2.15)：

根据以上离散策略，将式(2.1)～(2.8)形式的动态优化问题离散化为式(2.16)：

其中x为nu维变量，ff(x)和c(x)分别表示连续可微的目标函数和mu维约束方程；x中既包含离散化的待优控制变量u_i,q、终端时间t_f，也包含离散的各种状态变量，x^U和x^L分别表示变量的上下界约束；

中央控制单元中的非线性求解模块，采用如下方法对离散化的非线性优化问题进行求解：

步骤B1：对式(2.16)所示的优化问题，首先根据变量上下界约束和经验给x赋予初值；

步骤B2：在x第k次迭代，的迭代点x_k处对式(2.16)进行泰勒展开，k为大于等于零的整数，并忽略高次项和目标函数中的常数项，则原问题求解转为式(3.1)所示的QP子问题的迭代求解：

式中d_k为搜索方向，gg_k ^T和A_k ^T分别表示在x_k处目标函数的导数的转置和约束方程的雅克比矩阵，c_k表示在x_k处c(x_k)的值，W_k为拉格朗日函数的Hessian阵；其中拉格朗日函数为式(3.2)：

L(x,λ,v,π)＝ff(x)+λ^Tc(x)+v^T(x-x^U)-π^T(x-x^L) (3.2)；

λ、v和π分别表示与等式约束、上边界约束和下边界约束相关的拉格朗日乘子，λ^T、v^T和π^T分别表示相应乘子的转置；

步骤B3：在迭代过程中将搜索空间分解为两个子空间Y和Z，将式(3.1)表示的QP子问题转化为低维QP子问题；其中Z∈R^nu×(nu-mu)，由雅克比矩阵的零空间向量组成；Y∈R^nu×mu，由雅克比矩阵的值空间向量组成；在x_k处子空间Y和Z的值表示为Y_k和Z_k，Z_k满足式(3.3)：

A_k ^TZ_k＝0 (3.3)；

搜索方向d_k表示为式(3.4)：

d_k＝Y_kp_y+Z_kp_z (3.4)；

这里p_y和p_z表示值空间和零空间移动的矢量矩阵，且p_y∈R^mu，p_z∈R^nu-mu；将式(3.3)和式(3.4)带入到QP子问题(3.1)中的等式约束，得式(3.5)：

因此根据式(3.5)p_y被唯一确定为式(3.6)：

p_y＝-(A_k ^TY_k)^-1c_k (3.6)；

搜索方向为式(3.7)：

d_k＝-Y_k(A_k ^TY_k)^-1c_k+Z_kp_z (3.7)；

将式(3.7)中的搜索方向d_k代入到以上QP子问题(3.1)中，并去掉与变量p_z无关的常数项，则原QP子问题可表示为以p_z∈R^nu-mu为变量的以下QP子问题形式：

其中w_k为(nu-mu)×1矩阵

B_k为(nu-mu)×(nu-mu)矩阵

p_z为(nu-mu)维变量；

步骤B4：采用积极集方法求解式(3.8)获得p_z，然后根据式(3.4)获得搜索方向矢量的值d_k；

步骤B5：令x_k+1＝x_k+αd_k，得到下一个迭代点；其中α∈(0,1]，通过一维搜索方法求取；

步骤B6：求取搜索方向二范数norm(d_k,2)、拉格朗日函数梯度与搜索方向乘积的值和一阶优化条件值；如果一阶优化条件值小于容许误差ε，或者norm(d_k,2)的值和拉格朗日函数梯度与搜索方向乘积的值同时小于容许误差ε，则停止计算，得到最优目标函数值、最优发酵时间和在此时间段内控制变量值；然后将这些值发送给中央控制单元；否则令x_k＝x_k+1，转步骤B2，继续计算；

中央控制单元中的控制指令输出模块采用如下方法输出冷媒流量最优控制轨线：

步骤C1：中央控制单元利用非线性求解模块计算得到最优值信息，获得最优目标函数值、最优发酵时间和在此时间段内控制变量值，并与模型离散过程对应起来，得到[t₀,t_f]区间内t_i,j时刻对应的控制量u_i,j，u_i,j按照时间连接起来就是最优控制轨线值，其中i＝1,....ne；j＝1,…,K；

步骤C2：中央控制单元根据控制量u_i,j的最优控制轨线值得到冷媒流量的最优控制轨线值；控制量u_i,j表示冷媒对啤酒发酵罐的冷却速率，与冷媒的流量存在式(4.1)所示关系：

Q_i,j＝V_beeru_i,j/(ρ_coolC_pcool) (4.1)；

其中，Q_i,j表示冷媒在t_i,j时刻的流量，V_beer表示啤酒发酵罐内物料体积，ρ_cool表示冷媒密度，C_pcool表示冷媒热容；

中央控制单元根据式(4.1)得到不同时刻冷媒最优流量值Q_i,j，Q_i,j按照时间连接起来就是冷媒流量最优控制轨线值；

步骤C3：中央控制单元利用现场总线，通过A/D和D/A转换模块将得到的冷媒最优流量值Q_i,j和对应的时间发送给流量控制模块；

所述的流量控制模块接受到中央控制单元发送的冷媒流量最优控制轨线值，并以此作为对应时刻的设定值，以通过传感测量模块得到的冷媒实际流量作为反馈值，然后调用内部的PID控制算法对冷媒阀进行调节，使得冷媒流量实际值与设置值保持一致。

Images