Optimization of methanation with Peng-Robinson EoS and BARON solver

[Lukaboo]

Hello all,

For my current project I would like to implement a way to be able to run a process optimization for the methanation reaction.
For this I use the Generic Property Package Framework and want to use the Peng-Robinson Cubic Equations of State for non-ideal behavior.
I see from the documentation that an external function is called in the calculation of Z for this. For the final optimization I want to use the commercial solver BARON.
Unfortunately, it looks like it can't handle calling the external function. Before the optimization, BARON creates a .bar file based on all optimization variables and constraints. In this file the term _PR_proc_Z_vap appears, which leads to an ApplicationError in BARON during the optimization. I think _PR_proc_Z_vap is a reference to the external function mentioned.

Is there any way to work around this issue and be able to use BARON? Possibly by removing the function from the optimization problem and solving it beforehand?

Thank you very much for your help!

Kind regards
Luka

[andrewlee94]

Hi Luka,

The short answer is that we do not have an implementation of the Cubic EoS that does not use the external function. Whilst there are ways to implement a cubic root finder in an equation-oriented framework, they are rather temperamental and require a lot of effort to use. Further, the cubic root finder is the central part of any cubic EoS, thus it is not possible to remove it and solve beforehand - the result of this will vary as you change the conditions in the process.

So, some possible options are:

1. Use a solver that does support AMPL external functions. This will depend a bit on what you are trying to do - IPOPT is quite capable of optimizing these types of problems with external functions, but does only give a local optima.
2. Create a surrogate model for the property package you wish to use. IDAES has a number of tools for developing surrogate models, however we do currently lack a tool for automatically using these to make property packages (i.e. you would have to do that part yourself).
3. Create your own property package implementing a set of equations for the cubic root in place of the external functions. You should be able to build this off the existing property packages, but as I mentioned earlier it takes a bit of work to get this to solve.

[Lukaboo]

Thanks for your quick reply, it helped me a lot!
I have another question: I would like to calculate the pressure drop inside the PFR reactors and include it in the pressure balance. For this I need the calculated viscosity and density of the gas mixture inside the reactor. Is there a way to access this?
I can see from the documentation that there is at least the variable dens_mass_phase. But when I try to access this value at the output, I get the error message 'ScalarPort' object has no attribute 'dens_mass_phase'.
To be sure: If I set has_pressure_change=True in the Reactor model, there is no equation stored in the model to calculate the pressure drop, but I have to specify this calculation rule as a constraint for example, right?

[andrewlee94]

OK, so so the first thing to understand is that the inlet and outlets from a unit model do not actually contain any variables themselves, rather they just contain pointers to the state variables within the model. So whilst you can access things like flow rate, mole fractions, temperature and pressure through the Prots, you cannot access any other properties through them.

Instead, you need to go to the actual StateBlock in the model - for the PFR this is called pfr.control_volume.properties which is indexed by space and time. For example, the temeprature at the outlet at time t=0 is pfr.control_volume.properties[0, 1] where time is the first index and space is the second index with x=1 representing the outlet (note that the spatial domain is normailised). This is where you will be able to find things like density (note that we normally work is mole basis, so you can generally rely on a molar property being present but mass basis might be hit or miss).

Next, there is the question of what properties will actually be supported by the property package, which will depend on the specific property package you are using. The only package which supports viscosity is the generic property package framework (in which case you will need to provide a model for viscosity as we do not have one in our libraries yet) - molar density should be supported by most property packages however.

Finally, to confirm the behaviour of has_pressure_change=True - you are correct. When this option is set to True we add a new variable deltaP to the control volume (and a Reference to it in the unit model) and include this in the momentum balance. However, we do not add any constraint to define deltaP so it is up to the user to either fix this variable to a suitable value or add a constraint to describe this. N.B. It is important to note that for spatially discretized models, deltaP is defined in terms of pressure change per unit length (e.g. Pa/m) - make sure you account for this when defining deltaP (this is different from undiscretized models which define deltaP as a simple pressure difference).

[Lukaboo]

Thank you very much for your help!

[Lukaboo]

Hello, it's me again!
I'm using the adiabatic PFR reactor model and I'm having a lot of difficulties there to intialise the flowsheet reliably.
If I set a very small value for the reactor volume, there are no problems. However, if I set a very large value (--> reaching equilibrium --> high temperatures) the initialisation fails and IPOPT reports that the problem is infeasible.
I use the Generic Property Package Framework and select the PR-EOS. I wrote the reaction model using the existing examples ( see at the end).

Furthermore, I get AMPL error messages in the IPOPT log very often. For example like this:

Error evaluating constraint 60: can't evaluate pow(-8.92286e-10,0.304).
ERROR: Solver (ipopt) returned non-zero return code (1)
ERROR: See the solver log above for diagnostic information.

So the problems are in the equation for the reaction_rate. However, it is not necessarily always constraint 60, but can also affect the other self.pX expressions.
According to the error message, this expression takes a negative value. But the partial pressure cannot be negative. The state_bounds were defined like this:

"state_bounds": {#"flow_mol": (0, 100, 1000, pyunits.mol/pyunits.s),
                    "flow_mol_phase_comp": (1e-20, 30, 200, pyunits.mol/pyunits.s),
                     "temperature": (273.15, 570, 2000, pyunits.K),
                     "pressure": (5e4, 5e5, 31e5, pyunits.Pa)},
    "pressure_ref": (101325, pyunits.Pa),
    "temperature_ref": (298.15, pyunits.K),

That can't really be the problem either. So I suspect that this is a numerical error. And I suspect that this is amplified by the equation for self.keq2. For example, if I remove the equation and set self.k_eq2 = Expression(expr=1), the initialisation runs much more stable.

I suspect the problem is in the equation for self.k_eq2 as its values are strongly dependent on temperature and the methanation reaction is strongly exothermic. I was hoping to solve the problem by formulating k_rxn and k_eq as expressions, thus avoiding one additional variable and constraint each, as well as having to set 2 initial values less. But this did not achieve the desired success.

Sorry for the somewhat confused explanation and presentation, but I don't know how to explain it better.
Thank you very much for your help!

class HDAReactionParameterData(ReactionParameterBlock):
    """
    Reaction Parameter Block Class
    """

    def build(self):
        '''
        Callable method for Block construction.
        '''
        super(HDAReactionParameterData, self).build()

        self._reaction_block_class = HDAReactionBlock

        # Rate Reaction Index
        self.rate_reaction_idx = Set(initialize=["R1"])

        # Rate Reaction Stoichiometry
        self.rate_reaction_stoichiometry = {("R1", "Vap", "CO2"): -1,
                                            ("R1", "Vap", "H2"): -4,
                                            ("R1", "Vap", "CH4"): 1,
                                            ("R1", "Vap", "H2O"): 2}



        self.alpha = Param(default=0.91,
                           units=1/pyunits.Pa)
        self.gKat = Param(default=146.63E3,
                          doc="gKat")
        self.bulk_den = Param(default= 17.186E3) # 850 kg/m³; Paper: V=2.186e-5 m³, gCat = 0.375g -> 17.186kg/m³
        self.n = Param(default=0.152,
                       doc="Fit exponent")

        # Arrhenius Constant
        self.ko_ref = Param(
            default=95.43* self.bulk_den,
            doc="Rate constant ref",
            units=pyunits.mol/pyunits.Pa**(5*self.n)/pyunits.s/pyunits.m**3) #mol*bar^(-0.54)*s^(-1)*gcat*(-1)
        # Activation Energy
        self.energy_activation = Param(default=75.3e3,
                                       doc="Activation energy",
                                       units=pyunits.J/pyunits.mol)

        # Reference temperature
        self.T_ref = Param(default=(310+273.15),
                           doc="Reference temperature",
                           units=pyunits.K)

        # Heat of Reaction
        dh_rxn_dict = {"R1": -165e3}
        self.dh_rxn = Param(self.rate_reaction_idx,
                            initialize=dh_rxn_dict,
                            doc="Heat of reaction [J/mol]")

    @classmethod
    def define_metadata(cls, obj):
        obj.add_properties({
                # 'k_rxn': {'method': '_rate_constant', 'units': 'mol/Pa^0.54.s'}
                'reaction_rate': {'method':'_rxn_rate', 'units': 'mol/s'}
                # 'k_eq': {'method': '_keq_constant', 'units': '1/Pa^(2*0.076)'}
                })
        obj.add_default_units({'time': pyunits.s,
                               'length': pyunits.m,
                               'mass': pyunits.kg,
                               'amount': pyunits.mol,
                               'temperature': pyunits.K})


class _HDAReactionBlock(ReactionBlockBase):
    """
    This Class contains methods which should be applied to Reaction Blocks as a
    whole, rather than individual elements of indexed Reaction Blocks.
    """
    def initialize(blk, outlvl=idaeslog.NOTSET, **kwargs):
        '''
        Initialization routine for reaction package.

        Keyword Arguments:
            outlvl : sets output level of initialization routine

        Returns:
            None
        '''
        init_log = idaeslog.getInitLogger(blk.name, outlvl, tag="properties")
        init_log.info('Initialization Complete.')


@declare_process_block_class("HDAReactionBlock",
                             block_class=_HDAReactionBlock)
class HDAReactionBlockData(ReactionBlockDataBase):
    """
    An example reaction package for saponification of ethyl acetate
    """

    def build(self):
        """
        Callable method for Block construction
        """
        super(HDAReactionBlockData, self).build()

        # Heat of reaction - no _ref as this is the actual property
        add_object_reference(
                self,
                "dh_rxn",
                self.config.parameters.dh_rxn)


    def _rxn_rate(self):
        self.reaction_rate = Var(self.params.rate_reaction_idx,
                                 initialize=3,
                                 bounds=(0, 100),
                                 doc="Rate of reaction",
                                 units=pyunits.mol/pyunits.s/pyunits.m**3) ### ????
        try:
            def rate_rule(b, r):
                self.k_rxn2 = Expression(expr=self.params.ko_ref * exp(
                    (-self.params.energy_activation / const.gas_constant / self.state_ref.temperature)))

                self.k_eq2 = Expression(expr=exp(
                    (1 / 1.987) * (56000 / (self.state_ref.temperature ** 2) + (34633 / self.state_ref.temperature)
                                   - 16.4 * log(self.state_ref.temperature)
                                   + 0.00577 * self.state_ref.temperature) + 33.165))

                self.pH2O = Expression(expr=self.state_ref.pressure*1e-5/1.01325*(
                    self.state_ref.flow_mol_phase_comp["Vap", "H2O"] / (
                    self.state_ref.flow_mol_phase_comp["Vap", "H2O"] + self.state_ref.flow_mol_phase_comp["Vap", "CH4"] +
                    self.state_ref.flow_mol_phase_comp["Vap", "CO2"] + self.state_ref.flow_mol_phase_comp["Vap", "H2"])))

                self.pH2 = Expression(expr=self.state_ref.pressure*1e-5/1.01325*(
                    self.state_ref.flow_mol_phase_comp["Vap", "H2"] / (
                    self.state_ref.flow_mol_phase_comp["Vap", "H2O"] + self.state_ref.flow_mol_phase_comp["Vap", "CH4"] +
                    self.state_ref.flow_mol_phase_comp["Vap", "CO2"] + self.state_ref.flow_mol_phase_comp["Vap", "H2"])))

                self.pCH4 = Expression(expr=self.state_ref.pressure*1e-5/1.01325*(
                    self.state_ref.flow_mol_phase_comp["Vap", "CH4"] / (
                    self.state_ref.flow_mol_phase_comp["Vap", "H2O"] + self.state_ref.flow_mol_phase_comp["Vap", "CH4"] +
                    self.state_ref.flow_mol_phase_comp["Vap", "CO2"] + self.state_ref.flow_mol_phase_comp["Vap", "H2"])))

                self.pCO2 = Expression(expr=self.state_ref.pressure*1e-5/1.01325*(
                    self.state_ref.flow_mol_phase_comp["Vap", "CO2"] / (
                    self.state_ref.flow_mol_phase_comp["Vap", "H2O"] + self.state_ref.flow_mol_phase_comp["Vap", "CH4"] +
                    self.state_ref.flow_mol_phase_comp["Vap", "CO2"] + self.state_ref.flow_mol_phase_comp["Vap", "H2"])))

                return b.reaction_rate[r] == (
                            self.k_rxn2/(1+self.params.alpha*self.pH2O) *
                           ((self.pH2)**(4*self.params.n) *(self.pCO2)**self.params.n
                           -  (((self.pCH4)**(self.params.n)*(self.pH2O)**(2*self.params.n))/self.k_eq2**self.params.n)))
            self.rate_expression = Constraint(self.params.rate_reaction_idx,
                                              rule=rate_rule)
        except AttributeError:
            # If constraint fails, clean up so that DAE can try again later
            self.del_component(self.reaction_rate)
            self.del_component(self.rate_expression)
            raise







    def get_reaction_rate_basis(b):
        return MaterialFlowBasis.molar
    def model_check(blk):
        pass

[andrewlee94]

First, a fundamental concept; there is a reason that we treat equilibrium reactions differently from rate based reactions. Equilibrium reactions by nature have zero rate and run into numerical issues like I mentioned, thus we need to treat them differently in the formulation. Whilst on paper you can drive a rate based reaction to equilibrium in this fashion, it becomes numerically challenging, and this is why we generally separate those reactions that will go to equilibrium (or near enough) from those that don't. I.e., if you expect a reaction to go to equilibrium (o close to it), it is far better to model it using an equilibrium constant than to drive a rate based reaction to zero (one caveat however - equilibrium based reactions cannot be used to determine the required volume for the reactor).

On the subject of finite element, the ideal length of a finite element is completely problem dependent, and it related to the time constants of all the phenomena in the system and thus will vary with things like rate constant, temperature and concentration. This has two implications:

1. At the highest level, this means that if the time constants remain the same, then the length of the finite elements should remain constant. I.e. as you increase the length of the domain you should increase the number of finite elements proportionally.
2. However, in reactive systems the time constants are changing - notably as you rive the reaction closer to equilibrium the rate gets slower. This means that the length of the finite elements should ideally vary along the length of the reactor (smaller at the beginning where rates are fast and getting larger as you get closer to equilibrium. However, this has challenges - notably that if you were to do a dynamic study then the rate is changing in time (and thus the domain would ideally change in time). This quickly gets to be impractical to implement.

It should also be noted that there is generally an upper and lower limit on the length of a finite element that can be feasibly solved. If the finite elements are too small, then you will get inaccurate results (as the length of the finite element determines the accuracy of the approximation of the partial derivatives). However, if the finite elements are too small the problem becomes "stiff" as cannot be solved - basically the partial derivatives get too small and numerical errors begin to dominate and you can no longer converge the problem.

Adding bounds to the reaction rate will probably not help you here - adding bounds only reduces the feasible space of the problem, and you are already having issues with feasibility.

Finally, one thing I can tell from those solver failures, notably the first one. Those two errors generally do not mean the overall problem is infeasible, but rather that IPOPT went into restoration mode and could not get out of it. The number one cause of that problem is poor model scaling (which is not surprising given you are trying to go to equilibrium). You should take a look at the IDAES scaling tools and try scaling the constraint in your problem so that the body on each side of the constraints is between 1 and 100.

[Lukaboo]

Hi Andrew,

I have now had some time to look into the subject again. I checked the bounds once for the simple case of only one adiabatic PFR reactor using m.display(). I couldn't see any abnormalities for the variables. However, the mole_phase_comp_eq constraint seems to be violated with regard to CO2. Do you think this could be the reason for the infeasibility? That is only the case for the first element [0,0] in the reactor. In the following elements, values in the range of 1e-10 result.

Thanks again!

Kind regards
Luka

[Lukaboo]

After further testing with only one adiabatic reactor and variation of pressure, temperature, volume, length, finite elements, I can well understand your last answer. If I change the pressure and temperature at the inlet and adjust finite elements and volume/length accordingly, I can produce an optimal result for many cases. However, as soon as I change any of this only slightly, it quickly becomes infeasible.
Later, I would like to use my model for operational optimisation at different partial load points (0 - 10 MW CH4 ~ 0-40 mol/s H2, 0-10 mol/s CH4). Pressure and temperature at the inlet should be degrees of freedom. According to the latest findings, I wonder whether this is at all reliably possible without changing the finite elements each time by trial and error. The reactor geometry should be fixed at the end.

[andrewlee94]

This is a common problem with (P)DAE systems unfortunately, and is in part because you are pushing the system close to equilibrium. Due to the approach to equilibrium, the rate of reaction at the two ends of the reactor are very different - you have a large rate of reaction (and thus extent of reaction) at the inlet, but very small near the outlet. This in turns means that the time constants for the reaction of the problem are very different at the two ends as well.

The challenge with solving (P)DAE systems is that at each point in space and time there are normally a number of completing time constants (e.g. bulk flow, rate of reaction, etc.). These time constants tie into the choice of finite element size as they basically describe the magnitude of the accumulation terms (the numerical approximation to the integral of the partial derivatives). At a high level, you can thing of it as finite element size/time constant ~ accumulation. E.g. If you have a small time constant (fast phenomena) and large finite element the accumulation terms end up very big which results in very large Jacobian terms and thus a problem that is very sensitive to small changes. On the flip side, very large time constant and small finite elements results in very small accumulation terms that can be difficult to solve to within the desired numerical tolerance.

To make things even more complicated, in most cases the time constants of the problem are not well balanced - some might be very large and some might be very small, and thus you find yourself caught where there is no feasible size for the finite elements (too small and the large time constant part cannot be solved, too large and the large time constant part cannot be solved). In these cases it is generally necessary to try to simplify the problem and eliminate the least significant time constants (e.g. assuming fast phenomena (small time constant) terms are always at equilibrium - in this case that would mean replacing the rate based reaction terms with equilibrium ones).

So, there are a few things you can try. The first would be a non-uniform discretization across the reactor, with a smaller size near the inlet and a larger size near the outlet. You can use the length_domain_set to pass a set of custom points for the finite elements to the PFR to do this (note that these should be normalized to be between 0 and 1). You might find this helps by balancing things out better. The second approach is to ask yourself if you ever think the system will not reach equilibrium? If the answer is no, it will always get close to equilibrium, then an EquilibriumReactor (or StoichiometricReactor with conversion near 1) might be a better choice (although the volume of the reactor is not defined in those cases). However, it will be far easier to solve.

[Lukaboo]

Thank you again Andrew! This helps me a lot to understand the cause of my problems and to find possible solutions. I have a mechanical engineering background, so it is difficult for me to find such causes in the details of chemical reactions.

In fact, in my case, the reaction is rather slow. Especially in the range of 250-350°C. The reason I looked at equilibrium is to get a reference point for the design of the reactors. And also, I need 95% CH4 in the product, so I need to let the reaction progress very far.
If I use adiabatic reactors for this reaction, I have a small rate of reaction at the inlet (low temperature), a large rate of reaction in the middle (high temperature due to exothermic reaction) and a small rate of reaction again at the end because equilibrium is reached.
Based on your description, I assume that this case is very unpleasant.

The reason why I am trying to avoid an equilibrium reactor is that I hope for different plant behaviour in the partial load range. So that, for example, the cooling demand or power demand is not strictly linearly dependent on the amount of product gas and that, in terms of costs, there are operating points that are to be preferred.