Resolvendo o Set Cover com OpenSolver e Pyomo

A ideia deste exercício é resolver o problema da cobertura de conjuntos, utilizando o Pyomo (http://www.pyomo.org/) e o OpenSolver (www.opensolver.org).

Link para download: https://github.com/asgunzi/SetCover.

O Set cover é o problema de cobertura de conjuntos.

Tenho uma quantidade de elementos, digamos 50.

E tenho uma quantidade de conjuntos, digamos 15.

Cada conjunto tem um preço (abaixo está como peso), e cada conjunto contém os elementos.

Algumas aplicações: imagine que há várias doenças, e uma série de vacinas. Cada vacina cobre um conjunto diferente de doenças, e cada vacina tem um custo. Minimizar o custo total, cobrindo todas as doenças.

Open Solver:

Variável de decisão binária: utilizo ou não o conjunto.

Função objetivo: se utilizo o conjunto, é o preço deste.

Restrição: devo cobrir todos os elementos.

Resultado: FO de 223 e as escolhas pintadas em rosa.

Formulação Pyomo:

Leitura de dados utilizando o módulo pandas.

    df1 = pd.read_excel(io=’SetCover_Instancia01.xlsm’, sheet_name=’Base’, usecols = “B:P”, skiprows = 1, nrows = 1,header = None)

    df2 = pd.read_excel(io=’SetCover_Instancia01.xlsm’, sheet_name=’Base’, usecols = “B:P”, skiprows = 7, header = None)

    Nsets =  len(df1.columns)

    Nelementos =  len(df2)

Crio o model com o primeiro comando abaixo. Defino os índices (conjuntos e elementos).

   #—————- Modelo Pyomo ————#     

    model = ConcreteModel()

    model.S = list(range(0,Nsets)) #Num sets

    model.E = list(range(0,Nelementos)) #Num elementos

Defino a variável de decisão com o comando Var. O nome da variável é AtivaSet, ela é binária e inicializa com o valor zero.

    model.AtivaSet = Var(model.S, within=Binary, initialize=0)

Função objetivo.

Soma do produto entre a variável de decisão e o peso de cada conjunto.

def obj_rule(model):

        return sum(model.AtivaSet[s]*pesos[s] for s in model.S)

model.obj = Objective(rule = obj_rule, sense=minimize)

Restrição.

Para cada elemento, ele deve pertencer a um conjunto, ou seja, a soma do produto da matriz de alocação e a variável de decisão devem ser maior do que 1.

    def constrAlocacao(model, e):

        return sum(model.AtivaSet[s]*Malocacao[e][s] for s in model.S) >= 1

    model.restr1 = Constraint(model.E, rule = constrAlocacao)      

Algumas configurações do solver, e depois mando resolver

    SOLVER_NAME = ‘cbc’

    TIME_LIMIT = 60*1   #Em segundos

    solver = SolverFactory(SOLVER_NAME)

    solver.options[‘seconds’] = TIME_LIMIT

    results  = solver.solve(model)

    results.write() #para mostrar o status da solucao

O resultado obtido foi o mesmo do OpenSolver.

Além do CBC, outro solver free é o GLPK.

Todos os outros são pagos. O Pyomo pode chamar CPLEX, Gurobi, se tiver a licença destes.

Toda a manipulação de entrada de dados, matrizes, e saída de dados, é Python puro.

Vantagens: é free, e muito melhor que o Open Solver para manipular grande quantidade de dados.

Desvantagens: requer um grau maior de conhecimento que o OpenSolver. Não é tão fácil de manipular ou tão poderoso quanto um AIMMS + CPLEX.

Linear programming x Heuristics x Metaheuristics

A didactic example of some optimization techniques.

See file to play a bit.

Problem: there are 3 sources of material (blocks of wood, for example) and 3 destinations (mills, for example).

cid:image001.png@01D589BE.D5DD9D80

The distances are given below.

I want to do the distribution source->destination, minimizing the total distance.

cid:image002.png@01D589BE.D5DD9D80

Method 1: Solver (linear programming).

The formulation is straightforward and it will reach optimal solution.

In this case, FO is 106.700.

2 – Method 2: Heuristics.

A heuristic is a simple and easy rule of thumb.

For example, begin from the first source, distribute it to the first destination, till there is no more volume to distribute.

It is easy, computationally cheap (only O(N) interactions), but gives a bad solution: 113.400 vs 106.700.

(take a look on the VBA code).

3 – Method 3: Meta heuristics.

This method uses the heuristics as a basis, but the rule involves probabilities.

For example, I can toss a die to decide which source and which destination I want to distribute.

Do this probabilistic distribution for the whole sets, and calculate results.

Then, repeat it several times, say T, and keep the best solution.

The more trials, the better the solution (and worst the computational time ~ T*O(N)).

Each time I run, it will give a different result.

Ex. With 4 trials, it gave a FO of 115.200, worse than the pure heuristics.

With 10 trials, it gave 107.100, better than the heuristic, but still worse than the linear programming.

With 100 trials, it reaches the optimal solution.    

Some advantages of metaheuristics against linear programming:

Do not need license (Aimms can cost US$ 50.000)

It’s possible to model non-linear constraints  

Some disadvantages of metaheuristics against linear programming:

It needs people with knowledge to deal with code

Does not guarantee the global optimum, but a good solution    

As every other tool, the application depends on several constraints and the real process. There’s a best tool for each case.      


Ideias técnicas com uma pitada de filosofia: https://ideiasesquecidas.com



Shadow price

The excel solver has a feature of sensibility analysis.

For example, we want to maximize a linear programming like this.

Run solver, and select the reports you want to see:

The sensibility will be in another spreadsheet.

What does a shadow price means?

For example, the shadow price for constraint2 is equal to 4.

It means that if I increase the right hand size of constraint2 in one unit, the objective function will increase 4 units.

Lets test it:

If I increase constraint2, from 166 to 167, and run solver again, the results will be as below.

The FO was previously 709. Now, it is 713, because the shadow price for this constraint is 4.

If shadow price is zero, it means this constraint is not restricting the problem.

One detail. If I increase too much the value of the RHS of the restriction, it can give another optimal solution.

The maximum I can increase is 14, according to the column “Permitido aumentar” in the Excel report.

Take a look in the file.