Pyomo: многоцелевая оптимизация

Question

Я использую этот код для решения многоцелевой модели оптимизации (диспетчеризация мощности) и пытаюсь адаптировать пример в своем коде.

Пример:https://stackru.com/questions/50742999/multi-objective-optimization-example-pyomo.

И я пытаюсь пропустить часть "неэффективного фронта Парето" и напрямую построить "эффективный фронт Парето".

Первая вкладка может работать правильно и генерировать Cost_min, Cost_max, Emission_min, Emission_max.

from pyomo.environ import *
import matplotlib.pyplot as plt
import numpy as np
import matplotlib.pyplot as plt
import random

# create a model
model = AbstractModel()

# declare decision variables

model.N = Param(mutable=True)
model.J = RangeSet(model.N)
model.A = Param(model.J)
model.B = Param(model.J)
model.C = Param(model.J)
model.D = Param(model.J)
model.E = Param(model.J)
model.F = Param(model.J)
model.P_min = Param(model.J, within=PositiveReals)
model.P_max = Param(model.J, within=PositiveReals)
model.demand = Param(mutable=True)

# declare constraints

def Pbounds(model, j):
    return (model.P_min[j], model.P_max[j])
model.P = Var(model.J, bounds=Pbounds, domain=NonNegativeReals)

def P_LoadgenBalance(model):
    return sum(model.P[j] for j in model.J) >= model.demand
model.P_LoadgenBalance = Constraint(rule=P_LoadgenBalance)

# declare objective_cost

def obj_cost(model):
     return sum(model.A[j]* model.P[j] ** 2 + model.B[j] * model.P[j] + model.C[j] for j in model.J) 
model.cost= Objective(rule=obj_cost, sense=minimize)
# declare objective_emission
def obj_emission(model):
     return sum(model.E[j]* model.P[j] ** 2 + model.D[j] * model.P[j] + model.F[j] for j in model.J) 
model.emission= Objective(rule=obj_emission, sense=minimize)

# deactivate model.emission  calculate emission_max,cost_min

model.emission.deactivate()
instance = model.create_instance("E:\pycharm_project\PD\END-10units.dat")
opt = SolverFactory('Ipopt')
results = opt.solve(instance)
for i in instance.J:
    print(i,value(instance.P[i]))
print( 'cost = ' + str(value(instance.cost)) )
print( 'emission = ' + str(value(instance.emission)) )
emission_max = value(instance.emission)
cost_min = value(instance.cost)


# ## max emission  deactivate model.cost    calculate emission_min,cost_max

model.emission.activate()
model.cost.deactivate()
instance = model.create_instance("E:\pycharm_project\PD\END-10units.dat")
results = opt.solve(instance)


for i in instance.J:
    print(i,value(instance.P[i]))
print( 'cost = ' + str(value(instance.cost)) )
print( 'emission = ' + str(value(instance.emission)) )
emission_min = value(instance.emission)
cost_max = value(instance.cost)

После запуска кода на этой вкладке ошибок не было. Но при выводе фронта Парето на нем отображается только одна точка.

# ## apply normal $\epsilon$-Constraint

model.emission.deactivate()
model.cost.activate()
model.emission_value = Param(initialize=0, mutable=True)

def c_epsilon(model):
    return model.emission <= model.emission_value
model.C_epsilon = Constraint(rule=c_epsilon)
results = opt.solve(instance)

print('Each iteration will keep emission lower than some values between emission_min and emission_max, so ['       + str(emission_min) + ', ' + str(emission_max) + ']')

n = 5
step = int((emission_max - emission_min) / n)
steps = list(range(int(emission_min), int(emission_max), step)) + [emission_max]


# ## apply augmented $\epsilon$-Constraint

# max   emission + delta*epsilon <br>
#  s.t. emission - s = emission_value

model.del_component(model.cost)
model.del_component(model.emission)
model.del_component(model.C_epsilon)

model.delta = Param(initialize=0.00001)

model.s = Var(within=NonNegativeReals)

def obj_cost_1(model):
    return sum(model.cost+model.delta * model.s)
model.obj_cost_1 = Objective(rule=obj_cost_1, sense=maximize)

def C_e(model):
    return model.emission-model.s==model.emission_value
model.C_e= Constraint(rule=C_e)

cost_l = []
emission_l = []
for i in steps:
    model.emission_value = i
    results = opt.solve(instance)  
    cost_l.append(value(instance.cost))
    emission_l.append(value(instance.emission))
plt.plot(cost_l,emission_l,'o-.');
plt.title('efficient Pareto-front');
plt.grid(True);
plt.show()

Результат показан ниже. Я не знаю, почему это не может дать правильную диаграмму Парето. Я не знаю, какой шаг кода неправильный.

эффективный Парето-фронт Может ли кто-нибудь помочь мне с этим кодом? Спасибо, Виви

Answer 1

Пара вещей....:)

Что не так:

Внутри вашего цикла единственное, что влияет на модель, - это присвоение нового значения параметру model.e Это что? Я думаю, что это опечатка, и вы по ошибке просто объявляете новую и неиспользуемую переменную экземпляра компонента модели с именемe. Вот почему вы не получаете разных значений. Я думаю, ты хочешь сменить наmodel.emission.

Кроме того, я бы не стал пытаться сделать 1000 решений с первого раза, просто попробуйте 5.

Что нужно очистить:

вы создаете новый решатель в своем цикле. Не требуется. Вам не нужно 1000 разных решателей, просто решите заново. У вас уже есть заявленный ранее решатель.

добавление некоторых комментариев к вашему коду для ясности не вызовет у вас проблем и поможет в T/S, наряду с небольшой реорганизацией.

Дополнительно, model.A model.B model.C ...не очень информативно. Я бы предложил более четкие имена переменных, если это возможно.