diff --git a/idaes/core/util/model_diagnostics.py b/idaes/core/util/model_diagnostics.py index 67b7adc0f4..6ff4ee4def 100644 --- a/idaes/core/util/model_diagnostics.py +++ b/idaes/core/util/model_diagnostics.py @@ -286,7 +286,7 @@ def svd_sparse(jacobian, number_singular_values): CONFIG.declare( "parallel_component_tolerance", ConfigValue( - default=1e-4, + default=1e-8, domain=float, description="Tolerance for identifying near-parallel Jacobian rows/columns", ), @@ -1179,6 +1179,29 @@ def _collect_numerical_warnings(self, jac=None, nlp=None): self.display_constraints_with_extreme_jacobians.__name__ + "()" ) + # Parallel variables and constraints + partol = self.config.parallel_component_tolerance + par_cons = check_parallel_jacobian( + self._model, tolerance=partol, direction="row", jac=jac, nlp=nlp + ) + par_vars = check_parallel_jacobian( + self._model, tolerance=partol, direction="column", jac=jac, nlp=nlp + ) + if par_cons: + p = "pair" if len(par_cons) == 1 else "pairs" + warnings.append( + f"WARNING: {len(par_cons)} {p} of constraints are parallel" + f" (to tolerance {partol:.1E})" + ) + next_steps.append(self.display_near_parallel_constraints.__name__ + "()") + if par_vars: + p = "pair" if len(par_vars) == 1 else "pairs" + warnings.append( + f"WARNING: {len(par_vars)} {p} of variables are parallel" + f" (to tolerance {partol:.1E})" + ) + next_steps.append(self.display_near_parallel_variables.__name__ + "()") + return warnings, next_steps def _collect_numerical_cautions(self, jac=None, nlp=None): @@ -1441,7 +1464,6 @@ def report_numerical_issues(self, stream=None): lines_list=next_steps, title="Suggested next steps:", line_if_empty=f"If you still have issues converging your model consider:\n" - f"{TAB*2}display_near_parallel_constraints()\n{TAB*2}display_near_parallel_variables()" f"\n{TAB*2}prepare_degeneracy_hunter()\n{TAB*2}prepare_svd_toolbox()", footer="=", ) @@ -3619,7 +3641,13 @@ def ipopt_solve_halt_on_error(model, options=None): ) -def check_parallel_jacobian(model, tolerance: float = 1e-4, direction: str = "row"): +def check_parallel_jacobian( + model, + tolerance: float = 1e-4, + direction: str = "row", + jac=None, + nlp=None, +): """ Check for near-parallel rows or columns in the Jacobian. @@ -3627,6 +3655,11 @@ def check_parallel_jacobian(model, tolerance: float = 1e-4, direction: str = "ro as this means that the associated constraints or variables are (near) duplicates of each other. + For efficiency, the ``jac`` and ``nlp`` arguments may be provided if they are + already available. If these are provided, the provided model is not used. If + either ``jac`` or ``nlp`` is not provided, a Jacobian and ``PyomoNLP`` are + computed using the model. + This method is based on work published in: Klotz, E., Identification, Assessment, and Correction of Ill-Conditioning and @@ -3637,6 +3670,8 @@ def check_parallel_jacobian(model, tolerance: float = 1e-4, direction: str = "ro model: model to be analysed tolerance: tolerance to use to determine if constraints/variables are parallel direction: 'row' (default, constraints) or 'column' (variables) + jac: model Jacobian as a ``scipy.sparse.coo_matrix``, optional + nlp: ``PyomoNLP`` of model, optional Returns: list of 2-tuples containing parallel Pyomo components @@ -3651,7 +3686,8 @@ def check_parallel_jacobian(model, tolerance: float = 1e-4, direction: str = "ro "Must be 'row' or 'column'." ) - jac, nlp = get_jacobian(model, scaled=False) + if jac is None or nlp is None: + jac, nlp = get_jacobian(model, scaled=False) # Get vectors that we will check, and the Pyomo components # they correspond to. diff --git a/idaes/core/util/tests/test_model_diagnostics.py b/idaes/core/util/tests/test_model_diagnostics.py index e60c30e8b7..39dd09d561 100644 --- a/idaes/core/util/tests/test_model_diagnostics.py +++ b/idaes/core/util/tests/test_model_diagnostics.py @@ -975,10 +975,10 @@ def test_display_near_parallel_variables(self): model.v3 = Var() model.v4 = Var() - model.c1 = Constraint(expr=model.v1 == model.v2 - 0.99999 * model.v4) - model.c2 = Constraint(expr=model.v1 + 1.00001 * model.v4 == 1e-8 * model.v3) + model.c1 = Constraint(expr=1e-8 * model.v1 == 1e-8 * model.v2 - 1e-8 * model.v4) + model.c2 = Constraint(expr=1e-8 * model.v1 + 1e-8 * model.v4 == model.v3) model.c3 = Constraint( - expr=1e8 * (model.v1 + model.v4) + 1e10 * model.v2 == 1e-6 * model.v3 + expr=1e3 * (model.v1 + model.v4) + 1e3 * model.v2 == model.v3 ) dt = DiagnosticsToolbox(model=model) @@ -1121,7 +1121,7 @@ def test_collect_numerical_warnings_jacobian(self): warnings, next_steps = dt._collect_numerical_warnings() - assert len(warnings) == 3 + assert len(warnings) == 4 assert ( "WARNING: 2 Variables with extreme Jacobian values (<1.0E-08 or >1.0E+08)" in warnings @@ -1132,7 +1132,7 @@ def test_collect_numerical_warnings_jacobian(self): ) assert "WARNING: 1 Constraint with large residuals (>1.0E-05)" in warnings - assert len(next_steps) == 3 + assert len(next_steps) == 4 assert "display_variables_with_extreme_jacobians()" in next_steps assert "display_constraints_with_extreme_jacobians()" in next_steps assert "display_constraints_with_large_residuals()" in next_steps @@ -1338,8 +1338,7 @@ def test_report_numerical_issues_ok(self): Suggested next steps: If you still have issues converging your model consider: - display_near_parallel_constraints() - display_near_parallel_variables() + prepare_degeneracy_hunter() prepare_svd_toolbox() @@ -1369,9 +1368,11 @@ def test_report_numerical_issues_exactly_singular(self): Jacobian Condition Number: Undefined (Exactly Singular) ------------------------------------------------------------------------------------ -1 WARNINGS +3 WARNINGS WARNING: 2 Constraints with large residuals (>1.0E-05) + WARNING: 1 pair of constraints are parallel (to tolerance 1.0E-08) + WARNING: 1 pair of variables are parallel (to tolerance 1.0E-08) ------------------------------------------------------------------------------------ 0 Cautions @@ -1382,6 +1383,8 @@ def test_report_numerical_issues_exactly_singular(self): Suggested next steps: display_constraints_with_large_residuals() + display_near_parallel_constraints() + display_near_parallel_variables() ==================================================================================== """ @@ -1433,8 +1436,8 @@ def test_report_numerical_issues_jacobian(self): model.v2 = Var(initialize=0) model.v3 = Var(initialize=0) - model.c1 = Constraint(expr=model.v1 == model.v2) - model.c2 = Constraint(expr=model.v1 == 1e-8 * model.v3) + model.c1 = Constraint(expr=1e-2 * model.v1 == model.v2) + model.c2 = Constraint(expr=1e-2 * model.v1 == 1e-8 * model.v3) model.c3 = Constraint(expr=1e8 * model.v1 + 1e10 * model.v2 == 1e-6 * model.v3) dt = DiagnosticsToolbox(model=model) @@ -1445,14 +1448,15 @@ def test_report_numerical_issues_jacobian(self): expected = """==================================================================================== Model Statistics - Jacobian Condition Number: 1.407E+18 + Jacobian Condition Number: 1.118E+18 ------------------------------------------------------------------------------------ -3 WARNINGS +4 WARNINGS WARNING: 1 Constraint with large residuals (>1.0E-05) WARNING: 2 Variables with extreme Jacobian values (<1.0E-08 or >1.0E+08) WARNING: 1 Constraint with extreme Jacobian values (<1.0E-08 or >1.0E+08) + WARNING: 3 pairs of variables are parallel (to tolerance 1.0E-08) ------------------------------------------------------------------------------------ 4 Cautions @@ -1468,6 +1472,7 @@ def test_report_numerical_issues_jacobian(self): display_constraints_with_large_residuals() display_variables_with_extreme_jacobians() display_constraints_with_extreme_jacobians() + display_near_parallel_variables() ==================================================================================== """