This is a package for configuring Reaktoro as a gray box model in Pyomo, IDAES-PSE, and WaterTAP modeling libraries. This package is not meant to replace or act as a higher level API for Reaktoro - it is only meant to enable setting up Reaktoro equilibrium problems as blocks on Pyomo models and automate transferring Reaktoro data into Pyomo variables.
The user must familiarize thyself with Reaktoro options and usage, especially when selecting Reaktoro provided databases, database files, and activity models for aqueous, solid, and gas phases. This package does not automatically select any of these options, and will use ideal models by default.
- Please refer here for information on Reaktoro supported databases (all are supported)
- Please refer here for information on Reaktoro supported activity models (all are supported, included chain operations, or passing in pre-configured activity models)
By default the package will use:
- Database: PhreeqcDatabase
- Data file: pitzer.dat
- Aqueous activity models: ActivityModelIdealAqueous
- Gas activity models: ActivityModelIdealGas
- Solid activity models: ActivityModelIdealSolution
- Ion exchange activity modells: ActivityModelIonExchange
The Reaktoro blocks built by this package are designed to solve an equilibrium problem using user provided apparent species or true species, temperature, pressure, and pH, which are broken down to base elements and equilibrated within Reaktoro to provide exact speciation and equilibrium state. Using this state the block can return various information supported by Reaktoro:
- Converted Property Types - used by default as they provide exact derivatives for common properties:
- vaporPressure - vapor pressure for specie in Pa
- osmoticPressure - osmotic pressure for specie in Pa
- elementAmount - amount of element in system in mol
- charge - solution charge balance
- alkalinityAsCaCO3 - solution alkalinity as CaCO3
- scalingTendencySaturationIndex - Saturation index of a phase
- scalingTendency - scaling tendencies of a phase (10^saturation index)
- pH - solution pH
If a property is not available in converted property types, reaktoro-pse can access any property from reaktoro which wil use numerical derivatives when exact derivatives are not available:
- Chemical properties
- Aqueous properties
- Pyomo build properties, which are custom properties built in Pyomo that use chemical properties or aqueous properties as inputs
Note: Only Reaktoro properties that return a single floating point or real value are supported, as they will be directly matched to a Pyomo Var. Reaktoro functions that return arrays or strings are not supported.
Currently, repo includes several tutorials and examples.
Tutorials:
-
Demonstration of working with Reaktoro block that shows
- How to balance feed charge with ReaktoroBlock
- Provide exact speciation
- Build reaktoro block with speciation_block option
-
Demonstration add ReaktoroBlock to 1D Reverse Osmosis model
- How to add indexed ReaktoroBlocks WaterTAP RO1D model for calculation of osmotic pressure and scaling tendencies
-
Demonstration add ReaktoroBlock to Nano filtration model
- How to add indexed ReaktoroBlocks to WaterTAP NF ZO model for calculation of osmotic pressure and scaling tendencies
-
Treatment train with Softening -> acid -> desalination process modeling with Reaktoro-PSE
- How to model softening
- How to model acid addition
- How to model scaling tendencies
- How to optimize softening and acid addition to control scaling and minimize chemical dosing costs.
Examples:
- Example of adding ReaktoroBlock to basic desalination problem that demonstrates how to:
- Setup up basic ReaktoroBlock
- Calculate basic properties for Scaling Tendency, pH, and Osmotic pressure
- Optimize system pH for operation at target Scaling Tendency
- Example of thermal precipitation that demonstrates how to:
- Configure different database from default
- Get enthalpy and vapor pressure from Reaktoro
- Setup precipitation calculation
- Setup simulation for removal of Calcite over different temperatures and estimate required energy input
- Example of ion exchange calculations that demonstrates how to:
- Set up ReaktoroBlock for charge neutralizing the feed composition
- Use outputs from speciation block as inputs into a second property block
- Add Ion Exchange phase and species into ReaktoroBlock
- Optimize addition of acid and bases for maximizing Calcium removal selectivity over Magnesium
- Example of biogas combustion that demonstrates how to:
- Set up ReaktoroBlock for customized database
- Use Condensed phase
Comparisons of Reaktoro-pse to PhreeqC when simulating
These comparisons further demonstrate how to setup Reaktoro-pse for each type of calculation.
- Water removal from solution (e.g evaporative processes)
- Vapor pressure calculation at different temperatures
- Precipitation of mineral phases
- Mixing of two solution
To-date, this has only been tested with cyipopt - other solvers have not been tested.
- For accessing other HSL linear solvers beyond MUMPS (such as MA27, MA57, etc. solver common to WaterTAP and IDAES-PSE) follow these instructions for adding it to cyipopt.
In some cases Ipopt might struggle to close unscaled dual infeasibility even when primal is reduced to <1e-8. This will be typically seen in the solver trace, that shows inf_pr being reduced to <1e-8 and inf_du stops being reduced and solver does not terminate or terminates with Solved to Acceptable level. This is a result of using approximated hessian in the GrayBox model causing issues in calculations of dual infeasibility error.
A. Set Ipopt solver option "recalc_y" to "yes"
This option will force Ipopt to use least squares method to calculate dual infeasibility, potentially improving accuracy of its estimates and reduce it to below tolerance level. Details on the recalc_y and recalc_y_feas_tol.
cy_solver = get_solver(solver="cyipopt-watertap")
cy_solver.options['recalc_y']='yes'
cy_solver.options["recalc_y_feas_tol"] = 1e-2
B. Use exact derivatives instead of numeric
The numeric derivatives carry additional errors that reduce accuracy in estimates of dual infeasibility. You can check which outputs in your Reaktoro block are exact, calculated, or numeric by using your_reaktoro_block.display_jacobian_outputs().
If option "A" did not work, using exact derivatives can potentially solve this issue. This can be accomplished by using properties with exact derivatives listed in JacobianRows class. These properties can be used to write Pyomo constraints that calculate the desired property. These derivatvies can be used in two ways:
- through use of ConvertedPropTypes, here we apply chain rule to calculate exact derivatives for desired function
- through use of PyomoProperties, where we pass outputs with exact derivatives to a Pyomo constraint.
Examples and available properties can be found in /src/reaktoro_pse/core/reaktoro_outputs.py
In some cases you might experience a failed solve with error
EXIT: Iterates diverging; problem might be unbounded.
This, in general, is a result of bad scaling. If your model solves with out a ReaktoroBlock, then the likely culprit is autoscaling for Jacobian (and possibly variables as well). In this case its recommend to provide manual scaling for the jacobian via user scaling option when building Reaktoro, for example as follows (check thermal_precipitation.py example):
jacobian_options={
"user_scaling": {
("molarEnthalpy", None): 1,
("specificHeatCapacityConstP", None): 1,
},
}
You can check the jacobian scaling by calling:
your_reaktoro_block.display_jacobian_scaling()
An alternative is to increase the divergence tolerance:
solver.options["diverging_iterates_tol"] = 1e30
Please include a minimal example using Reaktoro for your specific feature or issue request if possible. Please also include full traceback for errors the occur.
Reaktoro-pse depends on the following packages and/or versions:
- Python 3.9 through 3.12
- Reaktoro>=2.12.3
- CyIpopt 1.4.1
- Pyomo>=6.8.0
- idaes-pse>=2.5.0
- A conda or miniforge distribution compatible with
conda-forge, e.g. Miniforge - Git (needed by setuptools_scm to set the version dynamically during installation of the Python package distribution)
conda activate $YOUR_ENV
conda install cyipopt reaktoro
pip install git+https://github.com/watertap-org/reaktoro-pse.gitgit clone https://github.com/watertap-org/reaktoro-pse.git
conda create --name reaktoro-pse-dev --yes python=3.12
conda activate reaktoro-pse-dev
conda install cyipopt reaktoro
install -e.conda activate reaktoro-pse-dev
pytest --pyargs reaktoro_pse --verboseblack .