 /*-----------------------------------------------------------------------------------------
 LIBRARY: REACTOR
 COMPONENT: reactorAB
 PARTITION: partition1
 EXPERIMENT: optimSnoptIdas
 TEMPLATE: TRANSIENT
 CREATION DATE: 06/03/2017
-----------------------------------------------------------------------------------------*/
USE OPTIM_METHODS

EIDAS calculoSens

CONST INTEGER numC = 3				-- número de restricciones del problema de optimización
CONST INTEGER numU = 2				-- número de variables de decisión

FUNCTION INTEGER coste_y_restricciones (IN REAL esnopt_x[], IN INTEGER needF, OUT REAL esnopt_F[], IN INTEGER explicit_derivatives, IN INTEGER needG, OUT REAL esnopt_G[])
	DECLS
		REAL arr_quad[ numC + 1 ]						-- array de integrales temporales
		REAL arr_quadsen[ (numC + 1) * numU ]		-- array de sensibilidad de las integrales
		REAL arr_quadsen_p[ (numC + 1) * numU ]	-- array de la derivada temporal de la sensibilidad de las integrales

	BODY
		FOR( j IN 1,numC + 1 )
			arr_quad[ j ] = 0
		END FOR
		FOR( j IN 1, (numC + 1) * numU )
			arr_quadsen[ j ] = 0
			arr_quadsen_p[ j ] = 0
		END FOR
		
		-- Actualizar las variables de decisión a los valores que propone el optimizador	
		calculoSens.setU( "q", esnopt_x[1] )
		calculoSens.setU( "Fr", esnopt_x[2] )
		
		-- Reiniciar el cálculo de sensibilidades con los nuevos valores de los parámetros	
		calculoSens.reiniciar()	
		
		calculoSens.calculateSens( TSTOP )
		calculoSens.getQuad( arr_quad, arr_quadsen )		-- cálculo de las cuadraturas y sus sensibilidades
		calculoSens.getQuadN( arr_quadsen_p, 1 )			-- cálculo de la derivada de las sensibilidades de las cuadraturas

		--	Introduccion de los gradientes de la funcion de costo y restricciones a SNOPT
		--	si la funcion es de camino, hay que introducir la derivada parcial de la cuadratura ( arr_quadsen )
		-- en caso contrario, si la restricción es de punto final, al usar arr_quadsen_p se usa la derivada parcial del valor de la función.
		IF ( explicit_derivatives == 1) THEN

			esnopt_G[1] = arr_quadsen_p[ 1 ]
			esnopt_G[2] = arr_quadsen_p[ 2 ]

			esnopt_G[3] = arr_quadsen_p[ 3 ]
			esnopt_G[4] = arr_quadsen_p[ 4 ]

			esnopt_G[5] = arr_quadsen[ 5 ]
			esnopt_G[6] = arr_quadsen[ 6 ]

			esnopt_G[7] = arr_quadsen[ 7 ]
			esnopt_G[8] = arr_quadsen[ 8 ]
		END IF

		-- Introducción de los valores de las funcines a SNOPT
		--	si la funcion es de camino, hay que introducir la cuadratura ( arr_quad )
		-- en caso contrario, si la restricción es de punto final, al usar F_optim se usa el valor de la función.
		esnopt_F[1] = F_optim[1]
		esnopt_F[2] = F_optim[2]
		esnopt_F[3] = arr_quad[3]
		esnopt_F[4] = arr_quad[4]

		RETURN 0
END FUNCTION

-- llamada a la función de residuos por parte de IDAS. No se debe modificar
-- IDAS propone unos valores de las variables y de las derivadas y pide calcular los residuos
FUNCTION NO_TYPE funcionResiduos( IN REAL tiempo, IN REAL var[], IN REAL der[], IN INTEGER numPar, IN STRING nomPar[], IN REAL par[], OUT REAL res[] )
	DECLS
	OBJECTS
	BODY
		FOR( i IN 1,numPar )
			setValueReal( nomPar[i], par[i] )
		END FOR

		FRES_ARGS( tiempo, var, der, res, FALSE )

END FUNCTION

-- llamada a las funciones de cuadraturas por parte de IDAS, para el cálculo de sensibilidades
-- utilizar el vector F_optim que se actualiza al llamar a los residuos
FUNCTION NO_TYPE funcionQuadraturas( IN REAL tiempo, IN REAL var[], IN REAL der[], IN INTEGER numPar, IN STRING nomPar[], IN REAL par[], OUT REAL quad[] )
	DECLS
	OBJECTS
	BODY
		FOR( i IN 1,numPar )
			setValueReal( nomPar[i], par[i] )
		END FOR

		FRES_ARGS( tiempo, var, der, der, FALSE )
		
		-- funciones de cuadraturas
		FOR( i IN 1,4 )
			quad[i] = F_optim[i]
		END FOR

END FUNCTION

FUNC_PTR<ptrFunRes> ptrRes = funcionResiduos
FUNC_PTR<ptrFunRes> ptrQuad = funcionQuadraturas

EXPERIMENT optimSnoptIdas ON reactorAB.partition1
   DECLS
		REAL dec_var[numU]      	-- valor inicial de las variables de decision, size numU. Necesario para las llamadas a ESnopt
		
	  	REAL xlow[numU]				-- valor inferior de las variables de decision, size numU. Necesario para las llamadas a ESnopt
   	REAL xupp[numU]				-- valor superior de las variables de decision, size numU. Necesario para las llamadas a ESnopt
  
   	REAL Flow[numC + 1]			-- valor inferior de la funcion objetivo y las restricciones, size (numC + función de coste). Necesario para las llamadas a ESnopt
   	REAL Fupp[numC + 1]			-- valor superior de la funcion objetivo y las restricciones, size (numC + función de coste). Necesario para las llamadas a ESnopt

		INTEGER calcularSens = 1
		INTEGER infoESnopt = 0
   OBJECTS
		VECTOR_STRING nombresX		-- variable auxiliar para la inicialización de las sensibilidades
--		VECTOR_STRING nombresU		-- variable auxiliar para la inicialización de las sensibilidades

   INIT
      -- initial values for state variables
      cA = 2.891
      Tr = 51.5
		T = T0

   BOUNDS
      -- Set equations for boundaries: boundVar = f(TIME;...)
      Fr = 50
      q = 2

   BODY
      -- creates an ASCII file with the results in table format
      -- REPORT_TABLE("results.rpt", "*",TRUE)

      -- set the debug level (valid range [0,4])
      DEBUG_LEVEL = 1
      -- select default integration solver. Valid methods are IDAS (_SPARSE), DASSL(_SPARSE), CVODE_BDF(_SPARSE), CVODE_AM, RK4, EULER, AM1, AM2 and AM4
      IMETHOD = IDAS -- default is DASSL, recommended is either IDAS or IDAS_SPARSE  
      -- Settings for different actions. Valid actions are: SEV_DISABLE, SEV_NONE, SEV_WARNING, SEV_ERROR, SEV_KILLPOINT, SEV_FATAL
      eSetErrorAction(ERR_BAD_OPER,SEV_KILLPOINT) -- Detect bad numerical operations (eg division by zero), default is SEV_NONE, recommended is SEV_KILLPOINT
      eSetErrorAction(ERR_NAN_INF, SEV_KILLPOINT) -- Detect NaN or Inf values, default is SEV_NONE, recommended is SEV_KILLPOINT
      eSetConfig(CFG_FORCE_STOP_CINT,TRUE) -- Force or not to stop the solver each CINT. Sometimes you can speed up the simulation if you select FALSE
      -- set tolerances and other important inputs
      REL_ERROR = 1e-06 -- transient solver relative tolerance
      ABS_ERROR = 1e-06 -- transient solver absolute tolerance
      TOLERANCE = 1e-06 -- steady solver relative tolerance
      INIT_INTEG_STEP = -1 -- initial integration step size (-1 means use the solver estimation) 
      MAX_INTEG_STEP = -1 -- maximum integration step size (-1 means use the solver estimation) 
      NSTEPS = 1 -- Only for explicit solvers use CINT/NSTEPS as integration step size
--      REPORT_MODE = IS_EVENT -- by default it reports results at every CINT and event detection. Other valid options are IS_STEP, IS_CINT and IS_MANUAL_REFRESH 

		-- Añadir las variables del jacobiano para el cálculo de sensibilidades
		-- NO MODIFICAR, LO REALIZA AUTOMÁTICAMENTE
		-- SINTAXIS:
		-- addX( nombre de la variable, valor inicial de la misma, booleano que indica si la variable es T_DYNAMIC )
		nombresX = getVariablesWithCategory( T_JACOBIAN_VARS )		
		FOR( i IN 1,nombresX.size() )
			calculoSens.addX( nombresX.at(i), getValueReal( nombresX.at(i) ), ( getVarCategory (nombresX.at(i)) == 2 ) )
		END FOR
		
		-- Iniciliazacion del algoritmos de optimización
		-- inicialización de las variables de decisión, y los límites
		-- SINTAXIS:
		-- addU( nombre de la variable, valor inicial de la misma, booleano que indica si el parámetro es un valor inicial )
		calculoSens.addU( "q", q, FALSE ) 
		calculoSens.addU( "Fr", Fr, FALSE ) 
		dec_var[1] = q
		dec_var[2] = Fr
		
		xlow[1] = 0.5
		xupp[1] = 5

		xlow[2] = 5
		xupp[2] = 90

		TIME = 0
      TSTOP = 15
      CINT = 0.1
		
		-- Inicializacion del algoritmo de integracion y cálculo de sensibilidades
		-- Inicializacion del algoritmos de optimización
		-- SINTAXIS:
		-- iniciarRes( tiempo inicial desde el que calcular las sensibilidades, puntero a la función de residuos definida previamente )
		-- iniciarQuad( número de cuadraturas a calcular, puntero al cálculo de las funciones de cuadratura )
		calculoSens.iniciarRes(TIME, ptrRes)
		calculoSens.iniciarQuad(numC + 1, ptrQuad)
		
		-- Configurar la tolerancia en el cálculo de sensibilidades (opcional)
		calculoSens.setTol( 1e-5 )
		
		-- inicialización de los límites de las restricciones y la función de coste
		Flow[1] = -1.0e6
		Fupp[1] =  1.0e6			
		
		Flow[2] = Liminfx
		Fupp[2] = 1			
		
		Flow[3] = 0
		Fupp[3] = 0

		Flow[4] = 0
		Fupp[4] = 0
					
		-- Representacion en el monitor de los valores previos a la optimizacion
		WRITE ("\n\n\t\t------------------- limits before optimization \n")
		FOR ( i IN 1,numU )
				WRITE ("\t\t var_dec[%d]:\t %g \t %g \t %g \n", i, xlow[i], dec_var[i], xupp[i])
		END FOR
		FOR ( i IN 1,numC + 1 )
				WRITE ("\t\t F_optim[%d]:\t %g \t %g \t %g \n", i, Flow[i], F_optim[i], Fupp[i])		
		END FOR
		WRITE ("\n\t\t------------------- END of limits  \n\n")

   	--Optimization extern routine call
		setSilentMode(TRUE)
		SET_REPORT_ACTIVE("#MONITOR",FALSE)
		
		esnopt_init (numU, numC)
		esnopt_set_variables_bounds_and_initial_values (xlow, xupp, dec_var)
		esnopt_set_constraints_bounds_and_initial_values (Flow, Fupp, F_optim)
		esnopt_set_cost_function_and_constraints (coste_y_restricciones)
		esnopt_set_explicit_derivatives ( calcularSens )
		esnopt_set_function_precision ( 1.0e-5 )
		esnopt_set_iterations_limit (200) 
		infoESnopt = esnopt ()
		
		-- Final de la optimización, obtención de los resultados para la simulación.
		setSilentMode(FALSE)
		SET_REPORT_ACTIVE("#MONITOR",TRUE)
		
		esnopt_print_data ()
		esnopt_get_variables_values(dec_var)
		esnopt_free ()	
		
		-- Llamada al integrador
		RESET()
--			SET_INIT_ACTIVE(TRUE)

		-- Modificación de los parámetros con los nuevos valores obtenidos		
		q = dec_var[1]
		Fr = dec_var[2]

      TIME = 0
		CINT = 0.1

  	   INTEG()
		
		-- Representacion en el monitor de los valores posteriores a la optimizacion
		WRITE ("\n\n\t\t------------------- The optimiser succeed, returned the value: %d. \n", infoESnopt)
		WRITE ("\n\n\t\t------------------- constraints and decision variables after optimization \n")
		FOR ( i IN 1,numU )
			WRITE ("\t\t var_dec[%d]:\t %g \t %g \t %g \n", i, xlow[i], dec_var[i], xupp[i])
		END FOR
		WRITE ("\n\n\t\t------------------- cost function and constraints after optimization \n")
		FOR ( i IN 1,numC + 1 )
			WRITE ("\t\t F_optim[%d]:\t %g \t %g \t %g \n", i, Flow[i], F_optim[i], Fupp[i])		
		END FOR
		WRITE ("\n\t\t------------------- end   \n\n")

END EXPERIMENT
