Skip to content

Revise bsm1 for metab #1597

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Draft
wants to merge 5 commits into
base: main
Choose a base branch
from

Conversation

adam-a-a
Copy link
Contributor

@adam-a-a adam-a-a commented Jun 24, 2025

Fixes/Resolves:

NA

Summary/Motivation:

Aid with model convergence when connecting METAB to activated sludge process.

Changes proposed in this PR:

  • revise asm1 flowsheet
  • on-the-fly code to get model to converge and draw quick comparisons

Legal Acknowledgement

By contributing to this software project, I agree to the following terms and conditions for my contribution:

  1. I agree my contributions are submitted under the license terms described in the LICENSE.txt file at the top level of this directory.
  2. I represent I am authorized to make the contributions and grant the license. If my employer has rights to intellectual property that includes these contributions, I represent that I have received permission to make contributions and grant the required license on behalf of that employer.

adam-a-a added 3 commits June 20, 2025 10:19
…o BSM1; only low flow rates cannot solve; need to increase initial flowrate by at least 2 orders of mag
…instead of 10/hr (can't solve at 10). Additionally, neither model is numerically stable; next, try remove P1 pump
@adam-a-a adam-a-a self-assigned this Jun 24, 2025
Comment on lines +348 to +364
# solve_stat = {}
# clones = {}
# for i, (k,v) in enumerate(new_comps.items()):
# clone_name = f"clone_change_{k}"
# m_clone = m.clone()
# old_val = pyo.value(m_clone.fs.feed.conc_mass_comp[0,k])
# m_clone.fs.feed.conc_mass_comp[0,k].fix(v)
# print(f"{k} was fixed.")
# res = solver.solve(m_clone, tee=True)
# if pyo.check_optimal_termination(res):
# solve_stat[k]=1
# else:
# solve_stat[k]=0

# clones[clone_name] = m_clone
# m_clone.fs.feed.conc_mass_comp[0,k].fix(old_val)
# del m_clone
Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@Morgan88888888 This bit here is for troubleshooting which component concentration might be problematic in terms of model convergence; it checks the converge by fixing concentration values one at a time and independently.

Uncomment this to (1) clone model, (2) fix concentration to new value for that one component, (3) attempt to solve that cloned model and store optimal termination as 1 and failure as 0 in solve_stat, (4) store model clone in clones dict for optional diagnostics later, (5) delete clone and restart, cycling to the next component to probe

Just seeing now that the steps where I record original default value and refix the original default value are not necessary anymore since I moved to just cloning the model in each iteration.

Comment on lines +367 to +369
for i, (k, v) in enumerate(new_comps.items()):
print(f"{k} was fixed.")
m.fs.feed.conc_mass_comp[0, k].fix(v)
Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

After ensuring that I can solve the model with each component's newly fixed concentration separately, I proceed to solve the model with ALL new concentration values fixed simultaneously to ensure that the model solves without issue.

Comment on lines +377 to +378
m.fs.feed.alkalinity.fix(47.7 * pyo.units.mol / pyo.units.m**3)
m.fs.feed.temperature.fix(308)
Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Being cautious about this, I then fix the new alkalinity and temp values to see if those cause issues.

Comment on lines +386 to +388
if "conc_mass_comp" in var.name:
print(var.name)
iscale.set_scaling_factor(var, 10 / (pyo.value(var) + 1e-6))
Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

scaling concentration is important

Comment on lines +391 to +393
m.fs.R3.KLa.fix(20)
m.fs.R4.KLa.fix(20)
m.fs.R5.KLa.fix(20)
Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

problems solving when Kla at 10

Comment on lines +394 to +409
# m.fs.feed.flow_vol.fix(0.00011 * 1e4
# # 134500221723699
# )

# for var in m.fs.component_data_objects(pyo.Var, descend_into=True):
# if "flow_vol" in var.name:
# iscale.set_scaling_factor(var, 1* pyo.value(var))
# m.fs.R1.volume.fix(1000 *5e-4 * pyo.units.m**3)
# m.fs.R2.volume.fix(1000 *5e-4 * pyo.units.m**3)
# m.fs.R3.volume.fix(1333 *5e-4* pyo.units.m**3)
# m.fs.R4.volume.fix(1333 *5e-4* pyo.units.m**3)
# m.fs.R5.volume.fix(1333 *5e-4* pyo.units.m**3)
# for i in range(5):
# blk = getattr(m.fs,f"R{i+1}")
# iscale.set_scaling_factor(blk.volume, 1e2)
# iscale.calculate_scaling_factors(m)
Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Tried lower influent flowrates but could not solve at 0.001 or 0.0001 (0.01 and higher did converge however). I also tried adjusting reactor volumes accordingly, with no success.

Comment on lines +416 to +418
COD_removal_w_metab = pyo.value(
1 - m.fs.Treated.properties[0].COD / m.fs.feed.properties[0].COD
)
Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

example recording COD removal with metab. Another line does so without metab. We can do this for other metrics too, but COD is likely most critical.

Comment on lines +451 to +530
brew_comp_vals = [
0.02,
5.6,
0.025,
1.156,
1.00e-10,
1.00e-10,
1.00e-10,
1.00e-10,
1.00e-10,
0.140067,
0.06,
0.05324,
]
m.fs.feed.alkalinity.fix(40.0 * pyo.units.mol / pyo.units.m**3)
m.fs.feed.temperature.fix(308)
for var in m.fs.component_data_objects(pyo.Var, descend_into=True):
# if "flow_vol" in var.name:
# iscale.set_scaling_factor(var, 1e1)
# if "temperature" in var.name:
# iscale.set_scaling_factor(var, 1e-1)
# if "pressure" in var.name:
# iscale.set_scaling_factor(var, 1e-3)
if "conc_mass_comp" in var.name:
# if pyo.value(var) <= 1e-10:
# iscale.set_scaling_factor(var, 1)
# else:
iscale.set_scaling_factor(var, 1e-3 / (pyo.value(var) + 1e-10))
iscale.calculate_scaling_factors(m.fs)
brew_comps = {comps[i]: brew_comp_vals[i] for i in range(len(comps))}
solve_stat = {}
clones = {}
for i, (k, v) in enumerate(brew_comps.items()):
clone_name = f"clone_change_{k}"
m_clone = m.clone()
old_val = pyo.value(m_clone.fs.feed.conc_mass_comp[0, k])
m_clone.fs.feed.conc_mass_comp[0, k].fix(v)
res = solver.solve(m_clone, tee=True)
if pyo.check_optimal_termination(res):
solve_stat[k] = 1
else:
solve_stat[k] = 0
print(f"{k} was fixed.")

clones[clone_name] = m_clone
m_clone.fs.feed.conc_mass_comp[0, k].fix(old_val)
del m_clone
print(solve_stat)
m.fs.R3.KLa.fix(20)
m.fs.R4.KLa.fix(20)
m.fs.R5.KLa.fix(20)
for i, (k, v) in enumerate(brew_comps.items()):
print(f"{k} was fixed.")
m.fs.feed.conc_mass_comp[0, k].fix(v)
# for var in m.fs.component_data_objects(pyo.Var, descend_into=True):
# # if "flow_vol" in var.name:
# # iscale.set_scaling_factor(var, 1e1)
# # if "temperature" in var.name:
# # iscale.set_scaling_factor(var, 1e-1)
# # if "pressure" in var.name:
# # iscale.set_scaling_factor(var, 1e-3)
# if "conc_mass_comp" in var.name:
# # if pyo.value(var) <= 1e-10:
# # iscale.set_scaling_factor(var, 1)
# # else:
# iscale.set_scaling_factor(var, 1e-1/(pyo.value(var)+1e-10))
# iscale.calculate_scaling_factors(m.fs)
res = solver.solve(m, tee=True)

from idaes.core.util.model_diagnostics import DiagnosticsToolbox

dt = DiagnosticsToolbox(m)
if pyo.check_optimal_termination(res):
print("Success!")
else:
raise RuntimeError("Failed to solve")

COD_removal_no_metab = pyo.value(
1 - m.fs.Treated.properties[0].COD / m.fs.feed.properties[0].COD
)
Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

without metab (using brewery compositions translated to ASM1)

@lbianchi-lbl lbianchi-lbl added the Priority:Normal Normal Priority Issue or PR label Jun 26, 2025
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
Priority:Normal Normal Priority Issue or PR
Projects
None yet
Development

Successfully merging this pull request may close these issues.

2 participants