Coverage for backend/ahuora-builder/src/ahuora_builder/flowsheet_manager.py: 84%
329 statements
« prev ^ index » next coverage.py v7.10.7, created at 2026-06-23 21:51 +0000
« prev ^ index » next coverage.py v7.10.7, created at 2026-06-23 21:51 +0000
1from collections import defaultdict, deque
2import re
3from pyomo.network import SequentialDecomposition, Port
4from pyomo.environ import (
5 Component,
6 ConcreteModel,
7 TransformationFactory,
8 SolverFactory,
9 Block,
10 Expression,
11 Constraint,
12 Objective,
13 minimize,
14 assert_optimal_termination,
15 ScalarVar,
16 units as pyunits,
17 TerminationCondition
18)
19from pyomo.environ import assert_optimal_termination, check_optimal_termination
20from pyomo.core.base.constraint import ScalarConstraint, IndexedConstraint
21from idaes.core import FlowsheetBlock
22from .model_statistics_optimizations import (
23 degrees_of_freedom,
24 install_model_statistics_optimizations,
25)
27install_model_statistics_optimizations()
29from idaes.core.util.model_statistics import report_statistics
30import idaes.logger as idaeslog
31from ahuora_builder_types import FlowsheetSchema
32from ahuora_builder_types.flowsheet_schema import SolvedFlowsheetSchema, UnitDiagnosticsFinding
33from .property_package_manager import PropertyPackageManager
34from .port_manager import PortManager
35from .arc_manager import ArcManager
36from .tear_manager import TearManager
37from .unit_model_manager import UnitModelManager
38from .methods.adapter_library import AdapterLibrary
39from .methods.adapter import (
40 serialize_properties_map,
41 deactivate_fixed_guesses,
42 deactivate_component,
43 deactivate_components,
44 add_corresponding_constraint,
45)
46from .methods.expression_parsing import parse_expression, ExpressionParsingError
47from .timing import start_timing
48from .methods.units_handler import get_value, get_attached_unit, check_units_equivalent
49from .methods.scaling_suffix import sanitize_scaling_suffix
50from .properties_manager import PropertiesManager
51from .custom.energy.power_property_package import PowerParameterBlock
52from .custom.energy.ac_property_package import acParameterBlock
53from .custom.energy.transformer_property_package import transformerParameterBlock
54from pyomo.core.base.units_container import units, _PyomoUnit
55from idaes.core.util.model_serializer import StoreSpec, from_json, to_json
56from .diagnostics.infeasibilities import get_top_infeasibilities
58# CSTR Imports from example flowsheet.
59# To be used as placeholders until bespoke functions can be developed.
60from .custom import hda_reaction as reaction_props
61from .custom.hda_ideal_VLE import HDAParameterBlock
62from ahuora_builder.properties_manager import PropertyComponent
63from ahuora_property_packages.build_package import build_package
64from .methods.change_detection import get_current_values, detect_changes, print_changes
65from .methods.pyomo_expressions import dimensionless_objective_expression, sum_pyomo_expressions
67_log = idaeslog.getLogger(__name__)
69# from amplpy import modules
70# Import required to allow the library to set the PATH and allow conopt to be found.
73def build_flowsheet(dynamic=False,time_set=[0]):
74 """
75 Builds a flowsheet block
76 """
77 # create the model and the flowsheet block
78 model = ConcreteModel()
79 model.fs = FlowsheetBlock(dynamic=dynamic, time_set=time_set, time_units=units.s)
80 model.fs.guess_vars = []
81 model.fs.controlled_vars = Block()
82 # Always add the power property package (as there's only one and there's no different property package types)
83 # TODO: should this be on-demand?
84 model.fs.power_pp = PowerParameterBlock()
85 model.fs.ac_pp = acParameterBlock()
86 model.fs.tr_pp = transformerParameterBlock()
87 # properties map: { id: pyomo variable or expression }
88 # used to map symbols for a sympy expression
89 model.fs.properties_map = PropertiesManager()
90 # list of component constraints to add later
91 model.fs.constraint_exprs = []
93 # Placholder property packages for CSTR reactor.
94 model.fs.BTHM_params = HDAParameterBlock()
95 #Hard-coded peng-robinson package for reactor:
96 model.fs.peng_robinson = build_package("peng-robinson",["benzene", "toluene", "hydrogen", "methane"], ["Liq","Vap"])
98 # Reaction package for the HDA reaction
99 model.fs.reaction_params = reaction_props.HDAReactionParameterBlock(
100 property_package=model.fs.peng_robinson
101 )
103 return model
106class FlowsheetManager:
107 """
108 Manages the flowsheet, including the property packages, unit models, ports, arcs, and tears
109 Includes methods to load, initialise, and solve the flowsheet
110 """
112 def __init__(self, schema: FlowsheetSchema) -> None:
113 """
114 Stores all relevant information about the flowsheet, without actually loading it
115 """
116 self.timing = start_timing()
117 self.timing.add_timing("initialise_flowsheet_manager")
119 self.model = build_flowsheet(dynamic=schema.dynamic, time_set=schema.time_set)
121 self.schema = schema
122 # Add property packages first, so that unit models can use them
123 self.property_packages = PropertyPackageManager(self)
124 # Add the port manager, so the unit models can register their ports
125 self.ports = PortManager()
126 # Add unit models
127 self.unit_models = UnitModelManager(self)
128 # Add arcs to connect the unit models together
129 self.arcs = ArcManager(self)
130 # set certain arcs as tears
131 self.tears = TearManager(self)
132 self.properties_map: PropertiesManager
134 def load(self) -> None:
135 """
136 Parses the schema and loads the model
137 """
138 self.timing.step_into("load_flowsheet")
139 # Load property packages first, so that unit models can use them
140 self.timing.add_timing("load_property_packages")
141 self.property_packages.load()
142 self.timing.step_into("load_unit_models")
143 self.unit_models.load()
144 self.timing.step_out()
145 # no need to load ports seperately, they are loaded by the unit models
146 # Load arcs to connect the unit models together
147 self.arcs.load()
148 # load any expressions
149 self.load_specs()
150 # if dynamics, apply finite difference transformation.
151 if self.schema.dynamic:
152 print("performing finite difference with", len(self.model.fs.time), "time steps")
153 TransformationFactory("dae.finite_difference").apply_to(
154 self.model.fs,
155 nfe=len(self.model.fs.time)-1, # Number of finite elements to use for discretization. We aren't adding any extra steps as our constraints dont work for that.
156 wrt=self.model.fs.time,
157 scheme="BACKWARD"
158 )
160 self.properties_map = self.model.fs.properties_map
161 self.timing.step_out()
163 def load_specs(self) -> None:
164 """
165 Loads expressions from the schema
166 """
167 fs = self.model.fs
168 specs = self.schema.expressions or []
170 ## Build dependency tree and sort specs before loading.
171 dependencies, result_ids = self.build_spec_dependency_tree(specs)
172 sorted_result_ids = self.topological_sort(dependencies, result_ids)
174 # load the specs (expressions within specifications tab)
175 for result_id in sorted_result_ids:
176 spec_config = next(spec for spec in specs if spec["id"] == result_id)
177 expression_str = spec_config["expression"]
178 try:
179 component_name = f"{spec_config['name']}_{spec_config['id']}"
180 def expression_rule(blk, time_index):
181 return parse_expression(expression_str, fs,time_index)
183 component = Expression(fs.time, rule=expression_rule)
184 fs.add_component(component_name, component)
185 fs.properties_map.add(
186 spec_config["id"], component, component.name, unknown_units=True
187 )
188 except ExpressionParsingError as e:
189 raise ExpressionParsingError(f"{e} when parsing expression '{expression_str}' for {component_name}: ")
191 # load constraints (expressions for specific property infos)
192 # can only handle equality constraints for now
193 for component, expr_str, id in fs.constraint_exprs:
194 # get the time index
195 #for time_index in fs.time:
196 if component.index_set().dimen != 0 and component.index_set() != fs.time: 196 ↛ 197line 196 didn't jump to line 197 because the condition on line 196 was never true
197 raise ExpressionParsingError(f"Cannot add constraint for {component}: only time-indexed components are supported.")
198 try:
199 def constraint_rule(blk, time_index):
200 expression = parse_expression(expr_str, fs,time_index)
201 # make sure the units of the expression are the same as the component
202 u1, u2 = get_attached_unit(component), get_attached_unit(expression)
203 if not check_units_equivalent(u1, u2): 203 ↛ 204line 203 didn't jump to line 204 because the condition on line 203 was never true
204 raise ValueError(
205 f"Failed to add constraint for {component}: units do not match (expected {u1}, got {u2})"
206 )
207 return pyunits.convert(component[time_index], to_units=u2) == expression
208 c = Constraint(component.index_set(), rule= constraint_rule)
209 name = f"equality_constraint_{id}"
210 fs.add_component(name, c)
211 add_corresponding_constraint(fs, c, id)
212 except ExpressionParsingError as e:
213 raise ExpressionParsingError(f"Failed to parse constraint expression '{expr_str}' for {component}: {expr_str}, error: {e}")
216 def build_spec_dependency_tree(self, specs) -> tuple:
217 """
218 Builds dependency tree for expressions based on the references in their respective expressions.
219 """
220 dependencies = defaultdict(
221 set
222 ) # Maps an expression's result_id to a set of result_ids it depends on
223 result_ids = set() # A set to track all result_ids
225 # Get a list of all result_id's.
226 for spec in specs:
227 result_id = spec["id"]
228 result_ids.add(result_id)
230 for spec in specs:
231 result_id = spec["id"]
232 expression = spec["expression"]
234 # Find the result_ids that this expression depends on
235 dependent_expressions = self.get_dependent_expressions(
236 expression, result_ids
237 )
239 if dependent_expressions:
240 # If the expression depends on another result_id, add dependency
241 for id in dependent_expressions:
242 dependencies[id].add(result_id)
244 return dependencies, result_ids
246 def get_dependent_expressions(self, expression: str, all_result_ids: set) -> list:
247 """
248 Gets all result_ids referenced in the expression.
249 """
250 # match result_ids starting with 'id_' followed by numbers
251 ids = re.findall(r"\b(id_\d+)\b", expression)
253 # Filter the referenced_ids to only include those that are in all_result_ids
254 valid_referenced_ids = []
255 for id in ids:
256 # get the numeric part after "id_" and check if it's in the all_result_ids
257 numeric_id = int(id[3:])
258 if numeric_id in all_result_ids:
259 valid_referenced_ids.append(numeric_id)
261 return valid_referenced_ids
263 def topological_sort(self, dependencies: dict, result_ids: set) -> list:
264 """
265 Performs topological sorting on the specification dependency tree.
266 """
267 # Track teh in-degree count for all expressions (edges coming into it)
268 in_degree = defaultdict(int)
270 # Count dependencies for each expression
271 for result_id in result_ids:
272 for dep in dependencies[result_id]:
273 in_degree[dep] += 1
275 # Initialise the queue with result_ids that have no dependencies (in-degree 0)
276 dequeue = deque(
277 [result_id for result_id in result_ids if in_degree[result_id] == 0]
278 )
280 sorted_result_ids = []
282 while dequeue:
283 result_id = dequeue.popleft()
284 sorted_result_ids.append(result_id)
286 # loop thtough and decrement each dependent expression's in_degree, and append it to the deque if it has an in-degree of 0.
287 for dependent_result_id in dependencies[result_id]:
288 in_degree[dependent_result_id] -= 1
289 if in_degree[dependent_result_id] == 0:
290 dequeue.append(dependent_result_id)
292 # If there are any result_ids left with non-zero in-degree, a cycle exists so error.
293 if len(sorted_result_ids) != len(result_ids): 293 ↛ 294line 293 didn't jump to line 294 because the condition on line 293 was never true
294 raise ValueError(
295 "Cycle detected in the dependency graph. Check an expression does not reference itself!"
296 )
298 return sorted_result_ids
300 def initialise(self) -> None:
301 """
302 Expands the arcs and initialises the model
303 """
304 # check if initialisation is disabled for this scenario
305 if getattr(self.schema, "disable_initialization", False):
306 # if disable initialisation is set to True, then we don't need to initialise the model
307 _log.info("Initialisation is disabled for this scenario.")
309 # We need to "expand_arcs" to make them a bidirection link that actually imposes constraints on the model.
310 TransformationFactory("network.expand_arcs").apply_to(self.model)
312 self.timing.step_into("initialise_model")
314 # load tear guesses (including var/constraint unfixing & deactivation and/or equality constraint deactivation)
315 self.tears.load()
317 tears = self.tears._tears
319 def init_unit(unit):
320 if getattr(self.schema, "disable_initialization", False): 320 ↛ 321line 320 didn't jump to line 321 because the condition on line 320 was never true
321 return
322 if (
323 getattr(
324 self.schema,
325 "skip_initialization_for_units_with_initial_values",
326 False,
327 )
328 and self.unit_models.has_initial_values(unit)
329 ):
330 _log.info(
331 f"Skipping initialization for unit {unit} because initial values exist."
332 )
333 return
335 _log.info("Initializing unit %s", unit)
336 self.timing.add_timing(f"init_{unit.name}")
337 try:
338 #unit.display()
339 unit.initialize(outlvl=idaeslog.INFO)
340 #unit.report()
341 except Exception as e:
342 _log.warning(
343 "Error occurred while initializing unit %s: %s",
344 unit,
345 e,
346 )
348 self.timing.add_timing("setup_sequential_decomposition")
349 # Use SequentialDecomposition to initialise the model
350 seq = SequentialDecomposition(
351 run_first_pass=True,
352 )
353 seq.set_tear_set(tears)
354 # use create_graph to get the order of sequential decomposition, and also to
355 # find any units that are not connected to the sequential decomposition
356 G = seq.create_graph(self.model)
357 order = seq.calculation_order(G)
358 seq_blocks = []
359 for o in order:
360 # o is a list of "levels" of the graph,
361 # all items in the first level have to no dependencies, the next level only depends on the previous levels.
362 for b in o:
363 seq_blocks.append(b)
364 _log.info("Order of initialisation: %s", [blk.name for blk in seq_blocks])
365 # set all the tear guesses before running the decomposition
366 for arc in tears:
367 port = arc.destination
368 # guesses used are initial values for each var
369 guesses = {}
370 guesses = {key: get_value(var) for key, var in port.vars.items()}
371 _log.info("Guess for %s: %s", port, guesses)
373 self.timing.step_into("run_sequential_decomposition")
375 # sequential decomposition completes when all vars across port
376 # equalities are within tol of each other
377 seq.options["tol"] = 1e-2
378 seq.options["solve_tears"] = False
379 seq.run(self.model, init_unit)
381 self.timing.step_out()
382 self.timing.add_timing("initialise_disconnected_units")
383 # Initialise any unit model that is not connected to the sequential decomposition
384 for blk in self.model.fs.component_data_objects(
385 Block, descend_into=False, active=True
386 ):
387 ports = list(blk.component_objects(Port, descend_into=False))
388 if len(ports) == 0:
389 continue # if the block has no ports, then it is not a unit model
390 if blk in seq_blocks:
391 continue # already initialised by sequential decomposition
392 init_unit(blk)
394 # unfix guess vars
395 deactivate_fixed_guesses(self.model.fs.guess_vars)
397 self.timing.step_out()
399 def serialise(self) -> SolvedFlowsheetSchema:
400 self.timing.add_timing("serialise_model")
402 initial_values = {}
403 for unit_model_id, unit_model in self.unit_models._unit_models.items():
404 initial_values[str(unit_model_id)] = to_json(unit_model, return_dict=True, wts=StoreSpec.value())
406 solved_flowsheet = SolvedFlowsheetSchema(
407 group_id=self.schema.group_id,
408 properties=serialize_properties_map(self.model.fs),
409 initial_values=initial_values
410 )
412 return solved_flowsheet
414 def report_statistics(self) -> None:
415 """
416 Reports statistics about the model
417 """
418 report_statistics(self.model)
419 # I think the diagnostics toolbox is taking too long to run, so commenting it out for now.
420 # dt = DiagnosticsToolbox(self.model)
421 # dt.report_structural_issues()
422 # dt.display_overconstrained_set()
423 # dt.display_underconstrained_set()
425 def diagnose_problems(self) -> None:
426 from idaes.core.util import DiagnosticsToolbox
428 _log.info("=== DIAGNOSTICS ===")
430 report_statistics(self.model)
431 dt = DiagnosticsToolbox(self.model)
432 dt.report_structural_issues()
433 dt.display_overconstrained_set()
434 dt.display_underconstrained_set()
435 #dt.display_components_with_inconsistent_units()
436 try:
437 dt.compute_infeasibility_explanation()
438 except Exception as e:
439 _log.warning(
440 "Error computing infeasibility explanation: %s",
441 e,
442 )
443 try:
444 dt.report_numerical_issues()
445 except Exception as e:
446 _log.warning("Error reporting numerical issues: %s", e)
447 try:
448 dt.display_near_parallel_constraints()
449 except Exception as e:
450 _log.warning("Error displaying near parallel constraints: %s", e)
452 dt.display_variables_at_or_outside_bounds()
453 dt.display_variables_near_bounds()
454 try:
455 dt.display_constraints_with_extreme_jacobians()
456 except Exception as e:
457 _log.warning("Error displaying constraints with extreme jacobians: %s", e)
458 dt.display_constraints_with_large_residuals()
459 _log.info("=== END DIAGNOSTICS ===")
461 def get_unit_diagnostics(self):
462 diagnostics: list[tuple[Component, str]] = []
463 for unit_model in self.unit_models._unit_models.values():
464 if hasattr(unit_model, "diagnose"):
465 try:
466 results : list[tuple[Component, str]] = unit_model.diagnose()
467 diagnostics.extend(results)
468 # The unit diagnostics methods may be unreliable, as each unit model has different methods.
469 # So we wrap in a try-except to avoid the whole diagnostics process failing if one unit model's diagnostics fails.
470 except Exception as e:
471 _log.warning("Error Diagnosing %s: %s", unit_model, e)
472 return diagnostics
474 def get_unit_diagnostics_report(self) -> list[UnitDiagnosticsFinding]:
475 """Gets the diagnostics, but finds the property value ID for each diagnostic based on the component it references."""
476 diagnostics = self.get_unit_diagnostics()
477 diagnostics_report = [
478 UnitDiagnosticsFinding(
479 severity="error",
480 property_value_id=self.properties_map.get_id_by_component(component),
481 message=message)
482 for component, message in diagnostics
483 ]
484 diagnostics_report.extend(
485 [UnitDiagnosticsFinding(property_value_id=id,
486 message=f"High Infeasibility: {infeasibility:.4f}")
487 for id, infeasibility in
488 get_top_infeasibilities(self.properties_map, threshold=1e-3)
489 ]
490 )
492 return diagnostics_report
495 def degrees_of_freedom(self) -> int:
496 """
497 Returns the degrees of freedom of the model
498 """
499 return int(degrees_of_freedom(self.model))
501 def check_model_valid(self) -> None:
502 """
503 Checks if the model is valid by checking the
504 degrees of freedom. Will raise an exception if
505 the model is not valid.
506 """
507 self.timing.add_timing("check_model_valid")
509 degrees_of_freedom = self.degrees_of_freedom()
510 if degrees_of_freedom != 0: 510 ↛ 512line 510 didn't jump to line 512 because the condition on line 510 was never true
511 #self.model.display() # prints the vars/constraints for debugging
512 raise Exception(
513 f"Degrees of freedom is not 0. Degrees of freedom: {degrees_of_freedom}"
514 )
516 def solve(self):
517 """
518 Solves the model. Uses petsc_dae_by_time_element for the petsc solver
519 (dynamic only), or a standard SolverFactory solver otherwise.
520 """
521 self.timing.add_timing("solve_model")
522 _log.info("=== Starting Solve ===")
524 if self.schema.solver_option == "petsc": 524 ↛ 525line 524 didn't jump to line 525 because the condition on line 524 was never true
525 from idaes.core.solvers import petsc as idaes_petsc
527 if not self.schema.dynamic:
528 raise ValueError("The PETSc solver is only supported for dynamic flowsheets.")
529 result = idaes_petsc.petsc_dae_by_time_element(
530 self.model,
531 time=self.model.fs.time,
532 ts_options={
533 "--ts_type": "beuler",
534 "--ts_dt": 1,
535 "--ts_monitor": "",
536 "--ts_save_trajectory": 1,
537 },
538 )
540 opt = SolverFactory(self.schema.solver_option)
541 # opt.options["max_iter"] = 5000
543 opt.options["max_iter"] = 1000
544 try:
545 starting_values = get_current_values(self.properties_map)
546 res = opt.solve(self.model, tee=True)
547 updated_values = get_current_values(self.properties_map)
548 changes = detect_changes(updated_values, starting_values)
549 print_changes(changes, self.properties_map)
551 success = check_optimal_termination(res)
552 return success
553 except ValueError as e:
554 if str(e).startswith("No variables appear"): 554 ↛ 558line 554 didn't jump to line 558 because the condition on line 554 was always true
555 # https://github.com/Pyomo/pyomo/pull/3445
556 return True
557 else:
558 raise e
560 def optimize(self) -> None:
561 if self.schema.optimizations is None or self.schema.optimizations == []:
562 return
564 self.timing.add_timing("optimize_model")
565 _log.info("=== Starting Optimization ===")
567 # ipopt doesn't support multiple objectives, so we need to create
568 # a single objective expression.
569 # this is done by summing all objectives in the model, adding
570 # or subtracting based on the sense (minimize or maximize)
571 objective_terms = []
572 for schema in self.schema.optimizations:
573 # get the expression component to optimize
574 # TODO: This is assuming optimisation is run on a steady-state simulation. We need to change how this works to handle dynamics.
575 # For now, just hardcoding time_index=0
576 objective_component = self.model.fs.properties_map.get_component(schema.objective)
577 # the objective component should be a scalar in non-dynamic models
578 # in a dynamic model it'll be indexed across all time steps
579 # We sum up all time steps so each time step is weighted equally.
580 objective = dimensionless_objective_expression(
581 sum_pyomo_expressions(objective_component.values())
582 )
583 # add or subtract the objective based on the sense
584 sense = schema.sense
585 if sense == "minimize":
586 objective_terms.append(objective)
587 else:
588 objective_terms.append(-objective)
592 # unfix relevant vars (add to degrees of freedom)
593 for dof_info in schema.unfixed_variables:
594 id = dof_info.id # Id of the propertyValue for this degree of freedom
595 var: PropertyComponent = self.model.fs.properties_map.get(id)
598 # TODO: may need to handle deactivating constraints,
599 # for expressions that are constrained (instead of state vars)
600 c = var.corresponding_constraint
601 if c is not None:
602 # TODO: better typing for constraints
603 if isinstance(c, ScalarConstraint) or isinstance(c, IndexedConstraint): 603 ↛ 604line 603 didn't jump to line 604 because the condition on line 603 was never true
604 c.deactivate()
605 else:
606 deactivate_components(c)
607 else:
608 for i in var.component.values():
609 if isinstance(i, ScalarVar): 609 ↛ 610line 609 didn't jump to line 610 because the condition on line 609 was never true
610 i.unfix()
611 # Because if not, it is ExpressionData, meaning it is already an expression and doesn't need to be unfixed. (we've already checked if there is a constraint for it above too.)
614 # TODO: set attributes for upper and lower bounds of property infos. i.e. use propertyinfo id.
615 # Var is either a Variable or Expression
616 # set the minimum or maximum bounds for this variable if they are enabled
617 #self.model.upper_bound_12 = Constraint(expr= var <= upper_bound_value )
619 upper_bound = dof_info.upper_bound
620 lower_bound = dof_info.lower_bound
622 c = var.component
624 if upper_bound is not None: 624 ↛ 630line 624 didn't jump to line 630 because the condition on line 624 was always true
625 def upper_bound_rule(model,index):
626 return c[index] <= upper_bound
627 upper_bound_constraint = Constraint(c.index_set(),rule=upper_bound_rule)
628 setattr(self.model,"upper_bound_" + str(id), upper_bound_constraint)
630 if lower_bound is not None: 630 ↛ 593line 630 didn't jump to line 593 because the condition on line 630 was always true
631 def lower_bound_rule(model,index):
632 return c[index] >= lower_bound
633 lower_bound_constraint = Constraint(c.index_set(),rule=lower_bound_rule)
634 setattr(self.model,"lower_bound_" + str(id), lower_bound_constraint)
637 # add the objective to the model
638 objective_expr = sum_pyomo_expressions(objective_terms)
639 self.model.objective = Objective(expr=objective_expr, sense=minimize)
641 # solve the model with the objective
642 opt = SolverFactory(self.schema.solver_option)
644 if self.schema.solver_option != "conopt": 644 ↛ 647line 644 didn't jump to line 647 because the condition on line 644 was always true
645 opt.options["max_iter"] = 1000
647 try:
648 res = opt.solve(self.model, tee=True)
649 assert_optimal_termination(res)
650 except ValueError as e:
651 if str(e).startswith("No variables appear"):
652 # https://github.com/Pyomo/pyomo/pull/3445
653 pass
654 else:
655 raise e