OiO.lk Community platform!

Oio.lk is an excellent forum for developers, providing a wide range of resources, discussions, and support for those in the developer community. Join oio.lk today to connect with like-minded professionals, share insights, and stay updated on the latest trends and technologies in the development field.
  You need to log in or register to access the solved answers to this problem.
  • You have reached the maximum number of guest views allowed
  • Please register below to remove this limitation

Adding conditional variable in gekko leads to no solution

  • Thread starter Thread starter jim
  • Start date Start date
J

jim

Guest
I'm using gekko to optimize a certain function. When I use a dummy objective function like 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[0]) for i in range(30)]

w_base = [m.Param(value=w) for i in range(w.shape[0])]
h_base = [m.Param(value=h) for i in range(h.shape[0])]
t_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]):
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[j] for i in range(30))
m.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_constants)

## 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(x[j]*h_base[j] for j in range(len(h_base))) +
h_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(x[j]*h_base[j] for j in range(len(h_base))) +
h_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(x[j]*h_base[j] for j in range(len(h_base))) +
h_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(x[j]*h_base[j] for j in range(len(h_base))) +
h_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(x[j]*h_base[j] for j in range(len(h_base))) +
h_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 x[j] is non-zero or not
## We will use the count of these b[j] values in our objective function

## 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[j] = m.if3(x[j] - epsilon, 0, 1)

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

Latest posts

I
Replies
0
Views
1
impact christian
I
Top