How to include penalty in a Objective Function with Python? GEKKO

Question

I'm trying to include a "great M" penalty in my objective function.

I want use the entry x vector values as entry values in a function. A fixed maximum value is took initially for the returned value of this function, and I want to avoid the solutions which, with x values given, returns values higher than the fixed one.

How is the best way to do it?

I wrote some code with an approach but I don't know if it is a good way to do that

#from calculator import calculateConcentration
'GEKKO MODELING'
from gekko import GEKKO
m = GEKKO()
m.options.SOLVER=1  # APOPT is an MINLP solver

a_max = 30

# Initialize variables
x = []
x1 = m.Var(value=20,lb=20, ub=6555)  #integer=True
x2 = m.Var(value=0,lb=0,ub=10000)  #integer=True
x3 = m.sos1([30, 42, 45, 55])

x = [x1, x2, x3]
# Equations
m.Equation((x1 * x2* x3) * 10 ** (-6)>=50)

def fun(x):
    return 44440 + ((np.pi * x[0] * x[1] * x[2]) * 10 ** (-4))**0.613 #+ penalty(x)

#def penalty(x):
#    a = calculateConcentration(x)
#    if (a>a_max):
#        return 10**10
#    else:
#        return 0

x = [400,300,19]

'GEKKO Optimization'
m.Obj(fun(x))

m.solve(disp=False) # Solve
print('Results')
print('x1: ' + str(x1.value))
print('x2: ' + str(x2.value))
print('x3: ' + str(x3.value))

print('Objective: ' + str(m.options.objfcnval))

Answer 1

I think you should stick with yyour current formulation. If you want to use a penalty, I would augment your function. Calculate your function itself, and then add a penalty of the form $$dJ = -(y - y_{max})^3$$. This will heavily punish any values of $y > y_{max}$. But this will tell your optimizer to seek low values of y. The best choice is if your optimizer offers bounds/constraints, enter $y < y_{max}$ as your constraint.

Answer 2

I made a mistake in the concept formulation. If you are using an optimization algorithm with constraints, then you just need to develop properly the constraints, as functions of the variables.

In this case, I was worried because some constraints are functions of the different problem variables, but another constrains were a result of properties values, which depend on the optimization problem variables. Everything should be in the constraints, and in this case, should be as follows:

'Generic packages'
import numpy as np
import time
inicio = time.time()

'GEKKO MODELING'
from gekko import GEKKO


A_max = 5

class PContr(object):
    def __init__(self, A, flA, cpB):
        self.A = A
        self.flA = flA
        self.cpB = cpB
    'Ps' 
    def PA(self, x):
        if self.A < 10:
            dPA = 8.2*10**(-3)*x[2] + 0.14
            return dPA
        else:
            dPA = 3.2*10**(-2)*x[2] + 0.14 
            return dPA

m = GEKKO()
x1 = m.Var(value=20,lb=20, ub=6555)  #integer=True
x2 = m.Var(value=1,lb=1,ub=100)  #integer=True
x3 = m.sos1([30, 42, 45, 55])
x3.value = 1.0

x = [x1, x2, x3]

obj = PContr(1.20, 300, 1.40)
dPA = obj.PA(x)

'Process Constraints'
m.Equation(dPA<=A_max)


def fun(x):
    return 22223 + (x[0] * x[1] * x[2])**0.83

m.Obj(fun(x))

# Change to True to initialize with IPOPT
init = False
if init:
    m.options.SOLVER=3  
    m.solve(disp=False) # Solve

m.options.SOLVER=1
m.solve(disp=True) # Solve

print('Results')
print('x1: ' + str(x1.value))
print('x2: ' + str(x2.value))
print('x3: ' + str(x3.value))
print('Objective: ' + str(m.options.objfcnval))

##Execution time
final = time.time()
print('Tiempo de ejecucion:', (final-inicio), 's')