Matcont with Ternary Reactive Distillation column

[Fatima]

Hello Everyone,

From the past few week I have been trying to learn Matcont for performing bifurcation analysis of Differential equations. But the examples are limited to only 3 or 4 equations. My column consists of 108 differential equations which in itself consists of many algebraic variables which needs to be computed beforehand using algebraic equations. How can I realize it using Matcont? I need to do bifurcation analysis in between VB (Vapour Boilup) and D (Distillate) to find out multiple steady states of the column using particular input values.

I will be very thankful, if any help is available.

below is the code of my model;

function xprime=rdcol_mod(~,X,U) 

% Inputs and disturbances
LT =    U(1);                            % Reflux
VB =    U(2);                            % Boilup
F(1) =  U(3);                         % Feedrate 1
F(2) =  U(4);                         % Feedrate 2
zF(1,1:4) = U(5:8);                                              % Feed composition
zF(2,1:4) = U(9:12);                                              % Feed composition

% Number of stages  (stages are counted from the bottom)  
    NS  = 5+1;               % number of stripping stages [---]  
    NRX = 15;                % number of reactive stages [---]
    NR  = 5;                 % number of rectifying stages [---]
    NT  = NS+NRX+NR+1 ;      % total number of stages  [---]
    NC  = 4;                 % number of components [---]

% Physical properties
Ef = 30000.;                             % activation energy forward [kcal/kmol]
Eb = 40000.;                             % activation energy backward [kcal/kmol]
Rc = 1.987;                              % gas law constant [kcal/kmolK]
Keq = 50.0;                              % equilibrium constant [---]
kf0 = 0.008;                             % rate coefficient / forward (for 366 K) [s^(-1)]
kr0 = kf0/Keq;                           % rate coefficient / backward (for 366 K) [s^(-1)]
af = kf0*exp(Ef/(Rc*366.));              % preexponential factor forward [1/s]
ab = kr0*exp(Eb/(Rc*366.));              % preexponential factor backward [1/s]
Lambda = -10.000;                        % heat of reaction [kcal/kmol]
Heatvap = 6.944;                         % heat of vaporization [kcal/kmol]
v = [-1 -1 1 0];                         % stoc. coeff. 
Avp=[12.34 11.65 10.96 12.34];
Bvp=[3862 3862 3862 3862];
P=8;
% alpha=[3.9749 1.9937 1 3.9749];
alpha = [4 2 1 4];                       % const. rel. volatilities
Mrx = 2000;                              % holdup of reactive trays [kmol]
M = 670;                                 % holdup of nonreactive trays [kmol]
MB0 = 23200;                             % holdup of reboiler [kmol]
MD0 = 24800;                             % holdup of condenser [kmol]

% Splitting the states
x(1,:)=X(1:NT)';
x(2,:)=X(NT+1:2*NT)';
x(3,:)=X(2*NT+1:3*NT)';
T=X(3*NT+1:4*NT)';

j=1:NT;
kf(j) = af*exp(-Ef./(Rc*T(j)));           % forward kinetic constant at mean temperature [1/s]
kb(j) = ab*exp(-Eb./(Rc*T(j)));           % backward kinetic constant at mean temperature [1/s]

%------------------------------------------------------------
% THE MODEL

for j=NS+1:NS+NRX
 rate(j)= Mrx*(kf(j)*x(1,j)*x(2,j)-kb(j)*x(3,j));
end

for i=1:NC-1 
 for j=1:NT-1
    y(i,j) = alpha(i)*x(i,j)/(1+(alpha(1)-1)*x(1,j)+(alpha(2)-1)*x(2,j)+(alpha(3)-1)*x(3,j)+(alpha(4)-1)*(1-x(1,j)-x(2,j)-x(3,j)));
 end
end

% Vapor flows assuming constant molar flows (from steady state energy balance) 
for j=1:NS
        V(j)=VB;
end
for j=NS+1:NS+NRX   
    V(j)=V(j-1) - rate(j)*(Lambda/Heatvap);
end
for  j=NS+NRX+1:NT-1   
     V(j)=V(j-1);
end

% Liquid flows (from steady state mat. bal. on liquid phase)
for j = NS+NRX+1:NT
      L(j) = LT;
end
      L(NS+NRX) = L(NS+NRX+1) + F(2) + rate(NS+NRX)*(Lambda/Heatvap-1);
for j = NS+NRX-1:-1:NS+2
      L(j) = L(j+1) + rate(j)*(Lambda/Heatvap-1);
end
      L(NS+1) = L(NS+2) + F(1) + rate(NS+1)*(Lambda/Heatvap-1);
for j = 2:NS
      L(j) = L(NS+1);
end

%%%%% Time derivatives %%%%%

for i =1:NC-1
% Reboiler 
B = L(2)    - VB;
dxdt(i,1)= (L(2)*x(i,2) - V(1)*y(i,1) - B*x(i,1))/MB0;

% stripping tray 
for j=2:NS
dxdt(i,j)= (L(j+1)*x(i,j+1) - L(j)*x(i,j) + V(j-1)*y(i,j-1) - V(j)*y(i,j))/M;
end

% lower feed tray 
dxdt(i,NS+1) = (L(NS+2)*x(i,NS+2) - L(NS+1)*x(i,NS+1)  + V(NS)*y(i,NS) ...
   - V(NS+1)*y(i,NS+1) + F(1)*zF(1,i) + v(i)*rate(NS+1))/Mrx;

% reaction trays  
for j= NS+2:NS+NRX-1
dxdt(i,j) = (L(j+1)*x(i,j+1) - L(j)*x(i,j)  + V(j-1)*y(i,j-1)   - V(j)*y(i,j) + v(i)*rate(j))/Mrx;
end

% upper feed tray 
dxdt(i,NS+NRX) = (L(NS+NRX+1)*x(i,NS+NRX+1) - L(NS+NRX)*x(i,NS+NRX)  + V(NS+NRX-1)*y(i,NS+NRX-1) ...
  - V(NS+NRX)*y(i,NS+NRX) + F(2)*zF(2,i) + v(i)*rate(NS+NRX))/Mrx;

% rectifying trays
for j=NS+NRX+1:NT-1
dxdt(i,j) = (L(j+1)*x(i,j+1) - L(j)*x(i,j)  + V(j-1)*y(i,j-1)  - V(j)*y(i,j))/M;
end

% Total condenser 
D = V(NT-1) - LT;
dxdt(i,NT)= (V(NT-1)*y(i,NT-1) - LT*x(i,NT) - D*x(i,NT))/MD0;
end
for j=1:NT
dTdt(j)=-(Bvp(1)*(alpha(1)*dxdt(1,j) + alpha(2)*dxdt(2,j) + alpha(3)*dxdt(3,j) - alpha(4)*(dxdt(1,j)+dxdt(2,j)+dxdt(3,j))))/(((Avp(1) - log((P*alpha(1))/(alpha(1)*x(1,j) + alpha(2)*x(2,j) + alpha(3)*x(3,j) + alpha(4)*(1-x(1,j)-x(2,j)-x(3,j)))))^2)*(alpha(1)*x(1,j) + alpha(2)*x(2,j) + alpha(3)*x(3,j) + alpha(4)*(1-x(1,j)-x(2,j)-x(3,j))));
end
xprime=[dxdt(1,:)';dxdt(2,:)';dxdt(3,:)';dTdt'];

[hilmeijer]

Studying Large and Complex Systems with (CL)Matcont#

I will first give some background to explain why the approach is complicated. You may directly skip to "implementing a new model" and the example below.

GENERAL STRUCTURE
The general structure of Matcont is as follows. The core is the CONTINUER. It needs some settings such as MaxStepsize, or the active parameters. Starting the GUI or the call "opt=conset" on the command line for the CL(Command Line)-version will set these variables. These can then be changed manually.
Next the CONTINUER needs a defining system. This defining system can involve equilibria, periodic and connecting orbits. In the defining system, AND its derivatives, one needs the right-hand side of the model. The continuer needs the derivatives for the continuation. IF they are not available, the CONTINUER uses finite differences for the Newton corrections.

IMPLEMENTING A NEW MODEL
What is needed now, are the model equations. There are two ways to do this:
1. Using the GUI. Following most tutorials, one will be able to specify a new system. Even if there are many variables, but that may be cumbersome to type. You can specify auxilary variables, but you cannot use other constructs such as for-loops or fsolve. Auxilary variables must all be specified before
NOTE that the GUI does not support indexed variables. Typically, this is a problem for specifying large systems.
2. By hand for the CL-version: The crucial point is to implement the general structure of the m-files. The easiest is to take an existing m-file from the Systems folder, copy and rename it, so that it has the name of your model. Then open the new m-file. The required parts are the functions "init" and "func". In principle all other others can be left empty.Then in the func ffunction you can specify your equations. In the call, function dxdt = func(t,kmrgd,par1,par2.....), make sure to specify all your parameters. You may change the name for the state variables "kmrgd" to, for instance, "X ". This contains all state variables, so you may access them via X(1), or X(5).
NOTE: you can only study the system now using the CL-version as some information required for the GUI is not specified.

EXAMPLE: Consider the Bratu-Gelfand PDE $u_t=u{xx}+\lambda exp(u)$, with boundary conditions $u(0)=u(1)=0$, and parameter $\lamdba$. More details can be found in LAB2.pdf in the documentation folder. The discretization leads to the following m-file:

function dxdt = func(t,u,lambda) 
%Auxilary variables and possibly complicated formulas
N=20;
h=1/(N+1);
u0=0;u1=0;
% Pre-allocation and fill-in
dxdt=nan(20,1);
dxdt(1)=(u0+u2-2*u(1))/h/h+lambda*exp(u(1));
for i=1:19
  dxdt(i)=(u(ii-1)+u(ii+1)-2*u(ii))/h/h+lambda*exp(u(ii));
end
dxdt(N)=(u(N-1)+u1-2*u(N))/h/h+lambda*exp(u(N));

TIPS:
1. If you are able to define the derivatives by hand/symbolically, It is possible to specify the jacobian.
2. If you have the parallel toolbox, and multiple cores, then change "for" into "parfor" on line XXX in cjac.m. This will speed up the computation of Jacobians.
3. There is a package CLMATCONT_L for large systems, see the paper: Numerical computation of bifurcations in large equilibrium systems in matlab (DOI: doi:10.1016/j.cam.2013.10.034). It supports continuation of equilibria and their bifurcations. Note that this is not part of the general MATCONT package. So ask support from the authors of that paper.
4. If continuation fails, the problem may be due to redundancy. That is, there may be a linear dependence between the variables.