NPV not changing even by changing product price

[Uzair1shah]

Hello,

I have been working on preprocessing area of biorefinery using BioSTEAM. I can run TEA and LCA successfully, however the results of TEA do not change when I change the price of product. In my case, S101 (feedstock) is the input of preprocessing and S114 (price defined by "dry feedstock" in price dictionary) is the product of preprocessing area. I am getting negative NPV and negative IRR, so I wanted to check at what price of dry feedstock, the preprocessing will be profitable. I used the following code "price = tea.solve_price(feed)" to get the breakeven price. But even by changing the price of output stream, the NPV and IRR value stays the same. I have been wondering if I made some error in the code or is there some bug? II went through documentation multiple times but still could not figure out this part. would really appreciate it if you can help me through it. The main code I have been working with is as follows:

import biosteam as bst
import thermosteam as tmo
import graphviz
from biosteam import main_flowsheet as F, units, Stream, settings
import math as math
import numpy as np
from biosteam.units import VibratingScreen, HammerMill, Mixer, Splitter

import functools

species = {}		#create an empty collection to hold your chemical species
  
#Biomass cellulose
species['cellulose'] = bst.Chemical('cellulose', phase='s', search_db=False, Cp=1.364, rho=1540, formula="C6H10O5", Hf=-1019000, MW =162.1424)            			#Cellulose, verified
species['cellA'] = bst.Chemical('cellA', search_db=False, default=True, phase='s', formula ="C6H10O5", Hf=-1019000, MW=162.1424)                                   	#Activated cellulose
  
# Biomass hemicellulose
species['gmsw'] = bst.Chemical('gmsw', search_db=False, default=True, phase='s', formula="C5H8O4",MW = 132.11612, Hf=-759200)                                      	#Verified
species['xyhw'] = bst.Chemical('xyhw', search_db=False, default=True, phase='s', formula="C5H8O4",MW = 132.11612, Hf=-759200)                                      	#Verified
species['hce1'] = bst.Chemical('hce1', search_db=False, phase='s', formula="C5H8O4",MW = 132.11612, Hf=-759200)                                      				#Verified
species['hce2'] = bst.Chemical('hce2', search_db=False, phase='s', formula="C5H8O4",MW = 132.11612, Hf=-759200)                                      				#Verified

# Biomass lignin
species['ligc'] = bst.Chemical('ligc', search_db=False, phase='s', formula="C15H14O4",MW=258.27376, Hf=-759390)                                    					#Verified
species['ligh'] = bst.Chemical('ligh', search_db=False, phase='s', formula="C22H28O9",MW=436.45892, Hf=-1722700)                                   					#Verified
species['ligO'] = bst.Chemical('ligO', search_db=False, phase='s', formula="C20H22O10",MW=422.38868, Hf=-1847500)                                  					#Verified
species['lig'] = bst.Chemical('lig', search_db=False, phase='s', formula="C11H12O4",MW=208.21388, Hf=-729310)                                     					#Secondary lignin intermediate, verified
species['ligcc'] = bst.Chemical('ligcc', search_db=False, phase='s', formula="C15H14O4",MW=258.27376, Hf=-759390)                                   				#Verified
species['ligoh'] = bst.Chemical('ligoh', search_db=False, phase='s', formula="C19H22O8",MW=378.37888, Hf=-1429200)                                  				#Verified

# Biomass extractives
species['tannin'] = bst.Chemical('tannin', search_db=False, phase='s', formula="C15H12O7", MW=304.25608,Hf=-1079700)                                 				#Tannic acid formula, verified
species['itann'] = bst.Chemical('itann', search_db=False, phase='s', formula="C8H4O4",MW=164.11736, Hf=-616980)                                     				#Intermediate tannin, verified
species['triglycerides'] = bst.Chemical('triglycerides', search_db=False, phase='s', formula="C76H52O46", MW=897.4168, Hf=-1546300)                          		#Triglycerides formula, update the triglycerides Hf

# Biomass inorganics
# This could be defined more explictly later
species['ash'] = bst.Chemical('CaO', phase='s', formula="CaO", Hf=5382027.245)                                         												#Ash, verify HF

# Pyrolysis solid products

species['char'] = bst.Chemical('carbon', phase='s', formula="C")        																							#Carbon, verified, properties from refs 1,2
species['hmwl'] = bst.Chemical('hmwl', search_db=False, phase='l', formula="C24H28O4", MW=380.48392, Hf=-958260)                                   					#Verified, High molecular weight lignin

# Fluidized bed media
species['sand'] = bst.Chemical('SiO2', formula="SiO2", Hf=-910860)                                            														#Sand, inert, properties from ref2

#metaplastic species
species['gco2'] = bst.Chemical('gco2',search_db=False, phase='s', formula="CO2", MW=44.01, Hf=-393510)                                            					#Trapped carbon dioxide, properties from ref2
species['gco'] = bst.Chemical('gco',search_db=False, phase='s', formula="CO",MW=28.01, Hf=-110520)                                               					#Trapped carbon monoxide, properties from ref2
species['gmeoh'] = bst.Chemical('gmeoh',search_db=False, phase='s', formula="CH3OH", MW=32.042, Hf=-238640)                                         				#Trapped methanol, verified, properties from ref2
species['gcoh2_l'] =bst.Chemical('gcoh2_l', search_db=False, phase='s', formula="CH2O",MW=30.026, Hf=-115890)                                        				#Trapped formaldehyde, properties from ref2, loose
species['gc2h6'] = bst.Chemical('gc2h6', search_db=False, phase='s', formula="C2H6", MW=30.07, Hf=-84667)                                           				#Trapped ethane, properties from ref2
species['gch4'] = bst.Chemical('gch4', search_db=False, phase='s', formula="CH4",MW=16.041, Hf=-74840)                                             					#Trapped methane, properties from ref2
species['gc2h4'] = bst.Chemical('gc2h4', search_db=False, phase='s', formula="C2H4", MW=28.052, Hf=52283)                                            				#Trapped ethylene, properties from ref2
species['gc6h5oh'] = bst.Chemical('gc6h5oh', search_db=False, phase='s', formula="C6H6O", MW=94.113, Hf=-96400)                                       				#Trapped phenol, properties from ref3
species['gcoh2_s'] = bst.Chemical('gcoh2_s', search_db=False, phase='s', formula="CH2O",MW=30.026, Hf=-115890)	                                      				#Trapped formaldehyde, properties from ref2, stiff
species['gh2'] = bst.Chemical('gh2', search_db=False, phase='s', formula="H2",MW=2.016, Hf=0)		                                              					#Trapped hydrogen, properties from ref2

#non-biomass species
species['argon'] = bst.Chemical('argon', formula="Ar")  
species['nitrogen'] = bst.Chemical('nitrogen', formula="N2")  
species['oxygen'] = bst.Chemical('oxygen', formula="O2")  

#Pyrolysis products
## Gas phase species
species['ethylene'] = bst.Chemical('ethylene', formula="C2H4", search_ID='ethylene', phase='g')                  
species['methane'] = bst.Chemical('methane', formula="CH4", search_ID='methane',  phase='g')                                   										#ref2
species['hydrogen'] = bst.Chemical('hydrogen', formula="H2", search_ID='hydrogen',  phase='g')                                     									#ref2
species['cacrbonmonoxide'] = bst.Chemical('carbon monoxide', formula="CO", search_ID='carbon monoxide', phase='g')                               					#ref2
species['carbondioxide'] = bst.Chemical('carbon dioxide', formula="CO2", search_ID='carbon dioxide', phase='g')                            							#ref2
species['formaldehyde'] = bst.Chemical('formaldehyde', formula="CH2O", search_ID='formaldehyde', phase='g')                            								#ref2

## Liquid phase species
species['acrolein'] = bst.Chemical('acrolein', formula="C3H4O", search_ID='acrolein', phase='l')      										#https://atct.anl.gov/Thermochemical%20Data/version%201.122e/species/?species_number=1209

species['sinapylaldehyde'] = bst.Chemical('sinapyaldehyde', search_db=False, phase='l', formula="C11H12O4", MW=208.21388, Hf=-483800)
species['methyllinoleate'] = bst.Chemical('methyllinoleate', phase='l', search_db=False, formula="C19H34O2", MW=294.4721, Hf=-604880)
species['glyoxal'] = bst.Chemical('glyoxal', search_ID='glyoxal', phase='l')							    														#NIST webbook
species['acetaldehyde'] = bst.Chemical('acetaldehyde', search_ID='acetaldehyde', phase='l')                       													#Acetaldehyde, properties from ref2 
#species['aceticacid'] = bst.Chemical('aceticacid', search_db=False, default=True, formula="CH3COOH", MW=60.05, Hf=-484930, Hvap=0, Psat=0, Tb=0)                   #Acetic acid, properties from ref2 
species['glycolaldehyde'] = bst.Chemical('glycolaldehyde', search_ID='glycolaldehyde', phase='l')                                    								#Source Pyrolysis of the Simplest Carbohydrate, Glycolaldehyde (CHO-CH2OH), and Glyoxal in a Heated Microreactor. Find additional sources
species['ethanol'] = bst.Chemical('ethanol', search_ID='ethanol', phase='l')                             															#Ethanol, properties from ref2 
species['propionaldehyde'] = bst.Chemical('propionaldehyde', search_ID='propionaldehyde', phase='l')                    							#Source: NIST webbook https://webbook.nist.gov/cgi/cbook.cgi?ID=C123386&Mask=2 
species['hydroxypropanal'] = bst.Chemical('hydroxypropanal', search_ID='hydroxypropanal',phase='l')  												#https://www.chemeo.com/cid/20-360-9/2%2C2-Dimethyl-3-hydroxypropionaldehyde
species['furfural'] = bst.Chemical('furfural', search_ID='furfural', phase='l')
species['xylosan'] = bst.Chemical('xylosan', phase='l', search_db=False, formula="C5H8O4", MW=132.115, Hf=-647000)                             				#University of Pittsburg Pitt Quantum Repository, DOI 10.17614/Q4SQ8QT7Q
species['levoglucosan'] = bst.Chemical('levoglucosan', search_ID='levoglucosan', phase='l')                      											#Gas phase Hf value, ref4.
species['phenol'] = bst.Chemical('phenol', search_ID='phenol', phase='l')                                 													#phenol, properties from ref5
species['hydroxymethylfurfural'] = bst.Chemical('hydroxymethylfurfural', search_db=False, formula="C6H6O3", MW=126.1100, Hf=-437400, phase='l')           	#Hydromethyl furfural DOI:10.1021/ie901012g
species['anisole'] = bst.Chemical('anisole', search_ID='anisole', phase='l')                                  												#DOI: 10.1021/jp507267f
species['pcoumarylalcohol'] = bst.Chemical('pcoumarylalcohol', search_db=False, formula="C9H10O2", MW=150.1772, Hf=-193500, phase='l')                 		#Unverified
species['formicacid'] = bst.Chemical('formicacid', search_db=False, formula="CH2O2", MW=46.0254, Hf=-431650, phase='l')                           			#ref: Heat of formation of liquid. Should it be of gas?
species['methanol'] = bst.Chemical('methanol', search_ID='methanol', phase='l')                               												#ref2
species['water'] = bst.Chemical('water', search_ID='water', phase='l')                                   													#ref2

species['propaniocacid'] = bst.Chemical('propaniocacid', search_db=False, formula="C3H6O2", Hf=-199481.24962000002, phase='l')
species['linalylpropioinate'] = bst.Chemical('linalylpropioinate', search_db=False, formula="C19H34O2", Hf=-194296.922694, phase='l')
species['vanillin'] = bst.Chemical('vanillin', search_ID='vanillin', phase='l')
species['cresol'] = bst.Chemical('cresol', search_db=False, formula="C7H8O", Hf=-34212.624394, phase='l')

# Create an empty Chemicals object and add each species to it
chemicals = tmo.Chemicals([])

for key, value in species.items():
    #print(key)
    value.default() # call the default method to populate all other properties not defined explicitly. Todo: some properties default to water's properties. Check individual chemical properties for accuracy once model is developed. 
    chemicals.append(value) # push a chemical object to back of the chemicals collection object

chemicals.compile()								# compile the loaded chemicals. Must do.
chemicals.set_synonym('water', 'moisture')		# set synonyms or common name for select species
chemicals.set_synonym('CaO', 'ash')

settings.set_thermo(chemicals, cache=False)				# Todo: read more set_thermo method.


# set raw material price (USD/kg)
price = {'biomass': 0.08, # @10% moisture content wet basis
         'water': 0.000353,
         'steam': 0.017,
         'waste': -0.33,
         'dry feedstock': 7.212,
         'gasoline': 0.756} # 2 USD/gal

#setting electructy price
settings.electricity_price = 0.1           # set electricity price to your desired rate

# set the name of your flowsheet
bst.main_flowsheet.set_flowsheet('fast_pyrolysis_biorefinery')

# define a chemical group with its composition in mass fraction here
chemicals.define_group('biomass',IDs = ['cellulose','xyhw','ligc','ligh','ligO','triglycerides','tannin','moisture','ash'], composition = [0.2347,0.1857,0.0915,0.0761,0.1290,0.0386,0.0176,0.10,0.0915], wt=True)


S101 = Stream('S101', biomass=333334.2, units='kg/hr', phase='s', price=price['biomass'])



# Vibrating Screen
S105 = Stream('S105', phase='s')
S106 = Stream('S106', phase='s')

UO104 = VibratingScreen('UO104', ins=S101, outs=(S105,S106), split=0.9)                     # Vibrating Screen


# Simulation of UO103 (Vibrating Screen)
UO104.simulate()
UO104.show(flow='kg/hr', composition = False, N=9)
r4 = UO104.results()
print(r4)
                                                                                           # Rotary Dryer
S107 = Stream('S107', phase='g') #Air
S108 = Stream('S108', phase='g') #Natural gas
S109 = Stream('S109', phase='s') #dry feedstock/solids
S110 = Stream('S110', phase='g') #hot air
S111 = Stream('S111', phase='g') #emissions
UO105 = units.DrumDryer('UO105', ins=(S105,S107,S108), outs=(S109,S110,S111),moisture_content = 0.01, split = (dict(moisture=0.9)) )
# Simulation of UO105 (Rotary dryer)
UO105.simulate()
UO105.show(flow='kg/hr', composition = False, N=10)
r5 = UO105.results()
print(r5)

S113 = Stream('S113', phase='s') #recycle stream
S114 = Stream('S114', phase='s', price=price['dry feedstock']) # Output stream of preprocessing area
UO106 = units.Mixer('UO106', ins=(S109, S113), outs=(S114))

# Hammermill
S112 = Stream('S112', phase='s') #crushed dry feedstock
S115 = Stream('S115', phase='s') 

UO107 = units.HammerMill('UO107', ins=(S114), outs=(S115))

#Vibrating Screen
S116 = Stream('S116', phase='s')
UO108 = VibratingScreen('UO108', ins=S115, outs=(S113,S116), split=0.1) 

Preprocessing = bst.main_flowsheet.create_system('Preprocessing')
Preprocessing.simulate()
Preprocessing.diagram('cluster', format='png')

UO106.show(flow='kg/hr', composition = False, N=9)
r6 = UO106.results()
print(r6)

UO107.show(flow='kg/hr', composition = False, N=9)
r7 = UO107.results()
print(r7)

UO108.show(flow='kg/hr', composition = False, N=9)
r8 = UO108.results()
print(r8)

class Preprocessing(bst.TEA):

    def __init__(self, system, IRR, duration, depreciation, income_tax,
                 operating_days, lang_factor, construction_schedule, WC_over_FCI,
                 labor_cost, fringe_benefits, property_tax,
                 property_insurance, supplies, maintenance, administration):
        # Huang et. al. does not take into account financing or startup
        # so these parameters are 0 by default
        super().__init__(system, IRR, duration, depreciation, income_tax,
                         operating_days, lang_factor, construction_schedule,
                         startup_months=0, startup_FOCfrac=0, startup_VOCfrac=0,
                         startup_salesfrac=0, finance_interest=0, finance_years=0,
                         finance_fraction=0, WC_over_FCI=WC_over_FCI)
        self.labor_cost = labor_cost
        self.fringe_benefits = fringe_benefits
        self.property_tax = property_tax
        self.property_insurance = property_insurance
        self.supplies= supplies
        self.maintenance = maintenance
        self.administration = administration


    def _DPI(self, installed_equipment_cost):
        return installed_equipment_cost


    def _TDC(self, DPI):
        return DPI

 
    def _FCI(self, TDC):
        return TDC

    
    def _FOC(self, FCI):
            return (FCI*(self.property_tax + self.property_insurance
                     + self.maintenance + self.administration)
                + self.labor_cost*(1+self.fringe_benefits+self.supplies))



sys = bst.main_flowsheet.create_system('sys')
tea =  Preprocessing(system=sys,
                  IRR=0.15,
                  duration=(2018, 2038),
                  depreciation='MACRS7',
                  income_tax=0.21, # Previously 35% in published study
                  operating_days=200,
                  lang_factor=3,
                  construction_schedule=(0.4, 0.6),
                  WC_over_FCI=0.05,
                  labor_cost=2.5e6,
                  fringe_benefits=0.4,
                  property_tax=0.11,
                  property_insurance=0.11,
                  supplies=0.20,
                  maintenance=0.01,
                  administration=0.005)

tea.show()
# tea.get_cashflow_table()

feed = bst.main_flowsheet('S101')
products = [bst.main_flowsheet('S114')]
costs = tea.production_costs(products)# USD/yr
np.round(costs / 1e6) # million USD / yr


price = tea.solve_price(feed) # USD/kg
round(price, 6)

tea.IRR = tea.solve_IRR()
tea.show()
tea.get_cashflow_table()

#LCA

system=sys
feedstock=UO104.ins[0]
F_mass=feedstock.F_mass
GWP = 'GWP 100yr'
bst.settings.define_impact_indicator(key=GWP, units='kg*CO2e')
S101.set_CF(GWP, 0.02931 * 0.30 / 0.25, basis='kg', units='kg*CO2e')
bst.settings.set_electricity_CF(GWP, 0.36, basis='kWhr', units='kg*CO2e')


GWP_preprocessing_displacement = sys.get_net_impact(key=GWP) / sys.get_mass_flow(S114)
#problem with this line
GWP_per_GGE = sys.get_property_allocated_impact(
    key=GWP, name='energy', basis='GGE', # Energy basis defaults to 'kJ'
    ignored={S106, S111}
)

GWP_preprocessing_energy = (
    GWP_per_GGE * S114.get_property('LHV', 'GGE/hr') / S114.F_mass
)

# GWP_per_USD = sys.get_property_allocated_impact(
#     key=GWP, name='revenue', basis='USD', # Revenue basis defaults to USD
# ) # kg-CO2e / USD

for s in sys.products: s.price = 50000000

GWP_per_USD = sys.get_property_allocated_impact(
    key=GWP, name='cost', basis='USD', # Revenue basis defaults to USD
) # kg-CO2e / USD

GWP_preprocessing_revenue = (
    GWP_per_USD * S114.price
) # kg-CO2e / kg-ethanol
print(
    "Cradle to biorefinery gate GWP of Preprocessing:\n "
    f"Displacement allocation: {GWP_preprocessing_displacement:.2f} [kg CO2e / kg]\n "
    f"Energy allocation: {GWP_preprocessing_energy:.2f} [kg CO2e / kg]\n "
    f"Economic allocation: {GWP_preprocessing_revenue:.2f} [kg CO2e / kg]"
)

[yoelcortes]

@Uzair1shah, I just did a quick check. What you have listed as products are not products:

assert products[0].isproduct() # This raises an AssertionError

I recommend looking at the diagram for using sys.products to identify them.

sys.products # -> [<Stream: S106>, <Stream: S110>, <Stream: S111>, <Stream: S116>]

Thanks,