Search code examples
pythonnonlinear-optimizationgekko

Discrete Optimization (SOS1 constraint) - GEKKO


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


Solution

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