Strongly enforce SLSQP Bounds

[andrewellis55]

Purpose

One of the key differences between the pyOptSparse SLSQP implementation and the SciPy SLSQP implementation is strong bound enforcement. This was added the the SciPy implementation in these two PRs

scipy/scipy#4502
scipy/scipy#13009

These changes are fairly small and easy enough to port over. I've proposed this with this PR. This I suppose should fix #301 and relate to #303.

Expected time until merged

This is not overly pressing, but my team uses OpenMDAO with the SciPy SLSQP rather than the pyOptSparse driver due to bound enforcement issues, so this will allow us to shift over to pyOptSparse once merged.

Type of change

• Bugfix (non-breaking change which fixes an issue)
• New feature (non-breaking change which adds functionality)
• Breaking change (non-backwards-compatible fix or feature)
• Code style update (formatting, renaming)
• Refactoring (no functional changes, no API changes)
• Documentation update
• Maintenance update
• Other (please describe)

Testing

I ran the hs015.py example and brought the initial guess outside to the bounds.

Checklist

• I have run flake8 and black to make sure the Python code adheres to PEP-8 and is consistently formatted
• I have formatted the Fortran code with fprettify or C/C++ code with clang-format as applicable
• I have run unit and regression tests which pass locally with my changes
• I have added new tests that prove my fix is effective or that my feature works
• I have added necessary documentation

[codecov]

Codecov Report

Attention: Patch coverage is 83.33333% with 1 line in your changes missing coverage. Please review.

Project coverage is 63.24%. Comparing base (bc021e4) to head (8f41509).
Report is 2 commits behind head on main.

Files with missing lines	Patch %	Lines
pyoptsparse/pySLSQP/pySLSQP.py	83.33%	1 Missing ⚠️

Additional details and impacted files

@@             Coverage Diff             @@
##             main     #413       +/-   ##
===========================================
- Coverage   74.92%   63.24%   -11.69%     
===========================================
  Files          22       22               
  Lines        3334     3338        +4     
===========================================
- Hits         2498     2111      -387     
- Misses        836     1227      +391     

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

[ewu63]

Hey, thanks for the PR and for fixing a long-standing issue with our SLSQP implementation. Could you clarify whether SLSQP will generate intermediate points which violate bounds, or if the final solution can violate the bounds also?

The high-level code looks fine but I haven't looked too closely yet. In the mean time, can I ask that you add a test to show the intended fix (i.e. a test that fails without the code changes here)? You can probably base it off of an existing test in test_h015.py or test_hs071.py

[andrewellis55]

Thanks for taking a look @ewu63, I've added a simple test that just raises a ValueError if the optimizer ever steps outside of the bounds. With the case I added, the main branch hits this ValueError, the new branch does not.

Also good point on ensure the clipping occurs to the final solution as well. I've added a few lines to ensure this and checked for it in the test I added.

[whophil]

Any opinions about making this behavior configurable with an option flag?

[marcomangano]

Thanks for the bugfix @andrewellis55! I will also take a closer look in the next few days, the implementation makes sense at a first pass. I believe the tests are currently failing because the code is missing a self.optName = "SLSQP" line in the test_slsqp_strong_bound_enforcement function.

@whophil I don't believe that there are scenarios where we want the optimizer to violate the design variables bounds, do you have an example in mind?
I think this bugfix is worth a new patch release anyway, we can highlight the fix / change in behavior in the release notes.

[akleb]

I am concerned about the behavior of clipping when you are trying to take steps with Hessian information on multi-dimensional problems.

If we have a slanted quadratic like the example unconstrained case in the MDO book (D.1.1), clipped Newton steps will not get the minimum within the design variable bounds if you have Hessian information. If you have the Hessian, you will always be taking steps exactly to the global optimum regardless of where you are in the design space. This means you will end up with the state pushed up against the bound aligned with the global optimum (purple in the attached image). If instead, you did not use Hessian information, the gradient would push you down to the correct bounded optimum where the gradient is perpendicular to the design space boundary (red in the attached figure).

If you are building an approximate Hessian I am not sure what would happen, but I am guessing you would have some sort of in between behavior. However, for a simple case like this, I think you would build the exact (linearized?) Hessian within one or two steps, so if you did not start exactly on the bounded optimum, you would not get the correct solution.

To solve this you have to either switch to a lower dimensional problem on bounds or include the design variable bounds in the Lagrangian. In my experience the better more robust solution is to include the design bounds in the Lagrangian, I think there are a lot of potential pitfalls in switching to a lower dimensional problem for problems with multimodality resulting in cycling behavior. I can describe my thoughts on this with more pictures if desired.

It could be that the design variables bounds are already included in the QP and some other weird behavior is causing the design variables to exceed their bounds. If this is the case, perhaps clipping is an okay solution. Regardless, I think that providing a multidimensional test case that shows we avoid this behavior would be beneficial. Also having this issue is much better than a solution that violates design variable bounds.

[kanekosh]

Thank you @andrewellis55. I wanted to ask you about the cases in which you observed SLSQP violating bounds. Does that happen during the optimization search or in the final optimized solution? If the latter, does it say optimization has failed, or does it falsely claim that it is successful?

In theory, SQP should never violate bound constraints (or generally, any linear constraints) unless:

• it takes alpha>1 in line search
• it treats bounds as nonlinear constraints
• it solves a relaxation QP problem (aka elastic problem), which may relax linear/bound constraints and allow bound violations.

I haven't checked SLSQP code or paper in detail yet, but if SLSQP does one of the above, that means the algorithm is designed to allow bound violation. If so, I'm a bit concerned with clipping because it could confuse the line search. In that case, I'd suggest making the clipping optional.

It is also possible that SLSQP doesn't do any of these (again I haven't checked the code) and violating the bounds is simply a bug. If so, I think clipping is a good, simple fix.

[andrewellis55]

@akleb You raise some good points. First of all, as other have mentioned, the optimizer generally should never step outside these bounds, meaning that as long as all goes well, the clipping should never occur and thus the hessian approximation should be unaffected. I've ran the example problem you've provided from the MDO book with the code as follows

from pyoptsparse import Optimization, SLSQP
import numpy as np


def slanted_quadratic(xdict):
    x1 = xdict['x1']
    x2 = xdict['x2']
    beta = 3/2  # beta = 3/2 as shown in the image
    f = x1**2 + x2**2 - beta * x1 * x2  # Function (D.1)
    funcs = {}
    funcs['f'] = f
    return funcs, False

def slanted_quadratic_grad(xdict, funcs):
    x1 = xdict['x1']
    x2 = xdict['x2']
    beta = 3/2
    grad_f = np.array([2*x1 - beta * x2, 2*x2 - beta * x1])
    grad = {}
    grad['f'] = {'x1': grad_f[0], 'x2': grad_f[1]}
    return grad, False

# Set up the optimization problem
optProb = Optimization('Slanted Quadratic Minimization', slanted_quadratic)
optProb.addVar('x1', 'c', lower=-10.0, upper=10.0, value=-9)
optProb.addVar('x2', 'c', lower=-10.0, upper=-5.0, value=9)
optProb.addObj('f')

# Define the optimizer
optimizer = SLSQP()

# Solve the optimization problem
sol = optimizer(optProb, sens=slanted_quadratic_grad)

# Print the results
print(sol)

I've added a print statement before the clipping at the python level in slfunc to plot the path of the variables

            # =================================================================
            # SLSQP - Objective/Constraint Values Function
            # =================================================================
            def slfunc(m, me, la, n, f, g, x):
                print(x)
                fobj, fcon, fail = self._masterFunc(np.clip(x, blx, bux), ["fobj", "fcon"])
                f = fobj
                g[0:m] = -fcon
                slsqp.pyflush(self.getOption("IOUT"))
                return f, g

Here is the trajectory

1: [-9. -5.]
2: [ 1.5 -8.5]
3: [-5.37931034 -6.20689655]
4: [-3.85207101 -5.        ]
5: [-3.75571076 -5.        ]
6: [-3.75 -5.  ]

Although it isn't clear in the image that this is the minima along the line, this can be easily verified

>>> def slanted_quadratic(x1, x2, beta=1.5):
...     return x1**2 + x2**2 - beta * x1 * x2
>>> slanted_quadratic(-5, -3.74, beta=1.5)
10.9376
>>> slanted_quadratic(-5, -3.75, beta=1.5)
10.9375
>>> slanted_quadratic(-5, -3.76, beta=1.5)
10.9376

As can be seen , in this example SLSQP actually never exceeds the bounds and behaves normally. There have been edge cases that I (as well as other users #301 ) have seen where this is not the case for a particular iteration. An example of this that I've personally run is a optimizer of a tube with design variables inner_diameter and thickness and a safety_margin constraint that is essentially a 1/x function of force. The optimizer has stepped the thickness to be negative (exceeding the lower bounds) and caused the optimization to fail as the gradients of the margin function point in the opposite direction on the left side of the asymtote.

[andrewellis55]

@kanekosh Both I and other users (#301) have observed SLSQP exceeding bounds. This was also observed in the SciPy implementation leading to this fix scipy/scipy#4502. My code is a pure copy paste of the fortran changes and a slight adaptation of the python ones that were implemented in SciPy.

As for why it's happening, I frankly havn't read through the details of the SLSQP code enough to say (I'm not much of a fortran guy), but if I were to hazard a guess based off my general knowledge of SQP algorithms, small bound violations could come from numerical issues and larger bound violations could come if the bounds are treated as general non-linear constraints and/or are not included in the active set.

The fix I've added is the one that SciPy has been running for the past ~9 years and so I think it's generally a fairly safe one.

[kanekosh]

@andrewellis55 Would you mind adding lines of code to raise a warning when the clipping occurs, like scipy does here?

There might be a scenario where this clipping would confuse the optimization algorithm, so raising a warning would be beneficial just in case. (The benefit of enforcing the bounds is much greater than the potential convergence issue by clipping, so I think this is a good fix and I appreciate it)

[whophil]

@whophil I don't believe that there are scenarios where we want the optimizer to violate the design variables bounds, do you have an example in mind?

@marcomangano my question was about enabling "fallback" to the previous implementation of the SLSQP algorithm without having to install an older version ofpyoptsparse.

[andrewellis55]

@andrewellis55 Would you mind adding lines of code to raise a warning when the clipping occurs, like scipy does here?

@kanekosh Done!

[andrewellis55]

@whophil I don't believe that there are scenarios where we want the optimizer to violate the design variables bounds, do you have an example in mind?

@marcomangano my question was about enabling "fallback" to the previous implementation of the SLSQP algorithm without having to install an older version ofpyoptsparse.

While that's probably simple enough on the python side, it would require a bit of modification on Fortran wrapper side too which would made this PR a bit more than a copy/paste job

[marcomangano]

Thanks again for the fix, I think this is a welcome improvement that addresses a long-standing bug. I don't see the need for further tweaks of the source code - if we had to include all the possible improvements we would probably end up with a substantially different code, or a different optimizer altogether. Moving forward, we should consider the discussion in #303 and switch to an external maintained version.