Related to the C-TAEA algorithm on Pymoo

[RafaHPSUnicamp]

Hi, Cyril. How was doing?

I am doing a multiobjective optimization about minimization and maximization. My algorithm does work well in NSGA-II and NSGA-III, but there is a problem with the C-TAEA.

As an example, the f1 results are supposed to be only negative (due to code), but there are also positive points (which does not appear in the NSGA-II and NSGA-III). In addition, the restrictions also does not to be obeyed. What is wrong to the C-TAEA code? Can you help me? I already send an e-mail to Julian Blank, but he recommends that I need to talk to you about this problem. The code is at the end of this email, writted using the Pymoo.

Here what Julian answers me about my C-TAEA problem.

Dear Rafael Henrique Pinto e Silva,

As much as I would like to give you a more concrete answer, I have to refer to you to the person who as provided the implementation for C-TAEA.

You can find his contact information on the very bottom of the tutorial page:

I am not sure if the inconsistency is due to the algorithm design or the concrete implementation. Please double check with him.

Code here:

import time
start_time = time.time()

pip install -U pymoo==0.4.2

import numpy as np

from pymoo.factory import get_problem
from pymoo.optimize import minimize

Valores iniciais

#valores iniciais (em R$/t)

S_preco = 2300

#Área (em ha)
#MA - 1
#TO - 2
#PI - 3
#BA - 4
#Nova terra MA - 5
#Nova terra TO - 6
#Nova terra PI - 7
#Nova terra BA - 8

AreaMAX = [600000, 650000, 500000, 1500000, 309853, 930190, 14197, 55525]

#produtividade das terras (em t/ha)

Produtividade = [3.206, 3.322, 3.377, 3.779, 3, 3.12, 3, 3.3]

#custos

#custo do plantio da soja (em R$/t)
S_custo = [983, 966.5, 1051.5, 1040.167, 983, 966.5, 1051.5, 1040.167]

#custo do uso da terra (em R$/t)

Custo_terra = [0, 0, 0, 0, 7500, 9000, 10000, 20500]

#demanda de soja a ser atendida

S_demanda = 5056921 + 2332769 - (5056921/117887685 + 2332769/117887685)*18344777

#Definindo os valores iniciais

#AreaSoja = [0, 0, 0, 0, 0, 0, 0, 0] - não utilizada no programa

#demais preços e custos (valores constantes)

S_preco = 2120
F_preco = 1800
O_preco = 4200
B_preco = 3800

F_custo = 36.72
O_custo = 814.5
B_custo = 1000

#demanda de farelo, óleo bruto e biodiesel

F_demanda = 1510028
O_demanda = 378913
B_demanda = 244031.02

#Valores iniciais da porcentagem de soja a ser vendida e da porcentagem de óleo a ser vendido

#S_per = 0
#O_per = 0

Modelo matemático

import numpy as np
from pymoo.model.problem import Problem

class MyProblem(Problem):

    def __init__(self):
        super().__init__(n_var=10,
                         n_obj=2,
                         n_constr=12,
                         xl=np.array([0,0,0,0,0,0,0,0,0,0]),
                         xu=np.array([AreaMAX[0],AreaMAX[1],AreaMAX[2],AreaMAX[3],AreaMAX[4],AreaMAX[5],AreaMAX[6],AreaMAX[7],1,1]))

    def _evaluate(self, X, out, *args, **kwargs):
        f1 = -1*((X[:,0]*Produtividade[0] + X[:,1]*Produtividade[1] + X[:,2]*Produtividade[2] + X[:,3]*Produtividade[3] + X[:,4]*Produtividade[4] + X[:,5]*Produtividade[5] + X[:,6]*Produtividade[6] + X[:,7]*Produtividade[7])*X[:,8]*S_preco + (X[:,0]*Produtividade[0] + X[:,1]*Produtividade[1] + X[:,2]*Produtividade[2] + X[:,3]*Produtividade[3] + X[:,4]*Produtividade[4] + X[:,5]*Produtividade[5] + X[:,6]*Produtividade[6] + X[:,7]*Produtividade[7])*(1 - X[:,8])*(F_preco*0.4 + 0.2*(X[:,9]*O_preco + (1 - X[:,9])*B_preco)) - (X[:,0]*(Produtividade[0]*S_custo[0] + Custo_terra[0]) + X[:,1]*(Produtividade[1]*S_custo[1] + Custo_terra[1]) + X[:,2]*(Produtividade[2]*S_custo[2] + Custo_terra[2]) + X[:,3]*(Produtividade[3]*S_custo[3] + Custo_terra[3]) + X[:,4]*(Produtividade[4]*S_custo[4] + Custo_terra[4]) + X[:,5]*(Produtividade[5]*S_custo[5] + Custo_terra[5]) + X[:,6]*(Produtividade[6]*S_custo[6] + Custo_terra[6]) + X[:,7]*(Produtividade[7]*S_custo[7] + Custo_terra[7])) - (X[:,0]*Produtividade[0] + X[:,1]*Produtividade[1] + X[:,2]*Produtividade[2] + X[:,3]*Produtividade[3] + X[:,4]*Produtividade[4] + X[:,5]*Produtividade[5] + X[:,6]*Produtividade[6] + X[:,7]*Produtividade[7])*(1 - X[:,8])*(F_custo + O_custo + 0.2*(1 - X[:,9])*(B_custo)))/1e10
        f2 = (X[:,0] + X[:,1] + X[:,2] + X[:,3] + X[:,4] + X[:,5] + X[:,6] + X[:,7])/1e7

        g1 = X[:,0] - AreaMAX[0]
        g2 = X[:,1] - AreaMAX[1]
        g3 = X[:,2] - AreaMAX[2]
        g4 = X[:,3] - AreaMAX[3]
        g5 = X[:,4] - AreaMAX[4]
        g6 = X[:,5] - AreaMAX[5]
        g7 = X[:,6] - AreaMAX[6]
        g8 = X[:,7] - AreaMAX[7]
        g9 = -1*((X[:,0]*Produtividade[0] + X[:,1]*Produtividade[1] + X[:,2]*Produtividade[2] + X[:,3]*Produtividade[3] + X[:,4]*Produtividade[4] + X[:,5]*Produtividade[5] + X[:,6]*Produtividade[6] + X[:,7]*Produtividade[7])*X[:,8] - S_demanda)
        g10 = -1*((X[:,0]*Produtividade[0] + X[:,1]*Produtividade[1] + X[:,2]*Produtividade[2] + X[:,3]*Produtividade[3] + X[:,4]*Produtividade[4] + X[:,5]*Produtividade[5] + X[:,6]*Produtividade[6] + X[:,7]*Produtividade[7])*(1 - X[:,8])*0.4 - F_demanda)
        g11 = -1*((X[:,0]*Produtividade[0] + X[:,1]*Produtividade[1] + X[:,2]*Produtividade[2] + X[:,3]*Produtividade[3] + X[:,4]*Produtividade[4] + X[:,5]*Produtividade[5] + X[:,6]*Produtividade[6] + X[:,7]*Produtividade[7])*(1 - X[:,8])*0.2*X[:,9] - O_demanda)
        g12 = -1*((X[:,0]*Produtividade[0] + X[:,1]*Produtividade[1] + X[:,2]*Produtividade[2] + X[:,3]*Produtividade[3] + X[:,4]*Produtividade[4] + X[:,5]*Produtividade[5] + X[:,6]*Produtividade[6] + X[:,7]*Produtividade[7])*(1 - X[:,8])*0.2*(1 - X[:,9]) - B_demanda)

        out["F"] = np.column_stack([f1, f2])
        out["G"] = np.column_stack([g1, g2, g3, g4, g5 ,g6, g7, g8, g9, g10, g11, g12])

vectorized_problem = MyProblem()

#X[:,0] - AreaSoja[0]
#X[:,1] - AreaSoja[1]
#X[:,2] - AreaSoja[2]
#X[:,3] - AreaSoja[3]
#X[:,4] - AreaSoja[4]
#X[:,5] - AreaSoja[5]
#X[:,6] - AreaSoja[6]
#X[:,7] - AreaSoja[7]
#X[:,8] - S_Per (porcentagem de soja vendido)
#X[:,9] - O_Per (porcentagem de oleo vendido)

Código de execução

import numpy as np
from pymoo.util.misc import stack
from pymoo.model.problem import Problem

class MyTestProblem(MyProblem):

    def _calc_pareto_front(self, *args, **kwargs):
        return func_pf(**kwargs)

    def _calc_pareto_set(self, *args, **kwargs):
        return func_ps(**kwargs)

problem = MyProblem()

from pymoo.algorithms.ctaea import CTAEA
from pymoo.factory import get_problem, get_reference_directions, get_sampling, get_crossover, get_mutation
from pymoo.optimize import minimize
from pymoo.visualization.scatter import Scatter

ref_dirs = get_reference_directions("das-dennis", 2, n_partitions=200)

algorithm = CTAEA(
    ref_dirs=ref_dirs)

# Nova seção

from pymoo.factory import get_termination

termination = get_termination("n_gen", 1000)

res = minimize(problem,
               algorithm,
               termination,
               seed=1,
               save_history=True,
               verbose=True)

opt = res.opt[0]
X, F, CV = opt.get("X", "__F__", "__CV__")

print("Best solution found: \nX = %s\nF = %s\nCV = %s\n" % (X, F, CV))

Grafico Pareto

from pymoo.visualization.scatter import Scatter

res.F[:,0] *= -1

plot = Scatter(title = "Espaço Objetivo")
plot.add(res.F, color="red")
plot.show()

X = res.pop.get("X")
np.savetxt("CTAEAPop.csv", X, delimiter=",")

F = res.pop.get("F")
np.savetxt("CTAEAResult.csv", F, delimiter=",")

Criando a fronteira de pareto

i = 0

x = []
y = []

while i < len(res.F):
  xa = res.F[i][0]
  ya = res.F[i][1]

  x.append(xa)
  y.append(ya)

  i += 1

j = 0

xall = []
yall = []

while j < len(F):
  xa = -F[j][0]
  ya = F[j][1]

  xall.append(xa)
  yall.append(ya)

  j += 1

from matplotlib import pyplot
import matplotlib.pyplot as plt

pyplot.plot(x, y,'o',color='red')
pyplot.plot(xall, yall, '.', color='blue')
plt.title("Espaço Objetivo")
plt.xlabel("Lucro obtido (em R$×10^10)")
plt.ylabel("Área utilizada (em ha×10^7)")

Indicadores de performance

from pymoo.factory import get_performance_indicator

gd = get_performance_indicator("gd", F)
print("GD", gd.calc(res.F))

hv = get_performance_indicator("hv", ref_point=np.array([1.2, 1.2]))
print("hv", hv.calc(res.F))

print("--- %s seconds ---" % (time.time() - start_time))

Kind regards,

Rafael

[cyrilpic]

Hi Rafael,

(I have removed the issue under the hvwfg repository and I will answer here. Actually, you could have opened a discussion on https://github.com/anyoptimization/pymoo)

It is hard to answer your question. If I run the code that you have provided and try to solve vectorized_problem, C-TAEA returns 19 valid results all with negative f1 and negative constraints (as it should since pymoo assumes g_i(x) <= 0). Also g1 to g8 seem redundant to setting the upper bounds (xu = ...).

Where do you see the positive value for f1? I run some tests and found some possible X that generate positive f1 values (e.g. X = array([[ 600000, 650000, 500000, 1500000, 309853, 930190, 14197, 55525, 0.1, 1]]).

Internally, C-TAEA will keep invalid solutions, but res.opt should contain only valid individuals.

I am also curious, what were the problems with NSGA-II and NSGA-III that you switched to C-TAEA?

[RafaHPSUnicamp]

So, Cyril. I wanted to analyze the results of NSGA-II, NSGA-III and C-TAEA with Pymoo, since they are the only evolutionary algorithms possible to use on my model. I recently did the same with the NSGA-II and NSGA-III and I found great results. The problems are the same, with the only difference being that I switched the C-TAEA part with the NSGA-II and NSGA-III part.

from pymoo.algorithms.nsga2 import NSGA2 from pymoo.factory import get_sampling, get_crossover, get_mutation

algorithm = NSGA2( pop_size=1100, n_offsprings=10, sampling=get_sampling("real_random"), crossover=get_crossover("real_sbx", prob=0.9, eta=15), mutation=get_mutation("real_pm", eta=20), eliminate_duplicates=True )

from pymoo.factory import get_termination

termination = get_termination("n_gen", 1000)

Also, another issue that I am seeing that the some C-TAEA points cannot obey the restrictions. It will more clear if I check any point, and see that restrictions related to g11 and g12 are not obeyed in some points. If you export the CSV archives, manually checking you can see that. In "Planilha 1".

So, it is common for C-TAEA to ignore some restrictions?

Kind regards,

Em ter., 3 de mai. de 2022 às 05:25, Cyril Picard @.***> escreveu:

Hi Rafael,

(I have removed the issue under the hvwfg repository and I will answer here. Actually, you could have opened a discussion on https://github.com/anyoptimization/pymoo)

It is hard to answer your question. If I run the code that you have provided and try to solve vectorized_problem, C-TAEA returns 19 valid results all with negative f1 and negative constraints (as it should since pymoo assumes g_i(x) <= 0). Also g1 to g8 seem redundant to setting the upper bounds (xu = ...).

Where do you see the positive value for f1? I run some tests and found some possible X that generate positive f1 values (e.g. X = array([[ 600000, 650000, 500000, 1500000, 309853, 930190, 14197, 55525, 0.1, 1]]).

Internally, C-TAEA will keep invalid solutions, but res.opt should contain only valid individuals.

I am also curious, what were the problems with NSGA-II and NSGA-III that you switched to C-TAEA?

— Reply to this email directly, view it on GitHub https://github.com/epfl-lamd/modact/issues/1#issuecomment-1115850690, or unsubscribe https://github.com/notifications/unsubscribe-auth/ARLXSODBORVFHDJU3V5MWL3VIDPGXANCNFSM5U5CMCTQ . You are receiving this because you authored the thread.Message ID: @.***>

-- Rafael Henrique Pinto e Silva Engenheiro de Energia pela Pontifícia Universidade Católica de Minas Gerais Mestre em Planejamento em Sistemas Energéticos pela Universidade Estadual de Campinas Administrador do Conexão Na7ural http://trakto.link/conexaona7

[cyrilpic]

Just to clarify, I am not the author of C-TAEA, but I have written the python implementation based on the C code from the author. This is the original paper:

K. Li, R. Chen, G. Fu, and X. Yao, “Two-Archive Evolutionary Algorithm for Constrained Multiobjective Optimization,” IEEE Transactions on Evolutionary Computation, vol. 23, no. 2, pp. 303–315, Apr. 2019, doi: 10.1109/TEVC.2018.2855411.

It uses two populations : the convergence archive and the diversity archive. The diversity archive never considers constraints. The convergence archive is updated normally only with individuals the respect the constraints, but, as an idea, to improve convergence, it can include individuals from the diversity archive. The paper gives the whole picture.

These are "internal" populations. At the end of your optimization, you should only consider res.opt.

I believe C-TAEA can suffer from normalization problems that lead to the behavior that you are seeing (only few individuals are in res.opt. This is expected based on the algorithm. Different algorithms work differently and it seems that C-TAEA is not adapted to your problem.