function xprime=colamod(t,X,U) %% colamod - This is a nonlinear model of a distillation column with%           NT-1 theoretical stages including a reboiler (stage 1) plus a%           total condenser ("stage" NT). The liquid flow dynamics are%           modelled by a simple linear relationship.%           Model assumptions: Two components (binary separation); constant%           relative volatility; no vapor holdup; one feed and two products;%           constant molar flows (same vapor flow on all stages); %           total condenser%%           The model is based on column A in Skogestad and Postlethwaite%           (1996). The model has 82 states.%% Inputs:    t    - time in [min].%            X    - State, the first 41 states are compositions of light%                   component A with reboiler/bottom stage as X(1) and %                   condenser as X(41). State X(42)is holdup in reboiler/%                   bottom stage and X(82) is hold-up in condenser. %            U(1) - reflux L,%            U(2) - boilup V,%            U(3) - top or distillate product flow D,%            U(4) - bottom product flow B,%            U(5) - feed rate F,%            U(6) - feed composition, zF.%            U(7) - feed liquid fraction, qF.%% Outputs:   xprime - vector with time derivative of all the states %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%------------------------------------------------------------% The following data need to be changed for a new column.% These data are for "colmn A".% Number of stages (including reboiler and total condenser:     NT=41; % Location of feed stage (stages are counted from the bottom):    NF=21;% Relative volatility    alpha=1.5;% Nominal liquid holdups    M0(1)=0.5;      	% Nominal reboiler holdup (kmol)    i=2:NT-1; M0(i)=0.5*ones(1,NT-2);% Nominal stage (tray) holdups (kmol)    M0(NT)=0.5;      	% Nominal condenser holdup (kmol)% Data for linearized liquid flow dynamics (does not apply to reboiler and condenser):    taul=0.063;     	% time constant for liquid dynamics (min)    F0=1;	 	% Nominal feed rate (kmol/min)     qF0 = 1; 		% Nominal fraction of liquid in feed     L0=2.70629;     	% Nominal reflux flow (from steady-state data)    L0b=L0 + qF0*F0;	% Nominal liquid flow below feed (kmol/min)    lambda=0;		% Effect of vapor flow on liquid flow ("K2-effect")    V0=3.20629;V0t=V0+(1-qF0)*F0;% Nominal vapor flows - only needed if lambda is nonzero % End data which need to be changed%------------------------------------------------------------% Splitting the statesx=X(1:NT)';                          % Liquid composition from btm to topM=X(NT+1:2*NT)';                     % Liquid hold up from btm to top% Inputs and disturbancesLT = U(1);                            % RefluxVB = U(2);                            % BoilupD  = U(3);                            % DistillateB  = U(4);                            % BottomsF  = U(5);                            % FeedratezF = U(6);                            % Feed compositionqF = U(7);                            % Feed liquid fraction% THE MODEL% Vapor-liquid equilibriai=1:NT-1;    y(i)=alpha*x(i)./(1+(alpha-1)*x(i));% Vapor Flows assuming constant molar flowsi=1:NT-1;    V(i)=VB*ones(1,NT-1);i=NF:NT-1;   V(i)=V(i) + (1-qF)*F;% Liquid flows assuming linearized tray hydraulics with time constant taul% Also includes coefficient lambda for effect of vapor flow ("K2-effect").i=2:NF;      L(i) = L0b + (M(i)-M0(i))./taul + lambda.*(V(i-1)-V0);i=NF+1:NT-1; L(i) = L0  + (M(i)-M0(i))./taul + lambda.*(V(i-1)-V0t);L(NT)=LT;% Time derivatives from  material balances for % 1) total holdup and 2) component holdup% Columni=2:NT-1;dMdt(i) = L(i+1)         - L(i)       + V(i-1)         - V(i);dMxdt(i)= L(i+1).*x(i+1) - L(i).*x(i) + V(i-1).*y(i-1) - V(i).*y(i);% Correction for feed at the feed stage% The feed is assumed to be mixed into the feed stagedMdt(NF) = dMdt(NF)  + F;dMxdt(NF)= dMxdt(NF) + F*zF;% Reboiler (assumed to be an equilibrium stage)dMdt(1) = L(2)      - V(1)      - B;dMxdt(1)= L(2)*x(2) - V(1)*y(1) - B*x(1);% Total condenser (no equilibrium stage)dMdt(NT) = V(NT-1)         - LT       - D;dMxdt(NT)= V(NT-1)*y(NT-1) - LT*x(NT) - D*x(NT);% Compute the derivative for the mole fractions from d(Mx) = x dM + M dxi=1:NT;dxdt(i) = (dMxdt(i) - x(i).*dMdt(i) )./M(i);% Outputxprime=[dxdt';dMdt'];