#*-------------------------------------------------------------------
* EMSO Model Library (EML) Copyright (C) 2004 - 2007 ALSOC.
*
* This LIBRARY is free software; you can distribute it and/or modify
* it under the therms of the ALSOC FREE LICENSE as available at
* http://www.enq.ufrgs.br/alsoc.
*
* EMSO Copyright (C) 2004 - 2007 ALSOC, original code
* from http://www.rps.eng.br Copyright (C) 2002-2004.
* All rights reserved.
*
* EMSO is distributed under the therms of the ALSOC LICENSE as
* available at http://www.enq.ufrgs.br/alsoc.
*
*--------------------------------------------------------------------
* Sample file for column model
*--------------------------------------------------------------------
* 
* This sample file needs VRTherm DEMO(www.vrtech.com.br) to run
* SectionColumn_Test and needs VRTherm full to run the distillation
* column flowsheet.
*
*----------------------------------------------------------------------
* Author: Paula B. Staudt
* $Id: sample_column.mso 500 2008-04-14 18:57:48Z paula $
*--------------------------------------------------------------------*#

using "stage_separators/column";
using "controllers/PIDIncr";

# column section with 2 trays
FlowSheet SectionColumn_Test_with2tray
	PARAMETERS
	PP	as Plugin(Brief="Physical Properties",
		Type="PP",
		Components = [ "isobutane", "benzene"],
		LiquidModel = "PR",
		VapourModel = "PR"
	);
	NComp	as Integer;

	SET
	NComp = PP.NumberOfComponents;
	
	DEVICES
	sec as Section_Column;
	feed as liquid_stream;
	reb as vapour_stream;
	cond as liquid_stream;
	zero as stream;
	
	CONNECTIONS
	feed to sec.trays(2).Inlet;
	zero to sec.trays(1).Inlet;
	reb to sec.trays(2).InletV;
	cond to sec.trays(1).InletL;
	
	SPECIFY
	feed.F = 113.4 * 'kmol/h';
	feed.T = 291 * 'K';
	feed.P = 168.3 * 'kPa';
	feed.z = [0.5, 0.5];

	zero.F = 0 * 'kmol/h';
	zero.T = 300 * 'K';
	zero.P = 1 * 'atm';
	zero.z = [0.5, 0.5];
	zero.v = 0;
	zero.h = 0 * 'J/mol';
	
	cond.F = 68 * 'kmol/h';
	cond.P = 150 * 'kPa';
	cond.T = 281.75 * 'K';
	cond.z = [0.6664, 0.3336];

	reb.P = 185 * 'kPa';
	reb.T = 328.12 * 'K';
	reb.z = [0.001848, 0.9982];
	
	sec.trays.Emv = 1;
	sec.trays(1).OutletV.F = 150 * 'kmol/h';

	SET
	sec.NTrays = 2;
	#COLUMN
	sec.trays.V = 4 * 'ft^3';
	sec.trays.Ah = 0.394 * 'ft^2';
	sec.trays.lw = 20.94 * 'in';
	sec.trays.hw = 0.125 * 'ft';
	sec.trays.Q = 0 * 'kW';
	sec.trays.beta = 0.6;
	sec.trays.alfa = 4;
	sec.trays.Ap = 3.94 * 'ft^2';

	INITIAL
	sec.trays.OutletL.T = 290 *'K';
	sec.trays.Level = 0.9 * sec.trays.hw;
	sec.trays.OutletL.z(1) = 0.5;
	
	OPTIONS
	TimeStep = 10;
	TimeEnd = 1000;
	
# After running few seconds of transient the steady-state
# can be obtained by using the results from that transient:

	#GuessFile="SectionColumn_Test_with2tray.rlt";
	#Dynamic = false;
end

# column section with 8 trays
FlowSheet SectionColumn_Test_with8tray
	PARAMETERS
	PP	as Plugin(Brief="Physical Properties",
		Type="PP",
		Components = [ "isobutane", "benzene"],
		LiquidModel = "PR",
		VapourModel = "PR"
	);
	NComp	as Integer;

	SET
	NComp = PP.NumberOfComponents;
	
	DEVICES
	sec as Section_Column;
	feed as liquid_stream;
	reb as vapour_stream;
	cond as liquid_stream;
	zero as stream;
	
	CONNECTIONS
	feed to sec.trays(5).Inlet;
	
	zero to sec.trays([1:4]).Inlet;
	zero to sec.trays([6:8]).Inlet;
	
	reb to sec.trays(8).InletV;
	cond to sec.trays(1).InletL;
	
	SPECIFY
	feed.F = 113.4 * 'kmol/h';
	feed.T = 291 * 'K';
	feed.P = 168.3 * 'kPa';
	feed.z = [0.5, 0.5];

	zero.F = 0 * 'kmol/h';
	zero.T = 300 * 'K';
	zero.P = 1 * 'atm';
	zero.z = [0.5, 0.5];
	zero.v = 0;
	zero.h = 0 * 'J/mol';
	
	cond.F = 68 * 'kmol/h';
	cond.P = 150 * 'kPa';
	cond.T = 281.75 * 'K';
	cond.z = [0.6664, 0.3336];

#	reb.F = 153 * 'kmol/h';
	reb.P = 185 * 'kPa';
	reb.T = 328.12 * 'K';
	reb.z = [0.001848, 0.9982];
	
	sec.trays.Emv = 1;
	sec.trays(1).OutletV.F = 150 * 'kmol/h';

	SET
	sec.NTrays = 8;
	#COLUMN
	sec.trays.V = 4 * 'ft^3';
	sec.trays.Ah = 0.394 * 'ft^2';
	sec.trays.lw = 20.94 * 'in';
	sec.trays.hw = 0.125 * 'ft';
	sec.trays.Q = 0 * 'kW';
	sec.trays.beta = 0.6;
	sec.trays.alfa = 4;
	sec.trays.Ap = 3.94 * 'ft^2';

	INITIAL
	sec.trays.OutletL.T = [290:(330-290)/(sec.NTrays-1):330] *'K';
	sec.trays.Level = 0.3 * sec.trays.hw;
	sec.trays.OutletL.z(1) = 0.5;
	
	OPTIONS
	TimeStep = 1;
	TimeEnd = 100;
	#GuessFile="SectionColumn_Test_with8tray.rlt";
	#Dynamic = false;
end


FlowSheet Distillation_kettle_cond_Test
	PARAMETERS
	PP	as Plugin(Brief="Physical Properties",
		Type="PP",
		Components = [ "isobutane", "n-pentane", "propylene",
		"benzene", "isobutene" ],
		LiquidModel = "PR",
		VapourModel = "PR"
	);
	NComp	as Integer;
	
	VARIABLES
	Qc as energy_source (Brief="Heat rate removed from condenser");
	Qr as energy_source (Brief="Heat rate supplied to reboiler");
	
	SET
	NComp = PP.NumberOfComponents;

	DEVICES
	col as Distillation_kettle_cond;
	feed as source;
	zero as stream;
	
	CONNECTIONS
	feed.Outlet to col.trays(5).Inlet;
	zero to col.reb.Inlet;
	zero to col.trays([1:4]).Inlet;
	zero to col.trays([6:col.NTrays]).Inlet;
	Qc.OutletQ to col.cond.InletQ;
	Qr.OutletQ to col.reb.InletQ;
	
	SPECIFY
	feed.Outlet.F = 113.4 * 'kmol/h';
	feed.Outlet.T = 291 * 'K';
	feed.Outlet.P = 168.3 * 'kPa';
	feed.Outlet.z = 1/NComp;
	
	zero.F = 0 * 'kmol/h';
	zero.T = 300 * 'K';
	zero.P = 1 * 'atm';
	zero.z = 1/NComp;
	zero.v = 0;
	zero.h = 0 * 'J/mol';
	
	col.sptop.Outlet2.F = 85 * 'kmol/h';
	col.reb.OutletL.F = 28.4 * 'kmol/h';
	col.sptop.frac = 0.444445;
	col.cond.OutletV.F = 0 * 'kmol/h';
	Qr.OutletQ.Q = 3.7743e6 * 'kJ/h';
	Qc.OutletQ.Q = -3.71e6 * 'kJ/h';
	col.pump1.dP = 16 * 'kPa';
	col.trays.Emv = 1;
	
	col.alfaTopo = 2;
	
	SET
	col.NTrays = 8;
	col.cond.V = 2 * 'm^3';
	col.cond.Across = 1 * 'm^2';
	col.trays.V = 4 * 'ft^3';
	col.trays.Ah = 0.394 * 'ft^2';
	col.trays.lw = 20.94 * 'in';
	col.trays.hw = 0.125 * 'ft';
	col.trays.Q = 0 * 'kW';
	col.trays.beta = 0.6;
	col.trays.alfa = 4;
	col.trays.Ap = 3.94 * 'ft^2';
	col.reb.V = 2 * 'm^3';
	col.reb.Across = 1 * 'm^2';
	
	INITIAL
	# condenser
	col.cond.OutletL.T = 260 *'K';
	col.cond.Level = 1 * 'm';
	col.cond.OutletL.z([1:4]) = [0.65, 0.05, 0.01, 0.01];

	# reboiler
	col.reb.OutletL.T = 330 *'K';
	col.reb.Level = 1 * 'm';
	col.reb.OutletL.z([1:4]) = [0.1, 0.7, 0.01, 0.01];

	# column trays
	col.trays.OutletL.T = [290:(330-290)/(col.NTrays-1):330] * 'K';
	col.trays.Level = 1.2 * col.trays.hw;
	col.trays.OutletL.z([1:4]) = [0.5, 0.05, 0.01, 0.01];

	OPTIONS
	TimeStep = 0.1;
	TimeEnd = 50;
	#GuessFile="Distillation_kettle_cond_Test.rlt";
	#Dynamic = false;	
end


FlowSheet Column_ctrl
	PARAMETERS
	PP	as Plugin(Brief="Physical Properties",
		Type="PP",
		Components = [ "isobutane", "n-pentane", "propylene",
			"benzene", "isobutene" ],
		LiquidModel = "PR",
		VapourModel = "PR"
	);
	NComp	as Integer;
	
	Qcmin as heat_rate (Brief="Minimum Condenser Heat supplied");
	Qcmax as heat_rate (Brief="Maximum Condenser Heat supplied");
	Qrmin as heat_rate (Brief="Minimum Reboiler Heat supplied");
	Qrmax as heat_rate (Brief="Maximum Reboiler Heat supplied");
	Frmin as flow_mol (Brief="Minimum bottom flow rate");
	Frmax as flow_mol (Brief="Maximum bottom flow rate");
	Fcmin as flow_mol (Brief="Minimum reflux flow rate");
	Fcmax as flow_mol (Brief="Maximum reflux flow rate");
	Hmint as length (Brief="Minimum liquid level in top tank");
	Hmaxt as length (Brief="Maximum liquid level in top tank");
    Hminb as length (Brief="Minimum liquid level in reboiler");
	Hmaxb as length (Brief="Maximum liquid level in reboiler");
	Pmax as pressure (Brief="Maximum column pressure");
	Pmin as pressure (Brief="Minimum column pressure");
	Tmax as temperature (Brief="Maximum column temperature");
	Tmin as temperature (Brief="Minimum column temperature");

	VARIABLES
	Qc as energy_source (Brief="Heat rate removed from condenser");
	Qr as energy_source (Brief="Heat rate supplied to reboiler");
	Had_top as Real (Brief="Dimensionless condenser level");
	Had_bot as Real (Brief="Dimensionless reboiler level");
	Pad as Real (Brief="Dimensionless pressure");
	Tad as Real (Brief="Dimensionless temperature");
	RR	as positive (Brief="Reflux ratio");

	SET
	NComp = PP.NumberOfComponents;

	DEVICES
	col as Distillation_kettle_cond;
	feed as source;
	zero as stream;
	TCcond as PIDIncr;
	LCtop as PIDIncr;
	LCbot as PIDIncr;
	PC as PIDIncr;

	CONNECTIONS
	feed.Outlet to col.trays(5).Inlet;
	zero to col.reb.Inlet;
	zero to col.trays([1:4]).Inlet;
	zero to col.trays([6:col.NTrays]).Inlet;
	Qc.OutletQ to col.cond.InletQ;
	Qr.OutletQ to col.reb.InletQ;

	EQUATIONS
   "Temperature Controller"
	TCcond.Parameters.tau = 0*'s';	
	TCcond.Parameters.tauSet = 0*'s';	
	TCcond.Parameters.alpha = 0.3;
	TCcond.Parameters.bias = 0.5;	
	TCcond.Parameters.gamma = 1;
	TCcond.Parameters.beta = 1;
	TCcond.Options.action = 1;
	TCcond.Options.clip = 1;
	TCcond.Options.autoMan = 0;
	TCcond.Parameters.intTime = 60*'s';
	TCcond.Parameters.gain = 0.6;
	TCcond.Parameters.derivTime = 1*'s';
	TCcond.Ports.setPoint = ((15+273.15) * 'K' - Tmin)/(Tmax-Tmin);
	TCcond.Ports.input = Tad;
	Tad = (col.cond.OutletL.T-Tmin)/(Tmax-Tmin);
	Qc.OutletQ.Q = Qcmin+(Qcmax-Qcmin)*TCcond.Ports.output;	

	"Pressure Controller"
	PC.Parameters.tau = 0*'s';	
	PC.Parameters.tauSet = 0*'s';	
	PC.Parameters.alpha = 0.3;
	PC.Parameters.bias = 0;	
	PC.Parameters.gamma = 1;
	PC.Parameters.beta = 1;
	PC.Options.action = -1;
	PC.Options.clip = 1;
	PC.Options.autoMan = 0;
	PC.Parameters.intTime = 50*'s';
	PC.Parameters.gain = 0.5;
	PC.Parameters.derivTime = 1*'s';
	PC.Ports.setPoint = (4.0*'bar'-Pmin)/(Pmax-Pmin);
	PC.Ports.input = Pad;
	Pad = (col.cond.OutletV.P-Pmin)/(Pmax-Pmin);
	col.cond.OutletV.F = (Fcmin+(Fcmax-Fcmin)*PC.Ports.output);	
	
	"Ttop Level Controller"
	LCtop.Parameters.tau = 0*'s';	
	LCtop.Parameters.tauSet = 0*'s';	
	LCtop.Parameters.alpha = 0.3;
	LCtop.Parameters.bias = 0.5;	
	LCtop.Parameters.gamma = 1;
	LCtop.Parameters.beta = 1;
	LCtop.Options.action = -1;
	LCtop.Options.clip = 1;
	LCtop.Options.autoMan = 0;
	LCtop.Parameters.intTime = 10*'s';
	LCtop.Parameters.gain = 1;
	LCtop.Parameters.derivTime = 0*'s';
	LCtop.Ports.setPoint = (1.0 * 'm' - Hmint)/(Hmaxt-Hmint);
	LCtop.Ports.input = Had_top;
	Had_top = (col.cond.Level-Hmint)/(Hmaxt-Hmint);
	col.sptop.Outlet1.F = Fcmin + (Fcmax-Fcmin) * LCtop.Ports.output;

	"Tbottom Level Controller"
	LCbot.Parameters.tau = 0*'s';	
	LCbot.Parameters.tauSet = 0*'s';	
	LCbot.Parameters.alpha = 0.3;
	LCbot.Parameters.bias = 0.5;	
	LCbot.Parameters.gamma = 1;
	LCbot.Parameters.beta = 1;
	LCbot.Options.action = -1;
	LCbot.Options.clip = 1;
	LCbot.Options.autoMan = 0;
	LCbot.Parameters.intTime = 100*'s';
	LCbot.Parameters.gain = 1;
	LCbot.Parameters.derivTime = 0*'s';
	LCbot.Ports.setPoint = (1.0 * 'm' - Hminb)/(Hmaxb-Hminb);
	LCbot.Ports.input = Had_bot;
	Had_bot = (col.reb.Level-Hminb)/(Hmaxb-Hminb);
	col.reb.OutletL.F = Frmin + (Frmax-Frmin) * LCbot.Ports.output;

	RR * (col.cond.OutletV.F + col.sptop.Outlet1.F) = col.sptop.Outlet2.F;
	
	if time < 1 * 'h' then
		col.sptop.Outlet2.F = 75 * 'kmol/h'; # reflux
	else
		col.sptop.Outlet2.F = 85 * 'kmol/h'; # reflux
	end

	SPECIFY
	feed.Outlet.F = 113.4 * 'kmol/h';
	feed.Outlet.T = 291 * 'K';
	feed.Outlet.P = 5 * 'bar';
	feed.Outlet.z = 1/NComp;
	
	zero.F = 0 * 'kmol/h';
	zero.T = 300 * 'K';
	zero.P = 1 * 'atm';
	zero.z = 1/NComp;
	zero.v = 0;
	zero.h = 0 * 'J/mol';
	
	Qr.OutletQ.Q = 4e6 * 'kJ/h';
	col.pump1.dP = 16 * 'kPa';
	col.trays.Emv = 1;
	col.alfaTopo = 2;

	SET
	col.NTrays = 8;
	col.cond.V = 2 * 'm^3';
	col.cond.Across = 1 * 'm^2';
	col.trays.V = 4 * 'ft^3';
	col.trays.Ah = 0.394 * 'ft^2';
	col.trays.lw = 20.94 * 'in';
	col.trays.hw = 0.125 * 'ft';
	col.trays.Q = 0 * 'kW';
	col.trays.beta = 0.6;
	col.trays.alfa = 4;
	col.trays.Ap = 3.94 * 'ft^2';
	col.reb.V = 2 * 'm^3';
	col.reb.Across = 1 * 'm^2';

	# Controllers type
	TCcond.PID_Select = "Ideal_AW";
	PC.PID_Select = "Ideal_AW";
	LCtop.PID_Select = "Ideal_AW";
	LCbot.PID_Select = "Ideal_AW";

	Qrmax = 5e6 * 'kJ/h';
	Qrmin = 1e6 * 'kJ/h';
	Frmin = 0 * 'kmol/h';
	Frmax = 60 * 'kmol/h';
	Fcmin = 0 * 'kmol/h';
	Fcmax = 120 * 'kmol/h';
	Hmint = 0 * 'm';
	Hmaxt = 2 * 'm';
	Hminb = 0 * 'm';
	Hmaxb = 2 * 'm';
	Pmin = 0.5 * 'bar';
	Pmax = 6 * 'bar';
	Qcmax = -5e5 * 'kJ/h';
	Qcmin = -5e6 * 'kJ/h';
	Tmax = (30+273.15) * 'K';
	Tmin = (-20+273.15) * 'K';
	
	INITIAL
	# condenser
	col.cond.OutletL.T = 260 *'K';
	col.cond.Level = 1 * 'm';
	col.cond.OutletL.z([1:4]) = [0.2, 0.2, 0.4, 0.05];

	# reboiler
	col.reb.OutletL.T = 350 *'K';
	col.reb.Level = 1 * 'm';
	col.reb.OutletL.z([1:4]) = [0.1, 0.4, 0.1, 0.3];

	# column trays
	col.trays.OutletL.T = [290:(330-290)/(col.NTrays-1):330] * 'K';
	col.trays.Level = 1.2 * col.trays.hw;
	col.trays.OutletL.z([1:4]) = [0.15, 0.3, 0.25, 0.2];

	OPTIONS
	TimeStep = 0.1;
	TimeEnd = 5;
	TimeUnit = 'h';
	InitialFile = "Column_ctrl.rlt";
	#GuessFile = "Column_ctrl.rlt";
	#Dynamic = false;
end

# packed column section with 8 trays
FlowSheet PackedSectionColumn
	PARAMETERS
	PP	as Plugin(Brief="Physical Properties",
		Type="PP",
		Components = [ "isobutane", "n-pentane", "propylene",
		"benzene", "isobutene" ],
		LiquidModel = "PR",
		VapourModel = "PR"
	);
	NComp	as Integer;

	SET
	NComp = PP.NumberOfComponents;
	
	DEVICES
	sec as Packed_Section_Column;
	feed as liquid_stream;
	reb as vapour_stream;
	cond as liquid_stream;
	zero as stream;
	
	CONNECTIONS
	feed to sec.stage(5).Inlet;
	zero to sec.stage([1:4]).Inlet;
	zero to sec.stage([6:sec.NStages]).Inlet;
	reb to sec.stage(8).InletV;
	cond to sec.stage(1).InletL;
	
	SPECIFY
	feed.F = 113.4 * 'kmol/h';
	feed.T = 291 * 'K';
	feed.P = 1.66 * 'atm';
	feed.z = [0.2, 0.2, 0.2, 0.2, 0.2];

	zero.F = 0 * 'kmol/h';
	zero.T = 300 * 'K';
	zero.P = 1 * 'atm';
	zero.z = [0.2, 0.2, 0.2, 0.2, 0.2];
	zero.v = 0;
	zero.h = 0 * 'J/mol';
	
	cond.F = 85 * 'kmol/h';
	cond.P = 2.33 * 'atm';
	cond.T = 283.5 * 'K';
	cond.z = [0.599, 0.044, 0.035, 0.007, 0.315];

	reb.F = 137.57 * 'kmol/h';
	reb.P = 2.46 * 'atm';
	reb.T = 325 * 'K';
	reb.z = [0.16, 0.542, 0.013, 0.008, 0.277];

	sec.dP = 0.008 * 'atm';

	SET
	sec.H = 4 * 'm';
	sec.NStages = 8;
	sec.stage.Q = 0 * 'kW';
	sec.stage.d = 1.009 * 'm';
	sec.stage.Cpo = 0.763;
	sec.stage.e = 0.951;
	sec.stage.a = 112.6 * 'm^2/m^3';
	sec.stage.Qsil = 0.1;

	INITIAL
	sec.stage.OutletL.T =[283:(325-283)/(sec.NStages-1):325] *'K';
	sec.stage.ML = 0.01 * 'kmol';
	sec.stage.OutletL.z([1:4]) = [0.3, 0.2, 0.002, 0.1];

	OPTIONS
	DAESolver(File="dassl");
	NLASolver(File="nlasolver");
	TimeStep = 10;
	TimeEnd = 100;
end

FlowSheet Packed_kettle_cond_Test
	PARAMETERS
	PP	as Plugin(Brief="Physical Properties",
		Type="PP",
		Components = [ "isobutane", "n-pentane", "propylene",
		"benzene", "isobutene" ],
		LiquidModel = "PR",
		VapourModel = "PR"
	);
	NComp	as Integer;
	
	VARIABLES
	Qc as energy_source (Brief="Heat rate removed from condenser");
	Qr as energy_source (Brief="Heat rate supplied to reboiler");
	
	SET
	NComp = PP.NumberOfComponents;

	DEVICES
	col as PackedDistillation_kettle_cond;
	feed as source;
	zero as stream;
	
	CONNECTIONS
	feed.Outlet to col.stage(5).Inlet;
	zero to col.reb.Inlet;
	zero to col.stage([1:4]).Inlet;
	zero to col.stage([6:col.NStages]).Inlet;
	Qc.OutletQ to col.cond.InletQ;
	Qr.OutletQ to col.reb.InletQ;
	
	SPECIFY
	feed.Outlet.F = 113.4 * 'kmol/h';
	feed.Outlet.T = 291 * 'K';
	feed.Outlet.P = 168.3 * 'kPa';
	feed.Outlet.z = 1/NComp;
	
	zero.F = 0 * 'kmol/h';
	zero.T = 300 * 'K';
	zero.P = 1 * 'atm';
	zero.z = 1/NComp;
	zero.v = 0;
	zero.h = 0 * 'J/mol';
	
	col.sptop.Outlet2.F = 85 * 'kmol/h';
	col.reb.OutletL.F = 28.4 * 'kmol/h';
	col.sptop.frac = 0.444445;
	col.cond.OutletV.F = 0 * 'kmol/h';
	Qr.OutletQ.Q = 3.7743e6 * 'kJ/h';
	Qc.OutletQ.Q = -3.71e6 * 'kJ/h';
	col.pump1.dP = 16 * 'kPa';

	col.dP = 0.024 * 'atm';

	SET
	col.H = 3.5 * 'm';
	col.NStages = 8;
	col.stage.Q = 0 * 'kW';
	col.stage.d = 2.24 * 'ft';
	col.stage.Cpo = 0.763;
	col.stage.e = 0.951;
	col.stage.a = 112.6 * 'm^2/m^3';

	col.cond.V = 2 * 'm^3';
	col.cond.Across = 1 * 'm^2';
	col.reb.V = 2 * 'm^3';
	col.reb.Across = 1 * 'm^2';
	col.stage.Qsil = 2;
	
	INITIAL
	# condenser
	col.cond.OutletL.T = 260 *'K';
	col.cond.Level = 1 * 'm';
	col.cond.OutletL.z([1:4]) = [0.65, 0.05, 0.01, 0.01];

	# reboiler
	col.reb.OutletL.T = 330 *'K';
	col.reb.Level = 1 * 'm';
	col.reb.OutletL.z([1:4]) = [0.1, 0.7, 0.01, 0.01];

	# column trays
	col.stage.OutletL.T = [290:(330-290)/(col.NStages-1):330] * 'K';
	col.stage.ML = 0.1 * 'kmol';
	col.stage.OutletL.z([1:4]) = [0.15, 0.5, 0.001, 0.1];

	OPTIONS
	TimeStep = 10;
	TimeEnd = 500;
end

FlowSheet PackedColumn_ctrl
	PARAMETERS
	PP	as Plugin(Brief="Physical Properties",
		Type="PP",
		Components = [ "isobutane", "n-pentane", "propylene",
			"benzene", "isobutene" ],
		LiquidModel = "PR",
		VapourModel = "PR"
	);
	NComp	as Integer;
	
	Qcmin as heat_rate (Brief="Minimum Condenser Heat supplied");
	Qcmax as heat_rate (Brief="Maximum Condenser Heat supplied");
	Qrmin as heat_rate (Brief="Minimum Reboiler Heat supplied");
	Qrmax as heat_rate (Brief="Maximum Reboiler Heat supplied");
	Frmin as flow_mol (Brief="Minimum bottom flow rate");
	Frmax as flow_mol (Brief="Maximum bottom flow rate");
	Fcmin as flow_mol (Brief="Minimum reflux flow rate");
	Fcmax as flow_mol (Brief="Maximum reflux flow rate");
	Hmint as length (Brief="Minimum liquid level in top tank");
	Hmaxt as length (Brief="Maximum liquid level in top tank");
    Hminb as length (Brief="Minimum liquid level in reboiler");
	Hmaxb as length (Brief="Maximum liquid level in reboiler");
	Pmax as pressure (Brief="Maximum column pressure");
	Pmin as pressure (Brief="Minimum column pressure");
	Tmax as temperature (Brief="Maximum column temperature");
	Tmin as temperature (Brief="Minimum column temperature");

	VARIABLES
	Qc as energy_source (Brief="Heat rate removed from condenser");
	Qr as energy_source (Brief="Heat rate supplied to reboiler");
	Had_top as Real (Brief="Dimensionless condenser level");
	Had_bot as Real (Brief="Dimensionless reboiler level");
	Pad as Real (Brief="Dimensionless pressure");
	Tad as Real (Brief="Dimensionless temperature");
	RR	as positive (Brief="Reflux ratio");

	SET
	NComp = PP.NumberOfComponents;

	DEVICES
	col as PackedDistillation_kettle_cond;
	feed as source;
	zero as stream;
	TCcond as PIDIncr;
	LCtop as PIDIncr;
	LCbot as PIDIncr;
	PC as PIDIncr;

	CONNECTIONS
	feed.Outlet to col.stage(5).Inlet;
	zero to col.reb.Inlet;
	zero to col.stage([1:4]).Inlet;
	zero to col.stage([6:col.NStages]).Inlet;
	Qc.OutletQ to col.cond.InletQ;
	Qr.OutletQ to col.reb.InletQ;

	EQUATIONS
   "Temperature Controller"
	TCcond.Parameters.tau = 0*'s';	
	TCcond.Parameters.tauSet = 0*'s';	
	TCcond.Parameters.alpha = 0.3;
	TCcond.Parameters.bias = 0.5;	
	TCcond.Parameters.gamma = 1;
	TCcond.Parameters.beta = 1;
	TCcond.Options.action = 1;
	TCcond.Options.clip = 1;
	TCcond.Options.autoMan = 0;
	TCcond.Parameters.intTime = 60*'s';
	TCcond.Parameters.gain = 0.6;
	TCcond.Parameters.derivTime = 1*'s';
	TCcond.Ports.setPoint = ((15+273.15) * 'K' - Tmin)/(Tmax-Tmin);
	TCcond.Ports.input = Tad;
	Tad = (col.cond.OutletL.T-Tmin)/(Tmax-Tmin);
	Qc.OutletQ.Q = Qcmin+(Qcmax-Qcmin)*TCcond.Ports.output;	

	"Pressure Controller"
	PC.Parameters.tau = 0*'s';	
	PC.Parameters.tauSet = 0*'s';	
	PC.Parameters.alpha = 0.3;
	PC.Parameters.bias = 0;	
	PC.Parameters.gamma = 1;
	PC.Parameters.beta = 1;
	PC.Options.action = -1;
	PC.Options.clip = 1;
	PC.Options.autoMan = 0;
	PC.Parameters.intTime = 50*'s';
	PC.Parameters.gain = 0.5;
	PC.Parameters.derivTime = 1*'s';
	PC.Ports.setPoint = (4.0*'bar'-Pmin)/(Pmax-Pmin);
	PC.Ports.input = Pad;
	Pad = (col.cond.OutletV.P-Pmin)/(Pmax-Pmin);
	col.cond.OutletV.F = (Fcmin+(Fcmax-Fcmin)*PC.Ports.output);	
	
	"Ttop Level Controller"
	LCtop.Parameters.tau = 0*'s';	
	LCtop.Parameters.tauSet = 0*'s';	
	LCtop.Parameters.alpha = 0.3;
	LCtop.Parameters.bias = 0.5;	
	LCtop.Parameters.gamma = 1;
	LCtop.Parameters.beta = 1;
	LCtop.Options.action = -1;
	LCtop.Options.clip = 1;
	LCtop.Options.autoMan = 0;
	LCtop.Parameters.intTime = 10*'s';
	LCtop.Parameters.gain = 1;
	LCtop.Parameters.derivTime = 0*'s';
	LCtop.Ports.setPoint = (1.0 * 'm' - Hmint)/(Hmaxt-Hmint);
	LCtop.Ports.input = Had_top;
	Had_top = (col.cond.Level-Hmint)/(Hmaxt-Hmint);
	col.sptop.Outlet1.F = Fcmin + (Fcmax-Fcmin) * LCtop.Ports.output;

	"Tbottom Level Controller"
	LCbot.Parameters.tau = 0*'s';	
	LCbot.Parameters.tauSet = 0*'s';	
	LCbot.Parameters.alpha = 0.3;
	LCbot.Parameters.bias = 0.5;	
	LCbot.Parameters.gamma = 1;
	LCbot.Parameters.beta = 1;
	LCbot.Options.action = -1;
	LCbot.Options.clip = 1;
	LCbot.Options.autoMan = 0;
	LCbot.Parameters.intTime = 100*'s';
	LCbot.Parameters.gain = 1;
	LCbot.Parameters.derivTime = 0*'s';
	LCbot.Ports.setPoint = (1.0 * 'm' - Hminb)/(Hmaxb-Hminb);
	LCbot.Ports.input = Had_bot;
	Had_bot = (col.reb.Level-Hminb)/(Hmaxb-Hminb);
	col.reb.OutletL.F = Frmin + (Frmax-Frmin) * LCbot.Ports.output;

	RR * (col.cond.OutletV.F + col.sptop.Outlet1.F) = col.sptop.Outlet2.F;
	
	if time < 1 * 'h' then
		col.sptop.Outlet2.F = 75 * 'kmol/h'; # reflux
	else
		col.sptop.Outlet2.F = 85 * 'kmol/h'; # reflux
	end

	SPECIFY
	feed.Outlet.F = 113.4 * 'kmol/h';
	feed.Outlet.T = 291 * 'K';
	feed.Outlet.P = 5 * 'bar';
	feed.Outlet.z = 1/NComp;
	
	zero.F = 0 * 'kmol/h';
	zero.T = 300 * 'K';
	zero.P = 1 * 'atm';
	zero.z = 1/NComp;
	zero.v = 0;
	zero.h = 0 * 'J/mol';
	
	Qr.OutletQ.Q = 4e6 * 'kJ/h';
	col.pump1.dP = 16 * 'kPa';
	
	col.dP = 0.024  * 'atm';

	SET
	col.H = 3.5 * 'm';
	col.NStages = 8;
	col.stage.Q = 0 * 'kW';
	col.stage.d = 2.24 * 'ft';
	col.stage.Cpo = 0.763;
	col.stage.e = 0.951;
	col.stage.a = 112.6 * 'm^2/m^3';

	col.cond.V = 2 * 'm^3';
	col.cond.Across = 1 * 'm^2';
	col.reb.V = 2 * 'm^3';
	col.reb.Across = 1 * 'm^2';

	col.stage.Qsil = 2;

	# Controllers type
	TCcond.PID_Select = "Ideal_AW";
	PC.PID_Select = "Ideal_AW";
	LCtop.PID_Select = "Ideal_AW";
	LCbot.PID_Select = "Ideal_AW";

	Qrmax = 5e6 * 'kJ/h';
	Qrmin = 1e6 * 'kJ/h';
	Frmin = 0 * 'kmol/h';
	Frmax = 60 * 'kmol/h';
	Fcmin = 0 * 'kmol/h';
	Fcmax = 120 * 'kmol/h';
	Hmint = 0 * 'm';
	Hmaxt = 2 * 'm';
	Hminb = 0 * 'm';
	Hmaxb = 2 * 'm';
	Pmin = 0.5 * 'bar';
	Pmax = 6 * 'bar';
	Qcmax = -5e5 * 'kJ/h';
	Qcmin = -5e6 * 'kJ/h';
	Tmax = (30+273.15) * 'K';
	Tmin = (-20+273.15) * 'K';
	
	INITIAL
	# condenser
	col.cond.OutletL.T = 260 *'K';
	col.cond.Level = 1 * 'm';
	col.cond.OutletL.z([1:4]) = [0.2, 0.2, 0.4, 0.05];

	# reboiler
	col.reb.OutletL.T = 350 *'K';
	col.reb.Level = 1 * 'm';
	col.reb.OutletL.z([1:4]) = [0.1, 0.4, 0.1, 0.3];

	# column trays
	col.stage.OutletL.T = [290:(330-290)/(col.NStages-1):330] * 'K';
	col.stage.ML = 0.1 * 'kmol';
	col.stage.OutletL.z([1:4]) = [0.15, 0.3, 0.25, 0.2];

	OPTIONS
	TimeStep = 0.1;
	TimeEnd = 5;
	TimeUnit = 'h';
	InitialFile = "Packed_Column_ctrl.rlt";
	#GuessFile = "Packed_Column_ctrl.rlt";
	#Dynamic = false;
end