From b0d39d8baddb41bcafc04a26e21dad6e904a9ff7 Mon Sep 17 00:00:00 2001 From: Bernard Knueven Date: Mon, 11 Dec 2023 16:18:41 -0700 Subject: [PATCH 1/2] build ill conditioning problem iterating over non-zeros --- idaes/core/util/model_diagnostics.py | 57 ++++++++++++---------------- 1 file changed, 25 insertions(+), 32 deletions(-) diff --git a/idaes/core/util/model_diagnostics.py b/idaes/core/util/model_diagnostics.py index cba946ee24..4df6e09fa8 100644 --- a/idaes/core/util/model_diagnostics.py +++ b/idaes/core/util/model_diagnostics.py @@ -39,6 +39,7 @@ Expression, Objective, Param, + RangeSet, Set, SolverFactory, value, @@ -3257,49 +3258,41 @@ def check_ill_conditioning( inverse_target_kappa = 1e-16 / target_feasibility_tol - # Get a mapping of row/column names to index in jac - con_dict = {i.name: idx for idx, i in enumerate(nlp.get_pyomo_constraints())} - var_dict = {i.name: idx for idx, i in enumerate(nlp.get_pyomo_variables())} + # Set up the components we will analyze, either row or column + if direction == "row": + components = nlp.get_pyomo_constraints() + components_set = RangeSet(0, len(components) - 1) + results_set = RangeSet(0, nlp.n_primals() - 1) + jac = jac.transpose().tocsr() + elif direction == "column": + components = nlp.get_pyomo_variables() + components_set = RangeSet(0, len(components) - 1) + results_set = RangeSet(0, nlp.n_constraints() - 1) + jac = jac.tocsr() # Build test problem inf_prob = ConcreteModel() - inf_prob.con_set = Set(initialize=con_dict.keys()) - inf_prob.var_set = Set(initialize=var_dict.keys()) - - if direction == "row": - set1 = inf_prob.con_set - set2 = inf_prob.var_set - else: - set2 = inf_prob.con_set - set1 = inf_prob.var_set - inf_prob.y_pos = Var(set1, initialize=0, bounds=(0, None)) - inf_prob.y_neg = Var(set1, initialize=0, bounds=(0, None)) - inf_prob.y = Expression(set1, rule=lambda m, i: m.y_pos[i] - m.y_neg[i]) + inf_prob.y_pos = Var(components_set, bounds=(0, None)) + inf_prob.y_neg = Var(components_set, bounds=(0, None)) + inf_prob.y = Expression(components_set, rule=lambda m, i: m.y_pos[i] - m.y_neg[i]) - inf_prob.res_pos = Var(set2, initialize=1e-2, bounds=(0, None)) - inf_prob.res_neg = Var(set2, initialize=1e-2, bounds=(0, None)) - inf_prob.res = Expression(set2, rule=lambda m, i: m.res_pos[i] - m.res_neg[i]) + inf_prob.res_pos = Var(results_set, bounds=(0, None)) + inf_prob.res_neg = Var(results_set, bounds=(0, None)) + inf_prob.res = Expression( + results_set, rule=lambda m, i: m.res_pos[i] - m.res_neg[i] + ) def b_rule(b, i): lhs = 0.0 - for j in set1: - if direction == "row": - iidx = var_dict[i] - jidx = con_dict[j] - jac_entry = jac[jidx, iidx] - else: - jidx = var_dict[j] - iidx = con_dict[i] - jac_entry = jac[iidx, jidx] - - if jac_entry != 0: - lhs += jac_entry * b.y[j] + row = jac.getrow(i) + for j, val in zip(row.indices, row.data): + lhs += val * b.y[j] return lhs == b.res[i] - inf_prob.by = Constraint(set2, rule=b_rule) + inf_prob.by = Constraint(results_set, rule=b_rule) # Normalization of y inf_prob.normalize = Constraint( @@ -3363,7 +3356,7 @@ def b_rule(b, i): val = value(inf_prob.y[i]) if abs(val) < cutoff: break - ill_cond.append((i, val)) + ill_cond.append((components[i].name, val)) return ill_cond From 5316a01bebc7f8ba491e098ea7613876656b04f8 Mon Sep 17 00:00:00 2001 From: Bernard Knueven Date: Tue, 12 Dec 2023 15:21:07 -0700 Subject: [PATCH 2/2] return component instead of its name --- idaes/core/util/model_diagnostics.py | 2 +- .../core/util/tests/test_model_diagnostics.py | 36 +++++++++---------- 2 files changed, 19 insertions(+), 19 deletions(-) diff --git a/idaes/core/util/model_diagnostics.py b/idaes/core/util/model_diagnostics.py index 4df6e09fa8..2912c475ba 100644 --- a/idaes/core/util/model_diagnostics.py +++ b/idaes/core/util/model_diagnostics.py @@ -3356,7 +3356,7 @@ def b_rule(b, i): val = value(inf_prob.y[i]) if abs(val) < cutoff: break - ill_cond.append((components[i].name, val)) + ill_cond.append((components[i], val)) return ill_cond diff --git a/idaes/core/util/tests/test_model_diagnostics.py b/idaes/core/util/tests/test_model_diagnostics.py index 8022b15292..713f128562 100644 --- a/idaes/core/util/tests/test_model_diagnostics.py +++ b/idaes/core/util/tests/test_model_diagnostics.py @@ -3236,34 +3236,34 @@ def test_beta(self, model, caplog): @pytest.mark.solver def test_rows(self, model): assert check_ill_conditioning(model, direction="row") == [ - ("c4", pytest.approx(0.50000002, rel=1e-5)), - ("c1", pytest.approx(0.49999998, rel=1e-5)), + (model.c4, pytest.approx(0.50000002, rel=1e-5)), + (model.c4, pytest.approx(0.49999998, rel=1e-5)), ] @pytest.mark.component @pytest.mark.solver def test_rows(self, afiro): assert check_ill_conditioning(afiro, direction="row") == [ - ("R09", pytest.approx(0.5, rel=1e-5)), - ("R09b", pytest.approx(0.5, rel=1e-5)), + (afiro.R09, pytest.approx(0.5, rel=1e-5)), + (afiro.R09b, pytest.approx(0.5, rel=1e-5)), ] @pytest.mark.component @pytest.mark.solver def test_columns(self, afiro): assert check_ill_conditioning(afiro, direction="column") == [ - ("X39", pytest.approx(1.1955465, rel=1e-5)), - ("X23", pytest.approx(1.0668697, rel=1e-5)), - ("X25", pytest.approx(-1.0668697, rel=1e-5)), - ("X09", pytest.approx(-0.95897123, rel=1e-5)), - ("X13", pytest.approx(-0.95897123, rel=1e-5)), - ("X06", pytest.approx(0.91651956, rel=1e-5)), - ("X10", pytest.approx(0.91651956, rel=1e-5)), - ("X36", pytest.approx(0.76204977, rel=1e-5)), - ("X31", pytest.approx(-0.39674454, rel=1e-5)), - ("X35", pytest.approx(-0.39674454, rel=1e-5)), - ("X16", pytest.approx(0.14679548, rel=1e-5)), - ("X38", pytest.approx(-0.14679548, rel=1e-5)), - ("X15", pytest.approx(-0.042451666, rel=1e-5)), - ("X37", pytest.approx(-0.036752232, rel=1e-5)), + (afiro.X39, pytest.approx(1.1955465, rel=1e-5)), + (afiro.X23, pytest.approx(1.0668697, rel=1e-5)), + (afiro.X25, pytest.approx(-1.0668697, rel=1e-5)), + (afiro.X09, pytest.approx(-0.95897123, rel=1e-5)), + (afiro.X13, pytest.approx(-0.95897123, rel=1e-5)), + (afiro.X06, pytest.approx(0.91651956, rel=1e-5)), + (afiro.X10, pytest.approx(0.91651956, rel=1e-5)), + (afiro.X36, pytest.approx(0.76204977, rel=1e-5)), + (afiro.X31, pytest.approx(-0.39674454, rel=1e-5)), + (afiro.X35, pytest.approx(-0.39674454, rel=1e-5)), + (afiro.X16, pytest.approx(0.14679548, rel=1e-5)), + (afiro.X38, pytest.approx(-0.14679548, rel=1e-5)), + (afiro.X15, pytest.approx(-0.042451666, rel=1e-5)), + (afiro.X37, pytest.approx(-0.036752232, rel=1e-5)), ]