Skip to content

Poorly formulated constraint in MCAS model #1547

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

Open
andrewlee94 opened this issue Jan 2, 2025 · 5 comments · May be fixed by #1553
Open

Poorly formulated constraint in MCAS model #1547

andrewlee94 opened this issue Jan 2, 2025 · 5 comments · May be fixed by #1553
Assignees
Labels
enhancement New feature or request Priority:Normal Normal Priority Issue or PR

Comments

@andrewlee94
Copy link
Collaborator

Description

The MCAS property model has the following constraint formulation (line 1464):

def rule_mass_frac_phase_comp(b, p, j):
    return b.mass_frac_phase_comp[p, j] == b.flow_mass_phase_comp[p, j] / sum(
        b.flow_mass_phase_comp[p, j] for j in self.params.component_list
    )

This leads to a potential division by zero as flow goes to zero and thus would be better reformulated as:

def rule_mass_frac_phase_comp(b, p, j):
    return b.flow_mass_phase_comp[p, j] == b.mass_frac_phase_comp[p, j]  * sum(
        b.flow_mass_phase_comp[p, j] for j in self.params.component_list
    )

However, doing this results in issues with the scaling methods, primarily due to the transform_property_constraints function which assumes that all constraints are written in the form of property = function(). I suspect that this might have further reaching effects than just this.

Motivation

Generally, constraints should be formulated in ways to avoid potential singularities, such as avoiding divisions where possible (unless the denominator is a simple constant). Doing so helps avoid numerical issues during solving and will generally make models more robust (and often better scaled).

Possible Implementation

This issue effectively has two parts:

  1. Looking for poorly formulated constraints and rewriting them (relatively easy, especially with the diagnostics tools), and
  2. Updating the scaling routines to accommodate these changes (somewhat more effort, especially with the above mentioned transform_property_constraints function). However, given the current push to using the new Scaler classes this hopefully becomes much easier (or at least can be done at the same time).

Additional Context

No response

@andrewlee94 andrewlee94 added the enhancement New feature or request label Jan 2, 2025
@MarcusHolly
Copy link
Contributor

@andrewlee94 I've started to look into this. What "issues with scaling methods" does addressing the potential division by zero bring up? I've made this change locally, but I can't seem to find any issues.

@andrewlee94
Copy link
Collaborator Author

@MarcusHolly The transform_property_constraints, which is used in many of the old-style scaling methods, breaks as it assumes the constraint has the form property = function.

If we look only at the new Scaling tools this issue goes away, but there is the need to keep backward compatibility and thus the old-style methods would need to be updated.

@MarcusHolly
Copy link
Contributor

And by breaking, do you mean that it's silently breaking? From what I can tell, all the tests and flowsheets appear to be running fine after reformulating the rule_mass_frac_phase_comp constraint.

Also, I'm not sure how well this utility function was really working in the first place. When I use report_scaling_factors on flowsheets with the GAC property model, the scaling factor set for eq_mass_frac_phase_comp is None before and after the reformulation. These flowsheets call iscale.calculate_scaling_factors, so it should be calling the transform_property_constraints which calls iscale.constraint_scaling_transform(c, sf). So there may actually some issue with the constraint_scaling_transform.

@andrewlee94
Copy link
Collaborator Author

Oh, no, it was very vocal for me. However, I might have done the following for the reformulation:

def rule_mass_frac_phase_comp(b, p, j):
    return (
        b.mass_frac_phase_comp[p, j]
        * sum(b.flow_mass_phase_comp[p, j] for j in self.params.component_list)) == b.flow_mass_phase_comp[p, j]
    )

Try doing this - if this does trigger a lot of warnings then it means there are two issues that need to be addressed. The first being that the function requires the form property = function() and the second being that it assumes the LHS is the property of interest without confirming it.

@MarcusHolly
Copy link
Contributor

The above formulation does not break any flowsheets or tests for me locally either. Perhaps I should create a fresh environment to see if I see the warnings then.

@ksbeattie ksbeattie added the Priority:Normal Normal Priority Issue or PR label Jan 23, 2025
@MarcusHolly MarcusHolly linked a pull request Jan 30, 2025 that will close this issue
@lbianchi-lbl lbianchi-lbl linked a pull request Feb 27, 2025 that will close this issue
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
enhancement New feature or request Priority:Normal Normal Priority Issue or PR
Projects
None yet
Development

Successfully merging a pull request may close this issue.

4 participants