Refactor ExternalPyomoModel to allow different implicit function solvers

[Robbybp]

Summary/Motivation:

I would like to be able to test the implicit function implementation in ExternalPyomoModel with different solvers for the square system. This PR allows this by defining an implicit function base class, a subclass of which can be passed to ExternalPyomoModel via the solver_class argument.

These updates allow the following performance comparisons to be made:

1. This compares the current default (calculate-variable-from-constraint for 1x1 systems and CyIpopt for higher dimension systems) with an implementation that uses the SciPy fsolve function for higher dimension systems:

Comparing timing profiles between the two implementations on the CLC dynamic optimization, we can see that the speedup is driven by the solve of the square NLPs.

With the original implementation:

Identifier                 ncalls   cumtime   percall      %
------------------------------------------------------------
root                            1     8.959     8.959  100.0
     -------------------------------------------------------
     __init__                  10     1.940     0.194   21.6
     hessian                   90     1.330     0.015   14.8
     jacobian                 100     0.565     0.006    6.3
     set_inputs               120     4.563     0.038   50.9
               ---------------------------------------------
               solve_nlp     2400     1.925     0.001   42.2
               other          n/a     2.637       n/a   57.8
               =============================================
     other                    n/a     0.562       n/a    6.3
     =======================================================
============================================================

and with fsolve for the higher dimension square systems:

Identifier                 ncalls   cumtime   percall      %
------------------------------------------------------------
root                            1     8.138     8.138  100.0
     -------------------------------------------------------
     __init__                  10     1.924     0.192   23.6
     hessian                  100     1.404     0.014   17.2
     jacobian                 110     0.625     0.006    7.7
     set_inputs               130     3.543     0.027   43.5
               ---------------------------------------------
               solve_nlp     2600     0.493     0.000   13.9
               other          n/a     3.050       n/a   86.1
               =============================================
     other                    n/a     0.643       n/a    7.9
     =======================================================
============================================================

2. We can also compare what happens when we switch out calculate_variable_from_constraint for a SciPy scalar Newton solver (here the secant-Newton hybrid). Both implementations here use fsolve for higher-dimension square systems:

Comparing timing profiles, the slowdown is due to the amount of time spent constructing extra PyomoNLPs in the latter case and the fact SciPy's scalar Newton/secant solver (with PyomoNLP's ASL interface for function/derivative evaluations) seems to be slower than calculate_variable_from_constraint.

With calculate_variable_from_constraint:

Identifier                 ncalls   cumtime   percall      %
------------------------------------------------------------
root                            1     8.010     8.010  100.0
     -------------------------------------------------------
     __init__                  10     1.915     0.191   23.9
               ---------------------------------------------
               PyomoNLP        10     0.461     0.046   24.1
               other          n/a     1.453       n/a   75.9
               =============================================
     hessian                  100     1.430     0.014   17.9
     jacobian                 110     0.619     0.006    7.7
     set_inputs               130     3.509     0.027   43.8
               ---------------------------------------------
               calc_var     45630     2.660     0.000   75.8
               solve_nlp     2600     0.470     0.000   13.4
               other          n/a     0.380       n/a   10.8
               =============================================
     other                    n/a     0.537       n/a    6.7
     =======================================================
============================================================

and with the secant-Newton wrapper around SciPy's scalar Newton solver:

Identifier                 ncalls   cumtime   percall      %
------------------------------------------------------------
root                            1    11.990    11.990  100.0
     -------------------------------------------------------
     __init__                  10     4.688     0.469   39.1
               ---------------------------------------------
               PyomoNLP        10     2.881     0.288   61.5
               other          n/a     1.807       n/a   38.5
               =============================================
     hessian                  100     1.447     0.014   12.1
     jacobian                 110     0.625     0.006    5.2
     set_inputs               130     4.572     0.035   38.1
               ---------------------------------------------
               solve_nlp    48230     4.049     0.000   88.6
               other          n/a     0.522       n/a   11.4
               =============================================
     other                    n/a     0.658       n/a    5.5
     =======================================================
============================================================

These time-profile comparisons use the utilities in #2651, but this branch does not depend on that PR.

Changes proposed in this PR:

• A base class, PyomoImplicitFunctionBase defining an API for implicit function solvers that can be used by ExternalPyomoModel
• Two subclasses, ImplicitFunctionSolver and SccImplicitFunctionSolver that implement this API with and without a strongly connected component decomposition
• A base class, NlpSolverBase for NLP solvers that can be used by the above implicit function solvers
• Two subclasses, CyIpoptSolverWrapper and ScipySolverWrapper, that implement this API with their respective solvers for the NLP
• Refactor ExternalPyomoModel to accept an implicit function solver via the solver_class (and solver_options) argument instead of interfacing directly to the NLP solver (and implementing the sequential SCC solve itself)

Note that this PR removes the use_cyipopt argument and also removes the ability to use a "Pyomo solver" (e.g. pyo.SolverFactory('ipopt')) within ExternalPyomoModel (unless such a solver is implemented as a subclass of PyomoImplicitFunctionBase, which I have not done). While this was useful functionality for testing in the case where CyIpopt was not available, the SciPy square solvers can now be used instead, which should always be available as SciPy is a hard dependency of most of this machinery. Also worth noting that working around the CyIpopt-not-available case was always a bit contrived, as CyIpopt is (currently) the only solver that ExternalPyomoModel works with!

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.

[jsiirola]

Why is this function defined here? It seems a bit unrelated to (and more general than) square solvers.

[Robbybp]

Good point. I think I originally had it here before I decided to make an implicit_functions.py file.

[jsiirola]

You are adding logging here, but removing it in solvers/scipy_solvers.py. Is this just leftover debugging?

[Robbybp]

This was intentional. I wanted to be able to distinguish between time spent in function versus derivative evaluation when switching between Newton and Secant methods. However, this is admittedly a very low-level and frequently called method to add a timer call to. Here is comparison showing the speedup when this call is removed:

I think this is a sufficient difference that these calls should be removed.

[Robbybp]

However, note that in these examples, Newton solves are performed very frequently on very small systems, which makes the function/derivative methods particularly performance-critical. I'm not sure this changes the calculus about whether to include timing calls in low-level methods.

[jsiirola]

Is this timing from debugging? Again, we seem to be inconsistent in pynumero where we are instrumenting timing and where we are removing it. I wonder if it should all be removed (pynumero is meant to be efficient, and timer adds (not completely insignificant overhead).

[Robbybp]

My motivation in the timing calls is so that timing statistics like this:

Identifier                 ncalls   cumtime   percall      %
------------------------------------------------------------
root                            1    11.990    11.990  100.0
     -------------------------------------------------------
     __init__                  10     4.688     0.469   39.1
               ---------------------------------------------
               PyomoNLP        10     2.881     0.288   61.5
               other          n/a     1.807       n/a   38.5
               =============================================
     hessian                  100     1.447     0.014   12.1
     jacobian                 110     0.625     0.006    5.2
     set_inputs               130     4.572     0.035   38.1
               ---------------------------------------------
               solve_nlp    48230     4.049     0.000   88.6
               other          n/a     0.522       n/a   11.4
               =============================================
     other                    n/a     0.658       n/a    5.5
     =======================================================
============================================================

can be generated without intrusive modification. I think for time-consuming methods that we might commonly want to profile, the timing calls are justified. See above for a performance diff when timing is introduced in small, fast methods that are called often, which suggests that they should not be a permanent feature of these methods.

[andrewlee94]

@Robbybp I just tried using "scipy.fsolve" in some IDAES Applications, and ran into the following error:

TypeError: solve() got an unexpected keyword argument 'tee'

At first glance it looks like the scipy solver interfaces are missing the tee argument - whilst there might not be much output from these solvers (if any), the argument should probably be supported for consistency. Otherwise, downstream users are going to need to determine whether or not to send the tee argument based on the solver being used.

[Robbybp]

@andrewlee94 You can try by branch in #2668