no solution for simple integer valued objective function

[josh beal]

Hi Everyone,

I was wondering if I could get some help with a problem (working code in python posted below). This is an integer valued problem so I'm using APOPT. Now, when I run gekko with a dummy objective, m.Obj(0), the solver is able to identify feasible solutions, usually within a few seconds, up to a minute or so. However, when I include a (relatively simple) objective function, denoted m.Obj(m.sum(mix_score)), the code runs for over an hour, and produces an error. This is surprising, because I would at least expect the solver to identify one feasible solution and return the objective function value (min sum of mix_score). Note that mix_score[i] is either 1 or 2 for this example. It is calculated from count[i] which is computed from the m.if3() a few lines above. I would just like to find the the optimal x[i][j] that lead to a min value for the sum of these mix_score[i] values. Any help would be much appreciated! Thank you in advance. Code below (currently m.Obj(0) is uncommented to identify feasible solutions. To run m.Obj(m.sum(mix_score)), simply comment out m.Obj(0) and uncomment the two lines where mix_score is defined as well as m.Obj(m.sum(mix_score)) ).

# ## Imports

from gekko import GEKKO

import numpy as np

## Set fixed variables

h_matrix = np.array(

[

[119.4, 119.4, 119.4, 119.4, 119.4, 119.4, 111.4, 111.4, 111.4, 111.4],

[119.4, 119.4, 119.4, 119.4, 119.4, 119.4, 111.4, 111.4, 111.4, 111.4],

[119.4, 119.4, 119.4, 119.4, 119.4, 111.4, 111.4, 111.4, 111.4, 111.4]

]

)

z_matrix = np.array(

[

[383.91, 383.91, 383.91, 383.91, 383.91, 383.91, 254.49, 254.49, 254.49, 254.49],

[383.91, 383.91, 383.91, 383.91, 383.91, 383.91, 254.49, 254.49, 254.49, 254.49],

[383.91, 383.91, 383.91, 383.91, 383.91, 254.49, 254.49, 254.49, 254.49, 254.49]

]

)

w = np.array([47.93, 66.37])

h = np.array([12.10, 8.6])

t = np.array([104, 48])

## Reshape the matrices to make calulations easier

h_matrix_reshaped = h_matrix.reshape(30,1)

z_matrix_reshaped = z_matrix.reshape(30,1)

# ## Initialize

# Initialize the model

m = GEKKO(remote=False)

## Fixed variables

h_constants = [m.Param(value=h_matrix_reshaped[i][0]) for i in range(30)]

z_constants = [m.Param(value=z_matrix_reshaped[i][0]) for i in range(30)]

w_base = [m.Param(value=w[i]) for i in range(w.shape[0])]

h_base = [m.Param(value=h[i]) for i in range(h.shape[0])]

t_base = [m.Param(value=t[i]) for i in range(t.shape[0])]

h_cm = m.Param(value=220)

ho_cm = m.Param(value=14.4)

w_kg = m.Param(value=21.2)

s_constraint = m.Param(value=40)

# ### Set up x var (main integer variable)

## Initialize x array

x = np.empty((30, t.shape[0]), dtype=object)

for i in range(30):

for j in range(t.shape[0]):

x[i][j] = m.Var(value=0, lb=0, integer=True)

# ### Set up Constraints

## Total constraint

for j in range(len(t_base)):

t_contraints = sum(x[i][j] for i in range(30))

m.Equation(t_contraints == t_base[j])

## Weight constraints

for i in range(30):

w_constraints = sum(x[i][j]*w_base[j] for j in range(len(w_base))) + w_kg

m.Equation(w_constraints <= z_constants[i])

## Height constraints

for i in range(30):

h_constraints = sum(x[i][j]*h_base[j] for j in range(len(h_base))) + ho_cm + h_constants[i]

m.Equation(h_constraints <= h_cm)

## Neighbor constraints

for i in range(9):

# set neighbor constraints horizontally over first row

neighbor_constraints_1 = m.abs3((sum(x[i][j]*h_base[j] for j in range(len(h_base))) +

h_constants[i]) - (sum(x[i+1][j]*h_base[j] for j in range(len(h_base))) +

h_constants[i+1]))

m.Equation(neighbor_constraints_1 <= s_constraint)

for i in range(10,19):

# set neighbor constraints horizontally over second row

neighbor_constraints_2 = m.abs3((sum(x[i][j]*h_base[j] for j in range(len(h_base))) +

h_constants[i]) - (sum(x[i+1][j]*h_base[j] for j in range(len(h_base))) +

h_constants[i+1]))

m.Equation(neighbor_constraints_2 <= s_constraint)

for i in range(20,29):

# set neighbor constraints horizontally over second row

neighbor_constraints_3 = m.abs3((sum(x[i][j]*h_base[j] for j in range(len(h_base))) +

h_constants[i]) - (sum(x[i+1][j]*h_base[j] for j in range(len(h_base))) +

h_constants[i+1]))

m.Equation(neighbor_constraints_3 <= s_constraint)

for i in range(10):

# set neighbor constrainst vertically A with B

neighbor_constraints_4 = m.abs3((sum(x[i][j]*h_base[j] for j in range(len(h_base))) +

h_constants[i]) - (sum(x[i+10][j]*h_base[j] for j in range(len(h_base))) +

h_constants[i+10]))

m.Equation(neighbor_constraints_4 <= s_constraint)

for i in range(10,20):

# set neighbor constrainst vertically B with C

neighbor_constraints_5 = m.abs3((sum(x[i][j]*h_base[j] for j in range(len(h_base))) +

h_constants[i]) - (sum(x[i+10][j]*h_base[j] for j in range(len(h_base))) +

h_constants[i+10]))

m.Equation(neighbor_constraints_5 <= s_constraint)

# ### Mix Score section below ##################

## Create a binary variable/array b that identifies if x[i][j] is non-zero or not

## We will use the count of these b[i][j] values in our objective function

## Constraint to set b[i][k] = 1 if x[i][k] > 0, and 0 otherwise

## Use if3 to set b directly based on x values

b = np.empty((30, len(t_base)), dtype=object)

epsilon = 1e-2 # Small margin allows floating-point considerations

for i in range(30):

for j in range(len(t_base)):

b[i][j] = m.if3(x[i][j] - epsilon, 0, 1)

# # Calculation of count(i) for each row

counts = [m.Intermediate(m.sum(b[i])) for i in range(30)]

# # x_sums for sum of each row in x

#x_sums = [m.Intermediate(m.sum(x[i])) for i in range(30)]

### Mix Score section above #############################

# ## Run Solver ##

## Set a dummy objective just to identify solutions that are feasible

m.Obj(0)

### Define the main objective function

#mix_score = [counts[i] for i in range(30)] # uncomment to run mix score objective

#m.Obj(m.sum(mix_score)) #uncomment to run mix score objective

# Set the solver options

m.options.SOLVER = 1 # APOPT solver for non-linear programs

# Increase max iterations because we don't care much about time

m.solver_options = ['minlp_gap_tol 1.0e-4',\

'minlp_maximum_iterations 50000',\

'minlp_max_iter_with_int_sol 40000']

## Solve

m.solve(disp=True)

[John Hedengren]

Josh,

 

Thanks for asking your question on StackOverflow. I’ve posted an answer to the question: https://stackoverflow.com/questions/78643391/adding-conditional-variable-in-gekko-leads-to-no-solution/78674665#78674665

 

The m.max3() function is much more reliable than the m.max2() function. Because you are using m.abs3() and m.if3() functions in your optimization problem, it isn’t a problem to add a few more binary variables by using the m.max3() function. The MPCC version of the function has troubles converging, especially when the solution is at the switching location where there is a saddle point that is a challenging for nonlinear optimizers.

 

-John Hedengren

 

 

 

--
--
APMonitor user's group e-mail list.
- To post a message, send email to apmo...@googlegroups.com
- To unsubscribe, send email to apmonitor+...@googlegroups.com
- Visit this group at http://groups.google.com/group/apmonitor
---
You received this message because you are subscribed to the Google Groups "apmonitor" group.
To unsubscribe from this group and stop receiving emails from it, send an email to apmonitor+...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/apmonitor/3c4f18df-8168-43a9-8d20-4ab01164dbcen%40googlegroups.com.

[josh beal]

This worked like a charm, John. Thanks a lot!

You received this message because you are subscribed to a topic in the Google Groups "apmonitor" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/apmonitor/jK3ikwJiek4/unsubscribe.
To unsubscribe from this group and all its topics, send an email to apmonitor+...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/apmonitor/CH3PR22MB5605CD9741CFCB62F818B8DAEED62%40CH3PR22MB5605.namprd22.prod.outlook.com.