This PR is motivated by discussion #1053 in which a situation the current implementation of find_eqpt does not find a good solution due to the problem being over/underconstrained. The main functional change is this PR is to help address that by adding a root_method method that allows a least squares solution to be obtained.

In addition, while looking through this (fairly old) code, it seemed worth making some additional changes:

• Renamed find_eqpt to find_operating_point to be more consistent with python-control naming conventions (and the fact that an operating point is not always an equilibrium point.
• Changed the return type for find_operating_point from a tuple of various lengths (depending on parameter settings) to a OperatingPoint object, which contains the states, inputs, and outputs corresponding to the operating point, along with the result of the call to the scipy.minimize.root function.
• Added root_method and root_kwargs parameters to allow customization of the root finding method.
• Maintained backward compatibility so that find_eqpt still works, as well as return_y (now return_outputs, with legacy processing of return_y) and return_result. In addition, the return object can be accessed as a tuple for consistency with the previous return signature.
• Updated docstrings and documentation to use the new conventions.
• Added unit tests to check all new functionality.

In the context of #1053, a solution to the issue there can be obtained by calling find_eqpt (now find_operating_point) with root_method='lm'.