Intermediate Arrays in GEKKO

Question

If I have an optimization problem with a grid, is it better to create arrays of variables or intermediates?. By "better" I am referring primarily to performance with a secondary interest in tractability.

For example, say I have some function H(x,s) that I need to integrate over x and s. If I use a Gauss quadrature I will evaluate H at each node on a large grid. What is the difference between the following approaches?

m = GEKKO()
x,s = get_quadrature_vectors()
nodes_x, nodes_s = np.meshgrid(x,s)

Approach 1:

H = m.Array(m.Var,(n,n))
for i in range(n):
    for j in range(n):
        m.Equation(H[i,j] == evaluate_weighted_H(nodes_x[i,j],nodes_s[i,j]))

Approach 2:

H = m.Array(m.Intermediate,(n,n))
for i in range(n):
    for j in range(n):
        m.Equation(H[i,j] == evaluate_weighted_H(nodes_x[i,j],nodes_s[i,j]))

Approach 3:

H = np.empty((n,n))
H = evaluate_weighted_H(nodes_x,nodes_s)

After evaluating the grid, of course, the weighted results will be summed to form the integral which will then be used in another equation or the objective.

Answer 1

Intermediates generally improve model performance for compile and solution time by explicitly calculating parts of the model. They are also useful to sub-divide equations if the total symbolic expression is over 15,000 characters (equation length limit in APMonitor). Intermediates also allow parts of the model to be evaluated and returned with the solution so that those values can be displayed as a post-processing step. Approach 3 is also fine, but the expression (not a single intermediate variable) is propagated to other parts of the model where elements of H are used.

One correction is that intermediates do not need to be explicitly defined.

H = [[None for _ in range(n)] for _ in range(n)]
for i in range(n):
    for j in range(n):
        H[i,j] == m.Intermediate(evaluate_weighted_H(nodes_x[i,j],nodes_s[i,j]))

Even better, H can be defined using list comprehensions.

H = [[m.Intermediate(evaluate_weighted_H(nodes_x[i,j],nodes_s[i,j])) 
                     for i in range(n)]
                     for j in range(n)]

To get timing values with different options, try using m.options.DIAGLEVEL=1 to display compile and solve time.