Fixed SNOPT workspace lengths (#197)

Co-authored-by: Neil Wu <[email protected]>

doc/optimizers/SNOPT.rst

@@ -40,6 +40,12 @@ The following are the options set in Python for the wrapper.
   - ``lambda``
   - ``condZHZ``
 - The option ``Start`` corresponds to the value directly passed to the SNOPT kernel, and will be overwritten if another option, e.g. ``Cold start`` is supplied.
+- The default value for ``Total character workspace`` is set internally to the minimum work array length of 500.
+  The default values for ``Total integer workspace`` and ``Total real workspace`` depend on the number of design variables and constraints.
+  They are computed based on recommendations in the SNOPT manual.
+- If SNOPT determines that the default values for ``Total character workspace``, ``Total integer workspace``, or ``Total real workspace`` are too small, the Python wrapper will overwrite the defaults with estimates for the required workspace lengths from SNOPT and initialize the optimizer for a second time.
+  SNOPT might still exit with ``82``, ``83``, or ``84``, but this should automate the storage allocation for most cases.
+  User-specified values are not overwritten.
 
 .. optimizertable:: pyoptsparse.pySNOPT.pySNOPT.SNOPT
    :type: options

examples/large_sparse.py

@@ -2,7 +2,6 @@
 # http://trac.openopt.org/openopt/browser/PythonPackages/FuncDesigner/FuncDesigner/examples/nlpSparse.py
 
 import numpy as np
-from numpy import arange
 from scipy import sparse
 import argparse
 from pyoptsparse import Optimization, OPT
@@ -58,16 +57,16 @@ def large_sparse(optimizer="SNOPT", optOptions=None):
 
     # Design Variables
     optProb.addVar("x", lower=-100, upper=150, value=0)
-    optProb.addVarGroup("y", N, lower=-10 - arange(N), upper=arange(N), value=0)
-    optProb.addVarGroup("z", 2 * N, upper=arange(2 * N), lower=-100 - arange(2 * N), value=0)
+    optProb.addVarGroup("y", N, lower=-10 - np.arange(N), upper=np.arange(N), value=0)
+    optProb.addVarGroup("z", 2 * N, upper=np.arange(2 * N), lower=-100 - np.arange(2 * N), value=0)
     # Constraints
     optProb.addCon("con1", upper=100, wrt=["x"])
     optProb.addCon("con2", upper=100)
     optProb.addCon("con3", lower=4, wrt=["x", "z"])
     optProb.addConGroup(
         "lincon",
         N,
-        lower=2 - 3 * arange(N),
+        lower=2 - 3 * np.arange(N),
         linear=True,
         wrt=["x", "y"],
         jac={"x": np.ones((N, 1)), "y": sparse.spdiags(np.ones(N), 0, N, N)},

examples/tp109.py

@@ -7,19 +7,19 @@
             -(x3-x4+0.550) <= 0
             -(2.25*10**(+6)-x1**2-x8**2) <= 0
             -(2.25*10**(+6)-x2**2-x9**2) <= 0
-            -(x5*x6*sin(-x3-0.250) + +x5*x7*sin(-x4-0.250)+2.0*x5**2*b)*ra+400.0-x1 = 0
-            -((x5*x6*sin(x3-0.250)+x6*x7*sin(x3-x4-0.250)+2.0*x6**2*b)*ra+400.0-x2) = 0
-            -((x5*x7*sin(x4-0.250)+x6*x7*sin(x4-x3-0.250)+2.0*x7**2*b)*ra+881.7790) = 0
-            -(x8+(x5*x6*cos(-x3-0.250)+x5*x7*cos(-x4-0.250)-2.0*x5**2*c)*ra+0.7533*10**(-3)*x5**2-200.00) = 0
-            -(x9+(x5*x6*cos(x3-0.250)+x7*x6*cos(x3-x4-0.250)-2.0*x6**2*c)*ra+0.7533*10**(-3)*x(6)**2-200.00) = 0
-            -((x5*x7*cos(x4-0.250)+x6*x7*cos(x4-x3-0.250)-2.0*x7**2*c)*ra+0.7533*10**(-3)*x7**2-22.9380) = 0
+            -(x5*x6*np.sin(-x3-0.250) + +x5*x7*np.sin(-x4-0.250)+2.0*x5**2*b)*ra+400.0-x1 = 0
+            -((x5*x6*np.sin(x3-0.250)+x6*x7*np.sin(x3-x4-0.250)+2.0*x6**2*b)*ra+400.0-x2) = 0
+            -((x5*x7*np.sin(x4-0.250)+x6*x7*np.sin(x4-x3-0.250)+2.0*x7**2*b)*ra+881.7790) = 0
+            -(x8+(x5*x6*np.cos(-x3-0.250)+x5*x7*np.cos(-x4-0.250)-2.0*x5**2*c)*ra+0.7533*10**(-3)*x5**2-200.00) = 0
+            -(x9+(x5*x6*np.cos(x3-0.250)+x7*x6*np.cos(x3-x4-0.250)-2.0*x6**2*c)*ra+0.7533*10**(-3)*x(6)**2-200.00) = 0
+            -((x5*x7*np.cos(x4-0.250)+x6*x7*np.cos(x4-x3-0.250)-2.0*x7**2*c)*ra+0.7533*10**(-3)*x7**2-22.9380) = 0
             0 <= xi, i = 1,2
             -0.55 <= xi <= 0.55, i = 3,4
             196.0 <= xi <= 252.0, i = 5,6,7
             -400.0 <= xi <= 800.0, i = 8,9
     where   a = 50.176
-            b = sin(0.25)
-            c = cos(0.25)
+            b = np.sin(0.25)
+            c = np.cos(0.25)
             ra = 1.0/50.176
 
 
@@ -29,7 +29,6 @@
 
 import numpy as np
 import argparse
-from numpy import sin, cos
 from pyoptsparse import Optimization, OPT
 
 parser = argparse.ArgumentParser()
@@ -43,37 +42,45 @@ def objfunc(xdict):
     x = xdict["xvars"]
 
     a = 50.1760
-    b = sin(0.250)
-    c = cos(0.250)
+    b = np.sin(0.250)
+    c = np.cos(0.250)
     funcs = {}
     funcs["obj"] = 3.0 * x[0] + (1e-6) * x[0] ** 3 + 0.522074e-6 * x[1] ** 3 + 2 * x[1]
     fcon = np.zeros(10, "D")
     fcon[0] = 2250000 - x[0] ** 2 - x[7] ** 2
     fcon[1] = 2250000 - x[1] ** 2 - x[8] ** 2
-    fcon[2] = x[4] * x[5] * sin(-x[2] - 0.25) + x[4] * x[6] * sin(-x[3] - 0.25) + 2 * b * x[4] ** 2 - a * x[0] + 400 * a
+    fcon[2] = (
+        x[4] * x[5] * np.sin(-x[2] - 0.25) + x[4] * x[6] * np.sin(-x[3] - 0.25) + 2 * b * x[4] ** 2 - a * x[0] + 400 * a
+    )
     fcon[3] = (
-        x[4] * x[5] * sin(x[2] - 0.25) + x[5] * x[6] * sin(x[2] - x[3] - 0.25) + 2 * b * x[5] ** 2 - a * x[1] + 400 * a
+        x[4] * x[5] * np.sin(x[2] - 0.25)
+        + x[5] * x[6] * np.sin(x[2] - x[3] - 0.25)
+        + 2 * b * x[5] ** 2
+        - a * x[1]
+        + 400 * a
+    )
+    fcon[4] = (
+        x[4] * x[6] * np.sin(x[3] - 0.25) + x[5] * x[6] * np.sin(x[3] - x[2] - 0.25) + 2 * b * x[6] ** 2 + 881.779 * a
     )
-    fcon[4] = x[4] * x[6] * sin(x[3] - 0.25) + x[5] * x[6] * sin(x[3] - x[2] - 0.25) + 2 * b * x[6] ** 2 + 881.779 * a
     fcon[5] = (
         a * x[7]
-        + x[4] * x[5] * cos(-x[2] - 0.25)
-        + x[4] * x[6] * cos(-x[3] - 0.25)
+        + x[4] * x[5] * np.cos(-x[2] - 0.25)
+        + x[4] * x[6] * np.cos(-x[3] - 0.25)
         - 200 * a
         - 2 * c * x[4] ** 2
         + 0.7533e-3 * a * x[4] ** 2
     )
     fcon[6] = (
         a * x[8]
-        + x[4] * x[5] * cos(x[2] - 0.25)
-        + x[5] * x[6] * cos(x[2] - x[3] - 0.25)
+        + x[4] * x[5] * np.cos(x[2] - 0.25)
+        + x[5] * x[6] * np.cos(x[2] - x[3] - 0.25)
         - 2 * c * x[5] ** 2
         + 0.7533e-3 * a * x[5] ** 2
         - 200 * a
     )
     fcon[7] = (
-        x[4] * x[6] * cos(x[3] - 0.25)
-        + x[5] * x[6] * cos(x[3] - x[2] - 0.25)
+        x[4] * x[6] * np.cos(x[3] - 0.25)
+        + x[5] * x[6] * np.cos(x[3] - x[2] - 0.25)
         - 2 * c * x[6] ** 2
         - 22.938 * a
         + 0.7533e-3 * a * x[6] ** 2

pyoptsparse/__init__.py

@@ -1,4 +1,4 @@
-__version__ = "2.4.0"
+__version__ = "2.4.1"
 
 from .pyOpt_history import History
 from .pyOpt_variable import Variable

pyoptsparse/pyNLPQLP/pyNLPQLP.py

@@ -260,7 +260,7 @@ def nlgrad(m, me, mmax, n, f, g, df, dg, x, active, wa):
                 self.hist.close()
 
             # Store Results
-            inform = np.asscalar(ifail)
+            inform = ifail.item()
             sol_inform = {}
             sol_inform["value"] = inform
             sol_inform["text"] = self.informs[inform]

pyoptsparse/pyOpt_optimizer.py

@@ -82,6 +82,9 @@ def __init__(
         # Store the Jacobian conversion maps
         self._jac_map_csr_to_csc = None
 
+        # Initialize metadata
+        self.metadata = {}
+
     def _clearTimings(self):
         """Clear timings and call counters"""
         self.userObjTime = 0.0

pyoptsparse/pyPSQP/pyPSQP.py

@@ -260,7 +260,7 @@ def pdcon(n, k, x, dg):
             self.optProb.comm.bcast(-1, root=0)
 
             # Store Results
-            inform = np.asscalar(iterm)
+            inform = iterm.item()
             if inform < 0 and inform not in self.informs:
                 inform = -10
             sol_inform = {}

pyoptsparse/pySLSQP/pySLSQP.py

@@ -257,7 +257,7 @@ def slgrad(m, me, la, n, f, g, df, dg, x):
             self.optProb.comm.bcast(-1, root=0)
 
             # Store Results
-            inform = np.asscalar(mode)
+            inform = mode.item()
             sol_inform = {}
             sol_inform["value"] = inform
             sol_inform["text"] = self.informs[inform]

pyoptsparse/pySNOPT/pySNOPT.py

@@ -54,6 +54,9 @@ def __init__(self, raiseError=True, options={}):
             "Start": [str, ["Cold", "Warm"]],  # used in call to snkerc, the option "Cold Start" etc are NOT allowed
             "Derivative level": [int, 3],
             "Proximal iterations limit": [int, 10000],  # very large # to solve proximal point problem to optimality
+            "Total character workspace": [int, None],  # lencw: 500
+            "Total integer workspace": [int, None],  # leniw: 500 + 100 * (ncon + nvar)
+            "Total real workspace": [int, None],  # lenrw: 500 + 200 * (ncon + nvar)
             "Save major iteration variables": [
                 list,
                 ["step", "merit", "feasibility", "optimality", "penalty"],
@@ -328,14 +331,28 @@ def __call__(
                     raise Error("Failed to properly open %s, ierror = %3d" % (SummFile, ierror))
 
             # Calculate the length of the work arrays
-            # --------------------------------------
+            # ---------------------------------------
             nvar = self.optProb.ndvs
-            lencw = 500
-            leniw = 500 + 100 * (ncon + nvar)
-            lenrw = 500 + 200 * (ncon + nvar)
-
-            self.setOption("Total integer workspace", leniw)
-            self.setOption("Total real workspace", lenrw)
+            lencw = self.getOption("Total character workspace")
+            leniw = self.getOption("Total integer workspace")
+            lenrw = self.getOption("Total real workspace")
+
+            # Set flags to avoid overwriting user-specified lengths
+            checkLencw = lencw is None
+            checkLeniw = leniw is None
+            checkLenrw = lenrw is None
+
+            # Set defaults
+            minWorkArrayLength = 500
+            if lencw is None:
+                lencw = minWorkArrayLength
+                self.setOption("Total character workspace", lencw)
+            if leniw is None:
+                leniw = minWorkArrayLength + 100 * (ncon + nvar)
+                self.setOption("Total integer workspace", leniw)
+            if lenrw is None:
+                lenrw = minWorkArrayLength + 200 * (ncon + nvar)
+                self.setOption("Total real workspace", lenrw)
 
             cw = np.empty((lencw, 8), "c")
             iw = np.zeros(leniw, np.intc)
@@ -353,20 +370,32 @@ def __call__(
             # Set the options into the SNOPT instance
             self._set_snopt_options(iPrint, iSumm, cw, iw, rw)
 
+            # Estimate workspace storage requirement
             mincw, miniw, minrw, cw = snopt.snmemb(iExit, ncon, nvar, neA, neGcon, nnCon, nnJac, nnObj, cw, iw, rw)
 
-            if (minrw > lenrw) or (miniw > leniw) or (mincw > lencw):
-                if mincw > lencw:
-                    lencw = mincw
-                    cw = np.array((lencw, 8), "c")
-                    cw[:] = " "
-                if miniw > leniw:
-                    leniw = miniw
-                    iw = np.zeros(leniw, np.intc)
-                if minrw > lenrw:
-                    lenrw = minrw
-                    rw = np.zeros(lenrw, np.float)
-
+            # This flag is set to True if any of the lengths are overwritten
+            lengthsChanged = False
+
+            # Overwrite lengths if the defaults are too small
+            if checkLencw and mincw > lencw:
+                lencw = mincw
+                self.setOption("Total character workspace", lencw)
+                cw = np.array((lencw, 8), "c")
+                cw[:] = " "
+                lengthsChanged = True
+            if checkLeniw and miniw > leniw:
+                leniw = miniw
+                self.setOption("Total integer workspace", leniw)
+                iw = np.zeros(leniw, np.intc)
+                lengthsChanged = True
+            if checkLenrw and minrw > lenrw:
+                lenrw = minrw
+                self.setOption("Total real workspace", lenrw)
+                rw = np.zeros(lenrw, np.float)
+                lengthsChanged = True
+
+            # Initialize SNOPT again if any of the lengths were overwritten
+            if lengthsChanged:
                 snopt.sninit(iPrint, iSumm, cw, iw, rw)
 
                 # snInit resets all the options to the defaults.
@@ -435,7 +464,7 @@ def __call__(
                 snopt.closeunit(self.getOption("iSumm"))
 
             # Store Results
-            inform = np.asscalar(inform)
+            inform = inform.item()
             sol_inform = {}
             sol_inform["value"] = inform
             sol_inform["text"] = self.informs[inform]