Skip to content

Commit

Permalink
Merge pull request #20 from bknueven/ill_cond_sparse
Browse files Browse the repository at this point in the history 
build ill conditioning problem iterating over non-zeros
  • Loading branch information
@andrewlee94
andrewlee94 authored Dec 13, 2023
2 parents 447055c + 5316a01 commit c4c9f5d
Show file tree
Hide file tree
Showing 2 changed files with 43 additions and 50 deletions.
57 changes: 25 additions & 32 deletions idaes/core/util/model_diagnostics.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -39,6 +39,7 @@
Expression,
Objective,
Param,
RangeSet,
Set,
SolverFactory,
value,
Expand Down Expand Up @@ -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(
Expand Down Expand Up @@ -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], val))

return ill_cond

Expand Down
36 changes: 18 additions & 18 deletions idaes/core/util/tests/test_model_diagnostics.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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)),
]

0 comments on commit c4c9f5d

Please sign in to comment.