4

I am using scipy.optimize.minimize for a small optimization problem with 9 free variables. My objective function is basically a wrapper around another function, and if I evaluate my objective function, the return type is 'numpy.float32'... which is a scalar? However, I am getting the following error when attempting to use the minimize function:

raise ValueError("Objective function must return a scalar")
ValueError: Objective function must return a scalar

Is it not possible to wrap the objective function around another function? The other function arguments are globally declared, but if this is not ok, I can hardcode them into the beam_shear function.

Relevant code snippets:

from numpy import array, shape, newaxis, isnan, arange, zeros, dot, linspace
from numpy import pi, cross, tile, arccos,sin, cos, sum, atleast_2d, asarray, float32, ones

from numpy import sum, reshape
from scipy.optimize import minimize

def normrow(A):
    A = atleast_2d(asarray(A, dtype=float32))
    return (sum(A ** 2, axis=1) ** 0.5).reshape((-1, 1))

def beam_shear(xyz, pt0, pt1, pt2, x):

    # will not work for overlapping nodes...
    s         = zeros((len(xyz), 3))
    xyz_pt0   = xyz[pt0, :]
    xyz_pt1   = xyz[pt1, :]
    xyz_pt2   = xyz[pt2, :]
    e01       = xyz_pt1 - xyz_pt0
    e12       = xyz_pt2 - xyz_pt1
    e02       = xyz_pt2 - xyz_pt0
    trip_norm = cross(e01, e12)
    mu        = 0.5 * (xyz_pt2 - xyz_pt1)
    l01       = normrow(e01)
    l12       = normrow(e12)
    l02       = normrow(e02)
    l_tn      = normrow(trip_norm)
    l_mu      = normrow(mu)
    a         = arccos((l01**2 + l12**2 - l02**2) / (2 * l01 * l12))
    k         = 2 * sin(a) / l02 # discrete curvature
    ex        = trip_norm / tile(l_tn, (1, 3))
    ez        = mu / tile(l_mu, (1, 3))
    ey        = cross(ez, ex)
    kb        = tile(k / l_tn, (1, 3)) * trip_norm
    kx        = tile(sum(kb * ex, 1)[:, newaxis], (1, 3)) * ex
    m         = x * kx
    cma       = cross(m, e01)
    cmb       = cross(m, e12)
    ua        = cma / tile(normrow(cma), (1, 3))
    ub        = cmb / tile(normrow(cmb), (1, 3))
    c1        = cross(e01, ua)
    c2        = cross(e12, ub)
    l_c1      = normrow(c1)
    l_c2      = normrow(c2)
    ms        = sum(m**2, 1)[:, newaxis]
    Sa        = ua * tile(ms * l_c1 / (l01 * sum(m * c1, 1)[:, newaxis]), (1, 3))
    Sb        = ub * tile(ms * l_c2 / (l12 * sum(m * c2, 1)[:, newaxis]), (1, 3))
    Sa[isnan(Sa)] = 0
    Sb[isnan(Sb)] = 0
    s[pt0, :] += Sa
    s[pt1, :] -= Sa + Sb
    s[pt2, :] += Sb
    return s

def cross_section_obj(x):
    s = beam_shear(xyz, pt0, pt1, pt2, x)
    l_s = normrow(s)
    val = sum(l_s)
    return val

xyz = array([[ 0, 0., 0.],
        [ 0.16179067,  0.24172157,  0.],
        [ 0.33933063,  0.47210142,  0.],
        [ 0.53460629,  0.68761389,  0.],
        [ 0.75000537,  0.88293512, 0.],
        [ 0.98816469,  1.04956383, 0.],
        [ 1.25096091,  1.17319961,  0.],
        [ 1.5352774,  1.22977204,  0.],
        [ 1.82109752,  1.18695051,  0.],
        [ 2.06513705, 1.03245579,  0.],
        [ 2.23725517,  0.79943842,  0.]])

pt0 = array([0, 1, 2, 3, 4, 5, 6, 7, 8])
pt1 = array([1, 2, 3, 4, 5, 6, 7, 8, 9])
pt2 = array([2, 3, 4, 5, 6, 7, 8, 9, 10])
EIx = (ones(len(pt1)) * 12.75).reshape(-1, 1)

bounds = []
for i in range(len(EIx)):
    bounds.append((EIx[i][0], EIx[i][0] * 100))


print(type(cross_section_obj(EIx)))
res = minimize(cross_section_obj, EIx, method='SLSQP', bounds=bounds)

As mentioned before:

print(type(cross_section_obj(EIx)))

returns:

<type 'numpy.float32'>

EIx is the set of initial values for the optimization, which is an array of shape (9, 1).

stagermane
  • 1,003
  • 2
  • 12
  • 29
  • there is a parenthesis missing at the start of `sum` which is assigned to `l_s`. additionally you don't need to reshape if you're just going to sum to a scalar – atomsmasher Oct 19 '17 at 15:05
  • 1
    Add an example we can run. Yes... the shapes might be a problem. For debugging: add ```print(type(val))``` in the line before ```return val```. – sascha Oct 19 '17 at 15:08
  • Also: ```10 free variables``` != ```EIx is the set of initial values for the optimization, which is an array of shape (9, 1).``` – sascha Oct 19 '17 at 15:14
  • Here's an example, hope this clarifies things somewhat! You were also right about the free variables (should be 9 for this case), and the sum parentheses have been fixed. – stagermane Oct 19 '17 at 16:15
  • That code has even an error in beam_shear (which we don't know). Impossible to work on the problems with scipy then. And you still did not do my debugging-recommendation (maybe you think it's covered by your print; it's not! at least not in general). – sascha Oct 19 '17 at 16:16
  • I tested your debugging suggestion, and it also returns `.` Beam shear works as I expected, I'm using this snippet in several other situations.... Is the point that scipy cannot work with a function inside a function? – stagermane Oct 19 '17 at 16:28
  • scipy does indeed work with a function inside a function just try this https://repl.it/MkJt/2 – atomsmasher Oct 19 '17 at 16:51
  • beam_shear throws an error (in the code above!), no discussion. That's enough for me that i can't debug it for you. I can imagine the overall problem, but touching it will (probably) break every other code of yours. – sascha Oct 19 '17 at 16:59
  • @sascha beam_shear does not throw an error for me, I have not edited it. version is python 3.6 – atomsmasher Oct 19 '17 at 17:02
  • numpy-version, scipy-version (1.13.1 / 0.19.0 here)? – sascha Oct 19 '17 at 17:03
  • scipy: 0.19.0, numpy: 1.12.1 – atomsmasher Oct 19 '17 at 17:05

2 Answers2

5

You may want to look at Utilizing scipy.optimize.minimize with multiple variables of different shapes. What is important to understand is that if you want to use minimize with arrays, you should pass in the flattened version and then reshape. For that reason I always include the desired shape as one of the arguments to the minimize function. In your case I would do something like so:

def cross_section_obj(x, *args):
    xyz, pt0, pt1, pt2, shape = args
    x = x.reshape(shape)
    s = beam_shear(xyz, pt0, pt1, pt2, x)
    l_s = normrow(s)
    val = sum(l_s)
    return val

Then your minimize call would change like so:

res = minimize(cross_section_obj, EIx.flatten(), method='SLSQP',
               bounds=bounds, args=(xyz, pt0, pt1, pt2, EIx.shape))
Grr
  • 15,553
  • 7
  • 65
  • 85
1

Your array EIx of parameter values is two-dimensional; it has shape (9,1). During the minimization process this array becomes one-dimensional after the first iteration. However, your function beam_shear does not work if x is one-dimensional.

You can fix this by changing cross_section to the following:

def cross_section_obj(x):
    x = x.reshape((-1,1))
    s = beam_shear(xyz, pt0, pt1, pt2, x)
    l_s = normrow(s)
    val = sum(l_s)
    return val

The code then runs, but of course you need to double-check if that is what you actually want to calculate.

user8153
  • 4,049
  • 1
  • 9
  • 18