Scipy.optimize.minimize objective function ValueEr

Question

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).

Answer 1

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))

Answer 2

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.