Un algoritmo para hallar los minimales de un conjunto

El algoritmo que muestro en esta sección, se basa en el trabajo original de P.L. Yu, “Cone convexity, cone extreme points, and nondominated solutions in decision problems with multiobjectives”, published by the Journal of Optimization Theory and Applications in September 1974, Volume 14, Issue 3, pp 319–377, (see: Yu, P.L. J Optim Theory Appl (1974) 14: 319. doi:10.1007/BF00932614)

El código que se añade, en lenguaje Phyton, obtiene las soluciones minimales de un conjunto poliédrico
Uno de los objetivos en la fase de investigación de mi tesis doctoral, consiste en modificar algunos algoritmos vectoriales existentes para su aplicación a la obtención de conjuntos l-minimales de una familia de conjuntos S, en el contexto de la optimización conjuntista.

Os muestro el código fuente del caso vectorial, que creo es bastante autoexplicativo. De todas formas, en mi próxima entrada tendréis la explicación teórica completa!

Espero que os resulte interesante.

# -*- coding: utf-8 -*-
"""
Created on Mon May 1 11:03:03 2017

@author: elena
"""

from numpy import array
import numpy as np
from scipy.spatial import ConvexHull
import matplotlib.pyplot as plt
from numpy import *
from scipy.linalg import *
from pylab import *
import random
from copy import deepcopy

import itertools

from scipy.optimize import linprog

#Este programa implementa el algoritmo de YU CONJUNTISTA

##############################################################################
#Definimos las funciones que usará el programa:
##############################################################################

#1.Esta funcion dado el numero de elemento p de la familia y sus dimensiones,
#recoge los valores que le da el usuario a su matriz A. 
#En el caso del programa vectorial, p=1 siempre. Notese que los conjuntos pueden
#ser definidos de este modo dada la limitacion de uso de conjuntos poliédricos
def Matrix_A(p,m,n): 
 i=0
 j=0
 aux1=0
 aux2=0
 A_aux = array([[j for j in range (n)] for i in range (m)])
 print ("\n Inserte los elementos de A por filas\n") 
 for i in range (0,m):
 for j in range (0,n):
 A_aux[i,j]=int(input("\n Elemento : "))
 print("\n Esta es su matriz A= \n", A_aux)
 return(A_aux)
# *****************************

#2.Esta funcion dada la dimesion m de la matriz A, recoge un vector b independiente
#de esa dimension que mete el usuario para ese conjunto 
def Array_b(m): 
 print("\n Define independent elements array\n")
 b = array([aux1 for aux1 in range (m)])
 for aux2 in range (0,m):
 b[aux2] = int (input ("Enter element of b: " ))
 print("\n Este es su vector b= \n", b)
 return(b)

# *******************************
#Esta clase devuelve una lista de dos elementos. El primero es una sublista con
#todas las A de todos los conjunto. El segundo una sublista con todas las b de
#todos los conjuntos

class Poly_Ineq:
 
 def __init__(self, f):
 m=0
 n=0
 p=0

 self. Polyhedros_Ineq=[]
 self.Polyhedros_b=[]
 
 for p in range (0,f):
 print("\nDefinimos la matriz que define el conjunto A",p) 
 m = int (input ("Enter m dimension of matrix A: " )) 
 n = int (input ("Enter n dimension of matrix A: " ))
 A=Matrix_A(p,m,n)
 b=Array_b(m) 
 self.Polyhedros_Ineq.insert(p,A)#Meto todas las A en una lista
 self.Polyhedros_b.insert(p,b)#Meto todas las B en una lista 
 self.Polyhedros_Ineq.insert(p+1,'0') #Marca de final
 self.Polyhedros_b.insert(p+1,'0') #marca de final
 
 def returned(self):
 self.Polyhedro_fin=[]
 self.Polyhedro_fin.insert(0,self.Polyhedros_Ineq)
 self.Polyhedro_fin.insert (1,self.Polyhedros_b)
 return self.Polyhedro_fin
 
def defineCone ():
 
 i=0
 j=0
 
 n=int(input("\nIntroduce la dimension de espacio del cono polar. Dimension= "))
 K = array([[j for j in range (n)] for i in range (n)]) #Crea la Matriz K con los valores que quiere 
 
 print ("\n Inserte los elementos del cono polar por filas consecutivas\n")
 for i in range (0,n):
 for j in range (0,n):
 print("\n Elemento en fila",i,"columna",j)
 K[i,j]=int(input("\n Elemento : "))
 return (K)

# ************************************************

def yu_set_step3_parte2 (K):
 
 N = len(K)
 step_size =[]
 
 print()
 
 for j in range (N):

print("\nDireccion h[",j, "]:")
 num=int(input("numero de pasos (num.conjuntos corte)?="))

step_size.append(num)
 
 return step_size
 
 
def yu_step3_parte3 (step_size, K, A, b):
 """
 r_list: vector que define los pasos r para cada una de las direcciones h del cono K entre sus valores mín y máx.
 K: lista que contiene los vectores del cono K
 A,b: matriz y vector que definen nuestro poliedro (Ax<=b)
 """

Y_efi = [] # lista de puntos no-dominados del conjunto Y definido por Ax<=b
 N = len(K) # número de vectores del cono K

# Calculamos los valores mínimos y máximos de h_j·x en el conjunto para 
 # todos los h_j de K y también todos los pasos r[j][step] en r_list
 r_min_list = []
 r_max_list = []
 r_list = []
 
 for j in range(N):
 
 # LP: Calculamos el valor mínimo r_min[j]=min h[j]·x s.t. Ax<=b
 c = [(+1)*(h_j_i) for h_j_i in K[j]]
 min_res = linprog(c, A, b)
 r = min_res.get('fun')
 r_min_list.append(r)
 
 # LP: Calculamos el valor máximo r_max[j]=max h[j]·x s.t. Ax<=b
 c = [(-1)*(h_j_i) for h_j_i in K[j]]
 max_res = linprog(c, A, b)
 r = -max_res.get('fun')
 r_max_list.append(r)
 
 # Calculamos el conjunto de pasos r_list[j] par el vector r[j]
 r_list.append(linspace(r_min_list[j],r_max_list[j],step_size[j]))
 print()
 print("r_min(",j,")=", r_min_list[j],"r_max(",j,")=", r_max_list[j])
 print("r(",j,")=",r_list[j])
 
 # Calculamos los conjuntos de corte para cada vector h_j:
 # j: índice al vector del cono para el que se calcula el conjunto de corte S_j(r_red)
 # K_red: conjunto de vectores del cono K excepto el vector h_j
 # R_red: vector de todos los pasos a dar en cada dirección de K* al que se quita precisamente la r_j
 for j in range(N):
 
 print ("\nVectores de K", K)
 # Calculamos para esta dirección h_j los K_red y R_red
 K_red = list(K); del K_red[j]; print (" Vectores de K_red", K_red)
 R_red = list(r_list); del R_red[j]; print (" Vectores de R_red", R_red)
 
 # Calculamos todas las posibles combinaciones de los pasos R_red para
 # las direcciones de los vectores que han quedado en K_red
 r_combinations = list(itertools.product(*R_red))
 print("\nConjuntos de corte (pasos)", r_combinations)
 
 # Cada vector r que es un elemento de r_combinations da lugar a un conjunto de corte a explorar:
 # Cada conjunto de corte se define sobre el propio poliedro Ax<=b, ampliando las matriz A y el
 # vector b:
 # - A se amplía añadiendo los vectores h que hay en (-1)*K_red cada uno en una fila
 # - b se amplía añadiendo cada componente de -r en una nueva fila
 
 K_red_inverted = -numpy.array(K_red)
 
 A_ampliada =numpy.vstack((A, K_red_inverted))
 
 for r in r_combinations:
 
 b_ampliada = numpy.concatenate((b,-numpy.array(r)))
 
 print("\nConjunto de corte (matrices A y b ampliadas):\n A'=", 
 A_ampliada, "\n b'=",b_ampliada)
 
 # Para cada S_j se debe resolver el programa lineal de optimización:
 # max <h_j,x> s.t. A_ampliada·x<=b_ampliada 
 # (se cambia el signo de la función objetivo para que mximice y no
 # minimice)
 c = [(-1)*(h_j_i) for h_j_i in K[j]]
 res = linprog(c, A_ampliada, b_ampliada)
 
 if (res.get('success') == True):
 
 y_efi = list(res.get('x'))
 print("y_efi=",y_efi)
 
 if (y_efi not in Y_efi):Y_efi.append(y_efi)
 
 else:
 
 print("no hay x_efi para este conjunto de corte")
 print(res.get('message'))


 return(Y_efi)

####################### Main del programa ##################################

##Esta primera parte esta dedicada a definir la familia de poliedros con la que
##se trabaja por medio de su H-representacion.
#

Family=[[array([[ 0, -1],
 [ 1, 0],
 [ 0, 1],
 [-1, -1]]), 
 array([[ 0, -1],
 [-1, 0],
 [-1, -1],
 [ 1, 1],
 [ 2, 1]]), '0', '0'],
 [array([0, 1, 1, -1]), array([ 0, 0, -1, 4, 5]), '0', '0']]

#Step1.Yu conjuntos: h1*(A1,A2...), h2*(A1,A2...)
Vertex_matrix_list=[array([[ 1., 0.],
 [ 0., 0.],
 [ 0., 1.],
 [ 1., 1.]])]
 
K=[[-1., 0.],[ 0., -1.]] 
 
# Se toman datos del usuario para preparar los cutting subsets. 
step_size = yu_set_step3_parte2 (K)

# Se calcula la solcuion vectorial eficiente para todos los conjuntos en Family
for num_conjunto in range(len(Family[1])-2):
 
 A = Family[0][num_conjunto]
 b = Family[1][num_conjunto]

print ("\nAPLICAMOS YU AL POLIEDRO DEFINIDO POR:")
 print ("\nA=\n",A,"\nb=\n",b)
 
 Y_efi = yu_step3_parte3 (step_size, K, A, b)

print("\n\nLos puntos eficientes son:\n")
 for y_efi in Y_efi: print(y_efi)