Factorable Programming for PWL Approximations

[michaelbynum]

Summary/Motivation:

The purpose of this module/transformation is to convert any nonlinear model
to the following form:

min/max e(x_i)*f(x_j) + b^T*x
s.t.
        g_j(x_i)*h_j(x_k) + a_j^T*x >/</== 0
        g_j(x_i)/h_j(x_k) + a_j^T*x >/</== 0
        g_j(x_i)**h_j(x_k) + a_j^T*x >/</== 0

By doing so, each nonlinear function is only a function of one or two variables.
If this transformation is used prior to the nonlinear_to_pwl transformation,
it can, in some cases, significantly reduce the complexity of the PWL approximation.

Legal Acknowledgement

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

1. I agree my contributions are submitted under the BSD license.
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.

[codecov]

Codecov Report

❌ Patch coverage is 93.89068% with 19 lines in your changes missing coverage. Please review.
✅ Project coverage is 90.11%. Comparing base (68cefde) to head (545c678).

Files with missing lines	Patch %	Lines
pyomo/contrib/piecewise/transform/factorable.py	93.79%	19 Missing ⚠️

Additional details and impacted files

@@            Coverage Diff             @@
##             main    #3821      +/-   ##
==========================================
+ Coverage   90.10%   90.11%   +0.01%     
==========================================
  Files         904      905       +1     
  Lines      107113   107424     +311     
==========================================
+ Hits        96512    96808     +296     
- Misses      10601    10616      +15     

Flag	Coverage Δ
builders	29.13% <22.82%> (-0.01%)	⬇️
default	86.44% <93.89%> (?)
expensive	35.53% <22.82%> (?)
linux	87.60% <93.89%> (-2.00%)	⬇️
linux_other	87.60% <93.89%> (+0.02%)	⬆️
oldsolvers	28.08% <22.18%> (-0.02%)	⬇️
osx	82.97% <93.89%> (+0.03%)	⬆️
win	86.04% <93.89%> (+0.01%)	⬆️
win_other	86.04% <93.89%> (+0.01%)	⬆️

Flags with carried forward coverage won't be shown. Click here to find out more.

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:

• ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

[jsiirola]

Some minor comments, but really a pretty big design question.

I think I will advocate refactoring this to avoid the node_to_var_map and degree_map entirely and restrict/repurpose the use of the substitution_map to only managing Expression objects (and maybe caching 1/x?). The issue is that as written, those maps will end up with one entry for every node of every expression of the entire model. Even for modest models, that will be huge. The other reason to not make the dict is that the only time things will be looked up in the maps will be when processing the parent node of the node in context. I think a simpler / more efficient approach will be to return a tuple of (node, degree, varset) from exitNode.

I am still working through verifying the logic in the addition/product/division/pow handlers.

[michaelbynum]

Thanks for the review, @jsiirola. I agree with most of your feedback. However, I am going to push back on the refactor for the following reasons:

1. I have a dream where we create a separate transformation that actually identifies common subexpressions. The current implementation would be able to exploit that. I would really like to avoid creating separate PWL approximations for xy in one constraint and yx in another constraint.
2. You did motivate me to do some profiling. The first time I profiled this, the time spent in FBBT was significant. I was able to cut the total time in the transformation by a factor of 2 by switching to compute_bounds_on_expr (I was originally being lazy - we don't actually need FBBT in most cases here). I'll push these changes shortly. After that fix, it does look like some improvement could be made by moving away from ComponentMap. However, I think there are much bigger fish to fry that I would prefer to focus on first.
3. These are implementation details that we can change at any time without impacting users, so there is no pressing reason to make these changes now.

I'll push fixes for your other comments shortly.

[jsiirola]

OK. I don't think I found any real problems with this code. I still feel strongly that the handlers should return tuples of (node, degree, varlist) and not rely on huge dictionaries to pass that information up the tree, but that can be implemented in a subsequent PR.

If you have a chance, it would be nice to resolve the unreachable code in the product handler before merging.

[jsiirola]

This makes sense, as the walker always rebuilds the expression, even if nothing is being substituted out. In the future, we should consider revisiting this design and borrow ideas from the ExpressionReplacementVisitor and only recreate the (parts of the) expressions that we have to in order to effect the substitutions.

[jsiirola]

This gets repeated a lot. Should it be a helper function?

[michaelbynum]

I'm not sure there is anything for a helper function to do here? This is really just creating local variables. All of the "work" is already in visitor.create_aux_var.