Skip to content

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

@Andreas-SM

Description

@Andreas-SM

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

Metadata

Metadata

Assignees

No one assigned

    Labels

    No labels
    No labels

    Type

    No type

    Projects

    No projects

    Milestone

    No milestone

    Relationships

    None yet

    Development

    No branches or pull requests

    Issue actions