#*-------------------------------------------------------------------
* 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 PID MIMO application
* (dual control of a distillation column)
*--------------------------------------------------------------------
* 
* This sample file needs VRTherm (www.vrtech.com.br) to run
* the distillation column flowsheet.
*
*----------------------------------------------------------------------
* Author: Argimiro R. Secchi
* $Id: sample_PID_MIMO.mso 313 2007-07-14 16:45:55Z arge $
*--------------------------------------------------------------------*#

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

FlowSheet Column_MPID
	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");
	RRmax as positive (Brief="Maximum reflux ratio");
	RRmin as positive (Brief="Minimum reflux ratio");
	
	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");
	P_sp as pressure (Brief="Pressure setpoint"); 
	Tad as Real (Brief="Dimensionless temperature");
	T_sp as temperature (Brief="Condenser temperature setpoint");
	RR	as positive (Brief="Reflux ratio");
	LC_sp as length (Brief="Condenser level setpoint");
	LR_sp as length (Brief="Reboiler level setpoint");
	xtop as fraction (Brief="Controlled top composition");
	xbot as fraction (Brief="Controlled bottom composition");
	xtop_sp as fraction (Brief="Controlled top composition setpoint");
	xbot_sp as fraction (Brief="Controlled bottom composition setpoint");

	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;
	PID11 as PID;
	PID22 as PID;
	PID12 as PID;
	PID21 as PID;

	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 = (T_sp - 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 = (P_sp-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 = 20*'s';
	LCtop.Parameters.gain = 1;
	LCtop.Parameters.derivTime = 0*'s';
	LCtop.Ports.setPoint = (LC_sp - 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 = (LR_sp - 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;

	"MIMO PID Controller"
	PID11.Parameters.tau = 0*'s';	
	PID11.Parameters.tauSet = 0*'s';	
	PID11.Parameters.alpha = 0.3;
	PID11.Parameters.bias = 0.1;	
	PID11.Parameters.gamma = 1;
	PID11.Parameters.beta = 1;
	PID11.Options.action = -1;
	PID11.Options.clip = 1;
	PID11.Options.autoMan = 0;
	PID11.Parameters.intTime = 60*'s';
	PID11.Parameters.gain = 0.6;
	PID11.Parameters.derivTime = 0*'s';
	PID11.Ports.setPoint = xtop_sp;
	PID11.Ports.input = xtop;

	PID12.Parameters.tau = 0*'s';	
	PID12.Parameters.tauSet = 0*'s';	
	PID12.Parameters.alpha = 0.3;
	PID12.Parameters.bias = 0.0;	
	PID12.Parameters.gamma = 1;
	PID12.Parameters.beta = 1;
	PID12.Options.action = -1;
	PID12.Options.clip = 1;
	PID12.Options.autoMan = 0;
	PID12.Parameters.intTime = 100*'s';
	PID12.Parameters.gain = 0.5;
	PID12.Parameters.derivTime = 0*'s';
	PID12.Ports.setPoint = xbot_sp;
	PID12.Ports.input = xbot;

	RR = RRmin + (RRmax-RRmin) * (PID11.Ports.output + PID12.Ports.output);

	PID21.Parameters.tau = 0*'s';	
	PID21.Parameters.tauSet = 0*'s';	
	PID21.Parameters.alpha = 0.3;
	PID21.Parameters.bias = 0.0;	
	PID21.Parameters.gamma = 1;
	PID21.Parameters.beta = 1;
	PID21.Options.action = 1;
	PID21.Options.clip = 1;
	PID21.Options.autoMan = 0;
	PID21.Parameters.intTime = 100*'s';
	PID21.Parameters.gain = 0.5;
	PID21.Parameters.derivTime = 0*'s';
	PID21.Ports.setPoint = xtop_sp;
	PID21.Ports.input = xtop;

	PID22.Parameters.tau = 0*'s';	
	PID22.Parameters.tauSet = 0*'s';	
	PID22.Parameters.alpha = 0.3;
	PID22.Parameters.bias = 0.5;	
	PID22.Parameters.gamma = 1;
	PID22.Parameters.beta = 1;
	PID22.Options.action = -1;
	PID22.Options.clip = 1;
	PID22.Options.autoMan = 0;
	PID22.Parameters.intTime = 60*'s';
	PID22.Parameters.gain = 0.6;
	PID22.Parameters.derivTime = 0*'s';
	PID22.Ports.setPoint = xbot_sp;
	PID22.Ports.input = xbot;
	
	Qr.OutletQ.Q = Qrmin + (Qrmax-Qrmin) * (PID21.Ports.output + PID22.Ports.output);

	RR * (col.cond.OutletV.F + col.sptop.Outlet1.F) = col.sptop.Outlet2.F;
	xtop = col.cond.OutletL.z(4); # benzene
	xbot = col.reb.OutletL.z(3); # propylene

	if time < 1 * 'h' then
		feed.Outlet.F = 113.4 * 'kmol/h';
	else
		feed.Outlet.F = 120 * 'kmol/h';
	end

	SPECIFY
	feed.Outlet.T = 291 * 'K';
	feed.Outlet.P = 568.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';
	
#	Qr.OutletQ.Q = 3e6 * 'kJ/h';
	col.pump1.dP = 16 * 'kPa';
	col.trays.Emv = 1;
#	col.sptop.Outlet2.F = 75 * 'kmol/h'; # reflux

	xtop_sp = 1e-3; # benzene
	xbot_sp = 1e-5; # propylene
	P_sp = 4.0*'bar'; 
	T_sp = (15+273.15) * 'K';
	LC_sp = 1.0 * 'm';
	LR_sp = 1.0 * 'm';

	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";
	PID11.PID_Select = "Ideal_AW";
	PID12.PID_Select = "Ideal_AW";
	PID21.PID_Select = "Ideal_AW";
	PID22.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 = 10 * 'bar';
	Qcmax = -5e5 * 'kJ/h';
	Qcmin = -5e6 * 'kJ/h';
	Tmax = (30+273.15) * 'K';
	Tmin = (-20+273.15) * 'K';
	RRmax = 5;
	RRmin = 0.1;

	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