USE MATH
--USE ESNOPT

COMPONENT estim_param_deposito (INTEGER nent =2, INTEGER nsal = 2)

-------------------------------------------------------------------------------
--
-- El sistema es un depósito al que llega un flujo de entrada (qe) y
-- del que sale un flujo (F) por gravedad a traves de una valvula con apertura (a). 
-- El flujo de salida (F) es proporcional a la raíz cuadrada del nivel en el depósito
-- y a la apertura de la valvula, siendo la constante de proporcionalidad (k) 
-- el factor de fricción. Además el depósito dispone de una resistencia eléctrica (R) calefactora.
--
--
--	Entradas:	qe	caudal de entrada al depósito
--			      V 	voltaje aplicado a la resistencia eléctrica
--					a apertura de la valvula
--
--	Perturbación no medida:
-- 			   Te	temperatura de la corriente de entrada al depósito
--
--	Salidas: (que coinciden con los estados del sistema)
--			      h	nivel de líquido en el depósito			
--			      T	temperatura en el depósito
--
--	Parámetros conocidos:
--			rho	densidad
--			Cp	calor específico
--			Tamb	temperatura del ambiente exterior al depósito
--
--	Parámetros a estimar:
--			k	factor de fricción de la tubería de salida
--			Uamb	coeficiente de pérdida de calor al ambiente
--			A	área del depósito
--			R	resistencia eléctrica
--
-------------------------------------------------------------------------------
-------------------------------------------------------------------------------
--
--	- Te :  perturbación no medida ni considerada en el modelo
--	 modelo: se considera Te constante
--	 datos reales: se generaron con Te siguiendo una
--	 variación senoidal de frecuencia diaria
--    				  	
--------------------------------------------------------------------------------

DATA

 		REAL rho = 1000	"densidad del líquido (kg/m3)"
		REAL cp = 4168		"calor específico (J/kgºC)"
		REAL Tamb = 20		"temperatura del ambiente exterior (ºC)"
 		REAL Te = 40		"temperatura de la corriente de entrada (ºC)"	

	-- Parámetros a estimar sin escalar, aqui se dan unos valores iniciales
    	REAL k = 0.18 			"factor de fricción (m 2.5/h)"
    	REAL Uamb = 62 		"coeficiente de pérdida de calor al ambiente (W/ºC)"
    	REAL A =  0.47			"área del depósito (m2)"
    	REAL R  = 52			"resistencia eléctrica (ohmios)"
		REAL a  = 100 			"% de apertura de la valvula de descarga (%)"

    	REAL tsamp = 0.2		-- Periodo de muestreo de recogida de datos (h) 

    	INTEGER N1 = 1		-- Inicio horizonte de estimación
    	INTEGER NE = 451	-- Número de muestras, datos reales tomados
	-- Pesos en el coste para las diferencias entre predicciones y datos reales
    	REAL pesos[2] = {1.0, 100.0}			-- pesos[nsal]

	-- Valores medios de los datos experimentales, factor de escala
    	REAL media[2] = {0.48, 45.4} 		-- Nivel, Temperatura

	REAL Liminfh = 0.3
	REAL Limsuph = 1
	REAL LiminfT = 20
	REAL LimsupT = 80


DECLS
	-- Variables manipuladas  
    	REAL V				"voltaje aplicado a la resistencia eléctrica (V)"
    	REAL qe			   "caudal de entrada al depósito (m**3/h)"	
		
	-- Salidas medidas del proceso
    	REAL h		 		"nivel de líquido en el depósito (m)"	    
    	REAL T				"temperatura del depósito (ºC)"
		REAL F				"caudal de salida (m3/h)"

   	REAL y_modelo[nsal] 		-- Salidas del modelo que se quiere ajustar
    	REAL y_real[nsal]		   -- Salidas muestreadas experimentales
    	REAL u_real[nent]		   -- Entradas experimentales

     	REAL coef = 0			-- Coeficiente que activa la función de coste
    	DISCR REAL J1				-- Subtotal del índice de coste
		DISCR REAL J2				-- Subtotal del índice de coste

    	BOOLEAN Sample = TRUE		-- Variable bucles

	-- Tablas correspondientes a los datos experimentales
   	 TABLE_1D 	tab1, tab2, tab3, tab4

	-- El índice que se minimiza
		DISCR REAL J_costo

	-- Vector con la funcion objetivo y las restricciones no lineales
		REAL F_optim[5]
		
-- Variable auxiliar para el calculo de maximos y minimos de nivel 
		REAL hmin, hmax
	-- Variable auxiliar para el calculo de maximos y minimos de temperatura
		REAL Tmin, Tmax


INIT

   	J_costo = 0.0
    	J1 = 0.0
    	J2 = 0.0

	-- Entradas
    	readTableCols1D ("datos_potencia.txt",1,2,tab1)
    	u_real[1] = linearInterp1D (tab1, TIME)
    	V = u_real[1]
       
   	readTableCols1D ("datos_caudal_entrada.txt",1,2,tab2)
   	u_real[2] = linearInterp1D (tab2, TIME)
    	qe = u_real[2]

	-- Salidas
   	readTableCols1D ("datos_nivel.txt",1,2,tab3)  
    	y_real[1] = linearInterp1D (tab3, TIME)		-- valor real, leído de planta
    	h = y_real[1]							-- valor inicial para el modelo
    	y_modelo[1] = y_real[1]				-- valor para el modelo

    	readTableCols1D ("datos_temperatura.txt",1,2,tab4)  
    	y_real[2] = linearInterp1D (tab4, TIME)		-- valor real, leído de planta
    	T = y_real[2]							-- valor inicial para el modelo
    	y_modelo[2] = y_real[2]				-- valor del modelo

					
	-- coeficiente que activa el índice de coste
		coef = 1
    	--coef = 1 AFTER N1*tsamp
		
   --  inicializacion de variables auxiliares
		hmin = h
     	hmax = h
	   Tmin = T
	   Tmax = T

	-- calculo de la funcion de costo a minimizar
	DISCRETE

      	WHEN (Sample) THEN

		-- Subtotales del índice de coste
		    	J1 += coef*((y_modelo[1] - y_real[1])/media[1])**2   
          	J2 += coef*((y_modelo[2] - y_real[2])/media[2])**2
          	J_costo = pesos[1]*J1 + pesos[2]*J2

          Sample = FALSE
          Sample = TRUE AFTER tsamp
      	END WHEN

CONTINUOUS

	-- Entradas medibles
   	u_real[1] = linearInterp1D (tab1, TIME)
    	V = u_real[1]				-- Voltaje aplicado a la resistencia

    	u_real[2] = linearInterp1D (tab2, TIME)
    	qe = u_real[2]				-- Caudal de entrada
    
	-- Salidas medibles
    	y_real[1] = linearInterp1D (tab3, TIME)	-- Nivel
    	y_real[2] = linearInterp1D (tab4, TIME)	-- Temperatura

-- ECUACIONES DIFERENCIALES DEL MODELO (mezcla perfecta)

	-- Balance de materia
    	A*h' = qe - F
		
		F = a*k*sqrt(h)/100								

	-- Balance de energía
    --	A*h*rho*cp*T' = qe*rho*cp*Te- F*rho*cp*T + (3600*V*V/R) - 3600*Uamb*(T-Tamb)
	 A*h*rho*cp*T' = qe*rho*cp*(Te-T)  + (3600*V*V/R) - 3600*Uamb*(T-Tamb)

	-- calcular las restricciones no lineales (expresiones  c(x) <= 0)

		F_optim[1] = J_costo
		F_optim[2] = minabs(h,hmin) 
	   F_optim[3] = maxabs(h,hmax) 
	   F_optim[4] = minabs(T,Tmin) 
	   F_optim[5] = maxabs(T,Tmax) 
		
--		F_optim[2] = MATH.min(h,Liminfh) - Liminfh
--	   F_optim[3] = MATH.max(h,Limsuph) - Limsuph
--	   F_optim[4] = MATH.min(T,LiminfT) - LiminfT
--	   F_optim[5] = MATH.max(T,LimsupT) - LimsupT




	-- Asignación de valores del modelo
	    y_modelo[1] = h
	    y_modelo[2] = T

END COMPONENT
