3

I'm trying to define an optimization problem with GEKKO in Python, and I want to use some design variables with predefined list of choices. Also, each choice has an associated cost and the constraint would be that the total cost should be under a specified limit.

Below is a common gekko example(found here) with the modification that x1 and x2 are sos1. Also with the index of the selected values of x1 and x2, I find their associated cost from another list and their sum should be less than a certain value(constraint). 

from gekko import GEKKO
def test(x1,x2,x3,x4):
    res = x1*x4*(x1+x2+x3)+x3
    return res

def check(x1,x2):
    tt = [1,2,3,4,5]
    cost = [10,10,10,2,1]
    if x1.value in tt:
        y1 = tt.index(x1.value)
        y2 = tt.index(x2.value)
        C = cost[y1]+cost[y2]
        return C
    return 10

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

# optional solver settings with APOPT
m.solver_options = ['minlp_maximum_iterations 500', \
                    # minlp iterations with integer solution
                    'minlp_max_iter_with_int_sol 10', \
                    # treat minlp as nlp
                    'minlp_as_nlp 0', \
                    # nlp sub-problem max iterations
                    'nlp_maximum_iterations 50', \
                    # 1 = depth first, 2 = breadth first
                    'minlp_branch_method 1', \
                    # maximum deviation from whole number
                    'minlp_integer_tol 0.05', \
                    # covergence tolerance
                    'minlp_gap_tol 0.01']

# Integer constraints for x3 and x4
x3 = m.Var(value=1,lb=1,ub=5,integer=True)
x4 = m.Var(value=2,lb=1,ub=5,integer=True)
x1 = m.sos1([1,2,3,4,5])
x2 = m.sos1([1,2,3,4,5])

# Equations
m.Equation(x1*x2*x3*x4>=25)
m.Equation(x1**2+x2**2+x3**2+x4**2==40)
m.Equation(check(x1,x2)<=5)
m.Obj(test(x1,x2,x3,x4)) # 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))

Note: I had to add an if block in the check function as the initial value of x1 and x2 seems to be zero.

This code doesn't work and I get the following error. 

> Exception has occurred: Exception
 @error: Equation Definition
 Equation without an equality (=) or inequality (>,<)
 true
 STOPPING...

I do not know what causes this error. How should I reformulate my model to get the desired result?

Edit: This example code is just my attempt to recreate the error. My actual application is designing an engineering system. For example, let's say the system has 2 components - battery and bulb. I have two options for the battery, Battery A weighs 10kg and its reliability is 0.97 and Battery B weighs 6kg and its reliability is 0.75. Similarly, there are different options for the bulb. I need to select a choice for the battery and the bulb such that the total system reliability is as high as possible(objective) and the total weight is less than 'x' kg(constraint). In the above code, think of x1 and x2 values as selected choices for the components and I find their index to get their associated weight/cost(if Battery A and Bulb B was chosen, I get their weights to check if the total weight is less than the allowed limit). Now my actual system has n components and m choices for each component. And each choice has associated weight, cost, reliability, etc. I'm trying to find the optimal combination to maximize the system reliability with constraints on system weight, cost, etc

user7440787
  • 831
  • 7
  • 22
Skyrider
  • 133
  • 6
  • 1
    The problem has to do with line `m.Equation(check(x1,x2)<=5)` as `gekko` doesn't seem to like using a defined function. I would be surprised if any solver would allow you to do that. Could you give a bit more context to the MINLP you're trying to solve? why are x values always integers between 1 and 5 and why do you have to use x1 and x2 to index a value in another array? – user7440787 Apr 23 '20 at 20:10
  • 2
    @user7440787 I have edited the post to include more information. – Skyrider Apr 23 '20 at 21:32

2 Answers2

2

I built a simple model based on your example description.

from gekko import GEKKO
import numpy as np

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

# optional solver settings with APOPT
m.solver_options = ['minlp_maximum_iterations 500', \
                    # minlp iterations with integer solution
                    'minlp_max_iter_with_int_sol 10', \
                    # treat minlp as nlp
                    'minlp_as_nlp 0', \
                    # nlp sub-problem max iterations
                    'nlp_maximum_iterations 50', \
                    # 1 = depth first, 2 = breadth first
                    'minlp_branch_method 1', \
                    # maximum deviation from whole number
                    'minlp_integer_tol 0.05', \
                    # covergence tolerance
                    'minlp_gap_tol 0.01']


x1 = m.Array(m.Var, 5, **{'value':0,'lb':0,'ub':1, 'integer': True}) # battery options
print(f'x1_initial: {x1}')
x2 = m.Array(m.Var, 5, **{'value':0,'lb':0,'ub':1, 'integer': True}) # bulb options
print(f'x2_initial: {x2}')
bat_cost = np.array([ 10, 2, 3, 4, 5])  # battery costs
bat_weigh = np.array([ 1, 25, 20, 19, 20])  # battery weighs
bulb_cost = np.array([ 2, 5, 33, 24, 5])  # bulb costs
bulb_weigh = np.array([ 6, 10, 2, 10, 20])  # bulb weighs
m.Equation( sum(bat_weigh * x1) + sum(bulb_weigh * x2) <= 25)  # limit total weigh 
m.Equation(m.sum(x1) == 1)  # restrict choice to a single battery 
m.Equation(m.sum(x2) == 1)  # restrict choice to a single bulb
m.Obj( sum(bat_cost * x1) + sum(bulb_cost * x2) ) # Objective

m.solve(disp=False) # Solve
print('Results:')
print(f'x1: {x1}')
print(f'x2: {x2}')
print(f'battery cost: {sum(np.array([i[0] for i in x1]) * bat_cost)}')
print(f'battery weigh: {sum(np.array([i[0] for i in x1]) * bat_weigh)}')
print(f'bulb cost: {sum(np.array([i[0] for i in x2]) * bulb_cost)}')
print(f'bulb weigh: {sum(np.array([i[0] for i in x2]) * bulb_weigh)}')
print('Objective value: ' + str(m.options.objfcnval))

The result is the following:

x1_initial: [0 0 0 0 0]
x2_initial: [0 0 0 0 0]
Results:
x1: [[0.0] [0.0] [0.0] [1.0] [0.0]]
x2: [[1.0] [0.0] [0.0] [0.0] [0.0]]
battery cost: 4.0
battery weigh: 19.0
bulb cost: 2.0
bulb weigh: 6.0
Objective value: 6.0

This is a very simple example to show how to represent the battery and bulb information. It can be made more complex but I would need more details and understand why you have the polynomial equations, what they represent.

And just to reiterate, the error you're getting, has to do with line:

m.Equation(check(x1,x2)<=5)
user7440787
  • 831
  • 7
  • 22
2

In addition to the good answer from user7440787, you need to lookup multiple values the from a pre-defined set of discrete design variables. Instead of using the predefined m.SOS1() function, you can use something like the following to tie one binary decision variable array to multiple correlations or variables.

from gekko import GEKKO
m = GEKKO(remote=False)
# design variable
y = m.Var(lb=1,ub=5)
# options
n = 4
# weight
weight=[19.05-y, 25.0-0.1*y**2, 29.3-0.02*y**3, 30.2]
# cost
cost = [3.2+y,2.4+0.01*y**2,1.6+y+0.001*y**3,5.2]
# SOS1 with binary variables
b = m.Array(m.Var,n,lb=0,ub=1,integer=True)
m.Equation(m.sum(b)==1) # only select one
# cost x weight
cxw = m.sum([b[i]*cost[i]*weight[i] for i in range(4)])
# minimize cost x weight
m.Minimize(cxw)
# change to APOPT solver
m.options.SOLVER = 1
m.solve(disp=False)
print('Design Variable: ' + str(y.value[0]))
print('Option: ' + str(b))

In this example, you have one design variable y and different equations for cost and weight that are based on the design variable. The overall objective is to minimize the product of cost and weight while adjusting y.

John Hedengren
  • 12,068
  • 1
  • 21
  • 25