Attempting to implement multiple objective optimisation lexicographic payoff table

[David]

Hi,

Mavrotas wrote a frequently cited paper on multiple objective optimisation using what he calls the AUGMECON method. This is available on the GAMS website:

https://www.gams.com/modlib/adddocs/epscm.pdf

I'm trying to write a Pyomo version of the technique.

As I understand it, lexicographic optimisation can be used to generate a payoff table which only contains  Pareto optimal solutions. The lexicographic payoff table is used to determine the nadir values used in later stages of the AUGMECON technique. (There's a flow chart of the technique in the presentation at: https://www.chemeng.ntua.gr/gm/gmsite_eng/index_files/mavrotas_MCDA64_2006.pdf.)

The presentation and some other sources I've read says that the lexicographic payoff table is generated by adding constraints in an iterative process as the multiple objective functions  are executed:

1. max f1(x) = z1 | x∈S 􀃆 z1*
2. max f2(x) = z2 | x∈S and f1(x) = z1* 􀃆 z2*1
3. max f3(x) = z3 | x∈S and f1(x) = z1* and f2(x) = z2*1 􀃆 z3*1,2

...

I've tried to do this in the Pyomo code below for the lexicographic payoff table but I'm just getting the same results for f1, f2 and f3 as I iterate through  the three objective functions. 

This is my first Pyomo project so I'm possibly using Pyomo incorrectly. And my first attempt at implementing a lexicographic table. And I haven't used GAMS before so that implementation is a bit opaque to me.

Very grateful for any suggestions on how to get this done. If anyone has implemented Mavrotas' generation example in Pyomo that would be great. Or another multiple objective optimisation using a lexicographic table.

Thanks

David

Current code:

from pyomo.environ import *
import matplotlib.pyplot as plt

#                                   Lignite     Oil     Natural Gas     RES
# Maximum production per year (GWh) 31000       15000   22000           10000
# Cost of production (€/MWh)        30          75      60              90
# CO2 emission coefficient (t/MWh)  1.44        0.72    0.45            0


# MIN 30 LIGN + 75 OIL + 60 NG + 90 RES
# MIN 1.44 LIGN + 0.72 OIL + 0.45 NG
# MIN OIL + NG
# ST
# LIGN – LIGN1 – LIGN2 = 0
# OIL – OIL2 – OIL3 = 0
# NG – NG1 – NG2 – NG3 = 0
# RES – RES1 – RES3 = 0
# LIGN <= 31000
# OIL <= 15000
# NG <= 22000
# RES <= 10000
# LIGN1 + NG1 + RES1 >= 38400
# LIGN2 + OIL2 + NG2 >= 19200
# OIL3 + NG3 + RES3 >= 6400

model = ConcreteModel()

model.LIGN = Var(within=NonNegativeReals)
model.LIGN1 = Var(within=NonNegativeReals)
model.LIGN2 = Var(within=NonNegativeReals)

model.OIL = Var(within=NonNegativeReals)
model.OIL2 = Var(within=NonNegativeReals)
model.OIL3 = Var(within=NonNegativeReals)

model.NG = Var(within=NonNegativeReals)
model.NG1 = Var(within=NonNegativeReals)
model.NG2 = Var(within=NonNegativeReals)
model.NG3 = Var(within=NonNegativeReals)

model.RES = Var(within=NonNegativeReals)
model.RES1 = Var(within=NonNegativeReals)
model.RES3 = Var(within=NonNegativeReals)

model.C1 = Constraint(expr = model.LIGN - model.LIGN1 - model.LIGN2 == 0)
model.C2 = Constraint(expr = model.OIL - model.OIL2 - model.OIL3 == 0)
model.C3 = Constraint(expr = model.NG - model.NG1 - model.NG2 - model.NG3 == 0)
model.C4 = Constraint(expr = model.RES - model.RES1 - model.RES3 == 0)

model.C5 = Constraint(expr = model.LIGN <= 31000)
model.C6 = Constraint(expr = model.OIL <= 15000)
model.C7 = Constraint(expr = model.NG <= 22000)
model.C8 = Constraint(expr = model.RES <= 10000)

model.C9 = Constraint(expr = model.LIGN1 + model.NG1 + model.RES1 >= 38400)
model.C10 = Constraint(expr = model.LIGN2 + model.OIL2 + model.NG2 >= 19200)
model.C11 = Constraint(expr = model.OIL3 + model.NG3 + model.RES3 >= 6400)

model.f1 = Var()
model.f2 = Var()
model.f3 = Var()

model.C_f1 = Constraint(expr= model.f1 == 30 * model.LIGN + 75 * model.OIL + 60 * model.NG + 90 * model.RES)
model.C_f2 = Constraint(expr= model.f2 == 1.44 * model.LIGN + 0.72 * model.OIL + 0.45 * model.NG)
model.C_f3 = Constraint(expr= model.f3 == model.OIL + model.NG)

#model.C_f3 = Constraint(expr= model.f3 == model.LIGN + model.RES)

model.O_f1 = Objective(expr= model.f1  , sense=minimize)
model.O_f2 = Objective(expr= model.f2  , sense=minimize)
model.O_f3 = Objective(expr= model.f3  , sense=minimize)

#model.O_f3 = Objective(expr= model.f3  , sense=maximize)

model.O_f2.deactivate()
model.O_f3.deactivate()

conventionalPayoffTable = []

solvername='cbc'
solverpath_exe='...'

solver=SolverFactory(solvername,executable=solverpath_exe)

#solver = SolverFactory('cbc')
solver.solve(model);

print('Optimising f1')
print( '( LIGN, LIGN1, LIGN2 ) = ( ' + str(value(model.LIGN)) + ' , ' + str(value(model.LIGN1)) + ' , ' + str(value(model.LIGN2)) + ' )')
print( '( OIL, OIL2, OIL3 ) = ( ' + str(value(model.OIL)) + ' , ' + str(value(model.OIL2)) + ' , ' + str(value(model.OIL3)) + ' )')
print( '( NG, NG1, NG2, NG3 ) = ( ' + str(value(model.NG)) + ' , ' + str(value(model.NG1)) + ' , ' + str(value(model.NG2)) + ' , ' + str(value(model.NG3)) + ' )')
print( '( RES, RES1, RES3 ) = ( ' + str(value(model.RES)) + ' , ' + str(value(model.RES1)) + ' , ' + str(value(model.RES3)) + ' )')

print( 'f1 = ' + str(value(model.f1)) )
print( 'f2 = ' + str(value(model.f2)) )
print( 'f3 = ' + str(value(model.f3)) )

f2_min = value(model.f2)
f3_min = value(model.f3)

conventionalPayoffTable.append({'Identifier': "Optimise f1", \
                    'f1':value(model.f1), \
                    'f2':value(model.f2), \
                    'f3':value(model.f3)})

model.display()

# f2

model.O_f2.activate()
model.O_f1.deactivate()


solver.solve(model);

print('Optimising f2')
print( '( LIGN, LIGN1, LIGN2 ) = ( ' + str(value(model.LIGN)) + ' , ' + str(value(model.LIGN1)) + ' , ' + str(value(model.LIGN2)) + ' )')
print( '( OIL, OIL2, OIL3 ) = ( ' + str(value(model.OIL)) + ' , ' + str(value(model.OIL2)) + ' , ' + str(value(model.OIL3)) + ' )')
print( '( NG, NG1, NG2, NG3 ) = ( ' + str(value(model.NG)) + ' , ' + str(value(model.NG1)) + ' , ' + str(value(model.NG2)) + ' , ' + str(value(model.NG3)) + ' )')
print( '( RES, RES1, RES3 ) = ( ' + str(value(model.RES)) + ' , ' + str(value(model.RES1)) + ' , ' + str(value(model.RES3)) + ' )')

print( 'f1 = ' + str(value(model.f1)) )
print( 'f2 = ' + str(value(model.f2)) )
print( 'f3 = ' + str(value(model.f3)) )

f2_max = value(model.f2)

conventionalPayoffTable.append({'Identifier': "Optimise f2", \
                    'f1':value(model.f1), \
                    'f2':value(model.f2), \
                    'f3':value(model.f3)})

model.display()

# f3

model.O_f3.activate()
model.O_f2.deactivate()

solver.solve(model);

print('Optimising f3')
print( '( LIGN, LIGN1, LIGN2 ) = ( ' + str(value(model.LIGN)) + ' , ' + str(value(model.LIGN1)) + ' , ' + str(value(model.LIGN2)) + ' )')
print( '( OIL, OIL2, OIL3 ) = ( ' + str(value(model.OIL)) + ' , ' + str(value(model.OIL2)) + ' , ' + str(value(model.OIL3)) + ' )')
print( '( NG, NG1, NG2, NG3 ) = ( ' + str(value(model.NG)) + ' , ' + str(value(model.NG1)) + ' , ' + str(value(model.NG2)) + ' , ' + str(value(model.NG3)) + ' )')
print( '( RES, RES1, RES3 ) = ( ' + str(value(model.RES)) + ' , ' + str(value(model.RES1)) + ' , ' + str(value(model.RES3)) + ' )')

print( 'f1 = ' + str(value(model.f1)) )
print( 'f2 = ' + str(value(model.f2)) )
print( 'f3 = ' + str(value(model.f3)) )

f3_min = value(model.f3)

conventionalPayoffTable.append({'Identifier': "Optimise f3", \
                    'f1':value(model.f1), \
                    'f2':value(model.f2), \
                    'f3':value(model.f3)})

print('conventionalPayoffTable', conventionalPayoffTable)

model.display()

#1. max f1(x) = z1 | x∈S 􀃆 z1*
#2. max f2(x) = z2 | x∈S and f1(x) = z1* 􀃆 z2*1
#3. max f3(x) = z3 | x∈S and f1(x) = z1* and f2(x) = z2*1 􀃆 z3*1,2

lexicographicPayoffTable = []

model.O_f1.activate()
model.O_f2.deactivate()
model.O_f3.deactivate()

solver.solve(model);

print( 'f1 = ' + str(value(model.f1)) )
print( 'f2 = ' + str(value(model.f2)) )
print( 'f3 = ' + str(value(model.f3)) )

lexicographicPayoffTable.append({'Identifier': "Optimise f1", \
                    'f1':value(model.f1), \
                    'f2':value(model.f2), \
                    'f3':value(model.f3)})



#Need f1 value as constraint on f2 optimisation


model.O_f1.deactivate()
model.O_f2.activate()

model.f1OptimisedValue = Param(initialize=value(model.f1), mutable=False)
model.C_f1LexValue = Constraint(expr = model.f1 <= model.f1OptimisedValue)

#f1OptimisedValue = value(model.f1)
#model.C_f1LexValue = Constraint(expr = model.f1 <= f1OptimisedValue)



solver.solve(model);

print( 'f1 = ' + str(value(model.f1)) )
print( 'f2 = ' + str(value(model.f2)) )
print( 'f3 = ' + str(value(model.f3)) )

lexicographicPayoffTable.append({'Identifier': "Optimise f2", \
                    'f1':value(model.f1), \
                    'f2':value(model.f2), \
                    'f3':value(model.f3)})


model.O_f2.deactivate()
model.O_f3.activate()

#Need f2 value as constraint on f3 optimisation
model.f2OptimisedValue = Param(initialize=value(model.f2), mutable=False)
model.C_f2LexValue = Constraint(expr = model.f2 <= model.f2OptimisedValue)

solver.solve(model);

print( 'f1 = ' + str(value(model.f1)) )
print( 'f2 = ' + str(value(model.f2)) )
print( 'f3 = ' + str(value(model.f3)) )


lexicographicPayoffTable.append({'Identifier': "Optimise f3", \
                    'f1':value(model.f1), \
                    'f2':value(model.f2), \
                    'f3':value(model.f3)})

print('lexicographicPayoffTable = ', lexicographicPayoffTable)

model.display()

[David]

Sorted. Wasn't a problem with how I was using Pyomo. To calculate the lexicographic table needed to calculate each objective function with all the other non-active objectives held at their optimised value. The GAMS code referenced above demonstrates  the approach.