Python GEKKO Objective Function and Integer Constraint Issues

[btda...@gmail.com]

I've just started using GEKKO in Python and like it a lot so far - I just have a couple of questions.

First, if I try to use any built in functions in my objective function declaration, I'm getting errors. For example:

from gekko import GEKKO

import numpy as np

m = GEKKO()

m.options.SOLVER=1

Variables = ['x1','x2','x3','x4']

# intiliaze GEKKO variables for NLP

Variable_Dict = {j : m.Var(value=1,lb=1,ub=5,integer=True) for j in Variables}

# per GEKKO documentation, if multiple objectives are given, they are summed

for i in Variables:

    test=np.exp(Variable_Dict[i])

    m.Obj(test)

This codes gives me "AttributeError: exp" which I can get around it by using 2.7182818284591**test instead. Everything runs fine from there on.

However, I'm also trying to use the scipy.stats CDF for the standard normal distribution in another problem and I'm getting "TypeError: object of type 'int' has no len()". Is there a way to use these built in functions in the objective function like the below code?

As for the integer constraints being ignored, my decision variables in a separate problem are declared using a dictionary, but from a Pandas DataFrame instead of a list like in my above example. Below is my code:

Vars = {i : m.Var(value=1,lb=1,ub=20,integer=True) for i in PandasDF.index}

I'm getting non-integer solutions for this example, but my objective function is being evaluated correctly and my constraints are satisfied, so I don't think these are causing any issues. Is there any concern with declaring the variables this way? This doesn't seem to be giving me any issues in small problems, but I thought I'd check to make sure.

Thanks!

-Brian

[John Hedengren]

Brian,

 

The only thing you need to do is replace np.exp with m.exp. You have to use gekko functions so that gradients can be computed and passed to the solver.

 

from gekko import GEKKO

import numpy as np

 

m = GEKKO()

m.options.SOLVER=1

 

Variables = ['x1','x2','x3','x4']

 

# intiliaze GEKKO variables for NLP

Variable_Dict = {j : m.Var(value=1,lb=1,ub=5,integer=True) for j in Variables}

 

 

# per GEKKO documentation, if multiple objectives are given, they are summed

for i in Variables:

    test=m.exp(Variable_Dict[i])

    m.Obj(test)

 

m.solve()

 

 

If there are functions that are not currently supported in GEKKO then I recommend one of three options:

 

1.       Request that the function be added as a feature request in Github. https://github.com/BYU-PRISM/GEKKO

2.       Use a cspline (1D) or bspline (2D) to approximate the function CDF. See https://apmonitor.com/wiki/index.php/Main/GekkoPythonOptimization (Problem #4) for an example of a cspline.

3.       Switch to a package such as Scipy.optimize.minimize that allows you to use Numpy or Scipy functions. Here are tutorials on both: http://apmonitor.com/che263/index.php/Main/PythonOptimization (See methods #2 and #3).

 

Best regards,

 

John Hedengren

--
--
APMonitor user's group e-mail list.
- To post a message, send email to apmo...@googlegroups.com
- To unsubscribe, send email to apmonitor+...@googlegroups.com
- Visit this group at http://groups.google.com/group/apmonitor

---
You received this message because you are subscribed to the Google Groups "apmonitor" group.
To unsubscribe from this group and stop receiving emails from it, send an email to apmonitor+...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[btda...@gmail.com]

Thanks for the quick response! I tried scipy.optimize.minimize but it didn't look like it could handle integer constraints. I will look into the approximations.

What about the second part? Sorry - I didn't make it very clear that it was two different issues.

As for the integer constraint issue, my decision variables in a separate problem are declared using a dictionary, but from a Pandas DataFrame instead of a list like in my above example. Below is my code:

[Logan Beal]

Defining the decision variables in a dictionary like that should work. If you're getting non-integer solutions double check that you're using the APOPT solver (m.options.SOLVER=1) and that the solver claims a successful solution. If both of these are true, please send more information so we can debug this. 

[btda...@gmail.com]

Thanks! I realized that the bounds on that variable were making the problem infeasible, so changing that ended up fixing the problem. The problem works on a small scale, but now that I have around 33,000 decision variables, I get the following output:

apm 151.181.70.82_gk_model4 <br><pre> ----------------------------------------------------------------

 APMonitor, Version 0.8.4

 APMonitor Optimization Suite

 ----------------------------------------------------------------

 

 

 --------- APM Model Size ------------

 Each time step contains

   Objects      :            0

   Constants    :            0

   Variables    :            0

   Intermediates:            0

   Connections  :            0

   Equations    :            0

   Residuals    :            0

 

 

 Model must contain at least 1 variable and 1 equation

 STOPPING...

When I print the file path and go there, the APM file appears to be generated correctly so I'm not sure why it's reporting 0 variables. It will run with about half the variables, but not for the full problem of 33k.

Thanks for all your help!

-Brian

[John Hedengren]

Brian,

In the APM file you should see the following sections:

Variables

Equations

Here is an example model and description of the different sections: http://apmonitor.com/do/index.php/Main/ModelFormulation

Could you check that GEKKO is correctly including these headers? I'm wondering if there is a file line issue or string length problem with Python. If this doesn't lead to the problem, could you send your APM file so that I can have a look at it? I've solved problems with 100,000+ variables in GEKKO with no issues. You may need to switch to the IPOPT or BPOPT (Interior point) solver if the number of degrees of freedom number is very large. There are also model reduction strategies that are listed in the GEKKO paper that can help for cases that have a lot of data and lead to a large model: https://www.mdpi.com/2227-9717/6/8/106

Best regards,

John Hedengren

To post to this group, send email to apmo...@googlegroups.com.
Visit this group at https://groups.google.com/group/apmonitor.

For more options, visit https://groups.google.com/d/optout.

--

Best regards,

John Hedengren

APMonitor Optimization Suite

http://apmonitor.com

[btda...@gmail.com]

Thanks for the quick reply!

I checked my APM file (gk_model4, attached), and the headers are there for the variables and equations. I've also attached the file for the model that worked with half as many variables (gk_model15).

I've already switched over to IPOPT since the model is no longer a MINLP and it has a large number of degrees of freedom. I tried BPOPT and got the same output.

-Brian

[John Hedengren]

Brian,

 

When I run your model with the local executable (apm.exe from https://github.com/APMonitor/apm_server/tree/master/apm_windows/online), I get the following for your case that originally failed:

 

C:\Users\johnh\Desktop\gk_model4>apm gk_model4.apm

----------------------------------------------------------------

APMonitor, Version 0.8.3

APMonitor Optimization Suite

----------------------------------------------------------------

 

 

--------- APM Model Size ------------

Each time step contains

   Objects      :  0

   Constants    :  0

   Variables    :  33194

   Intermediates:  0

   Connections  :  0

   Equations    :  33194

   Residuals    :  33194

 

My guess is that GEKKO is loading the large model file to the server (35 MB) and it is not able to even start loading the Equations section by the time the connection times out. Here are two suggestions:

 

·       If your model is scaling up with the number of data points, you should consider using IMODE = 2 instead. This way you only define one set of model equations and objective function and then the computational engine multiplies your model for you for each of the data sets. Your model will be very small but it is still very efficient for large amounts of data. Here is an example of IMODE = 2 for regression:

 

from gekko import GEKKO

import numpy as np

import matplotlib.pyplot as plt 

 

# measurements

xm = np.array([0,1,2,3,4,5])

ym = np.array([0.1,0.2,0.3,0.5,0.8,2.0])

 

# GEKKO model

m = GEKKO()

 

# parameters

x = m.Param(value=xm)

a = m.FV()

a.STATUS=1

 

# variables

y = m.CV(value=ym)

y.FSTATUS=1

 

# regression equation

m.Equation(y==0.1*m.exp(a*x))

 

# regression mode

m.options.IMODE = 2

 

# optimize

m.solve(disp=False)

 

# print parameters

print('Optimized, a = ' + str(a.value[0]))

 

plt.plot(xm,ym,'bo')

plt.plot(xm,y.value,'r-')

plt.show()

 

·       The other suggestion is to use the option remote=False so that you don’t need to send the problem to a remote server and where you run into server issues such as a timeout for large files. I’d recommend that you try the first suggestion initially and then use this suggestion if it is not possible.

 

-John Hedengren

[btda...@gmail.com]

Hi John,

Thanks again for your help on this. I'm not sure I understand how/if using IMODE = 2 can help with my problem, but I tried running locally and this allowed me to run larger versions of my model, but not as large as I need because of memory restrictions. I'm also wondering if I might be writing some of the code in a way that makes the APM file unnecessarily large. I've attached my code for reference - does it look like I could be writing this differently to achieve a smaller APM file?

Thanks,

Brian

[John Hedengren]

Brian,

It looks like your application would benefit from IMODE=2. You would define each equation or objective function once and then load your data into the .value for each of your parameters or variables. All of the .value arrays need to be the same length. If there is a parameter or variable that is common for all data points, such as a regression parameter, then you declare it as an FV (m.FV()) and turn the status on such as:

a = m.FV()

a.STATUS = 1

Let me know if you need any help with trying out IMODE 2.

Best regards,

John Hedengren

[btda...@gmail.com]

Thanks, John!

Here is what I have so far. Everything is much cleaner, but I'm getting a strange error where it's telling me I'm missing an opening parenthesis even though they all seem to be accounted for.

apm 151.181.70.82_gk_model31 <br><pre> ----------------------------------------------------------------

 APMonitor, Version 0.8.4

 APMonitor Optimization Suite

 ----------------------------------------------------------------

 

 

 --------- APM Model Size ------------

 Each time step contains

   Objects      :            0

   Constants    :            0

   Variables    :            6

   Intermediates:            0

   Connections  :            0

   Equations    :            2

   Residuals    :            2

 

 @error: Model Expression

 *** Error in syntax of function string: Missing opening parenthesis

[John Hedengren]

Brian,

 

You were naming one of your variables starting with a reserved keyword (cos()) and it was looking for the opening of the cosine operator. The original error message had a very long string with spaces that wrapped around onto multiple lines:

 

Position: 433                

 ((((((v1-1))*(forecast))-(short-(((((((((((v1-1))*(forecast)))/(sd)))*((1-((exp((((-(((((((v1-1))*(forecast)))/(sd)))^(2))))/(2))))/(((0.5569620253164557+((1.6)*((((((v1-1))*(forecast)))/(sd)))))+((0.8333333333333334)*(sqrt(((((((((v1-1))*(forecast)))/(sd)))^(2))+3))))))))))+((0.3989422804014327)*(exp((((((-(((((v1-1))*(forecast)))/(sd))))*((((((v1-1))*(forecast)))/(sd)))))/(2))))))-(((((v1-1))*(forecast)))/(sd))))*(sd)))))*(cost))

                                                                                                                                                                                                                                                                                                                                                                                                                                                 ?

 

 

If we break up the long expression into multiple smaller expressions with Intermediates then it is easier to see:

 

def Expected_Short(Factor,SD,Forecast):

    part1 = m.Intermediate(Service_Factor(Factor, SD, Forecast)*CDF(Service_Factor(Factor, SD, Forecast)))

    part2 = m.Intermediate(PDF(Service_Factor(Factor, SD, Forecast)))

    part3 = m.Intermediate(Service_Factor(Factor, SD, Forecast))

    return (part1 + part2 - part3) * SD

 

The new error points to the “cost” variable that starts with “cos”.

 

@error: Model Expression

*** Error in syntax of function string: Missing opening parenthesis

Position: 62                 

 ((((((v1-1))*(forecast))-(short-((((i3+i4)-i5))*(sd)))))*(cost))

                                                              ?

 

If you don’t name your variables then it isn’t a problem. When you use names, please check your “names” for starting with reserved keywords: http://apmonitor.com/wiki/index.php/Main/Equations  Perhaps we should automatically pre-pend the names in GEKKO with something like “_cost” to avoid that.

 

abs()	Absolute value	abs(x*y)=0
exp()	Exponentiation	exp(x*y)=0
log10	Base-10 Log	log10(x*y)=0
log	Natural Log	log(x*y)=0
sqrt()	Square Root	sqrt(x*y)=0
sinh()	Hyperbolic Sine	sinh(x*y)=0
cosh()	Hyperbolic Cosine	cosh(x*y)=0
tanh()	Hyperbolic Tangent	tanh(x*y)=0
sin()	Sine	sin(x*y)=0
cos()	Cosine	cos(x*y)=0
tan()	Tangent	tan(x*y)=0
asin()	Arc-sine	asin(x*y)=0
acos()	Arc-cos	acos(x*y)=0
atan()	Arc-tangent	atan(x*y)=0
erf()	Error function	erf(x*y)=0
erfc()	Complementary error function	erfc(x*y)=0

 

If you start your names with one of these, it will trigger an error that it cannot find an opening parenthesis.

 

Best regards,

[btda...@gmail.com]

That makes sense - thanks!

I hate to keep bugging you with questions, but now I'm having trouble with constraints. I'm trying to add a constraint to limit the total available product to allocate, but it's just printing "True" or "False" in the APM file. It seems like it's using the initial value for each variable and checking if the equation is true rather than using it as a constraint.

m.Equation(sum(factor.value * forecast.value) <= bound)

This just produces "True" in the APM file and gives me the following error:

 Equation without an equality (=) or inequality (>,<)

 true

 STOPPING...

I'm also trying to give the m.Var function an array for the upper and lower bounds without much success either. I was hoping to adjust the lower bounds by line for the input file by sending the "lbound" array to the m.Var function. I tried an array of 1's first to test it - here is my code/error for that:

starting_point = np.ones(len(data))

lbound = np.ones(len(data))

ubound = np.ones(len(data))*3

factor = m.Var(value = starting_point, lb=lbound, ub=ubound, integer=False)

 *** Error in syntax of function string: Missing operator

Position: 4                   

 1.1.1....1.1.1.

    ?

Thanks again for all your help, I really appreciate it!

-Brian

[John Hedengren]

Brian,

For the first problem, you can just drop the .value when you define the equation:

m.Equation(sum(factor * forecast) <= bound)

You may need to define an Intermediate variable array and then do the summation on that list such as:

p = [m.Intermediate(factor[i]*forecast[i]) for i in range(n)]

m.Equation(sum(p) <= bound)

For the second issue, I recommend that you assign the upper an lower bounds as follows if you are using IMODE = 2:

factor = m.Var(value = starting_point, integer=False)

lb = m.Param(value=lbound)

ub = m.Param(value=ubound)

m.Equation(factor>=lb)

m.Equation(factor<=ub)

Unfortunately, there is no way to feed in an array of lower and upper bounds to the data file right now so you'll need to define additional parameters and equations.

Best regards,

John Hedengren

[btda...@gmail.com]

Hi John,

I hope you had a nice Thanksgiving!

I'm still having some problems with constraints using IMODE=2. Here is a quick example I mocked up to test my understanding:

from gekko import GEKKO

import numpy as np

m = GEKKO(remote=False)

m.options.SOLVER = 3

m.options.IMODE = 2

xm = np.array([1,1,1,1,1,1])

ym = np.array([2,2,2,2,2,2])

x = m.Var(value=xm,lb=1,ub=5,integer=False)

y = m.Param(value=ym)

m.Obj(x*y)

test = [m.Intermediate(x[i]) for i in range(len(xm))]

m.Equation(sum(test) <= 7)

# Solve

m.solve(disp=True)

for i in x:

    print(i)

If I run this as is, it runs as expected with no errors and the solution is [1,1,1,1,1,1] with an objective function value of 12. If I change the constraint to greater than or equal to 7, then I get an error that the LP is infeasible. I tried changing the starting point to [2,1,1,1,1,1] with the greater than or equal to 7 constraint and that ran with no errors. However, it reported an optimal solution of [1,1,1,1,1,1] with objective function value of 12 which does not satisfy the constraint. 

After some more testing, it appears that the solver is only checking that the constraint is satisfied by the initial starting point. If the constraint is satisfied by the initial starting point, it proceeds while ignoring it. If it is not satisfied by the initial starting point, then it says it's infeasible. Am I doing something incorrectly?

-Brian

[John Hedengren]

Brian,

You can sum across variables that are in the same model instance but not across the data dimension. Here is an example of a summation across the model variables.

from gekko import GEKKO

import numpy as np

 

m = GEKKO(remote=False)

m.options.SOLVER = 3

m.options.IMODE = 2

 

xm = np.array([1,1,1,1,1,1])

ym = np.array([2,2,2,2,2,2])

 

# define x[0], x[1], x[2]

x = [m.Var(value=xm,lb=i+1,ub=5,integer=False) for i in range(3)]

 

y = m.Param(value=ym)

 

# define parameter that is the same across all data points

k = m.FV(lb=2,ub=4)

k.STATUS = 1

 

m.Obj(-k*sum(x)*y)

 

test = [m.Intermediate(x[i]) for i in range(3)]

m.Equation(sum(test) <= 7)

 

# Solve

m.solve(disp=True)

 

for i in x:

    print(i.value)

print('k: ' + str(k.value[0]))

 

If you do need to solve across the data dimension, I recommend that you use a dynamic mode such as IMODE = 6 and integrate the variable. Let me know if you need this and I can show you an example, similar to this one: https://apmonitor.com/do/index.php/Main/IntegralObjective  You wouldn't have a differential equation but this would allow you to integrate the across the data dimension.

Best regards,

John Hedengren

[btda...@gmail.com]

Sorry for the confusion - I thought the formulation with the intermediates was for IMODE=2. I was able to define constraints before I started using IMODE=2 (similar to the example in the body of your response), but had trouble once I switched. 

Do you have an example where constraints are defined across the data dimension using IMODE=6? That sounds like what I would have to do since I'm trying to relate subsets of the individual decision variables within the x variable in the example I posted (or the factor variable in the model I'm working on outside of that example).

Thanks,

Brian

[John Hedengren]

Brian,

As you probably discovered, the most flexible mode is IMODE=3 because it lets you define individual objectives, equations, variables and the relationship between them. The issue with your model is that you have a lot of data points and the model takes a long time to compile versus when you use IMODE=2 where it compiles the kernel model and then propagates it to all data points. I suppose that you want to use IMODE=3 if you need to define relationships in the data dimension that don't fit into the IMODE=2 or IMODE=6 structure. There are many examples of IMODE=6 at http://apmonitor.com/do/index.php/Main/DynamicOptimizationBenchmarks and http://apmonitor.com/do/index.php/Main/MoreDynamicOptimizationBenchmarks  Another option for IMODE=2 is to define data columns of parameters that include the subsets that you'd like to include such as:

p = [0,1,1,1,0,0,0,0]

y = sum([p[i]*data[i] for i in range(n)])

Best regards,

John Hedengren

[btda...@gmail.com]

Hi John,

Yes, IMODE=3 worked very well for small scale tests of my model. I will check out those other examples using IMODE=6.

Also, wouldn't the formulation using the p vector of 1's and 0's end up trying to sum across the data dimension? It looks like it would end up being similar to the example I posted, but multiplying the x variable by a constant within the intermediate to isolate a given subset of x unless I'm not completely understanding what you're doing with the p vector.

Thanks,

Brian

[John Hedengren]

​Brian,

You are correct - it won't work with IMODE=2. I think you could do this with IMODE=6 by integrating a parameter such as #11 here:

https://apmonitor.com/wiki/index.php/Main/GekkoPythonOptimization  There is also an example of connections with specific a data horizon point in #12.

Best regards,

John Hedengren

---

From: apmo...@googlegroups.com <apmo...@googlegroups.com> on behalf of btda...@gmail.com <btda...@gmail.com>
Sent: Thursday, November 29, 2018 8:59 PM
To: apmonitor