Paquete NonLinGloOpt

[Diana Edith Guere Lettich]

Estoy intentando realizar una optimización con restricciones utilizando este paquete pero no está considerando la restricción. Al parecer puede ser un problema el hecho de no definir la matriz gradiente.

Estoy adjuntando un archivo comprimido con los archivos necesarios. "OptSinRes" es la que realiza la optimización sin restricciones. Sin embargo al ejecutar "OptConRes" empieza a arrojar valores negativos en lugar de decir que no existe el valor óptimo.

Muchas gracias de antemano.

[Víctor de Buen Remiro]

Hola

Por lo que yo veo no hay ningún problema, pues la solución que encuentra

[NonLinGloOpt::@PipeLine::run] step 1 NLOPT 2.3.0.LN_COBYLA found optimal value =138171533.4993447 after 50 target evaluations in 45.84641500000001 seconds
Stop reason: MAXEVAL_REACHED

es perfectamente factible, tal y como se puede comprobar mediante la sentencia

Real opt_is_valid = pipe_line::problems::problem::is_valid(pipe_line::x.opt);

El hecho de que por el camino evalúe soluciones no factibles no tiene porqué ser ningún problema, eso depende del método empleado. Quizás en un método de punto interior podría ser válido tu argumento pero no en el algoritmo empleado: ​LN_COBYLA, que utiliza una sucesión de aproximaciones lineales para encontrar un mínimo local sin hacer uso del gradiente.

También se puede consultar los algoritmos disponibles y sus características en la página ​NLopt_Algorithms de la librería, desde el propio TOL tabulando el conjunto Set NonLinGloOpt::_.algorithmInfo , o en la hoja excel ​NLopt_algorithms.xls

Aún así si que hay un posible problema pues el criterio de parada es MAXEVAL_REACHED lo cual indica que no ha convergido todavía según los criterios establecidos.

 Stop criteria:
    Absolute tolerance for target: ?
    Relative tolerance for target: 1e-008
    Absolute tolerance for variables: 
    Relative tolerance for variables: 1e-008
    Maximum running time: ?
    Maximum number of evaluations: 50
    Sufficient target value: ?

El motivo es que sólo se le ha permitido efectuar maxEval=50 evaluaciones lo cual es realmente muy poco para un problema con 636 variables. En general es muy difícil saber qué número de evaluaciones debería ser suficiente. Si el problema es que tarda mucho tiempo es mejor poner un máximo de tiempo maxTime. Si nos conformamos con obtener cierta cantidad pondremos un valor de parada stopVal. Si se trata de una unidad discretizable como el dinero en el que nos basta con una aproximación de miles de euros podemos establecer una tolerancia absoluta de la función objetivo absoluteTolerance_target. Si no tenemos ni idea de cuál puede ser el nivel óptimo usaremos la tolerancia relativa de la función objetivo relativeTolerance_target. Las tolerancias absoluta y relativa de las variables hay que manejarlas con mucho cuidado y sólo usarlas cuando hay problemas graves de convergencia si se sabe muy bien lo que se hace. Para más detalles de los criterios deparada ver la página oficial de la librería ​NLOPT

Otras notas respecto al código utilizado

Las funciones Include admiten caminos relativos a la ubicación del fichero en el que se utilizan. Por ejemplo si el fichero "costos.bdt" está en el mismo directorio que donde se va a cargar no hace falta escribir el camino completo

   Set Coste    = IncludeBDT("C:\Users\ ... \Optimizacion\costos.bdt");

sino que bastaría poner

   Set Coste    = IncludeBDT("costos.bdt");

La ventaja es que si mueves todo el sistema a otro directorio o a otra máquina el código sigue funcionando mientras se respete la estructura interna de directorios.

[Víctor de Buen Remiro]

Poniéndole 5000 evaluaciones, mejora el valor objetivo de 138171533 a 144460819

[NonLinGloOpt::@PipeLine::run] step 1 NLOPT 2.3.0.LN_COBYLA found optimal value =144460819.4023014 after 5000 target evaluations in 180897.456103 seconds
Stop reason: MAXEVAL_REACHED

Tarda más de 50 horas lo cual me hace pensar que la programación es bastante ineficiente. Si me explicáis el problema podemos ver cómo hacerlo más eficaz.

[Víctor de Buen Remiro]

Comentarios sobre eficiencia de la programación

Vamos a estudiar cómo se puede mejorar el rendimiento de la siguiente función

////////////////////////////////////////////////////////////////////////////////
Real Benef(Matrix colGRPCat, Matrix colGRPCat.1)
////////////////////////////////////////////////////////////////////////////////
{
   Set G = For(1, Card(Traspose(Serie2Pairs(qt_GRP_4))[2]), Real(Real k){
        Real k
    });
   Set C = For(Card(Traspose(Serie2Pairs(qt_GRP_4))[2]) + 1, 
             Card(Traspose(Serie2Pairs(qt_GRP_4))[2])+
                    Card(Traspose(Serie2Pairs(catalogo))[2]), Real(Real k){
        Real k
    });
  
   Matrix ColGRP = SubRow(colGRPCat, G);
   Matrix ColCat   = SubRow(colGRPCat, C);

   Serie colGRP_ = Matrix2Serie(qt_GRP_4*0, Traspose(Serie2Pairs(qt_GRP_4))[1], ColGRP);
   Serie colCat_ = Matrix2Serie(catalogo*0, Traspose(Serie2Pairs(catalogo))[1], ColCat);

   Real Profit(colGRP_, colCat_)
};

Hay nada menos que cuatro llamadas a Serie2Pairs(qt_GRP_4) y 2 llamadas a Serie2Pairs(catalogo), lo cual es una forma tremendamente ineficiente de programar.

Lo lógico es crear las variables locales que sean necesarias para no repetir cálculos innecesariamente. Incluso aunque sólo se use una vez a veces es conveniente crear una variable auxiliar si con eso se mejora la legibilidad del código, o simplemente por simetría.

Por otra parte, el uso del For para crear conjuntos de secuencias numéricas se puede hacer de forma mucho más eficaz con la función

////////////////////////////////////////////////////////////////////////////////
Real Benef(Matrix colGRPCat, Matrix colGRPCat.1)
////////////////////////////////////////////////////////////////////////////////
{
   Set qt_GRP_4_pairs = Serie2Pairs(qt_GRP_4);
   Set catalogo_pairs = Serie2Pairs(catalogo);
   Set qt_GRP_4_pairs_1 = Traspose(qt_GRP_4_pairs)[1];
   Set catalogo_pairs_1 = Traspose(catalogo_pairs)[1];
   Set qt_GRP_4_pairs_2 = Traspose(qt_GRP_4_pairs)[2];
   Set catalogo_pairs_2 = Traspose(catalogo_pairs)[2];
   Real qt_GRP_4_pairs_n = Card(qt_GRP_4_pairs_2);
   Real catalogo_pairs_n = Card(catalogo_2);

   Set G = Range(1, qt_GRP_4_pairs_n, 1);
   Set C = Range(qt_GRP_4_pairs_n+1, qt_GRP_4_pairs_n+catalogo_pairs_n, 1);
  
   Matrix ColGRP = SubRow(colGRPCat, G);
   Matrix ColCat = SubRow(colGRPCat, C);

   Serie colGRP_ = Matrix2Serie(qt_GRP_4*0, qt_GRP_4_pairs_1, ColGRP);
   Serie colCat_ = Matrix2Serie(catalogo*0, catalogo_pairs_1, ColCat);

   Real Profit(colGRP_, colCat_)
};