CN104318090A - 基于最小二乘法支持向量机的溶菌酶发酵过程中的广义预测方法 - Google Patents

Abstract

本发明公开了一种基于最小二乘法支持向量机的溶菌酶发酵过程中的广义预测方法，该预测控制包括以下步骤：建立非线性预测模型，利用罐发酵被筛选过的产率较高的生产数据对最小二乘法支持向量机进行训练；将输入输出的非线性预测模型实时线性化，设定参考轨迹，滚动优化控制器设计，将基于LS-SVM的溶菌酶发酵过程中的广义预测控制算法智能的嵌入到上位机中；本发明采用最小二乘法支持向量机和广义预测控制相结合，使得模型在求解过程中避免了求解耗时的QP问题，运算简单，收敛速度快，精度高。将遗传算法和广义预测控制中的滚动优化相结合使得系统鲁棒性变强、能够有效克服系统滞后、干扰。

Description

基于最小二乘法支持向量机的溶菌酶发酵过程中的广义预测方法

技术领域

本发明涉及微生物发酵过程的先进控制领域，尤其是涉及基于最小二乘法支持向量机(LS-SVM)的溶菌酶发酵过程补料广义预测控制方法。

技术背景

微生物发酵过程中高度非线性、时变性和不确定性等多变量耦合系统，由于涉及生命体的生长和繁殖，机理十分复杂。发酵过程的建模是发酵工程的一个基本问题，它是为发酵过程的控制和优化服务的，准确的建模有利于得到更好的控制策略和优化方法。支持向量机(简称SVM)是一种比较好的实现结构风险最小化思想的方法，它是统计学习理论中最年轻的部分。支持向量机以训练误差作为优化问题的约束条件，以置信范围值最小化作为优化目标，通过某种非线性映射，将输入向量映射到一个高维的特征空间，通过解一个线性约束的二次规划问题得到全局最优解并且约束数量等于样本容量，因此当样本的容量很大时，训练消耗的时间很长。而最小二乘支持向量机用训练误差的平方代替了松弛变量，并且用等式约束代替不等式约束，训练过程只要求求解一个线性方程组，避免了求解耗时的QP问题，运算简单，收敛速度快，精度高，因而在模式识别和非线性过程建模等领域得到广泛应用。

广义预测控制(简称GPC)是一种鲁棒性强、能够有效克服系统滞后、应用于开环不稳定非最小相位系统的先进控制算法，它既吸收了自适应控制适用于随机系统、在线辨识等优点，又拥有预测控制算法中的滚动优化策略、对模型要求不高等优点，采用多步预测、动态优化和反馈校正等策略，控制效果相对而言更好，适用于发酵生产过程的控制。

发明内容

本发明的目的在于提供一种基于最小二乘法支持向量机的溶菌酶发酵过程中的广义预测方法及控制系统。

本发明解决其技术问题所采用的技术方案是：

权利要求书修改好后再添加技术方案：

基于最小二乘法支持向量机的溶菌酶发酵过程中的广义预测方法，包括以下步骤：

步骤一，利用罐发酵被筛选过的产率较高的生产数据对最小二乘法支持向量机进行训练，建立溶菌酶发酵过程中的非线性预测模型；

步骤二，将输入输出的非线性预测模型实时线性化，模型在采样点处利用泰勒公式展开得到线性预测模型，将模型经过转化成CARIMA模型，并通过引入丢番图方程得到溶菌酶菌发酵过程中的广义预测模型；

步骤三，在当前时刻k，广义预测模型利用历史和未来发酵罐的输入输出信息及未来的输入值来预测系统未来的输出状态菌体浓度产物浓度基质浓度

步骤四，设定参考轨迹，使发酵过程中在未来输出能沿着设定轨迹平稳的到达设定值，参考轨迹一般采用如下形式：

y_r(k+j)＝β^jy(k)+(1-β^j)y_s

其中y_r为菌体浓度参考轨迹，y_s为未来时刻的设定值，β为调节因子；

步骤五，滚动优化控制器设计，分别将实际输出y_b(k+j)与参考轨迹输出y_br(k+j)的误差、预测模型输出与实际输出y_b(k+j)的误差进行滚动优化，将滚动时域下的二次型目标函数变成适应值函数，通过遗传算法进行优化获得的控制量为全局最优；选取酸碱度pH为u₁(k)、温度T为u₂(k)、电机转速n为u₃(k)、溶氧量DO为u₄(k)、通气量Q为u₅(k)、压力为u₆(k)作为控制量并分别计算出当前时刻它们加入发酵罐的最优控制量；

步骤六，将基于最小二乘法支持向量机的溶菌酶发酵过程中的广义预测控制算法智能的嵌入到上位机中，根据上位机计算出来的各参数最优控制量通过PLC传输给执行机构从而实现对各参数的调节与控制。

进一步，所述步骤一在训练获得最小二乘法支持向量机的支持向量集时采用修剪算法。

进一步，所述步骤一采用径向基函数RBF作为核函数，通过最小二乘法支持向量机最终得到被控系统的非线性模型：

$y\left(x\right)=\underset{l=1}{\overset{N}{\mathrm{\Σ}}}{\mathrm{\α}}_l\exp \{-{\left|\right|x-x_l\left|\right|}_2^2/{\mathrm{\σ}}^2\}+\mathrm{\δ}$

其中为径向基函数，σ为核函数的核宽，α_l为拉尔朗日乘子，δ为常值偏差。

进一步，所述步骤五中没有给出反馈或闭环的表示。

进一步，所述步骤五中选取滚动时域下的二次型目标函数为：

$J\left(k\right)={\left|\right|Y\left(k\right|k)-Y_s\left(k\right)\left|\right|}_W^2+{\left|\right|\mathrm{\ΔU}\left(k\right|k)\left|\right|}_Q^2$

其中W和Q为对称矩阵，y_s(k+l)为输出在未来k+j时刻的参考值；Y_s(k)＝[y_s(k)^T]^T；y_s(k)＝[y_s(k+1)，y_s(k+2)，…，y_s(k+p)]^T。

进一步，所述步骤五中的遗传算法的具体步骤为：

步骤6.1，适应度函数的选择，其值在区间[0，1]中变化；

步骤6.2，遗传算法编码，在遗传算法优化中应用浮点数编码方法；

步骤6.3，生成初始种群；

步骤6.4，选择操作，选择操作的原则是按其适应函数值大小决定并采用比例选择的方法产生下一代种群；

步骤6.5，交叉操作和变异操作；

步骤6.6，终止判定，遗传算法终止计算是将两种方法相结合：一种是通过适度函数的值判断是否终止计算，另外一种是根据适度函数的收敛来判断。

进一步，所述步骤六具体为：

系统将采用PLC作为下位机，通过CAN总线与上位机通讯，各参数的检测仪器有温度传感器、pH传感器、DO传感器、编码器，检测元件的功能是把各参数转变成PLC可识别的变量，PLC获取各参数的测量后，值与设定值比较，输出各参数的最优控制量到各执行机构，从而实现对各参数的测量与控制。

本发明具有的有益效果是：根据发酵过程本身是一个复杂的、不确定的、非线性的时变动态过程，由于涉及生命体的生长和繁殖，机理十分复杂。本专利基于最小二乘法支持向量机的溶菌酶发酵过程的广义预测控制系统及方法，采用最小二乘法支持向量机和广义预测控制相结合，使得模型在求解过程中避免了求解耗时的QP问题，运算简单，收敛速度快，精度高。将遗传算法和广义预测控制中的滚动优化相结合使得系统鲁棒性变强、能够有效克服系统滞后、干扰。

附图说明

图1中基于最小二乘法支持向量机溶菌酶发酵过程中的广义预测控制的基本框图；

图2中基于遗传算法的流程图；

图3中pH、DO、T、n、P、Ω的最优控制量系统调节与控制框图；

图4中发酵罐控制系统的硬件连接图。

具体实施方式

下面详细描述本发明的实施例，所述实施例的示例在附图中示出，其中自始至终相同或类似的标号表示相同或类似的元件或具有相同或类似功能的元件。下面通过参考附图描述的实施例是示例性的，仅用于解释本发明，而不能理解为对本发明的限制。

下面按照图1所示的基本框架作详细说明。

(1)建立基于最小二乘支持向量机的非线性模型

在发酵过程中将溶菌酶发酵可以用下面的离散形式表示：

y(k+1)＝f[y_b(k)，y_b(k-1)，…，y_b(k-n+1)，u_i(k)，u_i(k-1)…，u_i(k-m+1)]

b＝1，2，3；i＝1，2，3，4，5，6；    (1)

其f是未知的一非线性函数，{y_b(k)}对应于发酵过程中的菌体浓度y₁(k)、产物浓度y₂(k)、基质浓度y₃(k)，控制量{u(k)}对应发酵过程中的酸碱度pH为u₁(k)、温度T为u₂(k)、电机转速n为u₃(k)、溶氧量DO为u₄(k)、通气量Q为u₅(k)、压力为u₆(k)；n，m分别为时间序列{y(k)}和{u(k)}的阶次。为了方便计算，现考虑SISO情况下的建模方法。

采用最小二乘支持向量机方法，该方法采用如下形式的函数对(1)进行估计：

y(x)＝w^TΦ(x)+δ  w，Φ(x)∈Rⁿ，δ∈R    (2)

式中，是核空间映射函数，将原始空间的输入数据映射到高维特征空间；w是参数列矢量，δ是常值偏差。

给定溶菌酶的发酵训练数据集得到LS-SVM回归的优化问题：

$\underset{w,b,e}{\min }J(w,e)=\frac{1}{2}{w}^{T}w+\frac{P}{2}\underset{l=1}{\overset{N}{\mathrm{\Σ}}}{e}_l^2,p>0---\left(3\right)$

约束条件

y_l＝w^Tl(x_l)+δ+e_l  l＝1，2，…，N    (4)

式中，权矢量误差变量e_l∈R；p是正则化参数。

根据优化函数，定义拉格朗日函数为：

$L(w,\mathrm{\δ},e;\mathrm{\α})=J(w,e)-\underset{l}{\overset{N}{\mathrm{\Σ}}}{\mathrm{\α}}_l[w^T\mathrm{\Φ}\left(x_l\right)+\mathrm{\δ}+e_l-y_l]---\left(5\right)$

式中，拉格朗日乘子α_l∈R，根据KKT条件得到：

$\left\{\begin{array}{c}\frac{\∂L}{\∂w}=0\→w=\underset{l=1}{\overset{N}{\mathrm{\Σ}}}{\mathrm{\α}}_l\mathrm{\φ}\left(x_l\right)=0\\ \frac{\∂L}{\∂e_l}=0\→{\mathrm{\α}}_i={\mathrm{pe}}_i\\ \frac{\∂L}{\∂{\mathrm{\α}}_l}=0\→w^T\mathrm{\φ}\left(x_l\right)+\mathrm{\δ}+e_l-y_l=0\\ \frac{\∂L}{\∂\mathrm{\δ}}=0\→\underset{l=1}{\overset{N}{\mathrm{\Σ}}}{\mathrm{\α}}_l=0\end{array}\right.---\left(6\right)$

得到矩阵方程：

$\left(\begin{array}{cc}0& 1^T\\ 1& \mathrm{\Γ}+P^{-1}I\end{array}\right)\left(\begin{array}{c}\mathrm{\δ}\\ \mathrm{\α}\end{array}\right)=\left(\begin{array}{c}0\\ y\end{array}\right)---\left(7\right)$

其中

y＝[y₁，y₂，…，y_N]^T，1＝[1，1，…，1]^T，α＝[α₁，α₂，…，α_N]^T，I为单位矩阵，

Γ_il＝φ(x_i)^Tφ(x_l)＝K(x_l，x_i)，K(·，·)为核函数。在式(7)中求出δ，α后可进一步得到w，得到LS-SVM的估计非线性模型如下：

$f\left(x\right)=\underset{l=1}{\overset{N}{\mathrm{\Σ}}}{\mathrm{\α}}_{l}K(x,x_l)+\mathrm{\δ}---\left(8\right)$

SVM的稀疏性具有很强的表示优势，对于新的输入，采用由少量的支持向量构成的决策或回归函数，能够大大地减少计算量。但在式(3)中的α的所有元素都不为零，使得LS-SVM的最大缺点恰恰是丧失了稀疏性，为了得到稀疏的支持向量集，采用一种修剪的算法，由于|α_l|的大小反映了每个训练样本在解向量构成中的相对重要性，则可将|α_l|降序排列后得到支持的向量集，可以删除那些|α_l|值较小的训练样本，再对新得到的训练集进行训练，就可以获得较好的稀疏性。

本专利使用如下RBF(径向基函数)作为核函数，最终LS-SVM的最终得到被控系统的非线性模型

$y\left(x\right)=\underset{l=1}{\overset{N}{\mathrm{\Σ}}}{\mathrm{\α}}_l\exp \{-{\left|\right|x-x_l\left|\right|}_2^2/{\mathrm{\σ}}^2\}+\mathrm{\δ}---\left(9\right)$

构造学习数据集(l＝1，2，…，N)。通过训练集中10批共152组生产数据训练学习，LS-SVM获得了溶菌酶发酵过程中的动态特性，及LS-SVM非线性模型，即求出δ＝-5.8878，p＝600，σ²＝10.24等特征参数；将测试采用其中7批进行训练，剩余3批数据作为检验样本用来测试得到非线性模型。

(2)实时线性化

令l＝1，2，3，…，N，将F_l(x)利用泰勒公式展开，将式(10)在点x₀处线性化，则有

$\begin{array}{c}{F}_l\left(x\right)=F_l\left(x_0\right)+\frac{\∂F_l}{\∂x\left(1\right)}{|}_{x=x_0}[x\left(1\right)-x_0\left(1\right)]+...+\frac{\∂F_l}{\∂x(m+n)}{|}_{x=x_0}[x(n+m)-x_0(n+m)]\\ =\frac{\∂F_l}{\∂x\left(1\right)}{|}_{x=x_0}x\left(1\right)+...+\frac{\∂F_l}{\∂x(m+n)}{|}_{x=x_0}x(n+m)+E_l\end{array}$

$E_l=F_l\left(x_0\right)-\frac{\∂F_l}{\∂x\left(1\right)}{|}_{x=x_0}{x}_0\left(1\right)-...-\frac{\∂F_l}{\∂x(n+m)}{|}_{x=x_0}{x}_0(n+m)---\left(10\right)$

其中E_l为当前采样周期某一常数，其值与x₀有关，通常取x₀＝0，这时E_l＝F_l(x0)。

将式(11)代入式(10)得到

$y\left(x\right)=\underset{l=1}{\overset{N}{\mathrm{\Σ}}}{\mathrm{\α}}_l\∂F_l^{T}x+\underset{l=1}{\overset{N}{\mathrm{\Σ}}}{\mathrm{\α}}_l{E}_l+\mathrm{\δ}---\left(11\right)$

其中， $\∂F_l={[\frac{\∂F_l}{\∂x\left(1\right)}{|}_{x=x_0},...,\frac{\∂F_l}{\∂x(n+m)}{|}_{x=x_0}]}^T.$

设x(k)是回归数据向量

x(k)＝[u(k)，u(k-1)，…，u(k-m+1)，y(k)，y(k-1)，…y(k-n+1)]    (12)

将式(12)写成离散差分方程形式

$y\left(k\right)=b_{1}u(k-1)+...+b_{m}u(k-m)-a_{1}y(k-1)-...a_{n}y(k-n)+\underset{l}{\overset{N}{\mathrm{\Σ}}}{\mathrm{\α}}_l{E}_l+\mathrm{\δ}---\left(13\right)$

其中， $\begin{array}{cc}{b}_1=\underset{l=1}{\overset{N}{\mathrm{\Σ}}}{\mathrm{\α}}_l\∂F_l\left(1\right),& a_1=\underset{l=1}{\overset{N}{\mathrm{\Σ}}}{\mathrm{\α}}_l\∂F_l(m+1)\\ .& .\\ .& .\\ .& .\\ b_m=\underset{l=1}{\overset{N}{\mathrm{\Σ}}}{\mathrm{\α}}_l\∂F_l\left(m\right),& a_n=\underset{l=1}{\overset{N}{\mathrm{\Σ}}}{\mathrm{\α}}_l\∂F_l(n+m)\end{array}$

经过简单的模型转换将常数消除得到

A(z^-1)y(k)＝B(z^-1)u(k-1)+δ    (14)

其中

A(z^-1)＝1+a₁z^-1+…+a_nz^-n

B(z^-1)＝1+b₁z^-1+…+b_mz^-m

(3)建立广义预测控制的线性模型

在GPC中采用了最小方差控制采用的CARIMA(受控自回归积分滑动平均模型)模型来描述受到随机干扰的装置。根据当前时刻的输入输出线性模型(14)考虑SISO的情况下的预测模型，将如下离散差分方程描绘系统特性：

A(z^-1)y(t)＝B(z^-1)u(t-1)+γ(t)    (15)

其中γ(t)是表示式(15)拟合当前时刻系统特性引起的误差及其他扰动，包含线性化过程中得到的常量δ；用对γ(t)建模，其中，Δ＝1-z^-1为差分算子；ε(t)是一个不相关的随机序列，表示一类随机噪声的影响。于是，系统当前时刻的输入输出方程为：

$A\left(z^{-1}\right)y\left(t\right)=B\left(z^{-1}\right)u(t-1)+\frac{\mathrm{\ϵ}\left(t\right)}{\mathrm{\Δ}}---\left(16\right)$

z^-1是后移算子，即z^-1y(k)＝y(k-1)，z^-1u(k)＝u(k-1)；设预测时域长度为P，则预测步数j＝1，2，…，P，控制时域为L(L≤P)；为了利用式(16)导出j步之后的输出y(t+j)的预测值，先要引入丢番图方程：

1＝E_j(z^-1)A(z^-1)Δ+z^-jF_j(z^-1)    (17)

其中，E_j(z^-1)，F_j(z^-1)是由A(z^-1)和预测长度j唯一确定的多项式，表达为

E_j(z^-1)＝e_j，0+e_j，1z^-1+…e_j，j-1z^-(j-1)

F_j(z^-1)＝f_j，0+f_j，1z^-1+…f_j，nz^-n

在式(16)两端乘以E_j(z^-1)Δz^j后得到

E_jAΔy(t+j)＝E_jBΔu(t+j-1)+E_jε(t+j)    (18)

并利用式(17)，可以写出k+j时刻的输出预测值

y(k+j|k)＝E_j(z^-1)B(z^-1)Δu(k+j-1|k)+F_j(z^-1)y(k)+E_j(z^-1)ξ(k+j)    (19)

由于在k时刻未来的噪声ξ(k+j)，i∈{1，…，j}都是未知的，所以对y(k+j)最合适的预测值可由下式得到

${\hat{y}}_m(k+j|k)=E_j\left(z^{-1}\right)B\left(z^{-1}\right)\mathrm{\Δu}(k+j-1|k)+F_j\left(z^{-1}\right)y\left(k\right)---\left(20\right)$

在式(20)中，记G_j(z^-1)＝E_j(z^-1)B(z^-1)。结合式(17)可得

$G_j\left(z^{-1}\right)=\frac{B\left(z^{-1}\right)}{A\left(z^{-1}\right)\mathrm{\Δ}}[1-z^{-j}{F}_j\left(z^{-1}\right)]---\left(21\right)$

再引入另一丢番图方程

$G_j\left(z^{-1}\right)=E_j\left(z^{-1}\right)B\left(z^{-1}\right)={\tilde{G}}_j\left(z^{-1}\right)+z^{-(j-1)}{H}_j\left(z^{-1}\right)---\left(22\right)$

其中

${\tilde{G}}_j\left(z^{-1}\right)=g_{j,0}+g_{j,1}{z}^{-1}+...+g_{j,j-1}{z}^{-(j-1)}$

H_j(z^-1)＝h_j，1z^-1+h_j，2z^-2+…+h_j，mz^-m

根据clarke给出的一个解丢番图方程采用递推方式求出E_j，F_j；则由式(21)和式(22)可以得到单输入预测以后的模型为：

${\tilde{y}}_m(k+j|k)={\tilde{G}}_j\mathrm{\Δu}(k+j-1|k)+H_j\mathrm{\Δu}\left(k\right)+F_{j}y\left(k\right)---\left(23\right)$

其中

$\hat{y}(k+j|k)={[\hat{y}(k+1),\hat{y}(k+2),...,\hat{y}(k+p)]}^T;$

Δu(k|k)＝[Δu(k|k)，…，Δu(k+p-1|k)]^T；

Δu(k)＝[Δu(k-1)，Δu(k-2)，…，Δu(k-m)]^T；

y(k)＝[y(k)，y(k-1)，…，y(k-n)]^T；

F＝[F₁(z^-1)F(z^-1)…F_p(z^-1)]^T；

$H=\left[\begin{array}{c}{G}_1-g_0\\ z(G_2-z^{-1}{g}_1-g_0)\\ .\\ .\\ .\\ z^{p-1}(G_p-z^{-p+1}-...-z^{-1}{g}_1-g_0)\end{array}\right];$

在当前时刻k，广义预测模型利用历史和未来发酵罐的输入输出信息及未来的输入值来预测系统未来k+j时刻预测输出的菌体浓度为产物浓度为基质浓度为在溶菌酶发酵过程中{y_b(k)}对应于发酵过程中k时刻的菌体浓度为y₁(k)、产物浓度为y₂(k)、基质浓度为y₃(k)，控制量{u(k)}对应发酵过程中的选取酸碱度pH为u₁(k)、温度T为u₂(k)、电机转速n为u₃(k)、溶氧量DO为u₄(k)、通气量Q为u₅(k)、压力为u₆(k)对于六输入三输出的预测控制系统将式(23)采用叠加的形式可以求出溶菌酶过程的GPC模型为：

Y(k|k)＝GΔU(k|k)+HΔU(k)+FY(k)    (24)

其中

$Y\left(k\right|k)={[{\hat{y}}_1{\left(k\right|k)}^T,{\hat{y}}_2\left(k\right|k),{\hat{y}}_3\left(k\right|k)]}^T;$

Y(k)＝[y₁(k)^T，y₂(k)^T，y₃(k)^T]^T； $G=\left[\begin{array}{cccccc}{\tilde{G}}_1& {\tilde{G}}_2& {\tilde{G}}_3& {\tilde{G}}_4& {\tilde{G}}_5& {\tilde{G}}_6\end{array}\right];$

ΔU(k|k)＝[Δu₁(k|k)^T，Δu₂(k|k)^T，Δu₃(k|k)^T，Δu₄(k|k)^T，Δu₅(k|k)，Δu₆(k|k)]^T；

F＝[F₁，F₂，F₃]；H＝[H₁ H₂ H₃ H₄ H₅ H₆]；

ΔU(k)＝[Δu₁(k)^T，Δu₂(k)^T，Δu₃(k)^T，Δu₄(k)，Δu₅(k|k)，Δu₆(k|k)]^T；

(4)设定参考轨迹

引入参考轨迹有利于发酵过程中在未来输出能沿着设定轨迹平稳的到达设定值。参考轨迹一般采用如下形式：

y_r(k+j)＝β^jy(k)+(1-β^j)y_s    (25)

其中y_r为菌体浓度的参考轨迹，y_s为未来时刻温度的设定值，β为调节因子，则通过上式可以引入参考输入轨迹菌体浓度y_1r(k+j)、产物浓度y_2r(k+j)、基质浓度y_3r(k+j)，上述参考轨迹的设定对闭环系统的动态特性和鲁棒性起很重要的作用。

(5)滚动优化控制器设计

所述滚动优化控制器设计中没有给出反馈或闭环的表示，因为它在滚动优化时，每一步都将实际输出与预测值进行比较来修正预测的不确定性，能够达到很好的预测效果。选取滚动时域下的二次型目标函数：

$J\left(k\right)={\left|\right|Y\left(k\right|k)-Y_s\left(k\right)\left|\right|}_w^2+{\left|\right|\mathrm{\ΔU}\left(k\right|k)\left|\right|}_Q^2---\left(26\right)$

其中W和Q为对称矩阵，y_s(k+l)为输出在未来k+j时刻的参考值(设定值)；Y_s(k)＝[y_s(k)^T]^T；y_s(k)＝[y_s(k+1)，y_s(k+2)，…，y_s(k+p)]^T。为了获取多步预测的控制律，提出了采用遗传算法求取全局最优解。

根据图2所示的具体步骤如下：

Step1：适应度函数的选择。式(26)是一个求最小值问题需要转化为遗传算法优化的最大值问题，采用遗传的适应度函数，其值在区间[0，1]中变化。这里将遗传算法的目标函数修改为：

$f(k,U)=\frac{1}{1+{\left|\right|Y\left(k\right|k)-Y_s\left(k\right)\left|\right|}_W^2+{\left|\right|\mathrm{\ΔU}\left(k\right|k)\left|\right|}_Q^2}$

约束条件为：

Y_min≤GΔU(k|k)+HΔU(k)+FY(k)≤Y_max

U_min≤U≤U_max，U＝[u₁，u₂，…，u_k+p1]          (27)

ΔU_min≤ΔU≤U_max，ΔU＝[Δu₁，Δu₂，…，Δu_k+p-1]

Step2：遗传算法编码。在遗传算法优化中应用浮点数编码方法，在GPC中每一个控制量可以用一个浮点数表示，每一个控制变量可用一个实数表示同时规定和加入约束，因此，每个染色体基因个数等于控制步长p来设定编码长度。个体的基因即为控制量Δu(k+j-1)，k时刻的编码结构如下：

Δu(k)

Δu(k+1)

…

Δu(k+p-2)

Δu(k+p-1)

Step3：生成初始种群。在生成初始种群中需要把约束通过空间限定法加入进来，通常约束采用的是上下限定法及仅在满足式(27)中约束条件随机选择n个个体来产生初始种群。当在k+1时刻初始种群生成的方法是将k时刻优化结果中的后p-1个基因左移一位，最后一个基因随机产生，构成k+1时刻的一个初始个体，该个体保存了前一个时刻K优化信息的个体。在初始种群中以一定数量比例λ放入该个体，为了保证初始种群的多样性，其余1-λ比例的个体随机产生。这种方法在保证了初始群体中优秀基因个体的同时保证群体的多样性，可以加速算法收敛的数度，减少计算时间。

Step4：选择操作。选择操作的原则是按其适应函数值大小决定并采用比例选择的方法产生下一代种群。具体操作如下：

(1)对群体的所有个体按其适应度大小进行降序排列；

(2)根据求解问题设计一个概率分配表，将各个概率值按上述排列次序分配给各个个体；

(3)以各个个体所分配的概率值作为遗传到一下代的概率，并用比例选择的方法产生下一代种群。

Step5：交叉操作和变异操作。交叉的染色体体数是由交叉概率p_c决定的，随机从初始种群中选中的染色体采用算法交叉的方法。变异操作按概率随机地在个体中选取一个体替代原来的个体。假设个体为和采用算法交叉则得到的新体为：

${\mathrm{\Ψ}}_1^{k+1}={\mathrm{\κ\Ψ}}_2^k+(1-\mathrm{\κ}){\mathrm{\Ψ}}_1^k$

${\mathrm{\Ψ}}_1^{k+1}={\mathrm{\κ\Ψ}}_2^k+(1-\mathrm{\κ}){\mathrm{\Ψ}}_1^k---\left(28\right)$

Step6：终止判定。遗传算法终止计算通常是将下面两种方法相结合。一种是通过适度函数的值判断是否终止计算，式(27)的适度函数其最大值为1，当遗传过程中出现适度函数为1的个体就是最优个体，停止计算。另外一种是根据适度函数的收敛来判断，给定一个代数c作为收敛程度判定。如果在取得的c中，每个群体中最优个体没有变化，则终止计算，选择该个体为最优个体。

上面2个条件中有一个满足时，便可停止计算选取选取最优个体作为当前最优解Δu*(k)；输出控制量，返回步骤3进行下次优化。其输出的最优控制量如下：

U′(k)＝U′(k-1)+ΔU^*(k)    (29)

实践表明，遗传算法对目标函数进行优化后，系统的鲁棒性得到了加强。

(6)控制系统的设计

如图4所示将基于最小二乘法溶菌酶发酵过程中的广义预测控制算法智能的嵌入到上位机中，采用PLC作为下位机，通过CAN总线与上位机通讯。根据基于遗传算法的滚动优化计算出的各参数的最优控制量U′₁(k)，U′₂(K)，U′₃(K)，U′₄(K)，U′₅(K)，U′6₍K)作为设定值通过CAN总线传输给检测元件。检测元件的功能是把各参数的最优控制量信号转变成PLC可识别的变量，PLC获取各参数的最优控制量，通过智能调节各执行机构包括固态继电器、电动调节阀、蠕动泵、变频器等，从而实现对各参数的调节控制量达到设定值，如图3所示给出了各个参数具体控制流程。本示例中采用的镇江日泰的RT-100L的发酵罐，具体的发酵罐配置为减速电机RF17DRE80S4、台达变频器VFD015M43B-A、pH变送器是PH300、DO变送器是DO300、蠕动泵TH15-100等。例如对pH的控制是通过检测元件把上位机传输过来的最优控制量作为设定值变成PLC可以识别的信号，pH传感器和变送器根据发酵罐中pH值输出4～20mA的电流信号，滤波后经通用模拟量模块实现A/D转换，转换成对应的数字量储存在PLC中并通过读取存储的值进行计算，便获得了发酵罐中当前的pH值，通过和设定值比较，执行预定的控制规律，输出数字量控制信号，数字量控制信号经放大后便可以控制蠕动泵的开和关，便可加入特定的物质使pH值稳定在设定值。

尽管上述已经示出和描述了本发明的实施例，本领域的普通技术人员可以理解：在不脱离本发明的原理和宗旨的情况下可以对这些实施例进行多种变化、修改、替换和变型，本发明的范围由权利要求及其等同物限定。

Claims (7)

1.基于最小二乘法支持向量机的溶菌酶发酵过程中的广义预测方法，其特征在于，包括以下步骤：

步骤一，利用罐发酵被筛选过的产率较高的生产数据对最小二乘法支持向量机进行训练，建立溶菌酶发酵过程中的非线性预测模型；

步骤二，将输入输出的非线性预测模型实时线性化，模型在采样点处利用泰勒公式展开得到线性预测模型，将模型经过转化成CARIMA模型，并通过引入丢番图方程得到溶菌酶菌发酵过程中的广义预测模型；

步骤三，在当前时刻k，广义预测模型利用历史和未来发酵罐的输入输出信息及未来的输入值来预测系统未来的输出状态菌体浓度产物浓度基质浓度

步骤四，设定参考轨迹，使发酵过程中在未来输出能沿着设定轨迹平稳的到达设定值，参考轨迹一般采用如下形式：

y_r(k+j)＝β^jy(k)+(1-β^j)y_s

其中y_r为菌体浓度参考轨迹，y_s为未来时刻的设定值，β为调节因子；

步骤五，滚动优化控制器设计，分别将实际输出y_b(k+j)与参考轨迹输出y_br(k+j)的误差、预测模型输出与实际输出y_b(k+j)的误差进行滚动优化，将滚动时域下的二次型目标函数变成适应值函数，通过遗传算法进行优化获得的控制量为全局最优；选取酸碱度pH为u₁(k)、温度T为u₂(k)、电机转速n为u₃(k)、溶氧量DO为u₄(k)、通气量Q为u₅(k)、压力为u₆(k)作为控制量并分别计算出当前时刻它们加入发酵罐的最优控制量；

步骤六，将基于最小二乘法支持向量机的溶菌酶发酵过程中的广义预测控制算法智能的嵌入到上位机中，根据上位机计算出来的各参数最优控制量通过PLC传输给执行机构从而实现对各参数的调节与控制。

2.根据权利1所述的基于最小二乘法支持向量机的溶菌酶发酵过程中的广义预测方法，其特征在于，所述步骤一在训练获得最小二乘法支持向量机的支持向量集时采用修剪算法。

3.根据权利要求1所述的基于最小二乘法支持向量机的溶菌酶发酵过程中的广义预测方法，其特征在于，所述步骤一采用径向基函数RBF作为核函数，通过最小二乘法支持向量机最终得到被控系统的非线性模型：

$y=\left(x\right)=\underset{l=1}{\overset{N}{\mathrm{\Σ}}}{\mathrm{\α}}_l\exp \{-{\left|\right|x-x_l\left|\right|}_2^2/{\mathrm{\σ}}^2\}+\mathrm{\δ}$

其中为径向基函数，σ为核函数的核宽，α_l为拉尔朗日乘子，δ为常值偏差。

4.根据权利要求1所述的基于最小二乘法支持向量机的溶菌酶发酵过程中的广义预测方法，其特征在于，所述步骤五中没有给出反馈或闭环的表示。

5.根据权利要求1所述的基于最小二乘法支持向量机的溶菌酶发酵过程中的广义预测方法，其特征在于，所述步骤五中选取滚动时域下的二次型目标函数为：

$J\left(k\right)={\left|\right|Y\left(k\right|k)-Y_s\left(k\right)\left|\right|}_W^2+{\left|\right|\mathrm{\ΔU}\left(k\right|k)}_Q^2$

其中W和Q为对称矩阵，y_s(k+l)为输出在未来k+j时刻的参考值；Y_s(k)＝[y_s(k)^T]^T；y_s(k)＝[y_s(k+1)，y_s(k+2)，…，y_s(k+p)]^T。

6.根据权利要求1所述的基于最小二乘法支持向量机的溶菌酶发酵过程中的广义预测方法，其特征在于，所述步骤五中的遗传算法的具体步骤为：

步骤6.1，适应度函数的选择，其值在区间[0，1]中变化；

步骤6.2，遗传算法编码，在遗传算法优化中应用浮点数编码方法；

步骤6.3，生成初始种群；

步骤6.4，选择操作，选择操作的原则是按其适应函数值大小决定并采用比例选择的方法产生下一代种群；

步骤6.5，交叉操作和变异操作；

步骤6.6，终止判定，遗传算法终止计算是将两种方法相结合：一种是通过适度函数的值判断是否终止计算，另外一种是根据适度函数的收敛来判断。

7.根据权利要求1所述的基于最小二乘法支持向量机的溶菌酶发酵过程中的广义预测方法，其特征在于，所述步骤六具体为：

系统将采用PLC作为下位机，通过CAN总线与上位机通讯，各参数的检测仪器有温度传感器、pH传感器、DO传感器、编码器，检测元件的功能是把各参数转变成PLC可识别的变量，PLC获取各参数的测量后，值与设定值比较，输出各参数的最优控制量到各执行机构，从而实现对各参数的测量与控制。