Coverage for backend/ahuora-builder/src/ahuora_builder/diagnostics/infeasibilities.py: 100%

1from ..properties_manager import PropertiesManager 

2from typing import Iterator, TextIO 

3import sys 

4import pyomo.environ as pyo 

5from pyomo.core.base.constraint import ScalarConstraint 

6from pyomo.contrib.incidence_analysis import IncidenceGraphInterface 

7import numpy as np 

8import idaes.logger as idaeslog 

9 

10 

11_log = idaeslog.getLogger(__name__) 

12 

13 

14 

15def compute_infeasibilities(properties_map: PropertiesManager) -> Iterator[tuple[any,float]]: 

16 """ 

17 Compute how far each constrained property is from its target value. 

18  

19 :param properties_map: PropertiesManager instance containing property components 

20 :type properties_map: PropertiesManager 

21 :return: A tuple list of property IDs and their infeasibility magnitudes 

22 :rtype: list[tuple[Any, float]] 

23 """ 

24 igraph = None 

25 for id, property_component in properties_map.items(): 

26 if property_component.corresponding_constraint is not None: 

27 # Calculate the constraint infeasibility squared error: 

28 var = next(iter(property_component.corresponding_constraint)) 

29 if not isinstance(var, pyo.Var): 

30 continue 

31 if igraph is None: 

32 igraph = IncidenceGraphInterface(var.model(), include_fixed=True) 

33 try: 

34 constraints = igraph.get_adjacent_to(var) # get all the constraints this variable is used in. 

35 # Calculate the Mean Squared Error across the constraints using the slack 

36 error = np.sqrt(np.mean(np.array([constraint.slack() for constraint in constraints])**2)) 

37 yield (id, error) 

38 except Exception as e: 

39 # handle properties that don't exist 

40 _log.warning("Error computing infeasibility for property %s: %s", id, e) 

41 continue 

42 

43def get_top_infeasibilities(properties_map: PropertiesManager, threshold: float = 1e-5) -> list[tuple[any, float]]: 

44 """ 

45 Get the top infeasibilities above a certain threshold. 

46  

47 :param properties_map: PropertiesManager instance containing property components 

48 :type properties_map: PropertiesManager 

49 :param threshold: Minimum infeasibility magnitude to report 

50 :type threshold: float 

51 :return: A list of property IDs and their infeasibility magnitudes above the threshold 

52 :rtype: list[tuple[Any, float]] 

53 """ 

54 return [(id, infeasibility) for id, infeasibility in compute_infeasibilities(properties_map) if infeasibility > threshold]