#*---------------------------------------------------------------------
* 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.
*
*----------------------------------------------------------------------
* 3. Vapor pressure data representation by polynomials and equations
*----------------------------------------------------------------------
*
*   Description:
*      This problem is part of a collection of 10 representative
*	problems in Chemical Engineering for solution by numerical methods
*	developed for Cutlip (1998).
*
*   Subject:
*       * Mathematical Methods
*		* Thermodynamics
*
*	Concepts utilized:
*		Use of polynomials, a modified Clausis-Clapeyron, and the
*	Antoine equations to model vapor pressure versus temperature data.
*
*	Numerical method:
*		* Nonlinear regression
*
*   Reference:
*       * CUTLIP et al. A collection of 10 numerical problems in 
*	chemical engineering solved by various mathematical software 
*	packages. Comp. Appl. in Eng. Education. v. 6, 169-180, 1998. 
*       * More informations and a detailed description of all problems
*	is available online in http://www.polymath-software.com/ASEE
* 
*----------------------------------------------------------------------
* Author: Rodolfo Rodrigues
* GIMSCOP/UFRGS - Group of Integration, Modeling, Simulation, 
*					Control, and Optimization of Processes
* $Id$
*--------------------------------------------------------------------*#
using "types";


#*---------------------------------------------------------------------
*	Simple polynomial equation
*--------------------------------------------------------------------*#

FlowSheet simple_polynomial
	PARAMETERS
#	N		as Integer		(Brief="Number of polynomial terms", Default=5);
#	a(N)	as Real			(Brief="Parameters (coefficients) to be determined", Default=1e-2);

	a1		as Real			(Default=100);
	a2		as Real			(Default=1);
	a3		as Real			(Default=1e-2);
	a4		as Real			(Default=1e-4);
	a5		as Real			(Default=1e-6);

	VARIABLES
	P		as pressure		(Brief="Vapor pressure", DisplayUnit='mmHg', Upper=1e4);
	T		as temperature	(Brief="Temperature"); # T(K)
	
	EQUATIONS
#	"Polynomial expression" # P(mmHg) and T(K)
#	P/'mmHg' = sum(a*(T/'K' - 273.15)^[0:N-1]);

	"Polynomial expression" # P(mmHg) and T(K)
#	P/'mmHg' = a1 + (T/'K' - 273.15)*(a2 + (T/'K' - 273.15)*(a3 +
#		(T/'K' - 273.15)*(a4 + a5*(T/'K' - 273.15))));
	P/'mmHg' = a1 + a2*(T/'K' - 273.15) + a3*(T/'K' - 273.15)^2 +
		a4*(T/'K' - 273.15)^3 + a5*(T/'K' - 273.15)^4;

#	SET
#	a = [100 1 1e-2 1e-4 1e-6];
	
	SPECIFY
	T = 300*'K';

	OPTIONS
	Dynamic = false;
end


#*---------------------------------------------------------------------
*	Clausius-Clapeyron equation
*--------------------------------------------------------------------*#

FlowSheet clausius_clapeyron
	PARAMETERS
	A		as Real			(Brief="Parameters to be determined", Default=10);
	B		as temperature	(Brief="Parameters to be determined", Default=1e3);
	
	VARIABLES
	P		as pressure		(Brief="Vapor pressure", DisplayUnit='mmHg', Upper=1e4);
	T		as temperature	(Brief="Temperature"); # T(K)
	
	EQUATIONS
	"Clausis-Clapeyron equation" # P(mmHg) and T(°C)
	P/'mmHg' = exp(A - B/T);
	
	SPECIFY
	T = 300*'K';

	OPTIONS
	Dynamic = false;
end


#*---------------------------------------------------------------------
*	Antoine equation
*--------------------------------------------------------------------*#

FlowSheet antoine
	PARAMETERS
	A		as Real			(Brief="Parameters to be determined", Default=10);
	B		as temperature	(Brief="Parameters to be determined", Default=1e3);
	C		as temperature	(Brief="Parameters to be determined", Default=20);

	VARIABLES
	P		as pressure		(Brief="Vapor pressure", DisplayUnit='mmHg', Upper=1e4);
	T		as temperature	(Brief="Temperature"); # T(K)
	
	EQUATIONS
	"Antoine equation" # P(mmHg) and T(K)
	P/'mmHg' = exp(A - B/(T + C));

	SPECIFY
	T = 300*'K';

	OPTIONS
	Dynamic = false;
end


#*---------------------------------------------------------------------
*	Parameter estimation to simple polynomial equation
*--------------------------------------------------------------------*#

Estimation est_pol as simple_polynomial
	ESTIMATE
	# PAR  	START	LOWER  UPPER	UNIT
#*	a(1)	100		0		1e3;
	a(2)	1		0		10;
	a(3)	1e-2	0		10;
	a(4)	1e-4	0		10;
	a(5)	1e-6	0		10;
*#
	a1		100		0		1e3;
	a2		1		0		10;
	a3		1e-2	0		10;
	a4		1e-4	0		10;
	a5		1e-6	0		10;

	EXPERIMENTS
	# FILE                WEIGTH
	"prob03.dat"      		1;
	# Data from Perry (1999), Table 2-8

	# Solution: a=[24.75, 1.61, 3.56e-2, 4.13e-4, 4.23e-6]

	OPTIONS
	NLPSolver(
#		File = "ipopt_emso",
		File = "complex",
		ObjTol = 1e-6);
	Dynamic = false;
end


#*---------------------------------------------------------------------
*	Parameter estimation to Clausius-Clapeyron equation
*--------------------------------------------------------------------*#

Estimation est_clau as clausius_clapeyron
	ESTIMATE
	# PAR	START	LOWER	UPPER	UNIT
	A		10		-3		100;
	B 		1e3		200		5e3		'K';

	EXPERIMENTS
	# FILE                WEIGTH
	"prob03.dat"      		1;
	# Data from Perry (1999), Table 2-8
	
	# Solution: A=8.75, B=2035.33K

	OPTIONS
	NLPSolver(
		File = "ipopt_emso",
#		File = "complex",
		ObjTol = 1e-6);
	Dynamic = false;
end


#*---------------------------------------------------------------------
*	Parameter estimation to Antoine equation
*--------------------------------------------------------------------*#

Estimation est_ant as antoine
	ESTIMATE
	# PAR	START	LOWER	UPPER	UNIT
	A		10		-3		100;
	B 		1e3		200		5e3		'K';
	C 		20		-200	200		'K';

	EXPERIMENTS
	# FILE                WEIGTH
	"prob03.dat"      		1;
	# Data from Perry (1999), Table 2-8
	
	# Solution: A=5.73, B=665.42, C=152.47-273.15

	OPTIONS
	NLPSolver(
		File = "ipopt_emso",
#		File = "complex",
		ObjTol = 1e-6);
	Dynamic = false;
end