Bug report: scipy>=1.5.0 breaks yaw optimization with SLSQP and equality constraints

[Bartdoekemeijer]

Bug description
When I use a recent version of Scipy (versions 1.5.0 or later), I cannot set an equality constraint in the SLSQP yaw optimization functions successfully. It seems that the optimization attempts to calculate a gradient for the equality-constrained turbines, thereby dividing some value by 0 and exiting the optimization without success after the first iteration.

To Reproduce
The following is a minimal example that reproduces the issue for me:

import os
import numpy as np

import floris.tools as wfct
from floris.tools.optimization.scipy.yaw import YawOptimization


# Instantiate the FLORIS object
file_dir = os.path.dirname(os.path.abspath(__file__))
fi = wfct.floris_interface.FlorisInterface(
    os.path.join(file_dir, "../../../example_input.json")
)

# Set turbine locations to 3 turbines in a row
D = fi.floris.farm.turbines[0].rotor_diameter
layout_x = [0, 7 * D, 14 * D]
layout_y = [0, 0, 0]
fi.reinitialize_flow_field(layout_array=(layout_x, layout_y))

# Optimize yaw under bounds
yaw_opt = YawOptimization(fi, bnds=[(0., 25.), (0., 25.), (0., 0.)])
yaw_angles = yaw_opt.optimize()

Expected behavior
The correct behavior is, when using scipy<1.5.0:

=====================================================
Optimizing wake redirection control...
Number of parameters to optimize =  3
=====================================================
  NIT    FC           OBJFUN            GNORM
    1     5    -1.000000E+00     4.618128E-02
    2    10    -1.001534E+00     6.848057E-02
    3    15    -1.041302E+00     1.736367E-01
    4    20    -1.144908E+00     4.190432E-02
    5    25    -1.147140E+00     3.612001E-02
    6    30    -1.150802E+00     1.643277E-02
    7    37    -1.150821E+00     1.476510E-02
    8    52    -1.150821E+00     1.476231E-02
    9    67    -1.150821E+00     1.475954E-02
   10    82    -1.150821E+00     1.475678E-02
   11    97    -1.150821E+00     1.475404E-02
   12   112    -1.150821E+00     1.475131E-02
   13   127    -1.150821E+00     1.474860E-02
   14   142    -1.150821E+00     1.474590E-02
   15   157    -1.150821E+00     1.474322E-02
   16   172    -1.150821E+00     1.474056E-02
   17   187    -1.150821E+00     1.473791E-02
   18   202    -1.150821E+00     1.473528E-02
   19   217    -1.150821E+00     1.473266E-02
   20   232    -1.150821E+00     1.473006E-02
   21   247    -1.150821E+00     1.472747E-02
   22   262    -1.150821E+00     1.472489E-02
   23   277    -1.150821E+00     1.472233E-02
   24   292    -1.150821E+00     1.471979E-02
   25   307    -1.150821E+00     1.471726E-02
   26   322    -1.150821E+00     1.471474E-02
   27   337    -1.150821E+00     1.471223E-02
   28   352    -1.150821E+00     1.470974E-02
   29   367    -1.150821E+00     1.470727E-02
   30   382    -1.150821E+00     1.470480E-02
   31   397    -1.150821E+00     1.470235E-02
   32   412    -1.150821E+00     1.469991E-02
   33   427    -1.150821E+00     1.469749E-02
   34   442    -1.150821E+00     1.469508E-02
   35   457    -1.150821E+00     1.469268E-02
   36   472    -1.150821E+00     1.469030E-02
   37   487    -1.150821E+00     1.468792E-02
   38   502    -1.150821E+00     1.468556E-02
   39   517    -1.150821E+00     1.468321E-02
   40   532    -1.150821E+00     1.468088E-02
   41   547    -1.150821E+00     1.467855E-02
   42   562    -1.150821E+00     1.467624E-02
   43   577    -1.150821E+00     1.467394E-02
   44   592    -1.150821E+00     1.467165E-02
   45   607    -1.150821E+00     1.466938E-02
   46   622    -1.150821E+00     1.466711E-02
   47   637    -1.150821E+00     1.466486E-02
   48   652    -1.150821E+00     1.466262E-02
   49   667    -1.150821E+00     1.466038E-02
   50   682    -1.150821E+00     1.465817E-02
   51   697    -1.150821E+00     1.465596E-02
Iteration limit exceeded    (Exit mode 9)
            Current function value: -1.1508212676552003
            Iterations: 51
            Function evaluations: 697
            Gradient evaluations: 51

Screenshots, if applicable
This is the error that is reported:

=====================================================
Optimizing wake redirection control...
Number of parameters to optimize =  3
=====================================================
/home/bdoekeme/.python_environments/.floris/lib/python3.8/site-packages/scipy/optimize/_numdiff.py:519: RuntimeWarning: invalid value encountered in true_divide
  J_transposed[i] = df / dx
  NIT    FC           OBJFUN            GNORM
    1     4    -1.000000E+00              NAN
Inequality constraints incompatible    (Exit mode 4)
            Current function value: -1.0
            Iterations: 1
            Function evaluations: 4
            Gradient evaluations: 1
No change in controls suggested for this inflow                    condition...

Floris Version
v2.3 (main branch as of to date)

System Information (please complete the following information):

• OS: Ubuntu 18.04 LTS through WSL 1 on a Windows 10 (20H2) machine
• Library versions
  • numpy=1.19.5
  • pytest=6.2.2
  • scipy=1.6.0

Additional context
A potential future-proof solution to this would be to condition what goes into the scipy.optimize.minimize function. Namely, the yaw optimization code could take the user-specified bounds on the floris-side, then determine what turbines have equality constraints, and finally only optimize for the remaining non-equality-constrained turbines. That way, we would present the user with the right value for the number of variables optimized (in this example, it reports 3 but it should be 2), and it provides greater freedom in the optimization methods. I haven't worked out the details really, but it could be something as simple as:

Adding this snippet to reinitialize_opt(...):

self.turbines_opt = np.where(np.diff(bnds) > 0.001)[0]

and then change the cost function to accept a subset of yaw angles, only the ones relevant for optimization:

def _yaw_power_opt(self, yaw_angles_subset):
    yaw_angles = [b[0] for b in self.bnds]
    yaw_angles[self.turbines_opt] = yaw_angles_subset
   
    yaw_angles = self._unnorm(
        np.array(yaw_angles), self.minimum_yaw_angle, self.maximum_yaw_angle
    )
    return (
        -1
        * self.fi.get_farm_power_for_yaw_angle(
            yaw_angles,
            include_unc=self.include_unc,
            unc_pmfs=self.unc_pmfs,
            unc_options=self.unc_options,
        )
        / self.initial_farm_power
    )

[Bartdoekemeijer]

Pull request #245 should fix this issue.