esa / pagmo2

A C++ platform to perform parallel computations of optimisation tasks (global and local) via the asynchronous generalized island model.
https://esa.github.io/pagmo2/
GNU General Public License v3.0
823 stars 161 forks source link

Archipelago does not seem to work with my objective function #514

Open Andreas-SM opened 2 years ago

Andreas-SM commented 2 years ago

Hello,

I am a a novice to pygmo and I have had a few issues with the archipelago class. I have tested archipelago so far with the toy_problem from the documentation and it worked fine, but when I try to use a custom udp the population within the archipelago will not evolve. I compared the final champions from the initial and final populations and they are exactly the same. The objective function that I am using is a power cycle configuration that I am bulding with Tespy. Here is the objective function that I have been trying to optimize:

`

from tespy.networks import Network
from tespy.components import (
    Compressor, Turbine, HeatExchanger, Source, Sink, CycleCloser
    )
from tespy.connections import Connection, Bus
import numpy as np

class SimpleBraytonCycle:

    def __init__(self):

        self.nw = Network(
            fluids=['HEOS::CO2', 'INCOMP::TVP1', 'Water'],
            p_unit='bar', T_unit='C', h_unit='kJ / kg',
            iterinfo=False
            )

        # Main cycle components

        cc = CycleCloser('cycle closer')
        turb = Turbine('turbine')
        cool = HeatExchanger('cooler')
        comp = Compressor('compressor')
        hx = HeatExchanger('main heater')

        # Thermal fluid components

        dti = Source('thermal fluid source')
        dto = Sink('thermal fluid sink')

        # Cooling water components

        wi = Source('cooling water source')
        wo = Sink('cooling water sink')

        # Connections

        # Main Cycle

        cc_turb = Connection(cc, 'out1', turb, 'in1', label='Point 3')
        turb_cool = Connection(turb, 'out1', cool, 'in1', label='Point 4')
        cool_comp = Connection(cool, 'out1', comp, 'in1', label='Point 1')
        comp_hx = Connection(comp, 'out1', hx, 'in2', label='Point 2')
        hx_cc = Connection(hx, 'out2', cc, 'in1', label='Point close')

        # Thermal fluid

        dti_hx = Connection(dti, 'out1', hx, 'in1', label='Point Q_in')
        hx_dto = Connection(hx, 'out1', dto, 'in1', label='Point Q_out')

        # Cooling water

        wi_cool = Connection(wi, 'out1', cool, 'in2', label='Point W_in')
        cool_wo = Connection(cool, 'out2', wo, 'in1', label='Point W_out')

        # Add connections

        self.nw.add_conns(cc_turb, turb_cool, cool_comp, comp_hx, hx_cc,
                          wi_cool, cool_wo,
                          dti_hx, hx_dto)

        # Busses

        # Power

        power = Bus('power')
        power.add_comps(
            {'comp': turb, 'char': -1}, {'comp': comp, 'char': -1}
            )

        # Heat 

        heat = Bus('heat')
        heat.add_comps(
            {'comp': cool, 'char': 1}, {'comp': hx, 'char': -1}
            )

        self.nw.add_busses(power, heat)

        # Set parameters

        hx.set_attr(pr1=1, pr2=1)
        cool.set_attr(pr1=1, pr2=1)

        comp.set_attr(eta_s=.85)
        turb.set_attr(eta_s=.85)

    def calculate_efficiency(self, x):

        # Set main cycle parameters

        # h_3 = cpp.PropsSI('Hmass', 'T', x[1] + 273.15, 'P', x[2] * 1e5, 'CO2')/1e3

        self.nw.get_conn('Point 3').set_attr(m=x[0],
                                              T=x[1],
                                              p=x[2],
                                              fluid={'CO2': 1, 'Water': 0, 'TVP1': 0})

        self.nw.get_conn('Point 4').set_attr(p=x[3])#, T=69.03)

        # Set thermal fluid parameters

        self.nw.get_conn('Point Q_in').set_attr(m=30.25,
                                                p=10.469,
                                                h=380.65,
                                                fluid={'CO2': 0, 'Water': 0, 'TVP1': 1})

        # self.nw.get_conn('Point Q_out').set_attr(h=210.2)

        # Set cooling water parameters

        self.nw.get_conn('Point W_in').set_attr(m=x[4],
                                                p=1.01325,
                                                h=146.72,
                                                fluid={'CO2': 0, 'Water': 1, 'TVP1': 0})

        self.nw.get_conn('Point W_out').set_attr(h=178.27)

        # Solve design

        self.nw.solve('design')

        for component in self.nw.comps['object']:
            if isinstance(component, HeatExchanger):
                if component.Q.val > 0:
                    return np.nan
            elif isinstance(component, Compressor):
                if component.P.val < 0:
                    return np.nan
            elif isinstance(component, Turbine):
                if component.P.val > 0:
                    return np.nan

        if self.nw.res[-1] > 1e-3 or self.nw.lin_dep:
            return np.nan
        else:
            return self.nw.busses['power'].P.val / (-1 * self.nw.get_comp('main heater').Q.val)

    def restrictions(self):

        restrictions = {
            'T_return': abs(self.nw.get_conn('Point Q_out').h.val - 215.31),
            'positive work': self.nw.busses['power'].P.val
            }

        return restrictions`

And here is the problem function plus archipelago code:

`

  from brayton_simple import SimpleBraytonCycle
  import pygmo as pg

  class brayton_optimization:

      def fitness(self, x):
          try:
              f1 = 1 / self.model.calculate_efficiency(x)
              ce1 = -1 * max(0, 1e-3 - self.model.restrictions()['T_return'])
              ci2 = -1 * min(0, self.model.restrictions()['positive work'])
          except:
              f1 = np.inf
              ce1 = np.inf
              ci2 = np.inf
          ci1 = -1 * min(0, x[2] - x[3])

          f = f1 + 1e3*ce1 + 1e2*ci1 + 1e2*ci2

          return [f]

      def get_bounds(self):
          return([5., 50., 100., 100., 5.],
                     [200., 221.85, 300., 300., 200.])

      #def get_nic(self):
      #    return 2

      #def get_nec(self):
      #    return 1

  if __name__ == 'main':
      brayton_cycle = brayton_optimization()
      brayton_cycle.model = SimpleBraytonCycle()

      algo = pg.algorithm(pg.de(gen=500))
      prob = pg.problem(brayton_cycle)

      archi = pg.archipelago(n=5,algo=algo, prob=prob, pop_size=50)
      pop_1 = archi.get_champions_f()

      archi.evolve()
      archi.wait()

      pop_2 = archi.get_champions_f()

      print(pop_1)
      print()
      print(pop_2)

`

P.S.: I have Windows in my computer