J

#### jim

##### Guest

`m.Ojb(0)`

just to test for feasibility, the solver is able to find a feasible solution. However, when I add a conditional variable `b`

(see Mix Score section below), that is not subject to any constraints, the solver fails to find a solution. Here's is some fully functioning test code below:
Code:

```
# ## Imports
from gekko import GEKKO
import numpy as np
## Set fixed variabls
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 contraints
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
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] / m.max2(x_sums[i], 1e-3) for i in range(30)]
# m.Obj(m.sum(mix_score))
# 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)
```

In the Mix score section, you'll see where I would like to add in conditional variables defined as the

`b[i][j]`

that are 1 or 0 depending on the value of `x[i][j]`

(1 if `x[i][j]`

is non-zero, and 0 otherwise). I would like to apply this binary variable to a count function I'm using in the objective function that is commented out`# diversity_score = [counts[i] / m.max2(x_sums[i], 1e-3) for i in range(30)]`

`# m.Obj(m.sum(diversity_score))`

I'm surprised when I implement the variable

`b`

(so simply un-comment the code denoted Mix Score section below and above), the solver fails to find a feasible solution, even when running the dummy objective `m.Obj(0)`

It is as if the way I'm defining `b`

is actually adding a constraint and reducing the feasibility space. Any help would be much appreciated. Thanks!<p>I'm using gekko to optimize a certain function. When I use a dummy objective function like <code>m.Ojb(0)</code> just to test for feasibility, the solver is able to find a feasible solution. However, when I add a conditional variable <code>b</code> (see Mix Score section below), that is not subject to any constraints, the solver fails to find a solution. Here's is some fully functioning test code below:</p>

<pre><code>

# ## Imports

from gekko import GEKKO

import numpy as np

## Set fixed variabls

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

*[0]) for i in range(30)]*

z_constants = [m.Param(value=z_matrix_reshaped

z_constants = [m.Param(value=z_matrix_reshaped

*[0]) for i in range(30)]*

w_base = [m.Param(value=ww_base = [m.Param(value=w

*) for i in range(w.shape[0])]*

h_base = [m.Param(value=hh_base = [m.Param(value=h

*) for i in range(h.shape[0])]*

t_base = [m.Param(value=tt_base = [m.Param(value=t

*) 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]):

xh_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

*[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# ### Set up Constraints

## Total constraint

for j in range(len(t_base)):

t_contraints = sum(x

*[j] for i in range(30))*

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

## Weight contraints

for i in range(30):

w_constraints = sum(xm.Equation(t_contraints == t_base[j])

## Weight contraints

for i in range(30):

w_constraints = sum(x

*[j]*w_base[j] for j in range(len(w_base))) + w_kg*

m.Equation(w_constraints <= z_constantsm.Equation(w_constraints <= z_constants

*)*

## Height constraints

for i in range(30):

h_constraints = sum(x## Height constraints

for i in range(30):

h_constraints = sum(x

*[j]*h_base[j] for j in range(len(h_base))) + ho_cm*

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(xm.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

*[j]*h_base[j] for j in range(len(h_base))) +*

h_constantsh_constants

*) - (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(xh_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

*[j]*h_base[j] for j in range(len(h_base))) +*

h_constantsh_constants

*) - (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(xh_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

*[j]*h_base[j] for j in range(len(h_base))) +*

h_constantsh_constants

*) - (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(xh_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

*[j]*h_base[j] for j in range(len(h_base))) +*

h_constantsh_constants

*) - (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(xh_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

*[j]*h_base[j] for j in range(len(h_base))) +*

h_constantsh_constants

*) - (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 xh_constants[i+10]))

m.Equation(neighbor_constraints_5 <= s_constraint)

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

## Create a binary variable/array b that identifies if x

*[j] is non-zero or not*

## We will use the count of these b## We will use the count of these b

*[j] values in our objective function*

## Constraint to set b## Constraint to set b

*[k] = 1 if x**[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## 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

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

# # Calculation of count(i) for each row

# counts = [m.Intermediate(m.sum(b# # Calculation of count(i) for each row

# counts = [m.Intermediate(m.sum(b

*)) for i in range(30)]*

# # x_sums for sum of each row in x

# x_sums = [m.Intermediate(m.sum(x# # x_sums for sum of each row in x

# x_sums = [m.Intermediate(m.sum(x

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

*/ m.max2(x_sums**, 1e-3) for i in range(30)]*

# m.Obj(m.sum(mix_score))

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

</code></pre>

<p>In the Mix score section, you'll see where I would like to add in conditional variables defined as the <code>b# m.Obj(m.sum(mix_score))

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

</code></pre>

<p>In the Mix score section, you'll see where I would like to add in conditional variables defined as the <code>b

*[j]</code> that are 1 or 0 depending on the value of <code>x**[j]</code> (1 if <code>x**[j]</code> is non-zero, and 0 otherwise). I would like to apply this binary variable to a count function I'm using in the objective function that is commented out</p>*

<p><code># diversity_score = [counts<p><code># diversity_score = [counts

*/ m.max2(x_sums**, 1e-3) for i in range(30)]</code></p>*

<p><code># m.Obj(m.sum(diversity_score))</code></p>

<p>I'm surprised when I implement the variable <code>b</code> (so simply un-comment the code denoted Mix Score section below and above), the solver fails to find a feasible solution, even when running the dummy objective <code>m.Obj(0)</code> It is as if the way I'm defining <code>b</code> is actually adding a constraint and reducing the feasibility space. Any help would be much appreciated. Thanks!</p><p><code># m.Obj(m.sum(diversity_score))</code></p>

<p>I'm surprised when I implement the variable <code>b</code> (so simply un-comment the code denoted Mix Score section below and above), the solver fails to find a feasible solution, even when running the dummy objective <code>m.Obj(0)</code> It is as if the way I'm defining <code>b</code> is actually adding a constraint and reducing the feasibility space. Any help would be much appreciated. Thanks!</p>