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

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) 

26 

27install_model_statistics_optimizations() 

28 

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 

57 

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 

66 

67_log = idaeslog.getLogger(__name__) 

68 

69# from amplpy import modules 

70# Import required to allow the library to set the PATH and allow conopt to be found. 

71 

72 

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 = [] 

92 

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"]) 

97 

98 # Reaction package for the HDA reaction 

99 model.fs.reaction_params = reaction_props.HDAReactionParameterBlock( 

100 property_package=model.fs.peng_robinson 

101 ) 

102 

103 return model 

104 

105 

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 """ 

111 

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") 

118 

119 self.model = build_flowsheet(dynamic=schema.dynamic, time_set=schema.time_set) 

120 

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 

133 

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 ) 

159 

160 self.properties_map = self.model.fs.properties_map 

161 self.timing.step_out() 

162 

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 [] 

169 

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) 

173 

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) 

182 

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}: ") 

190 

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}") 

214 

215 

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 

224 

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) 

229 

230 for spec in specs: 

231 result_id = spec["id"] 

232 expression = spec["expression"] 

233 

234 # Find the result_ids that this expression depends on 

235 dependent_expressions = self.get_dependent_expressions( 

236 expression, result_ids 

237 ) 

238 

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) 

243 

244 return dependencies, result_ids 

245 

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) 

252 

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) 

260 

261 return valid_referenced_ids 

262 

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) 

269 

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 

274 

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 ) 

279 

280 sorted_result_ids = [] 

281 

282 while dequeue: 

283 result_id = dequeue.popleft() 

284 sorted_result_ids.append(result_id) 

285 

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) 

291 

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 ) 

297 

298 return sorted_result_ids 

299 

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.") 

308 

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) 

311 

312 self.timing.step_into("initialise_model") 

313 

314 # load tear guesses (including var/constraint unfixing & deactivation and/or equality constraint deactivation) 

315 self.tears.load() 

316 

317 tears = self.tears._tears 

318 

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 

334 

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 ) 

347 

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) 

372 

373 self.timing.step_into("run_sequential_decomposition") 

374 

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) 

380 

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) 

393 

394 # unfix guess vars 

395 deactivate_fixed_guesses(self.model.fs.guess_vars) 

396 

397 self.timing.step_out() 

398 

399 def serialise(self) -> SolvedFlowsheetSchema: 

400 self.timing.add_timing("serialise_model") 

401 

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()) 

405 

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 ) 

411 

412 return solved_flowsheet 

413 

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() 

424 

425 def diagnose_problems(self) -> None: 

426 from idaes.core.util import DiagnosticsToolbox 

427 

428 _log.info("=== DIAGNOSTICS ===") 

429 

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) 

451 

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 ===") 

460 

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 

473 

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 ) 

491 

492 return diagnostics_report 

493 

494 

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)) 

500 

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") 

508 

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 ) 

515 

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 ===") 

523 

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 

526 

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 ) 

539 

540 opt = SolverFactory(self.schema.solver_option) 

541 # opt.options["max_iter"] = 5000 

542 

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) 

550 

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 

559 

560 def optimize(self) -> None: 

561 if self.schema.optimizations is None or self.schema.optimizations == []: 

562 return 

563 

564 self.timing.add_timing("optimize_model") 

565 _log.info("=== Starting Optimization ===") 

566 

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) 

589 

590 

591 

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) 

596 

597 

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.) 

612 

613 

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 ) 

618 

619 upper_bound = dof_info.upper_bound 

620 lower_bound = dof_info.lower_bound 

621 

622 c = var.component 

623 

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) 

629 

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) 

635 

636 

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) 

640 

641 # solve the model with the objective 

642 opt = SolverFactory(self.schema.solver_option) 

643 

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 

646 

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