Hide code cell content
###############################################################################
# The Institute for the Design of Advanced Energy Systems Integrated Platform
# Framework (IDAES IP) was produced under the DOE Institute for the
# Design of Advanced Energy Systems (IDAES).
#
# Copyright (c) 2018-2023 by the software owners: The Regents of the
# University of California, through Lawrence Berkeley National Laboratory,
# National Technology & Engineering Solutions of Sandia, LLC, Carnegie Mellon
# University, West Virginia University Research Corporation, et al.
# for full copyright and license information.
###############################################################################


# Heater Unit Model with Ideal Property Package#

Author: Anuja Deshpande
Maintainer: Brandon Paul
Updated: 2023-06-01

In this tutorial, we will heat a liquid mixture of benzene-toluene using a simple heater unit model, and an ideal property package with the phases specified to liquid apriori. The inlet specifications are as follows:

• Flow Rate = 1 kmol/hr

• Mole fraction (Benzene) = 0.4

• Mole fraction (Toluene) = 0.6

• Pressure = 101325 Pa

• Temperature = 353 K

In addition to the inlet specifications, there is one additional unit level specification that needs to be set:

• Option 1: Specify the outlet temperature

• Option 2: Specify the heat duty

Therefore, in this tutorial, we will simulate the following cases:

• Case 1: Compute the heat duty (J/s) required to heat the mixture to 363 K.

• Case 2: Compute the outlet temperature of the mixture when fixing the heat duty to 2 J/s.

IDAES documentation reference for heater model: https://idaes-pse.readthedocs.io/en/stable/reference_guides/model_libraries/generic/unit_models/heater.html

## Setting up the problem in IDAES#

# Import objects from pyomo package
from pyomo.environ import ConcreteModel, SolverFactory, value
from pyomo.opt import TerminationCondition, SolverStatus

# Import the main FlowsheetBlock from IDAES. The flowsheet block will contain the unit model
from idaes.core import FlowsheetBlock

# Import idaes logger to set output levels
import idaes.logger as idaeslog

# Create the ConcreteModel and the FlowsheetBlock, and attach the flowsheet block to it.
m = ConcreteModel()

m.fs = FlowsheetBlock(
dynamic=False
)  # dynamic or ss flowsheet needs to be specified here

# Import the BTX_ideal property package to create a properties block for the flowsheet
from idaes.models.properties.activity_coeff_models import BTX_activity_coeff_VLE

# Add properties parameter block to the flowsheet with specifications
m.fs.properties = BTX_activity_coeff_VLE.BTXParameterBlock(
valid_phase="Liq", activity_coeff_model="Ideal"
)

WARNING: Params with units must be mutable.  Converting Param
'fs.properties.pressure_critical' to mutable.

WARNING: Params with units must be mutable.  Converting Param
'fs.properties.temperature_critical' to mutable.

WARNING: Params with units must be mutable.  Converting Param
'fs.properties.mw_comp' to mutable.

WARNING: Params with units must be mutable.  Converting Param
'fs.properties.dh_form' to mutable.

WARNING: Params with units must be mutable.  Converting Param
'fs.properties.ds_form' to mutable.

# Import heater unit model from the model library
from idaes.models.unit_models.heater import Heater

# Create an instance of the heater unit, attaching it to the flowsheet
# Specify that the property package to be used with the heater is the one we created earlier.
m.fs.heater = Heater(property_package=m.fs.properties)

# Import the degrees_of_freedom function from the idaes.core.util.model_statistics package
# DOF = Number of Model Variables - Number of Model Constraints
from idaes.core.util.model_statistics import degrees_of_freedom

# Call the degrees_of_freedom function, get initial DOF
DOF_initial = degrees_of_freedom(m)
print("The initial DOF is {0}".format(DOF_initial))

The initial DOF is 6

# Fix the BT stream inlet conditions
m.fs.heater.inlet.flow_mol.fix(
1 * 1000 / 3600
)  # converting to mol/s as unit basis is mol/s
m.fs.heater.inlet.mole_frac_comp[0, "benzene"].fix(0.4)
m.fs.heater.inlet.mole_frac_comp[0, "toluene"].fix(0.6)
m.fs.heater.inlet.pressure.fix(101325)  # Pa
m.fs.heater.inlet.temperature.fix(353)  # K


## Case 1: Fix Outlet Temperature#

m.fs.heater.outlet.temperature.fix(363)
# Call the degrees_of_freedom function, get final DOF
DOF_final = degrees_of_freedom(m)
print("The final DOF is {0}".format(DOF_final))

The final DOF is 0


### Flowsheet Initialization#

IDAES includes pre-written initialization routines for all unit models. The output from initialization can be set to 7 different levels depending on the details required by the user. In general, when a particular output level is set, any information at that level and above gets picked up by logger. The default level taken by the logger is INFO. More information on these levels can be found in the IDAES documentation: https://idaes-pse.readthedocs.io/en/stable/reference_guides/logging.html

# Initialize the flowsheet, and set the output at WARNING
m.fs.heater.initialize(outlvl=idaeslog.WARNING)
# From the output it can be inferred that since there are no errors or warnings encountered during initialization, nothing is displayed

---------------------------------------------------------------------------
ApplicationError                          Traceback (most recent call last)
Cell In[6], line 2
1 # Initialize the flowsheet, and set the output at WARNING
----> 2 m.fs.heater.initialize(outlvl=idaeslog.WARNING)
3 # From the output it can be inferred that since there are no errors or warnings encountered during initialization, nothing is displayed

File ~/checkouts/readthedocs.org/user_builds/idaes-examples/envs/latest/lib/python3.8/site-packages/idaes/core/base/unit_model.py:540, in UnitModelBlockData.initialize(blk, *args, **kwargs)
537     c.deactivate()
539 # Remember to collect flags for fixed vars
--> 540 flags = blk.initialize_build(*args, **kwargs)
542 # If costing block exists, activate and initialize
543 for c in init_order:

File ~/checkouts/readthedocs.org/user_builds/idaes-examples/envs/latest/lib/python3.8/site-packages/idaes/core/base/unit_model.py:589, in UnitModelBlockData.initialize_build(blk, state_args, outlvl, solver, optarg)
585 opt = get_solver(solver, optarg)
587 # ---------------------------------------------------------------------
588 # Initialize control volume block
--> 589 flags = blk.control_volume.initialize(
590     outlvl=outlvl,
591     optarg=optarg,
592     solver=solver,
593     state_args=state_args,
594 )
596 init_log.info_high("Initialization Step 1 Complete.")
598 # ---------------------------------------------------------------------
599 # Solve unit

File ~/checkouts/readthedocs.org/user_builds/idaes-examples/envs/latest/lib/python3.8/site-packages/idaes/core/base/control_volume0d.py:1490, in ControlVolume0DBlockData.initialize(blk, state_args, outlvl, optarg, solver, hold_state)
1487 init_log = idaeslog.getInitLogger(blk.name, outlvl, tag="control_volume")
1489 # Initialize state blocks
-> 1490 in_flags = blk.properties_in.initialize(
1491     outlvl=outlvl,
1492     optarg=optarg,
1493     solver=solver,
1494     hold_state=hold_state,
1495     state_args=state_args,
1496 )
1498 if state_args is None:
1499     # If no initial guesses provided, estimate values for states
1500     blk.estimate_outlet_state(always_estimate=True)

File ~/checkouts/readthedocs.org/user_builds/idaes-examples/envs/latest/lib/python3.8/site-packages/idaes/models/properties/activity_coeff_models/activity_coeff_prop_pack.py:635, in _ActivityCoeffStateBlock.initialize(blk, state_args, hold_state, state_vars_fixed, outlvl, solver, optarg)
633 # Second solve for the active constraints
634 with idaeslog.solver_log(solve_log, idaeslog.DEBUG) as slc:
--> 635     res = solve_indexed_blocks(opt, [blk], tee=slc.tee)
636 init_log.info("Initialization Step 2 {}.".format(idaeslog.condition(res)))
638 # Activate activity coefficient specific constraints

File ~/checkouts/readthedocs.org/user_builds/idaes-examples/envs/latest/lib/python3.8/site-packages/idaes/core/util/initialization.py:285, in solve_indexed_blocks(solver, blocks, **kwds)
278     tmp._ctypes[Block] = [  # pylint: disable=protected-access
279         0,
280         nBlocks - 1,
281         nBlocks,
282     ]
284     # Solve temporary Block
--> 285     results = solver.solve(tmp, **kwds)
287 finally:
288     # Clean up temporary Block contents so they are not removed when Block
289     # is garbage collected.
290     tmp._decl = {}  # pylint: disable=protected-access

File ~/checkouts/readthedocs.org/user_builds/idaes-examples/envs/latest/lib/python3.8/site-packages/pyomo/opt/base/solvers.py:534, in OptSolver.solve(self, *args, **kwds)
531 def solve(self, *args, **kwds):
532     """Solve the problem"""
--> 534     self.available(exception_flag=True)
535     #
536     # If the inputs are models, then validate that they have been
537     # constructed! Collect suffix names to try and import from solution.
538     #
539     from pyomo.core.base.block import BlockData

File ~/checkouts/readthedocs.org/user_builds/idaes-examples/envs/latest/lib/python3.8/site-packages/pyomo/opt/solver/shellcmd.py:140, in SystemCallSolver.available(self, exception_flag)
138     if exception_flag:
139         msg = "No executable found for solver '%s'"
--> 140         raise ApplicationError(msg % self.name)
141     return False
142 return True

ApplicationError: No executable found for solver 'ipopt'

# Initialize the flowsheet, and set the output at INFO_HIGH
m.fs.heater.initialize(outlvl=idaeslog.INFO_HIGH)
# At INFO_HIGH level, details of all the initialization steps are displayed

2023-11-02 10:30:10 [INFO] idaes.init.fs.heater.control_volume.properties_in: Initialization Step 1 skipped.

2023-11-02 10:30:10 [INFO] idaes.init.fs.heater.control_volume.properties_in: Initialization Step 2 optimal - Optimal Solution Found.

2023-11-02 10:30:10 [INFO] idaes.init.fs.heater.control_volume.properties_in: Initialization Step 3 optimal - Optimal Solution Found.

2023-11-02 10:30:10 [INFO] idaes.init.fs.heater.control_volume.properties_in: Initialization Step 4 optimal - Optimal Solution Found.

2023-11-02 10:30:10 [INFO] idaes.init.fs.heater.control_volume.properties_in: Initialization Step 5 optimal - Optimal Solution Found.

2023-11-02 10:30:10 [INFO] idaes.init.fs.heater.control_volume.properties_out: Initialization Step 1 skipped.

2023-11-02 10:30:10 [INFO] idaes.init.fs.heater.control_volume.properties_out: Initialization Step 2 optimal - Optimal Solution Found.

2023-11-02 10:30:10 [INFO] idaes.init.fs.heater.control_volume.properties_out: Initialization Step 3 optimal - Optimal Solution Found.

2023-11-02 10:30:10 [INFO] idaes.init.fs.heater.control_volume.properties_out: Initialization Step 4 optimal - Optimal Solution Found.

2023-11-02 10:30:10 [INFO] idaes.init.fs.heater.control_volume.properties_out: Initialization Step 5 optimal - Optimal Solution Found.

2023-11-02 10:30:10 [INFO] idaes.init.fs.heater.control_volume.properties_out: State Released.

2023-11-02 10:30:10 [INFO] idaes.init.fs.heater.control_volume: Initialization Complete

2023-11-02 10:30:10 [INFO] idaes.init.fs.heater: Initialization Step 1 Complete.

2023-11-02 10:30:10 [INFO] idaes.init.fs.heater: Initialization Step 2 optimal - Optimal Solution Found.

2023-11-02 10:30:10 [INFO] idaes.init.fs.heater.control_volume.properties_in: State Released.

2023-11-02 10:30:10 [INFO] idaes.init.fs.heater: Initialization Complete: optimal - Optimal Solution Found


### Obtaining Simulation Results#

# Solve the simulation using ipopt
# Note: If the degrees of freedom = 0, we have a square problem
opt = SolverFactory("ipopt")
solve_status = opt.solve(m, tee=True)

Ipopt 3.13.2:

******************************************************************************
This program contains Ipopt, a library for large-scale nonlinear optimization.
Ipopt is released as open source code under the Eclipse Public License (EPL).

This version of Ipopt was compiled from source code available at
https://github.com/IDAES/Ipopt as part of the Institute for the Design of
Advanced Energy Systems Process Systems Engineering Framework (IDAES PSE
Framework) Copyright (c) 2018-2019. See https://github.com/IDAES/idaes-pse.

This version of Ipopt was compiled using HSL, a collection of Fortran codes
for large-scale scientific computation.  All technical papers, sales and
publicity material resulting from use of the HSL codes within IPOPT must
contain the following acknowledgement:
HSL, a collection of Fortran codes for large-scale scientific
computation. See http://www.hsl.rl.ac.uk.
******************************************************************************

This is Ipopt version 3.13.2, running with linear solver ma27.

Number of nonzeros in equality constraint Jacobian...:       45
Number of nonzeros in inequality constraint Jacobian.:        0
Number of nonzeros in Lagrangian Hessian.............:       20

Total number of variables............................:       17
variables with only lower bounds:        2
variables with lower and upper bounds:        6
variables with only upper bounds:        0
Total number of equality constraints.................:       17
Total number of inequality constraints...............:        0
inequality constraints with only lower bounds:        0
inequality constraints with lower and upper bounds:        0
inequality constraints with only upper bounds:        0

iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
0  0.0000000e+00 0.00e+00 1.00e+00  -1.0 0.00e+00    -  0.00e+00 0.00e+00   0

Number of Iterations....: 0

(scaled)                 (unscaled)
Objective...............:   0.0000000000000000e+00    0.0000000000000000e+00
Dual infeasibility......:   0.0000000000000000e+00    0.0000000000000000e+00
Constraint violation....:   0.0000000000000000e+00    0.0000000000000000e+00
Complementarity.........:   0.0000000000000000e+00    0.0000000000000000e+00
Overall NLP error.......:   0.0000000000000000e+00    0.0000000000000000e+00

Number of objective function evaluations             = 1
Number of objective gradient evaluations             = 1
Number of equality constraint evaluations            = 1
Number of inequality constraint evaluations          = 0
Number of equality constraint Jacobian evaluations   = 1
Number of inequality constraint Jacobian evaluations = 0
Number of Lagrangian Hessian evaluations             = 0
Total CPU secs in IPOPT (w/o function evaluations)   =      0.000
Total CPU secs in NLP function evaluations           =      0.000

EXIT: Optimal Solution Found.



### View Results#

# Display Heat Duty only
m.fs.heater.heat_duty.display()

heat_duty : Size=1, Index=fs._time, ReferenceTo=fs.heater.control_volume.heat
Key : Lower : Value              : Upper : Fixed : Stale : Domain
0.0 :  None : 459.10147722222246 :  None : False : False :  Reals

# Display a readable report
m.fs.heater.report()

====================================================================================
Unit : fs.heater                                                           Time: 0.0
------------------------------------------------------------------------------------
Unit Performance

Variables:

Key       : Value  : Units : Fixed : Bounds
Heat Duty : 459.10 :  watt : False : (None, None)

------------------------------------------------------------------------------------
Stream Table
Units         Inlet     Outlet
flow_mol                mole / second    0.27778    0.27778
mole_frac_comp benzene  dimensionless    0.40000    0.40000
mole_frac_comp toluene  dimensionless    0.60000    0.60000
temperature                    kelvin     353.00     363.00
pressure                       pascal 1.0132e+05 1.0132e+05
====================================================================================


## Case 2: Fix Heat Duty#

# Fix heat duty and solve the model
m.fs.heater.outlet.temperature.unfix()
m.fs.heater.heat_duty.fix(459.10147722222354)
solve_status = opt.solve(m, tee=True)

Ipopt 3.13.2:

******************************************************************************
This program contains Ipopt, a library for large-scale nonlinear optimization.
Ipopt is released as open source code under the Eclipse Public License (EPL).

This version of Ipopt was compiled from source code available at
https://github.com/IDAES/Ipopt as part of the Institute for the Design of
Advanced Energy Systems Process Systems Engineering Framework (IDAES PSE
Framework) Copyright (c) 2018-2019. See https://github.com/IDAES/idaes-pse.

This version of Ipopt was compiled using HSL, a collection of Fortran codes
for large-scale scientific computation.  All technical papers, sales and
publicity material resulting from use of the HSL codes within IPOPT must
contain the following acknowledgement:
HSL, a collection of Fortran codes for large-scale scientific
computation. See http://www.hsl.rl.ac.uk.
******************************************************************************

This is Ipopt version 3.13.2, running with linear solver ma27.

Number of nonzeros in equality constraint Jacobian...:       46
Number of nonzeros in inequality constraint Jacobian.:        0
Number of nonzeros in Lagrangian Hessian.............:       21

Total number of variables............................:       17
variables with only lower bounds:        3
variables with lower and upper bounds:        6
variables with only upper bounds:        0
Total number of equality constraints.................:       17
Total number of inequality constraints...............:        0
inequality constraints with only lower bounds:        0
inequality constraints with lower and upper bounds:        0
inequality constraints with only upper bounds:        0

iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
0  0.0000000e+00 7.28e-12 1.00e+00  -1.0 0.00e+00    -  0.00e+00 0.00e+00   0

Number of Iterations....: 0

(scaled)                 (unscaled)
Objective...............:   0.0000000000000000e+00    0.0000000000000000e+00
Dual infeasibility......:   0.0000000000000000e+00    0.0000000000000000e+00
Constraint violation....:   4.7656100143311205e-12    7.2759576141834259e-12
Complementarity.........:   0.0000000000000000e+00    0.0000000000000000e+00
Overall NLP error.......:   4.7656100143311205e-12    7.2759576141834259e-12

Number of objective function evaluations             = 1
Number of objective gradient evaluations             = 1
Number of equality constraint evaluations            = 1
Number of inequality constraint evaluations          = 0
Number of equality constraint Jacobian evaluations   = 1
Number of inequality constraint Jacobian evaluations = 0
Number of Lagrangian Hessian evaluations             = 0
Total CPU secs in IPOPT (w/o function evaluations)   =      0.000
Total CPU secs in NLP function evaluations           =      0.000

EXIT: Optimal Solution Found.



# Display outlet temperature only
m.fs.heater.outlet.temperature.display()

_temperature_outlet_ref : Size=1, Index=fs._time, ReferenceTo=fs.heater.control_volume.properties_out[...].component('temperature')
Key : Lower : Value : Upper : Fixed : Stale : Domain
0.0 :     0 : 363.0 :  None : False : False : NonNegativeReals

# Display a readable report
m.fs.heater.report()

====================================================================================
Unit : fs.heater                                                           Time: 0.0
------------------------------------------------------------------------------------
Unit Performance

Variables:

Key       : Value  : Units : Fixed : Bounds
Heat Duty : 459.10 :  watt :  True : (None, None)

------------------------------------------------------------------------------------
Stream Table
Units         Inlet     Outlet
flow_mol                mole / second    0.27778    0.27778
mole_frac_comp benzene  dimensionless    0.40000    0.40000
mole_frac_comp toluene  dimensionless    0.60000    0.60000
temperature                    kelvin     353.00     363.00
pressure                       pascal 1.0132e+05 1.0132e+05
====================================================================================


## Plotting Q vs. Outlet Temperature#

# Heat Duty vs Outlet Temperature
import matplotlib.pyplot as plt
import numpy as np

# Unfix the heat duty from case 2
m.fs.heater.heat_duty.unfix()

# Create a list of outlet temperatures for which corresponding heat duty values need to be obtained
outlet_temp_fixed = [
91.256405 + 273.15,
90.828456 + 273.15,
86.535145 + 273.15,
89.383218 + 273.15,
93.973657 + 273.15,
85.377274 + 273.15,
92.399101 + 273.15,
94.151562 + 273.15,
87.564579 + 273.15,
88.767855 + 273.15,
]

# Fix the outlet temperature values and solve the model to obtain the heat duties
heat_duty = []
for temp in outlet_temp_fixed:
m.fs.heater.outlet.temperature.fix(temp)
solve_status = opt.solve(m)
if solve_status.solver.termination_condition == TerminationCondition.optimal:
heat_duty.append(m.fs.heater.heat_duty[0].value)

# Plotting the results

plt.figure("Q vs. Temperature")
plt.plot(outlet_temp_fixed, heat_duty, "bo")
plt.xlim(358.15, 368.15)
plt.ylim(250, 700)
plt.xlabel("Outlet Temperature (K)")
plt.ylabel("Heat Duty (W)")
plt.grid()