Correct objective function for optimising multivariate logistic functions

Question

I have the following time series dataset:

Input1, Input2, Input3, Output

In the following plots, you can see each input separately plotted against the output.

Input1/Output

Input2/Output

Input3/Output

For each input column, I used Scipy Optimize - curve_fit function to fit the data with the following expression:

def sigmoid(x, a, b, L):
     y = L / (1 + np.exp(-b*(x-a)))
     return y

For each input column, I have the following constants (a, b, L):

Input1/Output: [2.86115648e+04 4.52333694e-05 6.49423842e-01]

Input2/Output: [6.15077795e+03 2.00771121e-04 6.02374706e-01]

Input3/Output: [3.90539815e+03, 7.80392947e-04, 5.77858431e-01]

Given an (arbitrary/example) constraint that: Input1 + Input2 + Input3 < 120,000

I want to maximise each input to gain the maximum output.

Can I express the objective function as simply the sum of each individual function? In other words, can my objective function be expressed as follows:

z = (L1 / (1 + np.exp(-b1*(x1-a1)))) + (L2 / (1 + np.exp(-b2*(x2-a2)))) +(L3 / (1 + np.exp(-b3*(x3-a3)))) + (L4 / (1 + np.exp(-b4*(x4-a4))))

Referring to the Gekko code:

# Equations
m.Equation(x1+x2+x3+x4<=120000)
m.Obj(L1 / (1 + np.exp(-b1*(x1-a1)))) + (L2 / (1 + np.exp(-b2*(x2-a2)))) + (L3 / (1 + np.exp(-b3*(x3-a3)))) + (L4 / (1 + np.exp(-b4*(x4-a4)))) # Objective
m.solve(disp=False) # Solve

Using the code below I was seeing the following error:

 m.options.SOLVER=1  # APOPT is an MINLP solver

# optional solver settings with APOPT

a1 = 1662.2548899281423
b1 = 0.0008133406683575547
L1 = 0.5456713295037908

a2 = 5273.826922188703 
b2 = 0.00024498561094814383 
L2 = 0.5730871268164875

a3 = 3836.9976232989725
b3 = 0.0007892890342618781
L3 = 0.5697505125863976

a4 = 27077.484569050346
b5 = 0.00004580885182095956
L6 = 0.6429567105559689

# Initialize variables
x1 = m.Var(value=1.00,lb=1.00)
x2 = m.Var(value=5.00,lb=1.00)
x3 = m.Var(value=5.00,lb=1.00)
x4 = m.Var(value=5.00,lb=1.00)

# Equations
m.Equation(x1+x2+x3+x4<=120000)
m.Maximize( (L1 / (1.00 + np.exp(-b1*(x1-a1)))) + (L2 / (1.00 + np.exp(-b2*(x2-a2)))) + (L3 / (1.00 + np.exp(-b3*(x3-a3)))) + (L4 / (1.00 + np.exp(-b4*(x4-a4))))) # Objective
m.solve(disp=False) # Solve

print('Results')
print('x1: ' + str(x1.value))
print('x2: ' + str(x2.value))
print('x3: ' + str(x3.value))
print('x4: ' + str(x4.value))
print('Objective: ' + str(m.options.objfcnval))

Error:

---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
/usr/local/lib/python3.7/site-packages/gekko/gk_operators.py in __getattr__(self, name)
     35         else:
---> 36             raise AttributeError(name)
     37     #%%Operator overloading for building functions

AttributeError: exp

The above exception was the direct cause of the following exception:

TypeError                                 Traceback (most recent call last)
<ipython-input-136-404bfd33128d> in <module>
     31 m.Equation(x1+x2+x3+x4<=1000000)
     32 
---> 33 m.Maximize( (L1 / (1.00 + np.exp(-b1*(x1-a1)))) + (L2 / (1.00 + np.exp(-b2*(x2-a2)))) + (L3 / (1.00 + np.exp(-b3*(x3-a3)))) + (L4 / (1.00 + np.exp(-b4*(x4-a4))))) # Objective
     34 m.solve(disp=False) # Solve
     35 print('Results')

TypeError: loop of ufunc does not support argument 0 of type GK_Operators which has no callable exp method

Rather than using the Numpy exponential function, using the Gekko exp function resolves the issue as per the below:

from gekko import GEKKO
m = GEKKO() # Initialize gekko
m.options.SOLVER=1  # APOPT is an MINLP solver

# optional solver settings with APOPT

a1 = 1662.2548899281423
b1 = 0.0008133406683575547
L1 = 0.5456713295037908

a2 = 5273.826922188703 
b2 = 0.00024498561094814383 
L2 = 0.5730871268164875

a3 = 3836.9976232989725
b3 = 0.0007892890342618781
L3 = 0.5697505125863976

a4 = 27077.484569050346
b5 = 0.00004580885182095956
L6 = 0.6429567105559689

# Initialize variables
x1 = m.Var(value=1.00,lb=1.00)
x2 = m.Var(value=5.00,lb=1.00)
x3 = m.Var(value=5.00,lb=1.00)
x4 = m.Var(value=5.00,lb=1.00)

# Equations
m.Equation(x1+x2+x3+x4>=1000)
m.Maximize( (L1 / (1.00 + m.exp(-b1*(x1.value-a1)))) + (L2 / (1.00 + m.exp(-b2*(x2.value-a2)))) + (L3 / (1.00 + m.exp(-b3*(x3.value-a3)))) + (L4 / (1.00 + m.exp(-b4*(x4.value-a4))))) # Objective
m.solve(disp=False) # Solve

print('Results')
print('x1: ' + str(x1.value))
print('x2: ' + str(x2.value))
print('x3: ' + str(x3.value))
print('x4: ' + str(x4.value))
print('Objective: ' + str(m.options.objfcnval))

Answer 1

Yes, you can write the objective function as you presented. That refers to a 'Weighted Sum' approach for multiobjective optimization.
One thing, if you want to "maximize" your objective function, you might want to put a negative sign on the whole objective function equation.

m.Obj(-(L1 / (1 + np.exp(-b1*(x1-a1)))) + (L2 / (1 + np.exp(-b2*(x2-a2)))) +(L3 / (1 + np.exp(-b3*(x3-a3))))) # Objective

Or, you can use a Gekko built-in function "maximize" instead.

m.maximize(L1 / (1 + np.exp(-b1*(x1-a1)))) + (L2 / (1 + np.exp(-b2*(x2-a2)))) +(L3 / (1 + np.exp(-b3*(x3-a3)))) # Objective