[[PageOutline]]

= Package NonLinGloOpt =

Este paquete TOL una API para usuario final al sistema de optimización global 
no lineal [http://ab-initio.mit.edu/nlopt NLopt] de Steven G. Johnson.

== Definición del problema ==

El tipo de problemas que resuelve este sistema se puede formular matemáticamente del siguiente modo generalizado:

Se quiere maximizar una función objetivo arbitraria:[[BR]] [[BR]]
[[LatexEquation( \underset{x\in\Omega}{\min}f\left(x\right)\wedge f:\Omega_f\subset\mathbb{R}^{n}\longrightarrow\mathbb{R} \wedge n\geq1 )]] [[BR]]

Sujeta a las siguientes
 * '''Restricciones de intervalo''': [[BR]] [[BR]]
   [[LatexEquation( l_{i}\leq x\leq u_{i}\wedge-\infty\leq l_{i}<u_{i}\leq\infty\wedge i=1\ldots n )]] [[BR]] [[BR]]
   Si una variable no tuviera cota inferior entonces sería [[LatexEquation( l_{i} = -\infty  )]] y si no tuviera cota superior entonces sería [[LatexEquation( u_{i} = \infty  )]] 
   A la región [[BR]] [[BR]]
   [[LatexEquation(\underset{i=1}{\overset{n}{\prod}}\left[l_{i},u_{i}\right] \subset\mathbb{R}^{n})]][[BR]] [[BR]] 
   se le llamará hiperrectángulo de búsqueda y en el caso de optimización global
   deberá ser siempre acotado en todas las dimensiones. [[BR]] [[BR]]
 * '''Restricciones de desigualdad''' arbitrarias: [[BR]] [[BR]]
   [[LatexEquation( g_{j}\left(x\right)\leq0 \wedge g_j:\Omega_{g_j}\subset\mathbb{R}^{n}\longrightarrow\mathbb{R} \forall j=1\ldots n_{D} \wedge n_{D} \ge 0)]]  [[BR]] [[BR]]
 * '''Restricciones de igualdad''' arbitrarias: [[BR]] [[BR]]
   [[LatexEquation( h_{k}\left(x\right)=0 \wedge h_k:\Omega_{h_k}\subset\mathbb{R}^{n}\longrightarrow\mathbb{R} \forall k=1\ldots n_{I} \wedge n_{I} \ge 0)]]  [[BR]] [[BR]]
   
A la región del espacio [[LatexEquation( \Omega\subset\mathbb{R}^{n} )]] de los puntos que cumplen todas las restricciones del problema le llamaremos región de factibilidad del problema. A los puntos de dicha región les llamaremos factibles o interiores y al resto no factibles o exteriores.

Hay una condición que se suele sobreentender por su obviedad, pero que es a veces la causa de que no funcionen los algoritmos: un punto factible debe pertenecer al dominio de la función objetivo, y al de las funciones de restricción, o sea: 

[[LatexEquation( \Omega\subset\Omega_{f}\cap\underset{j=1}{\overset{n_{D}}{\bigcap}}\Omega_{g_{j}}\cap\underset{k=1}{\overset{n_{I}}{\bigcap}}\Omega_{h_{k}} )]]

Por ejemplo, si en la función objetivo aparece el logaritmo de una variable se debe incluir alguna restricción o conjunto de ellas que obligue a que dicha variable sea positiva. Evidentemente, si no hay restricciones, pues todas ellas son en principio opcionales, el dominio de la función objetivo debe ser todo el espacio [[LatexEquation( \Omega_f = \Omega = \mathbb{R}^{n} )]]. 
  
Obsérvese que si lo que se quiere es minimizar la función objetivo basta con tomar la función opuesta [[LatexEquation( -f(x) )]].

Nótese que a ninguna de las funciones se le ha exigido ser diferenciable, y ni tan siquiera ser continua; aunque, como se verá más tarde, el usuario deberá informar de qué clase analitica es el problema para poder saber qué métodos son aplicables. 

En lo que afecta a este sistema, la clase analítica de una función ha de ser la 
máxima de las condiciones cumplidas por la función de entre las siguientes

 || '''Arbitraria''' || {{{Real ARBITRARY            = -1}}} || si la función no es continua. ||
 || '''Continua''' ||  {{{Real CONTINUOUS           =  0}}} || si la función es continua pero no diferenciable. ||
 || '''Diferenciable''' ||  {{{Real DIFFERENTIABLE       =  1}}} || si la función es continua y diferenciable sólo una vez. ||
 || '''Doblemente diferenciable''' ||  {{{Real TWICE_DIFFERENTIABLE =  2}}}  || si la función es continua y al menos doblemente diferenciable. ||

La clase analítica del problema será la mínima de entre las de su funciones objetivo 
y las de sus restricciones de igualdad y desigualdad. El usuario debe conocer esta 
característica del problema para informarlo en su definición.
 
== Guia del usuario ==

El objeto encargado de gestionar la optimización será una instancia de {{{Class @PipeLine}}} definida en [source:/tolp/OfficialTolArchiveNetwork/NonLinGloOpt/pipeLine.tol pipeLine.tol] que consiste en una
serie de etapas sucesivas en cada una de las cuáles actúa un método iterativo determinado de tal forma que cada uno comienza en el óptimo alcanzado por el anterior en un proceso de refinamiento de la solución (pipe-line)
 
El problema se define una sola vez como una instancia de {{{Class @Problem}}} definida en [source:/tolp/OfficialTolArchiveNetwork/NonLinGloOpt/problem.tol problem.tol]  

La cadena de métodos se define como un conjunto de métodos que será instancias de {{{Class @Method}}} definida en [source:/tolp/OfficialTolArchiveNetwork/NonLinGloOpt/method.tol method.tol] y que se irán añadiendo a la instancia {{{Class @PipeLine}}} mediante el método {{{AddMethod}}}. Si no se especifica ningún método se utilizará uno acorde con la definición del problema.

El algoritmo de optimización de cada método estará implementado como parte de un motor de optimización que estará definido como una instancia de la clase abstracta {{{Class @Engine}}} definida en [source:/tolp/OfficialTolArchiveNetwork/NonLinGloOpt/engine.tol engine.tol] y de la cual por el momento sólo estará disponible la herencia {{{@Engine_NLopt}}} definida en [source:/tolp/OfficialTolArchiveNetwork/NonLinGloOpt/engine_nlopt.tol engine_nlopt.tol] y que corresponde a la librería C++ [http://ab-initio.mit.edu/wiki/index.php/NLopt_C-plus-plus_Reference nlopt].

Cada método tendrá definido una serie de parámetros que controlen los criterios de parada y que se incorporarán como instancia de {{{Class @StopCriteria}}} definida en [source:/tolp/OfficialTolArchiveNetwork/NonLinGloOpt/stopCriteria.tol stopCriteria.tol]

=== Formulación de funciones objetivo y de restricción ===

Tanto si se calcula el gradiente como si no, las funciones objetivo y las de restricción deben tener la API TOL general común a todas ellas:

{{{
#!cpp
Real eval (Matrix x, Matrix gradient)
}}}

siendo tanto {{{x}}} como {{{gradient}}} matrices de {{{n}}} filas y 1 columna.

Aún en el caso de que se quiera usar un método que precise el gradiente, no  
necesariamente se exigirá calcular el mismo cada vez que se quiera evaluar 
la función. En tales casos se le pasará la matriz vacía para indicar que 
no se haga el cálculo. Cabría pensar si no sería mejor hacer una función 
para la evaluación de la función y otra para el gradiente, pero resulta que 
muy a menudo el cálculo del gradiente utiliza en parte las mismas componentes 
que la propia función, por lo que resulta más eficiente que se calculen de 
forma conjunta. La función de evaluación debería por tanto tener este aspecto:

{{{
#!cpp
Real eval (Matrix x, Matrix gradient)
{
  ...
  Real y = ...;
  Real If(Rows(gradient), {
    ...
    Matrix gradient := ...;
    True
  });
  y
};
}}}

Si el problema no es diferenciable o el cálculo del gradiente es muy costoso 
y se van a usar sólo algoritmos que no requieren el cálculo del gradiente, 
entonces se puede hacer caso omiso del argumento {{{ Matrix gradient }}} y 
calcular sólo el valor de la función.

A veces es conveniente que las funciones del problema se definan como métodos 
de una clase, bien porque manejan información auxiliar compleja, bien porque 
se trate de funciones pertenecientes a una familia paramétrica y sea más 
cómodo implementarlas de este modo. En esos casos, el método es un objeto 
TOL común a todas las instancias de la clase, por lo que no es suficiente 
por sí misma y habrá que pasar también la instancia de la clase para la que 
queremos evaluar, y lo haremos mediante un conjunto como éste
{{{
#!cpp
[[NameBlock owner, Code owner::method ]]
}}}

=== Características obligatorias en la definición del problema  ===

Las siguientes características del problema deben ser cumplimentadas por el usuario 
de forma obligatoria y deben ser congruentes entre sí.

{{{
#!cpp
//Set sign ==- 1 for minimization, sign == 1 for maximization
  Real sign;
  
//If true, a global optimal point will be searched. Else just a 
//local one.
  Real neededGlobal;
 
//Minimum analytical class of target and constraining functions must be one
//of elements of 
//Set NonLinGloOpt::AnalyticalClass =
//{[[
//  Real ARBITRARY            = -1,  // Arbitrary functions
//  Real CONTINUOUS           =  0,  // Class C_0: Continuous
//  Real DIFFERENTIABLE       =  1,  // Class C_1: Differentiable
//  Real TWICE_DIFFERENTIABLE =  2   // Class C_2: Twice differentiable
//]]};
  Real id_analyticalClass;

//Indicates if the gradient is known analytically or numerically. 
//Must be an element of 
//Set NonLinGloOpt::GradientKnowledge =
//[[
//Real NONE        = 0,  //Doesn't exists or it's too hard to calculate
//Real ANALITICAL  = 1,  //Use analitical gradient
//Real NUMERICAL   = 2   //Use numerical gradient
//]];
  Real id_useGradient;
  
//Number of variables
  Real n;

//Target function could be a simple Code or a class method of a
//particular instance.
// - Real target(Matrix x, Matrix gradient)
// - Set target = [[NameBlock owner, Code method]]
//When argument gradient is not empty, then the function, or the 
//method must update it. 
  Anything target;

}}}

=== Características opcionales en la definición del problema  ===

Hay otra serie de carcterísticas del problema que son opcionales pues o bien tienen 
un valor por defecto predefinido o bien existe un método interno que les da un valor 
adecuado.
{{{
#!cpp
//Lower bounds for the variables 
//These must be specified by an nx1 column matrix with the minimum value for each 
//variable, or -1/0 if the variable is not bounded below
//If there are no lower bounds at all then it's the empty matrix
  Matrix lower_bounds = Rand(0,0,0,0);

//Upper bounds for the variables 
//These must be specified by an nx1 column matrix with the maximum value for each 
//variable, or -1/0 if the variable is not bounded above
//If there are no upper bounds at all then it's the empty matrix
  Matrix upper_bounds = Rand(0,0,0,0);

//Set of inequality constraints g(x)<=0
//Set of Code elements or [[NameBlock owner, Code target]] for methods
//Constraining function could be a simple Code or a class method of a
//particular instance.
// - Real g (Matrix x, Matrix gradient)
// - Set g = [[NameBlock owner, Code method]]
//When argument gradient is not empty, then the function, or the 
//method must update it. 
  Set inequations = Copy(Empty);

//Set of equality constraints h(x)<=0
//Set of Code elements or [[NameBlock owner, Code target]] for methods
//Constraining function could be a simple Code or a class method of a
//particular instance.
// - Real h (Matrix x, Matrix gradient)
// - Set h = [[NameBlock owner, Code method]]
//When argument gradient is not empty, then the function, or the 
//method must update it. 
  Set equations = Copy(Empty);

//Tolerance for the evaluation of inequality constraints to avoid  
//rounding errors
  Real inequationTolerance = 1E-8;

//Tolerance for the evaluation of the equality constraints to avoid  
//rounding errors
  Real equationTolerance = 1E-8;

}}}


=== Creación y lanzamiento de la optimización ===

El esqueleto del código de uso del paquete sería el siguiente
{{{
#!cpp
//Defines the problem and the optimization pipe-line
NonLinGloOpt::@PipeLine pipe_line = [[
  //The problem (or problems) to be optimized. In this version just the first 
  //one will be solved
  Set problems = [[ NonLinGloOpt::@Problem problem = [[
    ...
  ]] ]];
  //The system will show a trace every this number of evaluations of target 
  //function 
  Real verboseEach = 100
]];

//Runs the optimization sequence
Real pipe_line::optimize(?);

//Access to optimal point
Matrix pipe_line::x.opt;
//Access to optimal value of target function
Real pipe_line::y.opt;

}}}


[[Image(source:/tolp/OfficialTolArchiveNetwork/NonLinGloOpt/doc/opt_methods.png)]]

== Tutorial == 

La mejor forma de aprender a usar cualquier cosa es empezar a usarla de forma
controlada, para lo cual se dispone de una batería de ejemplos a disposición de 
los usuarios.

=== Función objetivo bivariante no lineal con restricciones de desigualdad no lineal === 

En primer lugar se expondrá el 
[http://ab-initio.mit.edu/wiki/index.php/NLopt_Tutorial ejemplo] 
proporcionado en el tutorial de NLopt 

Se quiere encontrar un mínimo local de la función objetivo:[[BR]] [[BR]]
[[LatexEquation( \underset{x}{\min}f\left(x_1,x_2\right)=\sqrt{x_2} )]] [[BR]]

Sujeta a las 
 * restricciones de intervalo: [[BR]] [[BR]]
   [[LatexEquation(  -x_2\leq 0 )]] [[BR]] [[BR]]
 * y a las restricciones de desigualdad: [[BR]] [[BR]]
   [[LatexEquation( g_{1}\left(x\right)=\left(a_{1}x_{1}+b_{1}\right)^3 -x_2 \leq0 \wedge a_{1}=2\wedge b_{1}=0 )]]  [[BR]]
   [[LatexEquation( g_{2}\left(x\right)=\left(a_{2}x_{1}+b_{2}\right)^3 -x_2 \leq0\wedge a_{2}=-1\wedge b_{2}=1 )]]  [[BR]] [[BR]]

Nótese cómo se han reescrito las inecuaciones de restricción en la forma canónica [[LatexEquation( g\left(x\right) \leq0 )]] 

El punto inicial de la búsqueda es el punto factible

   [[LatexEquation(  x_1=1.234; x_2=5.678 )]] [[BR]]

Obsérvese que se trata de un problema al menos dos veces diferenciable.

El código TOL completo de este ejemplo está disponible 
[source:/tolp/OfficialTolArchiveNetwork/NonLinGloOpt/test/test_0001/test.tol aquí]

A continuación veremos ese mismo código por partes 

==== Función objetivo ====

Obérvese cómo la función objetivo incluye la evaluación condicional del gradiente
y cómo ésta hace uso de los cálculos empleados en evaluar la propia función.

{{{
#!cpp
////////////////////////////////////////////////////////////////////////////////
//Target function defined as simple code
Real my_function(Matrix x, Matrix grad)
////////////////////////////////////////////////////////////////////////////////
{
  Real y = Sqrt(MatDat(x,2,1));
  Real If(Rows(grad), {
    Matrix grad := Col(0.0,0.5/y);
    True
  });
  y
};
}}}

==== Restricciones de desigualdad ====
Como ambas restricciones de desigualdad responden a un mismo patrón es 
conveniente no repetir su formulación usando una clase para su definición.

{{{
#!cpp
////////////////////////////////////////////////////////////////////////////////
//When a set of constraining functions is part of a parametric family
//of functions we can define a class to facilitate the definition
Class @my_inequation 
////////////////////////////////////////////////////////////////////////////////
{
  Real a;
  Real b;
  Real constraint(Matrix x, Matrix grad)
  {
    Real x1 = MatDat(x,1,1);
    Real x2 = MatDat(x,2,1);
    Real ax1b = a*x1 + b;
    Real If(Rows(grad), {
      Matrix grad := Col(3 * a * ax1b^2,-1.0);
      True
    });
    Real g = (ax1b^3 - x2);
    g
  }
};

////////////////////////////////////////////////////////////////////////////////
//Inequality constraining functions defined as class instances
////////////////////////////////////////////////////////////////////////////////
@my_inequation ineq1 = [[Real a= 2,Real b=0]];
@my_inequation ineq2 = [[Real a=-1,Real b=1]];

}}}
 
=== Definición del problema y sus características ===

{{{
#!cpp
//Don't forget to load the package if it has not yet been loaded
#Require NonLinGloOpt;

//Defining the problem
NonLinGloOpt::@Problem problem = [[
//Optimization method. If not specified the system will choose one
//by calling method select_automatic_algorithm
//In this case we want to use MMA (Method of Moving Asymptotes)
// http://ab-initio.mit.edu/wiki/index.php/NLopt_Algorithms#MMA_.28Method_of_Moving_Asymptotes.29 
  Real id_algorithm = NonLinGloOpt::Algorithm::LD_MMA;  
  //We want to minimize the target function
  Real sign = -1;
  //We want just a local minimum
  Real neededGlobal = False;
  //The problem is almost twice differentiable
  Real id_analyticalClass = 
    NonLinGloOpt::AnalyticalClass::TWICE_DIFFERENTIABLE;
  //The gradient of all functions is known 
  Real id_useGradient = True;
  //Target given as simple Code
  Anything target = my_function;
  //Lower bounds 
  Matrix lower_bounds = Col(-1/0,  0); // -x2 <= 0
  //No upper bounds. You can skip next line
  Matrix upper_bounds = Col( 1/0, 1/0); 
  //There are two non linear constraining inequations 
  Set inequations = [[
   [[ineq1, ineq1::constraint]], //Inequation given by class method
   [[ineq2, ineq2::constraint]]  //Inequation given by class method
  ]];
  //Initial values
  Matrix x0 = Col(1.234, 5.678); 
  //Stopping criteria of relative tolerance for both x and y
  Real relativeTolerance = 1E-4; //
  //Stopping criteria of maximum running time
  Real maxTime = 60; 
  //Functions will be traced every "verboseEach" evaluations
  Real verboseEach = 1
]];

}}}

=== Creación y lanzamiento de la optimización ===

{{{
#!cpp

//Creating the optimizer instance
NonLinGloOpt::@Opt opt = problem::create_optimizer(?);

//Running the optimization
Real opt::optimize_problem(problem);

}}}

== Chequeo de robustez y eficiencia ==

Con el objetivo de chequear la validez de los resultados y la eficacia de la búsqueda se van 
a ir implementando una serie de tests estándar en el mundo de la optimización.
Para ello es preciso poder convertir los formatos .mod usados por AMPL y quizás también 
los .gms de GAMS.
También sería importante lo contrario, convertir los problemas de NonLinGloOpt a esos formatos.

Estas son algunas de las baterías de tests más importantes:

  * [http://www.mat.univie.ac.at/~neum/glopt/coconut/Benchmark/Benchmark.html COCONUT]: 
    A benchmark for global optimization and constraint satisfaction (Formatos: AMPL,GAMS,DAG)
  * [http://www.mcs.anl.gov/~more/cops/ COPS]: 
    Large-Scale Optimization Problems (Formato:AMPL)
  * [http://www.gamsworld.eu/performance/selconglobal/selcongloballib.htm SelConGlobal]:
    Selected Continuous Global Optimization Test Problems (Formato:GAMS)

== Guía del programador y del usuario avanzado == 

El sistema NLopt contiene 17 algoritmos principales con 41 variantes que intentan 
abarcar la mayor parte de los tipos de problemas tipo de la programación matemática 
no lineal y cuya descripción puede verse 
[http://ab-initio.mit.edu/wiki/index.php/NLopt_Algorithms#Augmented_Lagrangian_algorithm aquí] 
de forma detallada.

Existen métodos para optimización local y global, para funciones no lineales 
arbitrarias, continuas o diferenciables y con restricciones de igualdad y 
desigualdad no lineal o sin ellas.

El valor añadido del paquete TOL NonLinGloOpt pretende ser un sistema experto 
que, en el caso de que el usuario no proponga un método específico, busque uno
cuyas características sean compatibles con las del problema, a saber:

 || '''IsGlobal''' || 0,1 || Un método de optimización puede ser global o local, pero no ambas cosas a al vez, luego es determinante saber si se requiere un óptimo local o uno global para la selección del mismo. ||
 || '''AnalyticalClass''' || '''AR''': Arbitrary, [[BR]]'''C0''': Continuous, [[BR]]'''C1''': Differentiable, [[BR]]'''C2''': Twice differentiable || Es la clase analítica mínima capaz de tratar un algoritmo debe ser igual o mayor que la clase mínima de las funciones objetivo y de restricciones. Por ejemplo, si una de las funciones del problema es continua pero no diferenciable y todas las demás son diferenciables, la clase analítica del problema es C0 y sólo se podrán usar algoritmos e tipo AR ó C0, pero no los de tipo C1 ó C2. El usuario es el único que puede y debe ofrecer esta información acerca del problema.||
 || '''NeedsGradient''' || 0,1 || Los algoritmos de tipo C1 ó C2 pueden o no necesitar realmente calcular el gradiente. Por otra parte el que un problema sea clase C1 ó C2 no implica que se conozca analíticamente el gradiente, aunque siempre se puede calcular numéricamente, pero en cualquier caso puede que el coste computacional sea excesivo. Así pues, tanto si el gradiente no existe (AR,C0), como si existe pero no tenemos un método rápido de obtenerlo, se debe elegir un algoritmo que no utilice el gradiente. Es responsabilidad del usuario comunicar al sistema esta característica del problema.||
 || '''NeedsFeasiblePoint''' || 0,1 || Unos métodos precisan comenzar por un punto factible y otros no, y no siempre el usuario es capaz de proporcionar un punto adecuado. En estos casos el sistema debería utilizar un método que pueda comenzar por un punto exterior hasta encontrar uno interior y a partir de ahí se podría continuar con otro método cualquiera que fuera más eficiente. El propio sistema comprobará si la solución inicial proporcionada es factible. ||
 || '''NeedsBounds''' || 0,1 || En general todos los métodos globales precisan que se dé un intervalo de dominio a cada variable, es decir, que se defina un hiperrectángulo de búsqueda , mientras que los métodos locales no lo necesitan.||
 || '''SupportsBounds''' || 0,1 || La mayoría de los métodos soportan intervalos de dominio pero no siemrpe es así, por lo que hay que tenerlo en cuenta.||
 || '''SupportsInequalities''' || 0,1 || La mayoría de los métodos de optimización no soportan restricciones de desigualdad, aunque siempre existe la posibilidad de utilizar el método especial del Lagrangiano Aumentado para convertir el problema en otro sin restricciones que se pueda tratar con un método subsidiario adecuado.||
 || '''SupportsEqualities''' || 0,1 || Sólo unos pocos métodos soportan restricciones de igualdad, aunque lo mism o que con las desigualdades se puede aplicar el método especial del Lagrangiano Aumentado. Concretamente se puede aplicar sólo a las restricciones de igualdad manteniendo las de desigualdad si hay un método adecuado que lo pueda resolver.||
 || '''PreferredScale'''  || '''S''': Short [1...10], [[BR]]'''M''': Medium [11...100],[[BR]]'''L''': Large[101...1000][[BR]]'''E''': Extreme[1001...] || El número de variables con el que el algoritmo es capaz de tratar de forma eficiente ||

En la hoja de cálculo [https://www.tol-project.org/export/HEAD/tolp/OfficialTolArchiveNetwork/NonLinGloOpt/doc/NLopt_algorithms.xls NLopt_algorithms] se muestra un cuadro resumen de las características asociadas a cada variante algorítmica del paquete NLopt.

Se aconseja a los programadores y usuarios avanzados visitar de forma somera el [source:/tolp/OfficialTolArchiveNetwork/NonLinGloOpt código fuente] del paquete que a veces es la mejor forma de familiarizarse con un sistema. Estos son los archivos principales del sistema:

 * [source:/tolp/OfficialTolArchiveNetwork/NonLinGloOpt/NonLinGloOpt.tol NonLinGloOpt.tol]: Fichero raíz del paquete
 * [source:/tolp/OfficialTolArchiveNetwork/NonLinGloOpt/constant.tol constant.tol]: Constantes de codificación de algorimos, resultados, criterios de parada y otras características de la optimización y la definición de los problemas.
 * [source:/tolp/OfficialTolArchiveNetwork/NonLinGloOpt/problem.tol problem.tol]: Clase @Problem para la definición de problemas.
 * [source:/tolp/OfficialTolArchiveNetwork/NonLinGloOpt/opt.tol opt.tol]: Clase @Opt para la optimización de problemas mediante la API con NLOpt.
 * [source:/tolp/OfficialTolArchiveNetwork/NonLinGloOpt/CppTools/source/CppTools.cpp CppTools.cpp]: API C++ con NLOpt.