Scipy.optimize.minimize() two arrays for non-linear objective

Question

I have fit a lognormal model, and I would now like to use the coefficients from that model to find the combination of variable values that will provide the greatest response. I intend on introducing some constraints later, but for now I would just like to get the optimization running. The complicating issue is that the model was a mixed methods model, so I have coefficients for each individual in the model. As a result, I need to optimize the variable values for each individual, given their individual coefficients.

My example:

import numpy as np
from scipy.optimize import minimize

def objective(x, beta):
    x1, x2 = x
    beta1, beta2 = beta
    return -1 * np.sum(np.exp(3 + x1*beta1 + x2*beta2))

# initial guesses for variables x1 and x2
n = 2
x1 = np.zeros(n)
x1[0] = 1.0
x1[1] = 2.0

x2 = np.zeros(n)
x2[0] = 3.0
x2[1] = 4.0

x0 = np.vstack((x1,x2))

# the coefficients (weights) for each of n individuals, in each variable
beta1 = np.zeros(n)
beta1[0] = 1.1
beta1[1] = 1.01

beta2 = np.zeros(n)
beta2[0] = 1.1
beta2[1] = 1.01

beta0 = np.vstack((beta1, beta2))

# show initial objective
print('Initial SSE Objective: ' + str(objective(x0, beta0)))  # this works as intended

# but I'm not sure how to specify bounds given my shape
b = (1.0,5.0)  #min = 1, max = 5 on any one channel for any one individual
bnds = (b, b)

# running without bounds
solution = minimize(objective, x0, method='SLSQP',
                    args=beta0)

...gives the following error;


  File "", line 2, in objective
    x1, x2 = x

ValueError: too many values to unpack (expected 2)

Given that the objective function works if I pass it the args, is this a limitation of minimize() ? Is the function unable to take x in shape (2,2)?

Also, I'm having a hard time specifying the bounds correctly...

# running with bounds
solution = minimize(objective, x0, method='SLSQP',
                    args=beta0, bounds=bnds)

...gives the error;

ValueError: operands could not be broadcast together with shapes (4,) (2,) (2,)

I was able to do this with gekko DOCS HERE

    from gekko import GEKKO
    m = GEKKO() # Initialize gekko
    n = 2
    
    # Init the coefficients for each HCP
    alpha_list = np.random.normal(3, 0.1, n)
    beta1_list = np.random.normal(1.01, 0.1, n)
    beta2_list = np.random.normal(1.02, 0.1, n)
    beta3_list = np.random.normal(1.04, 0.1, n)
    
    # Initialize variables
    x1 = [m.Var(value=1,lb=0,ub=4) for i in range(n)]
    x2 = [m.Var(value=1,lb=0,ub=4) for i in range(n)]
    x3 = [m.Var(value=1,lb=0,ub=4) for i in range(n)]
    
    # Init the coefficients
    alpha = [m.Const(value=alpha_list[i]) for i in range(n)]
    beta1 = [m.Const(value=beta1_list[i]) for i in range(n)]
    beta2 = [m.Const(value=beta2_list[i]) for i in range(n)]
    beta3 = [m.Const(value=beta3_list[i]) for i in range(n)]
    
    # Inequality constraints
    m.Equation(m.sum(x1) + m.sum(x2) + m.sum(x3) <= n*10)
    
    m.Obj(-1 * m.sum([m.exp(alpha[i] + x1[i]*beta1[i] + x2[i]*beta2[i] + x3[i]*beta3[i]) for i in range(n)])) # Objective
    m.options.IMODE = 3 # Steady state optimization set to 3, change to 1 for integer
    m.options.MAX_ITER = 1000
    m.solve() # Solve
    print('Results')
    print('x1: ' + str(x1))
    print('x2: ' + str(x2))
    print('x3: ' + str(x3))

...but for my own learning I'd like to be able to do this with scipy too.

user_5422386 · Answer

I'm reporting back in case someone else ever runs in to the same problem. I think the limitation is that the input from minimize to the objective function must be a flat array, so my (2,2) shape wouldn't work (correct me if I'm wrong folks).

So, a workaround is to input a flat array, and then unpack it in the objective function itself. If you hand the function the number of individuals to expect, the function can process the arrays so that each channel gets a 1xN as input to the regression equations we're trying to optimize.

import numpy as np
from scipy.optimize import minimize

def objective(x, beta, n):
    x1, x2 = x.reshape(2,n)
    beta1, beta2 = beta.reshape(2,n)
    return -1 * np.sum(np.exp(3 + x1*beta1 + x2*beta2))

# initial guesses for variables x1 and x2
n = 2
x1 = np.zeros(n)
x1[0] = 1.0
x1[1] = 2.0

x2 = np.zeros(n)
x2[0] = 3.0
x2[1] = 4.0

x0 = np.concatenate((x1,x2))

# the coefficients (weights) for each of n individuals, in each variable
beta1 = np.zeros(n)
beta1[0] = 1.1
beta1[1] = 1.01

beta2 = np.zeros(n)
beta2[0] = 1.1
beta2[1] = 1.01

beta0 = np.concatenate((beta1, beta2))

# show initial objective
print('Initial SSE Objective: ' + str(objective(x0, beta0, 2)))  # this works as intended

# specifying bounds is much easier now we just have a flat array
b = (1.0,5.0)  #min = 1, max = 5 on any one channel for any one individual
bnds = (b , b)*n

# running with bounds
solution = minimize(objective, x0, method='SLSQP',
                    args=(beta0,n), bounds=bnds)
x = solution.x

# show final objective
print('Final SSE Objective: ' + str(objective(x, beta0, n)))

# print solution
print('Solution')
x_sol = x.reshape(2,n)
print('x1 = ' + str(x_sol[0]))
print('x2 = ' + str(x_sol[1]))

As expected, without constraints, the minimize call just maxes out the value of each variable in the equation. My next step will be to put constraints on this function.

Scipy.optimize.minimize() two arrays for non-linear objective

1 Answers