#*-------------------------------------------------------------------
* 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;
	
	Qrmin as heat_rate (Brief="Minimum Reboiler Heat supplied");
	Qrmax as heat_rate (Brief="Maximum Reboiler Heat supplied");
	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");
	P_sp as pressure (Brief="Pressure setpoint"); 
	T_sp as temperature (Brief="Condenser temperature setpoint");
	RR	as positive (Brief="Reflux ratio");
	LC_sp as fraction (Brief="Condenser level setpoint");
	LR_sp as fraction (Brief="Reboiler level setpoint");
out	xtop as fraction (Brief="Controlled top composition");
out	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 Section_Column;
	cond as condenser;
	reb as reboiler;
	sptop as splitter_column;
	pump1 as pump;
#	Valve_Distillate as valve_flow;
#	Valve_Bottom as valve_flow;
	feed as source;
	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.FeedTray;
	Qc.OutletQ to cond.InletQ;
	Qr.OutletQ to reb.InletQ;
#	sptop.Distillate to Valve_Distillate.Inlet;
	sptop.Reflux to pump1.Inlet;
	pump1.Outlet to col.LiquidInlet;
	col.VapourOutlet to cond.InletVapour;
	cond.OutletLiquid to sptop.Inlet;
	col.LiquidOutlet to reb.InletLiquid;
	reb.OutletVapour to col.VapourInlet;
#	reb.OutletLiquid to Valve_Bottom.Inlet;	
	cond.TI to TCcond.Input;
	cond.PI to PC.Input;
	cond.LI to LCtop.Input;
	reb.LI to LCbot.Input;
	xtop to PID11.Input;
	xbot to PID12.Input;
	xtop to PID21.Input;
	xbot to PID22.Input;

	EQUATIONS

	TCcond.SetPoint * 'K' = T_sp;
	Qc.OutletQ = TCcond.Output * 'kJ/h';

	PC.SetPoint * 'bar' = P_sp;
	cond.OutletVapour.F = PC.Output * 'kmol/h';

	LCtop.SetPoint = LC_sp;
	sptop.Distillate.F = LCtop.Output * 'kmol/h';

	LCbot.SetPoint = LR_sp;
	reb.OutletLiquid.F = LCbot.Output * 'kmol/h';

	"MIMO PID Controller"
	PID11.SetPoint = xtop_sp;
	PID12.SetPoint = xbot_sp;

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

	PID21.SetPoint = xtop_sp;
	PID22.SetPoint = xbot_sp;

	Qr.OutletQ = Qrmin + (Qrmax-Qrmin) * (PID21.Output + PID22.Output);

	RR * (cond.OutletVapour.F + sptop.Distillate.F) = sptop.Reflux.F;
	xtop = cond.OutletLiquid.z(4); # benzene
	xbot = reb.OutletLiquid.z(3); # propylene

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

	SPECIFY
	feed.T = 291 * 'K';
	feed.P = 568.3 * 'kPa';
	feed.Composition = 1/NComp;
	
#	Qr.OutletQ = 3e6 * 'kJ/h';
	pump1.Pincrease = 16 * 'kPa';
#	sptop.Reflux.F = 75 * 'kmol/h'; # reflux
	col.VapourDrawOffFlow = 0 * 'kmol/h' ;
	col.LiquidDrawOffFlow = 0 * 'kmol/h' ;

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

	SET
	col.FeedTrayLocation(1) = 5;
	col.NumberOfTrays = 8;
	col.LiquidSideStreamLocation = 2;
	col.VapourSideStreamLocation = 2;
	col.WeirLength = 20.94 * 'in';
	col.WeirHeight = 0.125 * 'ft';
	col.TrayLiquidPasses = 1;
	col.HeatSupply = 0 * 'kW' ;
	col.AerationFraction = 0.6;
	col.DryPdropCoeff = 4;
	col.MurphreeEff = 1;

	cond.Geometry.Lenght = 1 * 'm';
	cond.Geometry.Diameter = 1 * 'm';
	reb.Geometry.Lenght = 2 * 'm';
	reb.Geometry.Diameter = 1 * 'm';

	# 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";

	# Temperature Controller
	TCcond.tau = 1*'s';	
	TCcond.tauSet = 1*'s';	
	TCcond.alpha = 0.1;
	TCcond.bias = 0.5;	
	TCcond.gamma = 1;
	TCcond.beta = 1;
	TCcond.Action = "Reverse";
	TCcond.Clip = "Clipped";
	TCcond.Mode = "Automatic";
	TCcond.intTime = 60*'s';
	TCcond.gain = 0.6;
	TCcond.derivTime = 1*'s';

	# Pressure Controller
	PC.tau = 1*'s';	
	PC.tauSet = 1*'s';	
	PC.alpha = 0.1;
	PC.bias = 0;	
	PC.gamma = 1;
	PC.beta = 1;
	PC.Action = "Direct";
	PC.Clip = "Clipped";
	PC.Mode = "Automatic";
	PC.intTime = 50*'s';
	PC.gain = 0.5;
	PC.derivTime = 1*'s';

	# Ttop Level Controller
	LCtop.tau = 1*'s';	
	LCtop.tauSet = 1*'s';	
	LCtop.alpha = 0.1;
	LCtop.bias = 0.5;	
	LCtop.gamma = 1;
	LCtop.beta = 1;
	LCtop.Action = "Direct";
	LCtop.Clip = "Clipped";
	LCtop.Mode = "Automatic";
	LCtop.intTime = 20*'s';
	LCtop.gain = 1;
	LCtop.derivTime = 0*'s';

	# Tbottom Level Controller
	LCbot.tau = 1*'s';	
	LCbot.tauSet = 1*'s';	
	LCbot.alpha = 0.1;
	LCbot.bias = 0.5;	
	LCbot.gamma = 1;
	LCbot.beta = 1;
	LCbot.Action = "Direct";
	LCbot.Clip = "Clipped";
	LCbot.Mode = "Automatic";
	LCbot.intTime = 100*'s';
	LCbot.gain = 1.5;
	LCbot.derivTime = 0*'s';

	# MIMO PID Controller
	PID11.tau = 1*'s';	
	PID11.tauSet = 1*'s';	
	PID11.alpha = 0.1;
	PID11.bias = 0.1;	
	PID11.gamma = 1;
	PID11.beta = 1;
	PID11.Action = "Direct";
	PID11.Clip = "Clipped";
	PID11.Mode = "Automatic";
	PID11.intTime = 60*'s';
	PID11.gain = 0.6;
	PID11.derivTime = 0*'s';

	PID12.tau = 1*'s';	
	PID12.tauSet = 1*'s';	
	PID12.alpha = 0.1;
	PID12.bias = 0.0;	
	PID12.gamma = 1;
	PID12.beta = 1;
	PID12.Action = "Direct";
	PID12.Clip = "Clipped";
	PID12.Mode = "Automatic";
	PID12.intTime = 100*'s';
	PID12.gain = 0.5;
	PID12.derivTime = 0*'s';
	
	PID21.tau = 1*'s';	
	PID21.tauSet = 1*'s';	
	PID21.alpha = 0.1;
	PID21.bias = 0.0;	
	PID21.gamma = 1;
	PID21.beta = 1;
	PID21.Action = "Reverse";
	PID21.Clip = "Clipped";
	PID21.Mode = "Automatic";
	PID21.intTime = 100*'s';
	PID21.gain = 0.5;
	PID21.derivTime = 0*'s';

	PID22.tau = 1*'s';	
	PID22.tauSet = 1*'s';	
	PID22.alpha = 0.1;
	PID22.bias = 0.5;	
	PID22.gamma = 1;
	PID22.beta = 1;
	PID22.Action = "Direct";
	PID22.Clip = "Clipped";
	PID22.Mode = "Automatic";
	PID22.intTime = 60*'s';
	PID22.gain = 0.6;
	PID22.derivTime = 0*'s';

	PC.MinOutput = 0;
	PC.MaxOutput = 120;
	PC.MinInput = 0.5;
	PC.MaxInput = 10;
	LCtop.MinInput = 0;
	LCtop.MaxInput = 1;
	LCtop.MinOutput = 0;
	LCtop.MaxOutput = 120;
	LCbot.MinInput = 0;
	LCbot.MaxInput = 1;
	LCbot.MinOutput = 0;
	LCbot.MaxOutput = 60;
	TCcond.MaxOutput = -5e5;
	TCcond.MinOutput = -5e6;
	TCcond.MaxInput = (30+273.15);
	TCcond.MinInput = (-20+273.15);
	PID11.MinInput = 0;
	PID12.MinInput = 0;
	PID21.MinInput = 0;
	PID22.MinInput = 0;
	PID11.MaxInput = 1;
	PID12.MaxInput = 1;
	PID21.MaxInput = 1;
	PID22.MaxInput = 1;
	RRmax = 5;
	RRmin = 0.1;
	Qrmax = 5e6 * 'kJ/h';
	Qrmin = 1e6 * 'kJ/h';

	# Initial condition
	# condenser
	cond.Levelpercent_Initial = 0.5;
	cond.Initial_Temperature = 260 *'K';
	cond.Initial_Composition = [0.2, 0.2, 0.4, 0.05, 0.15];

	# reboiler
	reb.Levelpercent_Initial = 0.5;
	reb.Initial_Temperature = 350 *'K';
	reb.Initial_Composition = [0.1, 0.4, 0.1, 0.3, 0.1];

	# column trays
	col.INITIALIZATION.TopTemperature = 290 * 'K';
	col.INITIALIZATION.BottomTemperature = 330 * 'K';
	col.INITIALIZATION.LevelFraction = 0.1;
	col.INITIALIZATION.TopComposition = [0.2, 0.2, 0.4, 0.05, 0.15];
	col.INITIALIZATION.BottomComposition = [0.1, 0.4, 0.1, 0.3, 0.1];

	OPTIONS
	TimeStep = 0.1;
	TimeEnd = 5;
	TimeUnit = 'h';
	InitialFile = "Column_MPID.rlt";
	DAESolver(File="dassl");
	#GuessFile = "Column_MPID.rlt";
	#Dynamic = false;
end