Coverage for backend/ahuora-compounds/ahuora_property_packages/helmholtz/helmholtz_extended.py: 27%

1from pyomo.environ import ( 

2 Expression, 

3 SolverFactory, 

4 check_optimal_termination, 

5 units, 

6 Block, 

7 Constraint, 

8 Var 

9) 

10from idaes.models.properties.general_helmholtz.helmholtz_state import HelmholtzStateBlockData, _StateBlock 

11from idaes.models.properties.general_helmholtz.helmholtz_functions import HelmholtzParameterBlockData 

12from idaes.core import declare_process_block_class 

13from idaes.core.util.initialization import solve_indexed_blocks, revert_state_vars 

14from idaes.core.util.model_statistics import degrees_of_freedom 

15from idaes.core.util.exceptions import InitializationError 

16import idaes.logger as idaeslog 

17 

18from ahuora_property_packages.base.flexible_state_block import FlexibleStateBlockData, FlexibleStateBlock 

19from ahuora_property_packages.utils.fix_state_vars import fix_state_vars 

20 

21 

22class _ExtendedStateBlock(FlexibleStateBlock,_StateBlock): 

23 

24 def initialize(blk, *args, **kwargs): 

25 #blk.pre_initialize(*args, **kwargs) 

26 # In this case we actually don't actually need to do any of the FlexibleStateBlock pre and post initialisation steps, because we have already 

27 # written a custom initialize method that handles the extra constraints for helmholtz.  

28 outlvl = kwargs.get("outlvl", idaeslog.NOTSET) 

29 init_log = idaeslog.getInitLogger(blk.name, outlvl, tag="properties") 

30 solve_log = idaeslog.getSolveLogger(blk.name, outlvl, tag="properties") 

31 

32 set_vapor_frac_guesses(blk) 

33 

34 flag_dict = fix_state_vars(blk, kwargs.get("state_args", None)) 

35 

36 dof = degrees_of_freedom(blk) 

37 if dof != 0: 

38 raise InitializationError( 

39 f"{blk.name} Unexpected degrees of freedom during " 

40 f"initialization at property initialization step: {dof}." 

41 ) 

42 

43 res = None 

44 opt = SolverFactory('ipopt') 

45 with idaeslog.solver_log(solve_log, idaeslog.DEBUG) as slc: 

46 try: 

47 res = solve_indexed_blocks(opt, [blk], tee=slc.tee) 

48 except ValueError as e: 

49 if str(e).startswith("No variables appear"): 

50 # https://github.com/Pyomo/pyomo/pull/3445 

51 pass 

52 else: 

53 raise e 

54 

55 if res is not None and not check_optimal_termination(res): 

56 raise InitializationError( 

57 f"{blk.name} failed to initialize successfully. Please check " 

58 f"the output logs for more information." 

59 ) 

60 

61 if kwargs.get("hold_state") is True: 

62 return flag_dict 

63 else: 

64 blk.release_state(flag_dict) 

65 

66 init_log.info( 

67 "Property package initialization: {}.".format(idaeslog.condition(res)) 

68 ) 

69 return flag_dict 

70 

71 def release_state(blk, flags, outlvl=idaeslog.NOTSET): 

72 revert_state_vars(blk, flags) 

73 

74 

75 

76def set_vapor_frac_guesses(blk: Block) -> None: 

77 """ 

78 Vapor fraction specified 

79 valid combinations: {p, x}, {T, x} 

80 

81 We find a guess for h and p based on the specified vapor fraction and either p or T. 

82 

83 This is because the initial guess has to be within the vapor fraction region 

84 otherwise the solver has a hard time converging. 

85 """ 

86 blk[blk.index_set().first()].config.parameters 

87 

88 for sb in blk.values(): 

89 #vapor_frac_constraint is added by the platform when the user specifies a vapor fraction, 

90 # using the FlexibleStateBlock.constrain_component() method. this will run _convert_expression_to_var which will add 

91 # vapor_frac_constraint and vapor_frac_var to the state block. 

92 if not hasattr(sb, "vapor_frac_constraint"): 

93 continue 

94 

95 # we remove the vapor fraction constraint added by the platform, and add our own custom smooth constraint below.  

96 # This allows the vapor fraction to go below 0 and above 1, which helps with solver convergence as there is a smooth gradient to follow. 

97 del sb.vapor_frac_constraint # new smooth constraint will be added later instead 

98 

99 

100 # add custom vapor fraction constraint: h = h_sat_liq + x(h_sat_vap - h_sat_liq) 

101 # this is to make the vapor fraction continuous, rather than cutting off at 0 and 1. Helps with solving reliability as outside 

102 # the vapor fraction region, this will still be smooth and will not have a zero gradient. 

103 sb.add_component( 

104 "vapor_frac_constraint", 

105 Constraint(expr=sb.enth_mol == sb.enth_mol_sat_phase["Liq"] + sb.vapor_frac_var * (sb.enth_mol_sat_phase["Vap"] - sb.enth_mol_sat_phase["Liq"])) 

106 ) 

107 

108@declare_process_block_class("HelmholtzExtendedStateBlock", block_class=_ExtendedStateBlock) 

109class HelmholtzExtendedStateBlockData(HelmholtzStateBlockData, FlexibleStateBlockData): 

110 

111 def build(blk, *args): 

112 HelmholtzStateBlockData.build(blk, *args) 

113 FlexibleStateBlockData.build(blk, *args) 

114 

115 def constrain_component(blk, component: Var | Expression, value: float) -> Var | None: 

116 # We actually ignore mole_frac_comp constraints, because helmholtz is pure so they are not necessary. 

117 

118 if component.local_name.startswith("mole_frac_comp"): 

119 return None 

120 else: 

121 return super().constrain_component(component, value) 

122 

123 def add_extra_expressions(blk): 

124 super().add_extra_expressions() 

125 

126 # Generic property packages support temperature_bubble and temperature_dew, 

127 # Helmholtz is pure so these are the same, and so it only has 

128 # temperature_sat. https://idaes-pse.readthedocs.io/en/stable/reference_guides/model_libraries/generic/property_models/helmholtz.html#expressions 

129 # We add them here so that the property package interface is consistent. 

130 blk.add_component("temperature_bubble", Expression(expr=blk.temperature_sat)) 

131 blk.add_component("temperature_dew", Expression(expr=blk.temperature_sat)) 

132 

133 # For numerical stability when constraining temperature, we actually smooth out the temperature equilibrium curve.  

134 # Otherwise, the solver has difficulty solving through the vapor fraction region as there is zero gradient. 

135 

136 # rename temperature to old_temperature 

137 old_temperature = blk.temperature 

138 blk.del_component("temperature") 

139 blk.add_component("old_temperature", old_temperature) 

140 # create smooth temperature expression 

141 blk.add_component("temperature", Expression( 

142 expr=blk.old_temperature + (blk.enth_mol / (1 * units.J/units.mol))*0.000001 * units.K 

143 )) 

144 

145 

146@declare_process_block_class("HelmholtzExtendedParameterBlock") 

147class HelmholtzExtendedParameterBlockData(HelmholtzParameterBlockData): 

148 def build(self): 

149 super().build() 

150 self._state_block_class = HelmholtzExtendedStateBlock # noqa: F821