For general information on distillation: see here .
The model described below has been developed for a binary mixture. Corresponding files for the multicomponent case are found here
Contents:
The column is "column A" studied in several papers by Skogestad and Morari, e.g. see S. Skogestad and M. Morari, ``Understanding the Dynamic Behavior of Distillation Columns'', Ind. & Eng. Chem. Research, 27, 10, 1848-1862 (1988) and the book Multivariable feedback control (Wiley, 1996) by S. Skogestad and I. Postlethwaite. The model is the same as the one given in the book of Morari and Zafiriou (1989) - see their Appendix - except that we have here also included liquid flow dynamics, which is crucial if the model is used for feedback control studies. In addition, the model has recently been used in a tutorial paper: S. Skogestad, Dynamics and control of distillation columns - A tutorial introduction., Trans IChemE (UK), Vol. 75, Part A, Sept. 1997, 539-562 (Presented at Distillation and Absorbtion 97, Maastricht, Netherlands, 8-10 Sept. 1997).
The following assumptions are used: Binary mixture; constant pressure; constant relative volatility; equlibrium on all stages; total condenser; constant molar flows; no vapor holdup; linearized liquid dynamics, but effect of vapor flow ("K2"-effect) is included. These assumptions may seem restrictive, but they capture the main effects important for dynamics and control (except for the assumption about constant pressure).
The column has 40 theoretical stages and separates a binary mixture with relative volatility of 1.5 into products of 99% purity. It is relatively easy to change the model parameters and to simulate another column.
The MATLAB column model is given in the file colamod.m. Most of results given below can be generated from the file cola_test.m , and a collection of useful MATLAB-commands for linear analysis (poles, zeros, RGA, CLDG, singular values, etc.) are found in cola_commands.m . In addition, the file cola_paper.m (in the subdirectory paper) contains the files needed to generate the results in the tutorial paper (Skogestad, 1997).
This documentation was written by Sigurd Skogestad on 26 Nov 1996 and was last updated on 30 May 1997.
Thanks to Kjetil Havre, Magne Wiig Mathisen, Elisabeth Kjensjord and Atle Christiansen for their contributions.
However, to get other configurations, to lead you the through the examples below, and to provide you with easily accessible linear models (*.mat) - the following MATLAB (cola_) and SIMULINK (colas_) files are available:
G4.mat cola_init.m cola_test.m README.cola cola_init.mat colamod.m cola.dat cola_lb.m colas.m cola.html cola_lb_F1.m colas_PI.m cola4.m cola_linearize.m colas_lin.m cola4_F1.m cola_lv.m colas_lv_nonlin.m cola4_lin.m cola_lv_F1.m colas_nonlin.m cola_G4.m cola_lv_lin.m colas_test.m cola_G4_lin.mat cola_lv_lin.mat matlab_cola.tar.gz cola_G4u_lin.mat cola_lvu_lin.mat matlab_cola.zip cola_commands.m cola_rr.m numjac.m cola_db.m cola_rr_F1.m ode15s.m cola_db_lin.mat cola_rr_lin.m odeget.m cola_dv.m cola_rr_lin.mat odeset.m cola_dv_lin.mat cola_rru_lin.mat paper/ cola_init.dat cola_simulink_readme vrga.mIn addition, subdirectory paper contains the files (see in particular cola_paper.m) needed to generate the results in the paper: S. Skogestad, Dynamics and control of distillation columns - A tutorial introduction. Presented at Distillation and Absorbtion 97, Maastricht, Netherlands, 8-10 Sept. 1997.
All the files can be transferred one by one using your browser (e.g netscape), or you can transfer all the files to your local machines by transferring the file matlab_cola.tar.gz (unix) or matlab_cola.zip (PC). Afterwards you "unpack" the files by writing:
gunzip -c matlab_cola.tar.gz | tar xvf -
It is recommended that you run MATLAB in one window and have the m-file in another window, and transfer text using the mouse.
cola_init % Generates initial steady-state profile cola_simf % Simulates an increase in feed rate of 1% (using Matlab)If this works then you are in business. But: These two files are the only ready-made scipt files, so from now on you mustgo into the files and modify by yourself (as descibed in more detail below).
Note that we do not assume constant holdup on the stages, that is, we include liquid flow dynamics. This means that it takes some time (about (NT-2)*taul) from we change the liquid in the top of the column until the liquid flow into the reboiler changes. This is good for control as it means that the initial ("high-frequency") reponse is decoupled (if we have sufficiently fast control then we can avoid some of the strong interactions that exist at steady-state between the compositions at the top and bottom of the column).
Notation:
L_i and V_i - liquid and vapor flow from stage i [kmol/min],
x_i and y_i - liquid and vapor composition of light component on stage i [mole fraction],
M_i - liquid holdup on stage i [kmol],
D and B - distillate (top) and bottoms product flowrate [kmol/min],
L=LT and V=VB - reflux flow and boilup flow [kmol/min],
F, z_F - Feed rate [kmol/min] and feed composition [mole fraction],
q_F - fraction of liquid in feed
i - stage no. (1=bottom. NF=feed stage, NT=total condenser)
alpha - relative volatility between light and heavy component
taul - time constant [min] for liquid flow dynamics on each stage
lambda - constant for effect of vapor flow on liquid flow ("K2-effect")
We will write the model such that the states are x_i (i=1,NT) and M_i (i=1,NT) - a total of 2*NT states.
Model equations. The basic equations are (for more details see colamod.m):
1. Total material balance om stage i:
3. Algebraic equations. The vapor composition y_i is related to the liquid composition x_i on the same stage through the algebraic vapor-liquid equlibrium:
The above equations apply at all stages except in the top (condenser), feed stage and bottom (reboiler).
Feed stage, i=NF (we assume the feed is mixed directly into the liquid at the feed stage):
d(M_i x_i)/dt = L_{i+1} x_{i+1} + V_{i-1} y_{i-1} - L_i x_i - V_i y_i + F z_F
Total condenser, i=NT (M_NT = M_D, L_NT=L_T)
d(M_i x_i)/dt = V_{i-1} y_{i-1} - L_i x_i - D x_i
Reboiler, i=1 (M_i = M_B, V_i = V_B = V)
d(M_i x_i)/dt = L_{i+1} x_{i+1} - V_i y_i - B x_i
This results in a distillate product with D=0.5 [kmol/min] and composition y_D = x_NT = 0.99 [mole fraction units], and a bottoms product with B = 0.5 [kmol/min] and composition x_B=x_1 = 0.01 [mole fraction units].
For more complete steady-state data see the file cola.dat.
Remark. It is easy to change the column data (no. of stages, feed composition, flows, relative volatility, holdups) in the files colamod.m and cola_lv.m so that other columns can be simulated.
We may use the file cola4.m. to run this file using a MATLAB integration routine (we found ode15s to be efficient). It can also be run using SIMULINK; see below.
To find the steady-state we simulate the column for 20000 minutes starting from an inital state where all 41 liquid compositions are 0.5, and the 41 tray holdups also are 0.5 [kmol]:
[t,x]=ode15s('cola_lv',[0 20000],0.5*ones(1,82)');
lengthx=size(x); Xinit = x(lengthx(1),:)';
The resulting steady-state values of the states (Xinit) are given in
cola_init.dat and are saved
on MATLAB format in
cola_init.mat.
As expected,
the composition on the top stage is yD=x_41=0.9900
and at the bottom stage is xB=x_1=0.0100. The steady-state holdups are all 0.5
[kmol].
F=1.0+0.01; % Feedrate increase from 1 to 1.01
After saving the file we simulate using the command
[t,x]=ode15s('cola4_F1',[0 500],Xinit);
t0 = t; M1=x(:,42); xB = x(:,1); yD = x(:,41); % Save the data for plotting
We can plot the reboiler level using
plot(t0,M1)
and we see that (use the command: axis([0 5, 0.45 0.55]) ) after an initial ``delay'' of about 20*0.063 = 1.26 min (the time the liquid takes to go through 20 ``lags'' (stages) each with time constant 0.063 min), the reboiler level increases in a ramplike fashion. The liquid composition in the reboiler, xB, increases in a similar fashion for the first 100 minutes or so, but then settles at a new steady-state value of about 0.017 (with a first-order time constant of about 200 minutes). The top composition, yD, has more of a second-order response, but otherwise the dynamics are similar.
The above simulation was with the "open-loop" model with no control loops closed. Simulations with various configurations (LV, LB, DV, L/D-V/B) where two level loops have been closed are given next.
[t,x]=ode15s('cola4_F1',[0 500],Xinit);
t0 = t; M1=x(:,42); xB = x(:,1); yD = x(:,41); % Save the data for plot
Let us now compare the response when we use the LV-configuration
with tight
control of reboiler and condenser holdup. We must go into the file
cola_lv.m and change F to
1.01 (save the new file as cola_lv_F1.m),
and then use
[t,x]=ode15s('cola_lv_F1',[0 500],Xinit);
tlv=t; M1lv=x(:,42); xBlv = x(:,1); yDlv = x(:,41); %save data for plot
Comparing the simulations with those for the open-loop model
(e.g. use plot(t0,xB,t,xBlv)) )
we see that the reboiler holdup Mb
levels of very quickly at 0.501 because of the fast level control
where D is increased from 0.5 to 0.501. However,
there is essentially no difference in yD and only a small difference in xB (the
increase in xB is slightly larger in the latter case with fast level control).
Note that with
slow level control, the LV-model becomes identical to the open-loop model.
This demonstrates that the LV-configuration is rather insensitive to the
level tuning.
Note that it is only the LV-configuration which behaves similarly to the
open-loop model, and which is insensitive to the level tuning.
For example, consider the LB-configuration with tight level control, as given in the file cola_lb.m. We do a simulation of the same feedrate change
[t,x]=ode15s('cola_lb_F1',[0 500],Xinit);
tlb=t; M1lb=x(:,42); xBlb = x(:,1); yDlb = x(:,41);
Here the vapor flow is adjusted to keep the reboiler holdup constant, so
the feed increase results in an increased vapor flow up the column, so
the increased feed actually ends up leaving the top of the column. The
result is that xB and yD decrease rather than increasing
as they did for the LV-configuration. With a slowly tuned level controller,
e.g. KcB=0.1, this results in an inverse response for the compositions
xB and yD.
Finally, consider the response with the Double ratio (L/D-V/B) configuration given in the file cola_rr.m. In this case we set the ratios L/D and V/B externally, and let the levels be controlled using D and B (controlling the levels like this may not be the best solution; but for tight level control it makes no difference).
[t,x]=ode15s('cola_rr_F1',[0 500],Xinit);
trr=t; M1rr=x(:,42); xBrr = x(:,1); yDrr = x(:,41);
In this case, the system is almost ``self-regulating'' with respect
to the feed rate disturbance, because by the indirect action of the
level controllers all flows in the column are increased proportionally
to the feedrate change (which obviously is the right thing to do - at least
at steady state).
There is also a file for the DV-configuration in the file cola_dv.m, but for a feed rate disturbance it behaves as the LV-configuration (but otherwise the dynamic response and control properties are completely different).
To plot all the above results use:
plot(t0,M1,'-',tlv,M1lv,':',tlb,M1lb,'-.',trr,M1rr,'--');
title('MB: reboiler (bottom) holdup [mol]');
plot(t0,xB,'-',tlv,xBlv,':',tlb,xBlb,'-.',trr,xBrr,'--');
title('xB: reboiler (bottom) composition');
plot(t0,yD,'-',tlv,yDlv,':',tlb,yDlb,'-.',trr,yDrr,'--');
title('yD: distillate (top) composition');
Ls=2.706; Vs=3.206; Ds=0.5; Bs=0.5; Fs=1.0; zFs=0.5;
[A,B,C,D]=cola_linearize('cola4_lin',Xinit',[Ls Vs Ds Bs Fs zFs]);
G4u = pck(A,B,C,D);
The linear model has six inputs (the latter two are actually disturbances):
[LT VB D B F zF]
and four outputs
[yD xB MD MB]
To check the model we may compute its 82 eigenvalues, eig(A) , and we find that the 3 eigenvalues furthest to the right are 0, 0 and -5.1562e-03. The first two are from the integrating (uncontrolled) levels in the reboiler and the condenser, and the latter one, which corresponds to a time constant of 1/5.1562e-03 = 193.9 minutes is the dominant time constant for the composition response (it also behaves almost as an integrator over short time scales, so a ``stabilizing'' controller is also needed in practice for the compositions, e.g. by using a temperature controller).
G4u. The above linear model is called G4u - 4 denotes that there are 4 inputs (L,V,D,B) - in addition to the two disturbances - and the u denotes that the model is unscaled.
Linear simulation. To simulate a feed-rate increase we can use one of the commands from the control toolboxes, for example,
[Y,X,T]=step(A,B,C,D,5); plot(T,Y(:,1),T, Y(:,2))
which makes a unit step change, deltaF=1 [kmol/min], in input 5 (the feed rate). This increase is actually 100 times larger than in the nonlinear simulation above, but since we here use a linear model we need only reduce the output by a factor 100. Taking this into account we see that an increase in F of 0.01 to 1.01 [kmol/min] (1%) yields a steady-state increase in y_D of 0.0039 and in x_B of 0.0059.
Comparison with nonlinear simulation. Simulation with the nonlinear model gave a an increase in y_D of 0.0025 (to 0.9925) an in x_B of 0.0072 (to 0.0172). Thus, the increase in y_D is smaller with the nonlinear model, but larger in x_B. The reason is that y_D approaches higher purity (which is difficult to achieve) whereas x_B approaches lower purity (which is easy to achieve).
Remark. These nonlinear effects are to a large extent counteracted by using logarithmic compositions, X=ln(x_L/x_H) where L denotes light component and H heavy component. For the top and bottom composition of a binary mixture, we may use
X_B = ln (x_B/(1-x_B)); Y_D = ln (y_D/(1-y_D))
Make sure you get the sign for the bottom loop right if these logarithmic outputs are used for controller design.
Model reduction. The above linear model has 82 states, but using model reduction the order can easily be reduced to 10 states or less - without any noticable difference in the response. But beware - the open-loop ("uncontrolled") model is unstable as it has two integrators. The following MATLAB commands from the Mu-toolbox show how we may first decompose the model into a stable and unstable part, and then model reduce the stable part to 6 states using Hankel norm approximation (for more details see page 465-466 in Skogestad and Postlethwaite (1996)).
% First scale the model: % Scaling of the model is generally recommended before you do model reduction % The scaling should make all inputs of about equal importance, and all outputs of about equal importance. % For more details about the scaling; see below Du = diag([1 1 1 1]); % max inputs (scalings) Dd = diag([ .2 .1]); % max disturbances (scalings) De = diag([0.01 0.01 1 1]); % max output errors (scalings) Si = daug(Du,Dd); So = minv(De); % introduce scaling matrices G4 = mmult (So, G4u, Si); % scaled 82 state model % Now do model reduction [syss,sysu]=sdecomp(G4); [syssb,hsig]=sysbal(syss); sys1=hankmr(syssb,hsig,6,'d'); G4_8=madd(sys1,sysu); % scaled 8 state model G4u_8=mmult( minv(So), G4_8, minv(Si)); % "un"scaled 8 state modelThe resulting overall model G4u_8 has 8 states. Note that the model was scaled before the model reduction such that all the inputs and outputs are of equal importance (see below for the scalings used). We can compare the reponses of the 82-state and 8-state linear models to a feed rate change:
[A8,B8,C8,D8]=unpck(G4u_8); [Y8,X8]=step(A8,B8,C8,D8,5,T); plot(T,Y(:,1),T,Y8(:,1));As can be seen, there is no significant difference.
Analysis of linear model. The main advantage with a linear model is that it is suitable for analysis (RGA, RHP-zeros, CLDG, etc.) and for controller systhesis. Again, see Skogestad and Postlethwaite (1996) for more details, or see some useful commands in cola_commands.m
Linearized models for other configurations. To obtain linear models for other configurations, we can start from the ``open-loop'' linear model G4u as shown in Table 12.3 in the book by Skogestad and Postletwaite (1996). Alternatively, we can linearize the nonlinear model directly, for example,
For the LV-configuration (with 4 inputs and 2 outputs):
Ls=2.706; Vs=3.206, Fs=1.0; zFs=0.5;
[A,B,C,D]=cola_linearize('cola_lv_lin',Xinit',[Ls Vs Fs zFs]);
Glvu = pck(A,B,C,D);
For the double ratio (L/D-V/B) configuration (with 4 inputs and 2 outputs):
R1s = 2.70629/0.5; R2s=3.20629/0.5; Fs=1.0; zFs=0.5;
[A,B,C,D]=cola_linearize('cola_rr_lin',Xinit',[R1s R2s Fs zFs]);
Grru = pck(A,B,C,D);
u = U / Umax; y = Y / Ymax; d = D / Dmaxwhere U is the deviation in original units, Umax is the maximum allowed or expected deviation, and u is the scaled variable. For more details see Skogestad and Postlethwaite (1996). NOTE: "Unscaled" linear models are here denoted with "u", e.g. G4u.
To get the scaled model G4 referred to in Table 12.3 on page 491 in the book of Skogestad and Postethwaite (1996), we need to scale the above linear model G4u as follows:
% The following max. changes are used (for scaling the model):
Du = diag([1 1 1 1]); % max inputs (scalings)
Dd = diag([ .2 .1]); % max disturbances (scalings)
De = diag([0.01 0.01 1 1]); % max output errors (scalings)
% This implies the folling in terms of the scaled model G4:
% Units for inputs (L,V,D,B): 1 = 1 kmol/min = F (the feed rate)
% Comment: For simplicity we here used the same value for all inputs
% Normally we select the nominal flow as the max input, i.e. Du = diag([2.7 3.2 0.5 0.5])
% Units for disturbance 1 (F): 1 = 0.2 kmol/min (20% change)
% Units for disturbance 2 (z_f): 1 = 0.1 mole fraction units (20% change)
% Units for outputs 1 and 2 (y_d and x_b): 1 = 0.01 mole fraction units
% Units for outputs 3 and 4 (M_d and M_b): 1 = 1 kmol
% The scaled model is then G4:
Si = daug(Du,Dd); So = minv(De); % introduce scaling matrices
G4 = mmult(So, G4u, Si);
All the command needed to generate this model ``from scratch'' are
given in cola_G4.m.
See also work by Mejdell et al. for use of column temperatures to estimate product compositions.
The SIMULINK interface to colamod.m is colas.m ; the latter file simply calls the first file to get derivatives of the states for a given state X and input U. The nominal initial steady-state conditions, Xinit are loaded from the file cola_init.mat (you may easily generate Xinit by running the file cola_init.m ).
Some SIMULINK files:
load cola_init
[A,B,C,D] = linmod('colas',Xinit,Uinit);
colas_lin
The linear SIMULINK simulation is faster than with the nonlinear
model, but still
slow compared to using MATLAB (see above).
In any case you it is recommended to use dynamic simulation and run it to steady-state (e.g. using a modefified cola_init.m ) to obtain steady-state data which can be used as initial data for later simulations (saved in cola_init.mat for use with Matlab simulations or with Simulink simulations).
You then need to change the steady-state data in colamod.m (most of the changes are needed because of the linearized flow dynamics). For simulations in Matlab you also need to change the "Inputs and disturbances" for the individual configurations, such as cola4.m or cola_lv.m . For simulations with Simulink you need to change NT (no. of stages) in colas.m .
Good Luck.