15.1 工業乙醇生產與計算機仿真
乙醇作為可再生清潔能源不僅可以代替四乙基鉛作為汽油的防爆劑，還可以制造汽油醇。這一巨大的潛在需求促使人們去尋找提高乙醇工業生產率的途徑，使人們着手於發酵工程的研究。微生物學、發酵研究的發展，使微生物反應過程的種類和規模不斷擴大，其應用也深入到多個工業領域，然而由於反應涉及活細胞的參與，菌體生長及產物生成等機理復雜多變，目前尚難為人們所了解和把握，更難以進行統一的描述。人們一般通過實驗的方法尋找微生物的生長規律，通過數值分析和擬合了解發酵過程的規律。這樣的做法周期性長，需要消耗的資源多，且無法在短時間內對工業流程中出現的問題提出及時的應對方案。
計算機仿真技術作為分析和研究系統運行行為、揭示系統動態過程和運動規律的一種重要手段和方法，隨着系統科學研究的深入，以及控制理論、計算技術、計算機科學與技術的發展而形成一門新興的學科。近年來，隨着信息處理技術的突飛猛進，仿真技術得到迅速發展。“仿真是一種基於模型的活動”，設計多學科、多領域的知識和經驗。成功進行仿真研究的關鍵是有機、協調地組織實施仿真全生命周期的各類活動。這里的“各類活動”，就是指“系統建模”、“仿真建模”、“仿真實驗”，聯系這些活動的要素是“系統”、“模型”、“計算機”。其中，系統是研究的對象，模型是系統的抽象，仿真是為了通過對模型的實驗達到跟真實系統對比並接近於研究其特性的目的。
15.2 工業乙醇發酵流程
A：預發酵罐
B~G：6個發酵罐
H：離心分離機
1、5、7、9：干法糖液
2：營養液
3：酵母
15：產品，含11%~12%的乙醇
14：廢酵母，其中一部分循環回去
發酵溫度31~33℃。
15.3 乙醇發酵動力學方程
- 葡萄糖液在酵母作用下分解產生乙醇；
- 酵母細胞本身吸收葡萄糖液來促進自身生長。
\[\begin{array}{l}
\frac{{{\rm{d}}x}}{{{\rm{d}}t}} = \mu x\\
\frac{{{\rm{d}}p}}{{{\rm{d}}t}} = \upsilon x\\
\frac{{{\rm{d}}s}}{{{\rm{d}}t}} = - (\alpha + \beta )x
\end{array}\]
- 對每個發酵罐，發酵液流入量=流出量；
- 對每種物質而言，其質量變化=流入/流出質量變化+反應生成/消耗質量變化。
\[\begin{array}{l}
{F_{in}} = {F_{out}}\\
\frac{{{\rm{d}}x}}{{{\rm{d}}t}} = \mu x + D({x_0} - x)\\
\frac{{{\rm{d}}p}}{{{\rm{d}}t}} = \upsilon x + D({p_0} - p)\\
\frac{{{\rm{d}}s}}{{{\rm{d}}t}} = - (\frac{\mu }{{{Y_{x/s}}}} + \frac{\upsilon }{{{Y_{p/s}}}})x + D({p_0} - p)
\end{array}\]
因為實際生產過程中基質濃度較高，需要考慮底物和產物對細胞生長的抑制作用。
\[\begin{array}{l}
\mu = {\mu _{\max }}\frac{S}{{{K_s} + S + {s^2}/{K_{sl}}}} \cdot \frac{{{K_p}}}{{{K_p} + p + {p^2}/{K_{pl}}}}\\
\upsilon = {\upsilon _{\max }}\frac{S}{{{K_{sp}} + S + {s^2}/{K_{spl}}}} \cdot \frac{{{K_{pp}}}}{{{K_{pp}} + p + {p^2}/{K_{ppl}}}}
\end{array}\]
發酵過程模型流程圖：
15.4 MATLAB求解
- [T,Y]=solver(odefun,tspan,y0,options,p1,p2,...)
參考：https://ww2.mathworks.cn/help/matlab/ref/ode45.html
solver為命令ode45、ode23、ode113、ode15s、ode23s、ode23t、ode23tb之一，odefun為顯式常微分方程dy/dt=f(t,y)或包含混合矩陣的方程M(t,y)×dy/dt=f(t,y)。命令ode23只能求解常數混合矩陣的問題；命令ode23t與ode15s可以求解奇異矩陣的問題。tspan為積分區間的向量，tspan=[t0,tf]。要獲得問題在指定時間點t0、t1、t2、...上的解，可以令tspan=[t0,t1,t2,...,tf]，但必須是單調的。y0為包含初值條件的向量，options為用命令odeset設置的可選積分參數。p1、p2、...為傳遞給函數odefun的可選參數。
function f=func(t,n,V) %發酵罐的動力學函數,輸入參數為時間,輸入矩陣,體積 global x0 p0 s0 Yxs Yps umax vmax Ks Kps Kp Kpp Ksi Kpsi Kpi Kppi T F Vy f=zeros(3,1); %按照三行擴展列來保存每一步的三種物料的濃度 u=umax*n(3)/(Ks+n(3)+n(3)^2/Ksi)*Kp/(Kp+n(2)+n(2)^2/Kpi); v=vmax*n(3)/(Kps+n(3)+n(3)^2/Kpsi)*Kpp/(Kpp+n(2)+n(2)^2/Kppi); D=Vy/V; f(1)=u*n(1)+D*(x0-n(1)); %Dx n(1)-x,n(2)-p,n(3)-s f(2)=v*n(1)+D*(p0-n(2)); %Dp f(3)=-(u/Yxs+v/Yps)*n(1)+D*(s0-n(3)); %Ds
也可以用龍格-庫塔法求解，四階示例如下。
function [x,y]=runge_kutta1(ufunc,y0,h,a,b,Vg) %參數表順序依次是微分方程組的函數名稱，初始值向量，步長，時間起點，時間終點,發酵罐體積 n=floor((b-a)/h); %求步數 x(1)=a; %時間起點 y(:,1)=y0; %賦初值，可以是向量，但是要注意維數 for ii=1:n x(ii+1)=x(ii)+h; k1=ufunc(x(ii),y(:,ii),Vg); k2=ufunc(x(ii)+h/2,y(:,ii)+h*k1/2,Vg); k3=ufunc(x(ii)+h/2,y(:,ii)+h*k2/2,Vg); k4=ufunc(x(ii)+h,y(:,ii)+h*k3,Vg); y(:,ii+1)=y(:,ii)+h*(k1+2*k2+2*k3+k4)/6; %按照龍格庫塔方法進行數值求解 end
>> [x,y]=runge_kutta1(@func,[10;0;200],0.1,0,15,500);
>> plot(x,y);
>> legend({'酵母菌體','乙醇','葡萄糖液'});
15.5 Simulink求解
Simulink常用的微分方程解法：
- 根據方程搭建模型，較為復雜，通用；
- 用S函數，需要求解狀態方程；
- 用傳輸函數，對時域方程進行拉普拉斯變換；
- 向量法，用Fcn模塊表示一個微分方程中的一位輸入，三個Fcn模塊並聯。
15.6 連續流程仿真
15.7 GUIDE設計
15.7.1 GUIDE介紹
在命令窗口輸入guide。
例：設計按鈕控制圖片顯示的界面。
拖出axes、pushbutton控件，右擊pushbutton，選擇View Callbacks→Callback。
% --- Executes on button press in pushbutton1.
function pushbutton1_Callback(hObject, eventdata, handles)
% hObject handle to pushbutton1 (see GCBO)
% eventdata reserved - to be defined in a future version of MATLAB
% handles structure with handles and user data (see GUIDATA)
axes(handles.axes1);
imshow('flower.jpg');
15.7.2 乙醇仿真界面
